cancel
Turn on suggestions
Auto-suggest helps you quickly narrow down your search results by suggesting possible matches as you type.
Showing results for 
Search instead for 
Did you mean: 
BookmarkSubscribeRSS Feed
JerryLee
Calcite | Level 5 JerryLee
Calcite | Level 5

call NLPNRA and NLPFDD in sas 9.2 and 9.3

Posted 06-04-2012 10:47 AM (873 views)

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;

  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;

     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;

     end;

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

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

   lyd=lyd+log(fyd);

   end;

  end;

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;

end;

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;

enable2fact.txt.zip
0 Likes
Reply
1 REPLY 1
Rick_SAS
SAS Super FREQ Rick_SAS
SAS Super FREQ

Re: call NLPNRA and NLPFDD in sas 9.2 and 9.3

Posted 06-12-2012 10:28 AM (759 views)  |   In reply to JerryLee

I suggest improving the efficiency of your code by using the ideas in the blog post "Eight tips to make your simulation run faster." All those loops over "i"

I'll even add a few more tips:
Tip #9: Don't read the data in the F_GLOBAL function. This is KILLING your performance because you are reading the data at EVERY STEP of the optimization!  Instead, read the data once and then include those vectors as global variables by using the GLOBAL statement

Tip #10: det(inv(X) = 1/det(X) so there is no need to take the inverse

0 Likes
Reply

sas-innovate-2024.png

Don't miss out on SAS Innovate - Register now for the FREE Livestream!

Can't make it to Vegas? No problem! Watch our general sessions LIVE or on-demand starting April 17th. Hear from SAS execs, best-selling author Adam Grant, Hot Ones host Sean Evans, top tech journalist Kara Swisher, AI expert Cassie Kozyrkov, and the mind-blowing dance crew iLuminate! Plus, get access to over 20 breakout sessions.

 

Register now!

Multiple Linear Regression in SAS

Learn how to run multiple linear regression models with and without interactions, presented by SAS user Alex Chaplin.

Watch Multiple Linear Regression in SAS on YouTube

Find more tutorials on the SAS Users YouTube channel.

From The DO Loop
Want more? Visit our blog for more articles like these.
Discussion stats
  • 1 reply
  • ‎06-04-2012 10:47 AM
  • 874 views
  • 0 likes
  • 2 in conversation