call NLPNRA and NLPFDD in sas 9.2 and 9.3

JerryLee — Mon, 04 Jun 2012 14:47:43 GMT

I wrote a simulation program in IML to obtain Maximum likelihood estimator and its variance. There results from SAS 9.2 and SAS 9.3 are different for some samples. SAS 9.2 gave negative variance but SAS 9.3 gave positive variance. I need to decide which version I should use since I am running a massive simulation. I can run the simulation on a server (SAS 9.2) or on my own computer (SAS 9.3) which is very slow. The following is the code. A data set (enable2fact.txt) based on which 9.2 and 9.3 gave different results is attached. Thanks a lot for any help on this problem.

proc sort data=enable2fact;

by ID;

run;

proc means data=enable2fact;

by id;

var tij;

output out=no1;

run;

data no2;

set no1;

rename tij=N;

if _STAT_='N';

keep tij;

run;

proc iml;

start F_GLOBAL(x);

c=13;bc=6;p=2;

lambda1t=x[1];lambda2t=x[2];lambda1=x[3];lambda2=x[4];

phi1=x[5]*x[5];phi2=x[6]*x[6];phi3=x[7]*x[7];phi=x[8];

beta=x[,9:14];

use no2;

read all into n;

use enable2fact;

read all var {tij_b} into tijb;

read all var {tij} into tij;

read all var {trt} into trt;

read all var {factpal} into y;

read all var {intercept} into intcep;

read all var {duration_month} into s;

read all var {death} into D;

m=0;lyd=0;tune=0;

do i=1 to nrow(n);

m=m+n;

si=s;

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

rando=j(n,n,1);

epsi=I(n);

Ri=j(n,n,1);

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

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

ybeta=J(n,1,.);

betaA=J(n,1,.);

fyd1=J(n+1,1,.);

fyd2=J(n+1,1,.);

fyd3=J(n+1,1,.);

onej1=J(n+1,1,1);

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

do j=1 to n;

do k=j to n;

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

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

end;

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

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

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

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

ind1=J(n,1,.);

do l=1 to n;

if tib<=bc then ind1=1;else ind1=0;

end;

ind2=J(n,1,1)-ind1;

tind1=tib#ind1;

tind2=tib#ind2;

xib=incepti||trti||ind2||tind1||tind2;

beta1=J(5,1,1);

beta1[1]=beta[1];beta1[2]=beta[2];beta1[3]=bc*(beta[3]+trt*beta[5]-beta[4]-trt*beta[6]);

beta1[4]=beta[3]+trt*beta[5];beta1[5]=beta[4]+trt*beta[6];

if D=1 then do;

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

*************************************************************************log likelihood for observed events;

if trt=1 then do;

if si<=c then ld=log(lambda1t)-lambda1t*si;else ld=log(lambda2t)-c*lambda1t-lambda2t*(si-c);

end;

if trt=0 then do;

if si<=c then ld=log(lambda1)-lambda1*si;else ld=log(lambda2)-c*lambda1-lambda2*(si-c);

end;

lyd=lyd+ly+ld;

end;

*************************************************************************log likelihood for censored events;

if D=0 then do;

lam1=lambda1t;lam2=lambda2t;

if trt=0 then do; lam1=lambda1;lam2=lambda2;end;

if si>c then do;

***************************************************************;

do j=1 to n+1;

if j=1 then do;

ybeta=yi+J(n,1,1)*(-beta[1]-beta[2]*trt)+(beta[3]+beta[5]*trt)*ti;

betaA=J(n,1,1)*(beta[3]+beta[5]*trt);

mu22=(t(ybeta)*inv(Vi)*betaA-lam2)/(t(betaA)*inv(Vi)*betaA);

sigma2=1/(t(betaA)*inv(Vi)*betaA);

if si>ti[1]+bc then fyd1[1]=0;

if si<=ti[1]+bc then fyd1[1]=sqrt(det(inv(Vi)))

*exp(-tune+0.5*(t(ybeta)*inv(Vi)*betaA-lam2)**2/(t(betaA)*inv(Vi)*betaA)-0.5*t(ybeta)*inv(Vi)*ybeta)

