topic 'nlpnra' command in SAS/IML Software and Matrix Computations

'nlpnra' command

telescopic — Sat, 29 Oct 2011 20:51:14 GMT

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);

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

'nlpnra' command

Rick_SAS — Mon, 31 Oct 2011 14:56:19 GMT

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};

Re: 'nlpnra' command

telescopic — Mon, 31 Oct 2011 16:53:16 GMT

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;

Re: 'nlpnra' command

Rick_SAS — Mon, 31 Oct 2011 17:36:50 GMT

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;