Statistical programming, matrix languages, and more

call NLPNRA and NLPFDD in sas 9.2 and 9.3

Reply
Occasional Contributor
Posts: 10

call NLPNRA and NLPFDD in sas 9.2 and 9.3

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;

Attachment
SAS Super FREQ
Posts: 3,222

Re: call NLPNRA and NLPFDD in sas 9.2 and 9.3

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

Post a Question
Discussion Stats
  • 1 reply
  • 259 views
  • 0 likes
  • 2 in conversation