*lam2*exp(-c*(lam1-lam2))*sqrt(sigma2)

*(sdf('NORMAL',(si-mu22)/sqrt(sigma2))-sdf('NORMAL',(ti[1]+bc-mu22)/sqrt(sigma2)));

end;

if j>1 & j<n+1 then do;

ybeta[1:j-1,]=yi[1:j-1,]+J(j-1,1,1)*(-beta[1]-beta[2]*trt-bc*(beta[3]+beta[5]*trt-beta[4]-beta[6]*trt))+ti[1:j-1,]

*(beta[4]+beta[6]*trt);

betaA[1:j-1,]=J(j-1,1,1)*(beta[4]+beta[6]*trt);

ybeta[j:n,]=yi[j:n,]+J(n-j+1,1,1)*(-beta[1]-beta[2]*trt)+ti[j:n,]*(beta[3]+beta[5]*trt);

betaA[j:n,]=J(n-j+1,1,1)*(beta[3]+beta[5]*trt);

mu22=(t(ybeta)*inv(Vi)*betaA-lam2)/(t(betaA)*inv(Vi)*betaA);

sigma2=1/(t(betaA)*inv(Vi)*betaA);

if si>ti+bc then fyd1=0;

if si<=ti[j-1]+bc then

fyd1=sqrt(det(inv(Vi)))*exp(-tune+0.5*(t(ybeta)*inv(Vi)*betaA-lam2)**2/(t(betaA)*inv(Vi)*betaA)-0.5*t(ybeta)*inv(Vi)*ybeta)

*lam2*exp(-c*(lam1-lam2))*sqrt(sigma2)

*(sdf('NORMAL',(ti[j-1]+bc-mu22)/sqrt(sigma2))-sdf('NORMAL',(ti+bc-mu22)/sqrt(sigma2)));

if si>ti[j-1]+bc & si<=ti+bc then do;

fyd1=sqrt(det(inv(Vi)))*exp(-tune+0.5*(t(ybeta)*inv(Vi)*betaA-lam2)**2/(t(betaA)*inv(Vi)*betaA)-0.5*t(ybeta)*inv(Vi)*ybeta)

*lam2*exp(-c*(lam1-lam2))*sqrt(sigma2)

*(sdf('NORMAL',(si-mu22)/sqrt(sigma2))-sdf('NORMAL',(ti+bc-mu22)/sqrt(sigma2)));end;

end;

if j=n+1 then do;

ybeta=yi+J(n,1,1)*(-beta[1]-beta[2]*trt-bc*(beta[3]+beta[5]*trt-beta[4]-beta[6]*trt))+(beta[4]+beta[6]*trt)*ti;

betaA=J(n,1,1)*(beta[4]+beta[6]*trt);

mu22=(t(ybeta)*inv(Vi)*betaA-lam2)/(t(betaA)*inv(Vi)*betaA);

sigma2=1/(t(betaA)*inv(Vi)*betaA);

if si>ti[n]+bc then fyd1=sqrt(det(inv(Vi)))

*exp(-tune+0.5*(t(ybeta)*inv(Vi)*betaA-lam2)**2/(t(betaA)*inv(Vi)*betaA)-0.5*t(ybeta)*inv(Vi)*ybeta

+logsdf('NORMAL',(si-mu22)/sqrt(sigma2)))

*lam2*exp(-c*(lam1-lam2))*sqrt(sigma2);

if si<=ti[n]+bc then fyd1=sqrt(det(inv(Vi)))

*exp(-tune+0.5*(t(ybeta)*inv(Vi)*betaA-lam2)**2/(t(betaA)*inv(Vi)*betaA)-0.5*t(ybeta)*inv(Vi)*ybeta

+logsdf('NORMAL',(ti[n]+bc-mu22)/sqrt(sigma2)))

*lam2*exp(-c*(lam1-lam2))*sqrt(sigma2);

end;

********************************************************;

fyd=t(onej1)*fyd1;

lyd=lyd+log(fyd);

end;

if si<c then do;

***************************************************************;

do j1=1 to n+1;

if j1=1 then do;

ybeta=yi+J(n,1,1)*(-beta[1]-beta[2]*trt)+(beta[3]+beta[5]*trt)*ti;

