ERROR: call NLPNRA

JerryLee · Posted 09-26-2011 10:55 AM

I am trying to use call NLPNRA to generate a maximum likelihood estimate. IML gives an ERROR message saying I am taking log of zero. But I checked the arguments in all the log commands, none of them is zero. Below is the output of the program and error message. Thanks for help.

11108 %macro infere(rep1,rep2,mcno);

11109

11110 proc datasets library=JM;

11111 delete parajm&rep1;

11112 run;

11113

11114 proc iml;

11115 u=J(1,2,.);

11116 uu=J(1,2,.);

11117 zero=J(1,2,0);

11118 var=I(2);

11119 call randseed(1);

11120 do h1=1 to &mcno;

11121 norma=randnormal(1,zero,var);

11122 do g1=1 to 2;

11123 u[g1]=cdf("NORMAL",norma[g1]);

11124 end;

11125 uu=uu//u;

11126 end;

11127

11128 uu1=uu[2:nrow(uu),];

11129 create unifor from uu1;

11130 append from uu1;

11131 quit;

11132

11133 %do i=&rep1 %to &rep2;

11134

11135 data enable2fact;

11136 set JM.finaldata1;

11137 if rep=&i;

11138 run;

11139

11140 proc sort data=enable2fact;

11141 by ID;

11142 run;

11143

11144 proc means data=enable2fact;

11145 by id;

11146 var tij;

11147 output out=no1;

11148 run;

11149

11150 data no2;

11151 set no1;

11152 rename tij=N;

11153 if _STAT_='N';

11154 keep tij;

11155 run;

11156

11157 proc iml;

11158 *******************************************define log likelihood;

11159 start F_GLOBAL(x);

11160 p=x[1];

11161 lambda1t=x[2];lambda1=x[3];

11162 phi1=x[4];phi2=x[5];phi3=x[6];phi=x[7];

11163 beta=x[,8:11];

11164

11165 use no2;

11166 read all into n;

11167 use enable2fact;

11168 read all var {tij_b} into tijb;

11169 read all var {tij} into tij;

11170 read all var {trt} into trt;

11171 read all var {factpal} into y;

11172 read all var {intercept} into intcep;

11173 read all var {duration_month} into s;

11174 read all var {death} into D;

11175 use unifor;

11176 read all var {col1 col2} into uu;

11177

11178 m=0;lyd=0;

11179

11180 do i=1 to nrow(n);

11181 m=m+n;

11182 si=s;

11183 **********************************log likelihood of mixed model part;

11184 rando=j(n,n,1);

11185 epsi=I(n);

11186 Ri=j(n,n,1);

11187 tib=tijb[m-n+1:m];

11188 ti=tij[m-n+1:m];

11189 *******************************generate Ri;

11190 do j=1 to n;

11191 do k=j to n;

11192 Ri[k,j]=exp(-phi*((abs(tib-tib))**p));

11193 Ri[j,k]=exp(-phi*((abs(tib-tib))**p));

11194 end;

11195 end;

11196

11197 Vi=phi1*rando+phi2*epsi+phi3*Ri;

11198 yi=y[m-n+1:m];

11199 incepti=intcep[m-n+1:m];

11200 trti=trt[m-n+1:m];

11201 one=J(n,1,1);

11202

11203 xib=incepti||trti||tib;

11204 beta1=J(3,1,1);

11205 beta1[1]=beta[1];beta1[2]=trt*beta[2];beta1[3]=beta[3]+trt*beta[4];

11206 if D=1 then do;

11207

11208 ly=-0.5*(log(det(Vi))+(t(yi-xib*beta1))*inv(Vi)*(yi-xib*beta1));

11209 *************************************************************************log likelihood

11209! for observed events;

11210 if trt=1 then ld=log(lambda1t)-lambda1t*si;

11211 if trt=0 then ld=log(lambda1)-lambda1*si;

11212 lyd=lyd+ly+ld;

11213 end;

11214 *************************************************************************log likelihood for

11214! censored events;

11215 if D=0 then do;

11216 ff=0;

11217 do l=1 to nrow(uu);

11218 if trt=0 then surv=(1/lambda1)*log(1/(1-uu[l,1]));

11219 if trt=1 then surv=(1/lambda1t)*log(1/(1-uu[l,2]));

11220 xbeta=(beta[1]+beta[2]*trt)*one+(beta[3]+beta[4]*trt)*one*surv-(beta[3]+beta[4]*t

11220! rt)*ti;

11221 fyd=exp(-0.5*(t(yi-xbeta))*inv(Vi)*(yi-xbeta));

11222 if surv<si then fyd=0;

11223 ff=ff+fyd;

11224 end;

11225

