'nlpnra' command

[telescopic]

Hi, I encounter problem again. I can't debug the error for 'nlpnra' command. I am trying to max the loglikelihood function. It returns the following error(see below). Is there something wrong with the argument used in 'nlpnra'?

Thanx

Elvin

Error message-----------------------------------------------------

                         Newton-Raphson Optimization with Line Search

                   Minimum Iterations                                     0

                   Maximum Iterations                                   200

                   Maximum Function Calls                               500

                   ABSGCONV Gradient Criterion                      0.00001

                   GCONV Gradient Criterion                            1E-8

                   GCONV2 Gradient Criterion                              0

                   ABSFCONV Function Criterion                            0

                   FCONV Function Criterion                    2.220446E-16

                   FCONV2 Function Criterion                              0

                   FSIZE Parameter                                        0

                   ABSXCONV Parameter Change Criterion                    0

                   XCONV Parameter Change Criterion                       0

                   XSIZE Parameter                                        0

                   ABSCONV Function Criterion                  1.340781E154

                   Line Search Method                                     2

                   Starting Alpha for Line Search                         1

                   Line Search Precision LSPRECISION                    0.9

                   DAMPSTEP Parameter for Line Search                     .

                   MAXSTEP Parameter for Line Search                      0

                   FD Derivatives: Accurate Digits in Obj.F    15.653559775

                   Singularity Tolerance (SINGULAR)                    1E-8

                   Constraint Precision (LCEPS)                        1E-8

                   Linearly Dependent Constraints (LCSING)             1E-8

                   Releasing Active Constraints (LCDEACT)                 .

                         Newton-Raphson Optimization with Line Search

                                   Without Parameter Scaling

                            Gradient Computed by Finite Differences

                          CRP Jacobian Computed by Finite Differences

                              Parameter Estimates               2

                              Lower Bounds                      2

                              Upper Bounds                      2

                                      Optimization Start

Active Constraints                           0  Objective Function                -92.77402163

Max Abs Gradient Element          37.072886626

ERROR: (execution) Invalid argument to function.

operation : LOG at line 4298 column 30

operands  : lamda

lamda      1 row       1 col     (numeric)

-1.38E-17

statement : ASSIGN at line 4298 column 5

traceback : module LOGLIK at line 4298 column 5

                                                       Objective   Max Abs            Slope of

                 Function       Active      Objective   Function  Gradient     Step     Search

Iter   Restarts     Calls  Constraints       Function     Change   Element     Size  Direction

   0*         0        29            0      -92.77402          0   37.0729        0    -237202

                                                X1                X2

                           Parms               0.2               0.3

---Code------------------------------------------------

Data radiotherapy;

   input y age gender $ censor_ind;

datalines;

7 68 f 1

9 69 f 1

12 68 f    1

12 71 f    1

19 77 m    1

23 70 f    1

24 67 f    1

24 68 m    1

24 88 m    1

24 89 m    1

29 28 m 0

34 73 m 1

41 60 f 1

54 60 f 1

72 44 m 0

78 82 f 1

80 62 f 0

83 53 f 0

92 66 f 0

139 63 f 0

139 68 m 0

;

RUN;

PROC IML;

RESET LOG PRINT;

USE radiotherapy VAR _ALL_;

READ all VAR{y} INTO Y;

READ all VAR{age} INTO AGE;

READ all VAR{gender} INTO GENDER;

READ all VAR{censor_ind} INTO CENSOR_IND;

CLOSE radiotherapy;

start LogLik(param) global (Y,CENSOR_IND);

   lamda = param[1];

   gamma = param[2];

   n = nrow(Y);

   f=0;

    DO i=1 to nrow(Y);

    a =CENSOR_IND* (   log(lamda)*gamma+log(gamma)+(gamma-1)*log(Y)     )    ;

    b=(-(lamda*Y)##gamma) ;

    f=a+b+f;

    END;

    return ( f );

finish;

p={0.2 0.3};

ans=LogLik(p);

con={0 0, 1 1};

optn ={1 4};

call nlpnra(rc, xres, "LogLik", p, optn, con);

---------------------------------------------------------------------------------------------

[Rick_SAS]

The error is in your constraints. When lamda=0, your log-likelihood function is undefined. Even when lamda gets very small (-1.38E-17) in your example), log(lamda) numerically overflows.

You can bound the constraint for lada away from 0.  Also, rewrite the LogLik function to eliminate the loop over the elements of  and Censor_Ind.  It might look something like this (untested):

start LogLik(param) global (Y,CENSOR_IND);

   lamda = param[1];

   gamma = param[2];

   a = CENSOR_IND # ( log(lamda)*gamma+log(gamma)+(gamma-1)*log(Y) );

   b = -(lamda*Y)##gamma ;

   f = sum(a + b);

   return ( f );

finish;

con={1e-8 0, 1 1};

[telescopic]

Thank you. This work is about survival analysis. The boldedpart(see Code) below is the K-M estimator graph created by PROC lifetest.

I am fitting a parametric Weibull survival function (S_y(t)=exp{-(lamda*t)^gamma})onthe data. Now having the estimated gamma and lamda.

To plot the Weibull survival function, I will create thefunction and call the function for different time t, then plot using IML thevector.

My question: How to combine the plot with the graph createdby PROC lifetest? Any suggestion?

Code------------------------------------------------

Data radiotherapy;

   input x age gender $ censor_ind;

datalines;

7 68 f 1

9 69 f 1

12 68 f    1

12 71 f    1

19 77 m    1

23 70 f    1

24 67 f    1

24 68 m    1

24 88 m    1

24 89 m    1

29 28 m 0

34 73 m 1

41 60 f 1

54 60 f 1

72 44 m 0

78 82 f 1

80 62 f 0

83 53 f 0

92 66 f 0

139 63 f 0

139 68 m 0

;

RUN;

PROC PRINT data=radiotherapy;

RUN;

ODS GRAPHICS ON;

PROC LIFETEST DATA=    radiotherapy PLOTS=s(ATRISK );

TIME x*censor_ind(0);

RUN;

ODS GRAPHICS OFF;

PROC LIFETEST DATA=    radiotherapy OUTSURV=a;

TIME x*censor_ind(0);

RUN;

PROC PRINT data=a;

RUN;

PROC IML;

RESET LOG PRINT;

USE radiotherapy VAR _ALL_;

READ all VAR{x} INTO X;

READ all VAR{age} INTO AGE;

READ all VAR{gender} INTO GENDER;

READ all VAR{censor_ind} INTO CENSOR_IND;

CLOSE radiotherapy;

start LogLik(param) global (X,CENSOR_IND);

   lamda = param[1];

   gamma = param[2];

   n = nrow(X);

    a =CENSOR_IND# (   log(lamda)*gamma+log(gamma)+(gamma-1)*log(X)     )    ;

    b=-(lamda*X)##gamma ;

    f=SUM(a+b);

    return ( f );

finish;

p={0.2 0.3};

con={1e-4 1e-4, 4 4};

opt ={1,4};

call nlpnra(rc, xres, "LogLik", p, opt, con);

print xres;

[Rick_SAS]

Well, you obviously know that you can get various output data sets from the procedure.

You can merge that data with your Weibull data, and then use SGPLOT (or the GTL, if it is really complicated) to overlay the curves on the same plot by using multiple SERIES statements.

IF, for some reason, there is not an output data set that contains all of the LIFETEST results that you need, you can always use ODS OUTPUT to generate the data that underlies the K-M plot. For example,

   ods output SurvivalPlot=plot;