betaA=J(n,1,1)*(beta[3]+beta[5]*trt);

mu=(t(ybeta)*inv(Vi)*betaA-lam1)/(t(betaA)*inv(Vi)*betaA);

mu2=(t(ybeta)*inv(Vi)*betaA-lam2)/(t(betaA)*inv(Vi)*betaA);

sigma2=1/(t(betaA)*inv(Vi)*betaA);

if si>ti[1]+bc then fyd2[1]=0;

if si<=ti[1]+bc & c>ti[1]+bc then fyd2[1]=sqrt(det(inv(Vi)))

*exp(-tune+0.5*(t(ybeta)*inv(Vi)*betaA-lam1)**2/(t(betaA)*inv(Vi)*betaA)-0.5*t(ybeta)*inv(Vi)*ybeta)

*lam1*sqrt(sigma2)*(sdf('NORMAL',(si-mu)/sqrt(sigma2))-sdf('NORMAL',(ti[1]+bc-mu)/sqrt(sigma2)));

if c<=ti[1]+bc then fyd2[1]=sqrt(det(inv(Vi)))

*exp(-tune+0.5*(t(ybeta)*inv(Vi)*betaA-lam1)**2/(t(betaA)*inv(Vi)*betaA)-0.5*t(ybeta)*inv(Vi)*ybeta)

*lam1*sqrt(sigma2)*(sdf('NORMAL',(si-mu)/sqrt(sigma2))-sdf('NORMAL',(c-mu)/sqrt(sigma2)));

if c>ti[1]+bc then fyd3[1]=0;

if c<=ti[1]+bc then fyd3[1]=sqrt(det(inv(Vi)))

*exp(-tune+0.5*(t(ybeta)*inv(Vi)*betaA-lam2)**2/(t(betaA)*inv(Vi)*betaA)-0.5*t(ybeta)*inv(Vi)*ybeta)

*lam2*exp(-c*(lam1-lam2))*sqrt(sigma2)

*(sdf('NORMAL',(c-mu2)/sqrt(sigma2))-sdf('NORMAL',(ti[1]+bc-mu2)/sqrt(sigma2)));

end;

if j1>1 & j1<n+1 then do;

ybeta[1:j1-1,]=yi[1:j1-1,]+J(j1-1,1,1)*(-beta[1]-beta[2]*trt-bc*(beta[3]+beta[5]*trt-beta[4]-beta[6]*trt))+ti[1:j1-1,]

*(beta[4]+beta[6]*trt);

betaA[1:j1-1,]=J(j1-1,1,1)*(beta[4]+beta[6]*trt);

ybeta[j1:n,]=yi[j1:n,]+J(n-j1+1,1,1)*(-beta[1]-beta[2]*trt)+ti[j1:n,]*(beta[3]+beta[5]*trt);

betaA[j1:n,]=J(n-j1+1,1,1)*(beta[3]+beta[5]*trt);

mu=(t(ybeta)*inv(Vi)*betaA-lam1)/(t(betaA)*inv(Vi)*betaA);

mu2=(t(ybeta)*inv(Vi)*betaA-lam2)/(t(betaA)*inv(Vi)*betaA);

sigma2=1/(t(betaA)*inv(Vi)*betaA);

if si>ti[j1]+bc then fyd2[j1]=0;

if c<=ti[j1-1]+bc then fyd2[j1]=0;

if si<=ti[j1-1]+bc & c>ti[j1]+bc then fyd2[j1]=sqrt(det(inv(Vi)))

*exp(-tune+0.5*(t(ybeta)*inv(Vi)*betaA-lam1)**2/(t(betaA)*inv(Vi)*betaA)-0.5*t(ybeta)*inv(Vi)*ybeta)

*lam1*sqrt(sigma2)*(sdf('NORMAL',(ti[j1-1]+bc-mu)/sqrt(sigma2))-sdf('NORMAL',(ti[j1]+bc-mu)/sqrt(sigma2)));

if si>ti[j1-1]+bc & si<=ti[j1]+bc & c>ti[j1]+bc then fyd2[j1]=sqrt(det(inv(Vi)))