11226 lyd=lyd+log(ff/nrow(uu))-0.5*log(det(Vi));

11227 end;

11228 end;

11229

11230 f=lyd;

11231

11232 return(f);

11233 finish F_GLOBAL;

11234

11235 ******************Making upper and lower bounds for the parameters;

11236 con=J(3,13,.);

11237 boundint=500;

11238 boundbeta=500;

11239 positive=10**(-10);

11240 varianlower=3;

11241 varianupper=500;

11242 boundphi=10;

11243 boundlam=1;

11244 do i=1 to 2;

11245 do j=2 to 11;

11246 if i=1 then do;

11247 if j<=3 then con[i,j]=positive;

11248 if j>=4 & j<=6 then con[i,j]=varianlower;

11249 if j=7 then con[i,j]=positive;

11250 if j=8 then con[i,j]=positive;

11251 if j>=9 then con[i,j]=-boundbeta;

11252 end;

11253 if i=2 then do;

11254 if j<=3 then con[i,j]=boundlam;

11255 if j>=4 & j<=6 then con[i,j]=varianupper;

11256 if j=7 then con[i,j]=boundphi;

11257 if j=8 then con[i,j]=boundint;

11258 if j>=9 then con[i,j]=boundbeta;

11259 end;

11260 end;

11261 end;

11262

11263 con[3,1]=1;con[3,12]=0;con[3,13]=2;

11264

11265 print con;

11266

11267 ********initial value;

11268 x={2 0.04954 0.07649 270.560 118.793 151.389 0.0284 108.702 12.560 1 -0.5};

11269

11270 ********threshold parameters;

11271 opt=J(1,9,.);

11272 opt[1]=1;

11273 opt[2]=5;

11274 tc=J(1,11,.);

11275 ********maximization;

11276 call nlpnra(rc,xres,"F_GLOBAL",x,opt,con,tc);

11277 create solution from xres;

11278 append from xres;

11279 quit;

11280

11281 proc datasets;

11282 append base=JM.parajm&rep1 data=solution;

11283 run;

11284

11285 %end;

11286

11287 %mend;

11288

11289 %infere(rep1=2,rep2=2,mcno=10000);

NOTE: PROCEDURE DATASETS used (Total process time):

real time 1:50.81

cpu time 7.31 seconds

NOTE: Writing HTML Body file: sashtml132.htm

Directory

Libref WORK

Engine V9

Physical Name C:\Users\Zhigang\AppData\Local\Temp\SAS Temporary Files\_TD2780

Filename C:\Users\Zhigang\AppData\Local\Temp\SAS Temporary Files\_TD2780

Member File

# Name Type Size Last Modified

1 AAA DATA 13312 26Sep11:10:12:46

2 ENABLE2FACT DATA 25600 26Sep11:10:49:31

3 FF1 DATA 5120 26Sep11:10:12:46

4 FF11 DATA 5120 26Sep11:10:36:27

5 FYD DATA 46080 26Sep11:09:56:52

6 FYD11 DATA 46080 26Sep11:09:57:17

7 NO1 DATA 9216 26Sep11:10:49:31

8 NO2 DATA 5120 26Sep11:10:49:31

9 SASMACR CATALOG 5120 26Sep11:09:15:51

10 SOLUTION DATA 9216 26Sep11:09:44:16

11 UNIFOR DATA 164864 26Sep11:10:49:31

12 _FMTDESC DATA 9216 26Sep11:10:01:38

13 _QWFFMT DATA 9216 26Sep11:10:01:38

NOTE: Appending WORK.SOLUTION to JM.PARAJM2.

NOTE: BASE data set does not exist. DATA file is being copied to BASE file.

NOTE: There were 1 observations read from the data set WORK.SOLUTION.

NOTE: The data set JM.PARAJM2 has 1 observations and 11 variables.

Rick_SAS · Posted 09-26-2011 11:32 AM

Because of the macro call, the ERROR message doesn't give any useful information about which LOG call is failing.Two suggestions:

1) For debugging, define %LET i=2; %LET mcno = 10000; and get rid of the other macro statements.

2) You have many LOG statements, but we don't know which temporary argument is going negative. Go through the code and assign real variables to the arguments:

temp1 = ff/nrow(uu);

temp2 = det(Vi);

etc.

Then run the program. When the program fails you'll get a real line number that points to the error, and you'll get the actual name of a variable that is negative. For example, if temp1 is negative, you'll get the error message:

ERROR: (execution) Invalid argument to function.

operation : LOG at line XXXX column CCC

operands : TEMP1

TEMP1 1 row 1 col (numeric)

0

Rick

SAS/IML blog: http://blogs.sas.com/content/iml

