Date: Mon, 11 Apr 2011 10:43:05 -0700
Reply-To: Dale McLerran <stringplayer_2@YAHOO.COM>
Sender: "SAS(r) Discussion" <SAS-L@LISTSERV.UGA.EDU>
From: Dale McLerran <stringplayer_2@YAHOO.COM>
Subject: Re: NLMIXED versus PROC MIXED
In-Reply-To: <1d3b88ebb42b8ead09bee6909eef5e5c.squirrel@webmail.pitt.edu>
Content-Type: text/plain; charset=iso-8859-1
Francis,
When I encounter problems estimating a particular model, I always
try to build up the final model from simple to complex. Rather
than going directly for the joint model, I would suggest estimating
the parameters for models 1 and 2 separately to make sure that
the code for the separate models converges OK. Thus, I would
estimate the following models:
proc mixed data=mysample(where=(model=1)) covtest method=ml CL;
class stnum model;
model Y = day3apache newsex newrace charlson daidx /s chisq;
random intercept daidx /type=UN subject=stnum g v;
ods output covparms=cp_m1;
ods output parameterestimates=parms_m1;
run;
proc mixed data=mysample(where=(model=2)) covtest method=ml CL;
class stnum model;
model Y = day3apache newsex newrace charlson daidx /s chisq;
random intercept daidx /type=UN subject=stnum g v;
ods output covparms=cp_m2;
ods output parameterestimates=parms_m2;
run;
proc nlmixed data=mysample(where=(model=1))
tech=TRUREG qpoints=1 absgconv=1E-5 optcheck
ITDETAILS ;
parms
/*Long Model 1*/
bL01=3.27 bL11=0.02 bL21=-0.03 bL31= 0.6 bL41= -0.17 bL51= -0.33
se1=1.98
/*covariances*/
s1=3.15 s2=0.72 s3=0.74;
bounds s1 s3 se1 >= 0 ;
/*longitudinal Tobit Model part*/
eta = bL01 +
bL11*day3apache +
bL21*newsex*(newsex=1) +
bL31*newrace*(newrace=1) +
bL41*charlson +
bL51*daidx +
a1 +
b1*daidx; /*linear predictor for Y1*/
se=se1;
pi=CONSTANT('PI');
lllong=log(1/sqrt(2*pi*se))*exp(-(Y-eta)**2/(2*se));
/*non censored values pick up this piece*/;
model last ~ general(lllong);
random a1 b1 ~ normal([0,0],
[s1,
s2, s3])
subject=stnum;
run;
proc nlmixed data=mysample(where=(model=2))
tech=TRUREG qpoints=1 absgconv=1E-5 optcheck
ITDETAILS ;
parms
/*long Model 2*/
bL02=1.84 bL12=0.007 bL22=0.004 bL32=-0.02 bL42= 0.09 bL52= -0.01
se2=0.44
/*covariances*/
s6=0.11 s9=0.04 s10=0.02;
bounds s6 s10 se2 >= 0 ;
eta = bL02 +
bL12*day3apache +
bL22*newsex*(newsex=1) +
bL32*charlson +
bL42*newrace*(newrace=1) +
bL52*daidx +
a2 +
b2*daidx; /*linear predictor for Y2*/
se=se2;
pi=CONSTANT('PI');
lllong=log(1/sqrt(2*pi*se))*exp(-(Y-eta)**2/(2*se));
/*non censored values pick up this piece*/;
model last ~ general(lllong);
random a2 b2 ~ normal([0,0],
[s6,
s9, s10])
subject=stnum;
run;
Do these two NLMIXED models converge with parameter estimates
which are similar to the parameter estimates from the two
invocations of PROC MIXED? Note that the parameter estimates
are not expected to be identical between the two procedures
because of differences in estimation methods. But for a
model with normal residuals and normal random effects, the
parameter estimates should be similar.
In looking at your code, I am curious about the reference to
a variable LAST on your NLMIXED MODEL statement. This variable
does not appear anywhere in your PROC MIXED code. What is the
purpose of this variable? Is it nonmissing for all observations
which are employed by your MIXED code? Could this variable be
the source of some of your problems?
Dale
---------------------------------------
Dale McLerran
Fred Hutchinson Cancer Research Center
mailto: dmclerra@nFoHsCpRaCm.org
Ph: (206) 667-2926
Fax: (206) 667-5977
---------------------------------------
--- On Mon, 4/11/11, SUBSCRIBE SAS-L Ameribrit <FRP3@PITT.EDU> wrote:
> From: SUBSCRIBE SAS-L Ameribrit <FRP3@PITT.EDU>
> Subject: NLMIXED versus PROC MIXED
> To: SAS-L@LISTSERV.UGA.EDU
> Date: Monday, April 11, 2011, 9:10 AM
> Hi Guys I am trying to reproduce the
> results from Proc Mixed in NLMIXED.
>
> Here is the Proc Mixed code for a bivariate Model (ignoring
> censoring)
>
>
> proc mixed data=mysample covtest method=ml CL;
> class stnum model;
>
> model Y = model day3apache newsex newrace charlson daidx
> model*day3apache
> model*newsex model*charlson model*newrace model*daidx /s
> chisq;
> random model model*daidx /type=UN subject=stnum g v;
> repeated /type=VC grp=MODEL subject=stnum;
> ods output covparms=cp;
> ods output parameterestimates;
> run;
>
> /*here is the NLMIXED code I am trying to use to reproduce
> the results*/
>
> proc nlmixed data=mysample tech=TRUREG qpoints=1
> absgconv=1E-5 optcheck
> ITDETAILS ;
>
> parms
> /*Long Model 1*/
> bL01=3.27 bL11=0.02 bL21=-0.03 bL31= 0.6
> bL41= -0.17 bL51= -0.33
> se1=1.98
> /*long Model 2*/
> bL02=1.84 bL12=0.007 bL22=0.004 bL32=
> -0.02 bL42= 0.09 bL52= -0.01
> se2=0.44
> /*covariances*/
> s1=3.15 s2=0.72 s3=0.74 s4=-0.50 s5=-0.14 s6=0.11 s7=-0.15
> s8=-0.12
> s9=0.04 s10=0.02;
>
> bounds s1 s3 s6 s10 se1 se2 >= 0 ;
>
> /*longitudinal Tobit Model part*/
> if Model=1 then do;
> eta=bL01+bL11*day3apache+bL21*newsex*(newsex=1)+bL31*newrace*(newrace=1)+bL41*charlson+bL51*daidx+a1+b1*daidx;
> /*linear predictor for Y1*/
> se=se1;
> end;
>
> if Model=2 then do;
> eta=bL02+bL12*day3apache+bL22*newsex*(newsex=1)+bL32*charlson+bL42*newrace*(newrace=1)+bL52*daidx+a2+b2*daidx;
> /*linear predictor for Y2*/
> se=se2;
> end;
>
> pi=CONSTANT('PI');
>
> lllong=log(1/sqrt(2*pi*se))*exp(-(Y-eta)**2/(2*se));
> /*non censored
> values pick up this piece*/;
>
> model last ~ general(lllong);
> random a1 b1 a2 b2 ~ normal([0,0,0,0], [s1, s2,
> s3,s4,s5,s6,s7,s8,s9,s10])
> subject=stnum;
>
> run;
>
> This just gives no valid parameter points found. The
> estimates are from
> the proc mixed run which is almost instantaneous. If i
> remove all the
> parms it just runs forever and eventually says it needs
> more iterations.
>
> Any ideas about what I am doing wrong would be
> appreciated.
> Francis
>