*exp(-tune+0.5*(t(ybeta)*inv(Vi)*betaA-lam1)**2/(t(betaA)*inv(Vi)*betaA)-0.5*t(ybeta)*inv(Vi)*ybeta)

*lam1*sqrt(sigma2)*(sdf('NORMAL',(si-mu)/sqrt(sigma2))-sdf('NORMAL',(ti[j1]+bc-mu)/sqrt(sigma2)));

if si<=ti[j1-1]+bc & c<=ti[j1]+bc & c>ti[j1-1]+bc then fyd2[j1]=sqrt(det(inv(Vi)))

*exp(-tune+0.5*(t(ybeta)*inv(Vi)*betaA-lam1)**2/(t(betaA)*inv(Vi)*betaA)-0.5*t(ybeta)*inv(Vi)*ybeta)

*lam1*sqrt(sigma2)*(sdf('NORMAL',(ti[j1-1]+bc-mu)/sqrt(sigma2))-sdf('NORMAL',(c-mu)/sqrt(sigma2)));

if si>ti[j1-1]+bc & c<=ti[j1]+bc then fyd2[j1]=sqrt(det(inv(Vi)))

*exp(-tune+0.5*(t(ybeta)*inv(Vi)*betaA-lam1)**2/(t(betaA)*inv(Vi)*betaA)-0.5*t(ybeta)*inv(Vi)*ybeta)

*lam1*sqrt(sigma2)*(sdf('NORMAL',(si-mu)/sqrt(sigma2))-sdf('NORMAL',(c-mu)/sqrt(sigma2)));

if c>ti[j1]+bc then fyd3[j1]=0;

if c<=ti[j1-1]+bc then fyd3[j1]=sqrt(det(inv(Vi)))

*exp(-tune+0.5*(t(ybeta)*inv(Vi)*betaA-lam2)**2/(t(betaA)*inv(Vi)*betaA)-0.5*t(ybeta)*inv(Vi)*ybeta)

*lam2*exp(-c*(lam1-lam2))*sqrt(sigma2)

*(sdf('NORMAL',(ti[j1-1]+bc-mu2)/sqrt(sigma2))-sdf('NORMAL',(ti[j1]+bc-mu2)/sqrt(sigma2)));

if c>ti[j1-1]+bc & c<=ti[j1]+bc

then fyd3[j1]=sqrt(det(inv(Vi)))

*exp(-tune+0.5*(t(ybeta)*inv(Vi)*betaA-lam2)**2/(t(betaA)*inv(Vi)*betaA)-0.5*t(ybeta)*inv(Vi)*ybeta)

*lam2*exp(-c*(lam1-lam2))*sqrt(sigma2)

*(sdf('NORMAL',(c-mu2)/sqrt(sigma2))-sdf('NORMAL',(ti[j1]+bc-mu2)/sqrt(sigma2)));

end;

if j1=n+1 then do;

ybeta=yi+J(n,1,1)*(-beta[1]-beta[2]*trt-bc*(beta[3]+beta[5]*trt-beta[4]-beta[6]*trt))+(beta[4]+beta[6]*trt)*ti;

betaA=J(n,1,1)*(beta[4]+beta[6]*trt);

mu=(t(ybeta)*inv(Vi)*betaA-lam1)/(t(betaA)*inv(Vi)*betaA);

mu2=(t(ybeta)*inv(Vi)*betaA-lam2)/(t(betaA)*inv(Vi)*betaA);

sigma2=1/(t(betaA)*inv(Vi)*betaA);

if c<=ti[n]+bc then fyd2[j1]=0;

if si>ti[n]+bc then fyd2[j1]=sqrt(det(inv(Vi)))

*exp(-tune+0.5*(t(ybeta)*inv(Vi)*betaA-lam1)**2/(t(betaA)*inv(Vi)*betaA)-0.5*t(ybeta)*inv(Vi)*ybeta

+logsdf('NORMAL',(si-mu)/sqrt(sigma2)))

*lam1*sqrt(sigma2)

-sqrt(det(inv(Vi)))