PS. Unless I'm mistaken, your first call to PROC IML just fills an (&mcno x 2) matrix with uniform random values. You can just use

uu1=J(&mcno,2,.);

call randgen(uu1, "Uniform");

JerryLee · Posted 09-26-2011 01:00 PM

Hi Rick,

Thanks a lot for the prompt advice. I tried what you suggested and found that the error message refered to temp1=ff/nrow(uu). To find out which step it is equal to zero, I printed out all the values of temp1. The strange thing is that none of them is zero. And the objective function F_GLOBAL(x) can be calculated the at the initial value. It is equal to -431.3269953. Do you know why IML give this error message? Thanks a lot. Below has the new output of the error message.

12092 proc iml;

NOTE: IML Ready

12093 *******************************************define log likelihood;

12094 start F_GLOBAL(x);

12095 p=x[1];

12096 lambda1t=x[2];

12096! lambda1=x[3];

12097 phi1=x[4];

12097! phi2=x[5];

12097! phi3=x[6];

12097! phi=x[7];

12098 beta=x[,8:11];

12099

12100 use no2;

12101 read all into n;

12102 use enable2fact;

12103 read all var {tij_b} into tijb;

12104 read all var {tij} into tij;

12105 read all var {trt} into trt;

12106 read all var {factpal} into y;

12107 read all var {intercept} into intcep;

12108 read all var {duration_month} into s;

12109 read all var {death} into D;

12110 use unifor;

12111 read all var {col1 col2} into uu;

12112

12113 m=0;

12113! lyd=0;

12114

12115 do i=1 to nrow(n);

12116 m=m+n;

12117 si=s;

12118 **********************************log likelihood of mixed model part;

12119 rando=j(n,n,1);

12120 epsi=I(n);

12121 Ri=j(n,n,1);

12122 tib=tijb[m-n+1:m];

12123 ti=tij[m-n+1:m];

12124 *******************************generate Ri;

12125 do j=1 to n;

12126 do k=j to n;

12127 Ri[k,j]=exp(-phi*((abs(tib-tib))**p));

12128 Ri[j,k]=exp(-phi*((abs(tib-tib))**p));

12129 end;

12130 end;

12131

12132 Vi=phi1*rando+phi2*epsi+phi3*Ri;

12133 yi=y[m-n+1:m];

12134 incepti=intcep[m-n+1:m];

12135 trti=trt[m-n+1:m];

12136 one=J(n,1,1);

12137

12138 xib=incepti||trti||tib;

12139 beta1=J(3,1,1);

12140 beta1[1]=beta[1];

12140! beta1[2]=trt*beta[2];

12140! beta1[3]=beta[3]+trt*beta[4];

12141 if D=1 then do;

12142

12143 ly=-0.5*(log(det(Vi))+(t(yi-xib*beta1))*inv(Vi)*(yi-xib*beta1));

12144 *************************************************************************log likelihood

12144! for observed events;

12145 if trt=1 then ld=log(lambda1t)-lambda1t*si;

12146 if trt=0 then ld=log(lambda1)-lambda1*si;

12147 lyd=lyd+ly+ld;

12148 end;

12149 *************************************************************************log likelihood for

12149! censored events;

12150 if D=0 then do;

12151 ff=0;

12152 do l=1 to nrow(uu);

12153 if trt=0 then surv=(1/lambda1)*log(1/(1-uu[l,1]));

12154 if trt=1 then surv=(1/lambda1t)*log(1/(1-uu[l,2]));

12155 xbeta=(beta[1]+beta[2]*trt)*one+(beta[3]+beta[4]*trt)*one*surv-(beta[3]+beta[4]*t

12155! rt)*ti;

12156 fyd=exp(-0.5*(t(yi-xbeta))*inv(Vi)*(yi-xbeta));

12157 if surv<si then fyd=0;

12158 ff=ff+fyd;

12159 end;

12160 temp1=ff/nrow(uu);

12160! temp2=det(Vi);

12161 lyd=lyd+log(temp1)-0.5*log(temp2);

12162 print temp1;

12163

12164 end;

12165 end;

12166

12167 f=lyd;

12168

12169 return(f);

12170 finish F_GLOBAL;

NOTE: Module F_GLOBAL defined.

12171

12172 ******************Making upper and lower bounds for the parameters;

12173 con=J(3,13,.);

12174 boundint=500;

12175 boundbeta=500;

12176 positive=10**(-10);

12177 varianlower=3;

12178 varianupper=500;

12179 boundphi=10;

12180 boundlam=1;

12181 do i=1 to 2;

12182 do j=2 to 11;

12183 if i=1 then do;

12184 if j<=3 then con[i,j]=positive;

12185 if j>=4 & j<=6 then con[i,j]=varianlower;

12186 if j=7 then con[i,j]=positive;

12187 if j=8 then con[i,j]=positive;

12188 if j>=9 then con[i,j]=-boundbeta;

12189 end;

12190 if i=2 then do;

12191 if j<=3 then con[i,j]=boundlam;

12192 if j>=4 & j<=6 then con[i,j]=varianupper;

12193 if j=7 then con[i,j]=boundphi;

12194 if j=8 then con[i,j]=boundint;

12195 if j>=9 then con[i,j]=boundbeta;

12196 end;

12197 end;

12198 end;

12199

12200 con[3,1]=1;

12200! con[3,12]=0;

12200! con[3,13]=2;

12201

12202 print con;

NOTE: Writing HTML Body file: sashtml147.htm

12203

12204 ********initial value;

12205 *x={2 0.04954 0.07649 270.560 118.793 151.389 0.0284 108.702 12.560 1 -0.5};

12206 x={2 0.0584913446 0.0810734988 270.5599998 118.79300677 151.38885874 0.0378699563

12206! 108.70197893 12.56070818 0.9993222093 -0.500002042};

12207 xxx=F_GLOBAL(x);

12208 print xxx;

12209 ********threshold parameters;

12210 opt=J(1,9,.);

12211 opt[1]=1;

12212 opt[2]=5;

12213 tc=J(1,11,.);

12214 ********maximization;

12215 call nlpnra(rc,xres,"F_GLOBAL",x,opt,con,tc);

ERROR: (execution) Invalid argument to function.

operation : LOG at line 12161 column 16

operands : temp1

temp1 1 row 1 col (numeric)

0

statement : ASSIGN at line 12161 column 5

traceback : module F_GLOBAL at line 12161 column 5

ERROR: Execution error as noted previously. (rc=100)

operation : NLPNRA at line 12215 column 1

operands : *LIT1116, x, opt, con, tc

*LIT1116 1 row 1 col (character, size 😎

F_GLOBAL

x 1 row 11 cols (numeric)

opt 1 row 9 cols (numeric)

1 5 . . . . . . .

con 3 rows 13 cols (numeric)

tc 1 row 11 cols (numeric)

statement : CALL at line 12215 column 1

12216 create solution from xres;

ERROR: Matrix xres has not been set to a value.

statement : CREATE at line 12216 column 1

12217 append from xres;

ERROR: No data set is currently open for output.

statement : APPEND at line 12217 column 1

12218 quit;

NOTE: Exiting IML.

NOTE: 190554 workspace compresses.

NOTE: The SAS System stopped processing this step because of errors.

NOTE: PROCEDURE IML used (Total process time):

real time 1:37.33

cpu time 1:35.27

Rick_SAS · Posted 09-26-2011 01:43 PM

I think you have to print temp1 BEFORE you call the LOG function. You are only seeing the result of previous iterations.

I think I might see what is happening. Your loop is essentially as follows:

ff=0;

do l=1 to nrow(uu);

surv= ...;

if surv<si then fyd=0;

ff=ff+fyd;

end;

I think si is large compared to surv so that surv < si is always true.

Of course, the real question is what to do about it, and hopefully you can figure that out. If the value of Duration_Month (which determines si) is bad, you can fix it. Otherwise, it could be that lambda1t and lambda1 are getting small during the optimization. This would case surv to get large. If this is the problem, you might want to make sure they remain bounded away from zero.

By the way, if you haven't read this paper by Charlie Hallahan on MLE estimation, I recommend it: http://www2.sas.com/proceedings/sugi22/STATS/PAPER293.PDF

JerryLee · Posted 10-04-2011 04:01 PM

Hi Rick,

Thank you so much! This is really helpful. I got the problem resolved. I was just wondering if IML or "call NLPNRA" also offer the calculation of Fisher information matrix which is a very natural thing associated with maximum likelihood estiamtors (MLEs). Thanks again.

Rick_SAS · Posted 10-05-2011 09:27 AM

Glad you resolved the problem.

I am not very familiar with the Fisher information matrix, so I'll let someone else offer advice on that. As far as I know, there are no "built in" routines for computing it, since the NLP functions are designed for general optimization, not for the specific case of MLE.

JerryLee · Posted 11-08-2011 04:57 PM

Hi Rick,

Do you know if IML has a built in function or command to calculate the second order derivatives of a function? Thanks a lot.

Rick_SAS · Posted 11-09-2011 09:55 AM

The NLPFDD function uses computes numerical first and second derivatives:

http://blogs.sas.com/content/iml/2011/10/14/hints-for-derivatives/

SAS Innovate 2025: Save the Date