*exp(-tune+0.5*(t(ybeta)*inv(Vi)*betaA-lam1)**2/(t(betaA)*inv(Vi)*betaA)-0.5*t(ybeta)*inv(Vi)*ybeta

+logsdf('NORMAL',(c-mu)/sqrt(sigma2)))

*lam1*sqrt(sigma2);

if c>ti[n]+bc & si<=ti[n]+bc then fyd2[j1]=sqrt(det(inv(Vi)))

*exp(-tune+0.5*(t(ybeta)*inv(Vi)*betaA-lam1)**2/(t(betaA)*inv(Vi)*betaA)-0.5*t(ybeta)*inv(Vi)*ybeta

+logsdf('NORMAL',(ti[n]+bc-mu)/sqrt(sigma2)))

*lam1*sqrt(sigma2)

-sqrt(det(inv(Vi)))

*exp(-tune+0.5*(t(ybeta)*inv(Vi)*betaA-lam1)**2/(t(betaA)*inv(Vi)*betaA)-0.5*t(ybeta)*inv(Vi)*ybeta

+logsdf('NORMAL',(c-mu)/sqrt(sigma2)))

*lam1*sqrt(sigma2);

if c>ti[n]+bc then fyd3[j1]=sqrt(det(inv(Vi)))

*exp(-tune+0.5*(t(ybeta)*inv(Vi)*betaA-lam2)**2/(t(betaA)*inv(Vi)*betaA)-0.5*t(ybeta)*inv(Vi)*ybeta

+logsdf('NORMAL',(c-mu2)/sqrt(sigma2)))

*lam2*exp(-c*(lam1-lam2))*sqrt(sigma2);

if c<=ti[n]+bc then fyd3[j1]=sqrt(det(inv(Vi)))

*exp(-tune+0.5*(t(ybeta)*inv(Vi)*betaA-lam2)**2/(t(betaA)*inv(Vi)*betaA)-0.5*t(ybeta)*inv(Vi)*ybeta

+logsdf('NORMAL',(ti[n]+bc-mu2)/sqrt(sigma2)))

*lam2*exp(-c*(lam1-lam2))*sqrt(sigma2);

end;

********************************************************;

fyd=t(onej1)*(fyd2+fyd3);

lyd=lyd+log(fyd);

end;

f=lyd;

return(f);

finish F_GLOBAL;

x={0.05234 0.03309 0.07719 0.018472 18.22067 11.35843 9.95227 0.018541 108.4382 12.02538 3.994826 0.088183 -1.36527

-0.06034};

opt=J(1,9,.);

opt[1]=1;

opt[2]=3;

tc=j(1,11,.);

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

con=J(2,14,.);

boundbeta=500;

positive=10**(-10);

varianlower=3;

varianupper=500;

boundphi=5;

boundlam=1;

do i=1 to 2;

do j=1 to 14;

if i=1 then do;

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

if j>=5 & j<=7 then con[i,j]=varianlower;

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

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

end;

if i=2 then do;

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

if j>=5 & j<=7 then con[i,j]=varianupper;

if j=8 then con[i,j]=boundphi;

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

end;

print con;

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

CALL NLPFDD( f1, g1, h1, "F_GLOBAL", xr);

print f1,g1,h1;

B=I(14);

V=inv(-B*h1*(t(B)));

print V;

var=VECDIAG(V);

do s=1 to nrow(var);

if var~~<0 then var~~=0;~~~~

end;

std=t(sqrt(var));

para=xr[,1:ncol(xr)];

print para std, V;

truex={0.05234 0.03309 0.07719 0.018472 18.22067 11.35843 9.95227 0.018541 108.4382 12.02538 3.994826 0.088183 -1.36527

-0.06034};

CP=J(1,ncol(std),.);

do r=1 to ncol(std);

if truex<=para+1.96*std & truex>=para-1.96*std then CP=1;else CP=0;

end;

resultpara=para||std||CP;

create result from resultpara;

append from resultpara;

quit;

topic call NLPNRA and NLPFDD in sas 9.2 and 9.3 in SAS/IML Software and Matrix Computations

call NLPNRA and NLPFDD in sas 9.2 and 9.3

Re: call NLPNRA and NLPFDD in sas 9.2 and 9.3