Date: Mon, 15 Dec 2003 12:42:06 -0500
Reply-To: Susan Stewart <sstewart@NBER.ORG>
Sender: "SAS(r) Discussion" <SAS-L@LISTSERV.UGA.EDU>
From: Susan Stewart <sstewart@NBER.ORG>
Subject: NLMIXED model unstable?
In-Reply-To: <20031006224203.39602.qmail@web21107.mail.yahoo.com>
Content-Type: TEXT/PLAIN; charset=US-ASCII
We have a complicated ordered probit model that we are running using the
PROC NLMIXED code that you generously provided to us in October. However,
we are concerned about the stability of the model.
Our model includes 21 interaction terms between sets of variables, for
example we estimate a single b for the interaction between 6 variables
representing mental health and 16 variables representing physical function
(these include both main effects and interaction terms). This is something
we came up with to avoid having to test hundreds of smaller interactions.
Do you think it is OK to do this?
Here is an excerpt of that code (the full model is farther below):
b178* ((b4*sfcn1 + b5*pac1n1 + b6*pac41 + b7*pac51 + b8*pac61 + b9*pac14 +
b10*pac15 + b11*pac16 + b12*drvbus + b14*pac45 + b15*pac46 + b16*pac56 +
b18*dbpac11 + b19*dbpac41 + b20*dbpac51 + b21*dbpac61) * (b40*depres1
+b41*anxious1 +b42*sleep1 +b43*depanx1 +b44*depslp1 +b45*anxsl1))
(All of these are 0/1 variables. Our sample size is 1420, so we should
have plenty of DF. The sample sizes for all terms are reasonable; the
lowest are above 75 and most are much higher.)
The version of our model that we originally wanted includes 58 main
effects and interactions, and 21 interaction terms between sets of
variables. It runs fairly well, with reasonably low standard errors (above
1.0 for some terms but not above 3.0). However, if one or more other such
terms are added or removed, the model often crashes, with very high SEs
for most terms (higher than 400 for some).
As another way of testing the stability of the original model, we tried
specifying different start values, for example half to 4 and half to -4.
The results are almost the same when we do this, as long as I specify a
starting value for the intercept as well (for example 6). If I don't
specify a starting value for the intercept, the model runs terribly, with
no SE calculated for most of the terms involving pain.
When I run it with starting values of 10 and -10, it seems to
automatically switch to a different convergence criterion (FCONV rather
than GCONV) and give crazy output, with SEs of zero for about half of the
terms, only 150 iterations, and a -2LL of 66590 (The original model with
GCONV did 244 iterations and had -2LL of 3073).
The model also yields slightly different results depending on how the
data is sorted, but hopefully this is just rounding error.
We realize that this is a rather large question but appreciate any bits of
advice you might have on the validity and stability of our model, possible
explanations for our strange output, and other ways of testing whether or
not we have a good model.
Thank you,
Susan
Here is our code. Let me know if sample output or any other
information would also be helpful.
proc sort; by record;
proc nlmixed maxiter = 1000;
parms
a1=1 a2=1 a3=1
b0 = 6
b200=4 b201=-4 b202=4 b203=-4 b204=4 b205=-4 b206=4 b207=-4 b208=4 b209=-4
b210=4 b211=-4 b212=4 b213=-4 b214=4
b1=4 b2=-4 b3=4 b4=-4 b5=4 b6=-4 b7=4 b8=-4 b9=4 b10=-4 b11=4 b12=-4 b14=4
b15=-4 b16=4 b18=-4 b19=4 b20=-4 b21=4 b22=-4 b23=4 b24=-4 b25=4 b26=-4
b35=4 b36=-4 b37=4 b40=-4 b41=4 b42=-4 b43=4 b44=-4 b45=4 b50=-4 b60=4
b61=-4 b62=4 b63=-4 b64=4 b65=-4 b70=4 b71=-4 b72=4
b170=-4 b171=4 b172=-4 b173=4 b174=-4 b175=4 b176=-4 b177=4 b178=-4 b179=4
b180=-4 b181=4 b182=-4 b183=4 b184=-4 b185=4 b186=-4 b187=4 b188=-4 b189=4
b190=-4
/*
b0=10
b200=10 b201=-10 b202=10 b203=-10 b204=10 b205=-10 b206=10 b207=-10
b208=10 b209=-10
b210=10 b211=-10 b212=10 b213=-10 b214=10
b1=10 b2=-10 b3=10 b4=-10 b5=10 b6=-10 b7=10 b8=-10 b9=10 b10=-10 b11=10
b12=-10 b14=10
b15=-10 b16=10 b18=-10 b19=10 b20=-10 b21=10 b22=-10 b23=10 b24=-10 b25=10
b26=-10
b35=10 b36=-10 b37=10 b40=-10 b41=10 b42=-10 b43=10 b44=-10 b45=10 b50=-10
b60=10
b61=-10 b62=10 b63=-10 b64=10 b65=-10 b70=10 b71=-10 b72=10
b170=-10 b171=10 b172=-10 b173=10 b174=-10 b175=10 b176=-10 b177=10
b178=-10 b179=10
b180=-10 b181=10 b182=-10 b183=10 b184=-10 b185=10 b186=-10 b187=10
b188=-10 b189=10
b190=-10
*/
;
array eta {4};
eta1 =
-(b0 + b200*age50_54 + b201*age55_59 + b202*age60_64 + b203*age65_69 +
b204*age70_74 + b205*age75_79 + b206*age80up + b207*sex + b208*sex5054 +
b209*sex5559 + b210*sex6064 + b211*sex6569 + b212*sex7074 + b213*sex7579 +
b214*sex80up
+ b1*sac0or11 + b2*sac21 + b3*s0or1_21
+ b4*sfcn1 + b5*pac1n1 + b6*pac41 + b7*pac51 + b8*pac61 + b9*pac14 +
b10*pac15 + b11*pac16 + b12*drvbus + b14*pac45 + b15*pac46 +
b16*pac56 + b18*dbpac11 + b19*dbpac41 + b20*dbpac51 + b21*dbpac61
+b22*pain1 +b23*rash1 +b24*ubsex1 +b25*painus1 +b26*painrh1
+b35*sick1 +b36*cough1 +b37*sickcgh1
+b40*depres1 +b41*anxious1 +b42*sleep1 +b43*depanx1 +b44*depslp1
+b45*anxsl1
+b50*tired1
+b60*sexprob1 +b61*ent1 +b62*appear1 +b63*meddiet1 +b64*limb1 +b65*hedach1
+b70*eyepain1 +b71*glasses1 +b72*talk1
+ b170*((b1*sac0or11 + b2*sac21 + b3*s0or1_21)* (b4*sfcn1 + b5*pac1n1 +
b6*pac41 + b7*pac51 + b8*pac61 + b9*pac14 + b10*pac15 + b11*pac16 +
b12*drvbus + b14*pac45 + b15*pac46 + b16*pac56
+ b18*dbpac11 + b19*dbpac41 + b20*dbpac51 + b21*dbpac61))
+ b171*((b1*sac0or11 + b2*sac21 + b3*s0or1_21)*(b22*pain1
+b23*rash1 +b24*ubsex1 +b25*painus1 +b26*painrh1))
+ b172*((b1*sac0or11 + b2*sac21 + b3*s0or1_21)*
(b35*sick1 +b36*cough1 +b37*sickcgh1))
+ b173*((b1*sac0or11 + b2*sac21 + b3*s0or1_21)* (b40*depres1
+b41*anxious1 +b42*sleep1 +b43*depanx1 +b44*depslp1 +b45*anxsl1))
+ b174*((b1*sac0or11 + b2*sac21 + b3*s0or1_21)* tired1)
+ b175*((b1*sac0or11 + b2*sac21 + b3*s0or1_21)*
(b70*eyepain1 +b71*glasses1 +b72*talk1))
+ b176* ((b4*sfcn1 + b5*pac1n1 + b6*pac41 + b7*pac51 + b8*pac61 + b9*pac14
+ b10*pac15 + b11*pac16 + b12*drvbus + b14*pac45 + b15*pac46
+ b16*pac56 + b18*dbpac11 + b19*dbpac41 + b20*dbpac51 +
b21*dbpac61) * (b22*pain1 +b23*rash1 +b24*ubsex1 +b25*painus1
+b26*painrh1))
+ b177* ((b4*sfcn1 + b5*pac1n1 + b6*pac41 + b7*pac51 + b8*pac61 + b9*pac14
+ b10*pac15 + b11*pac16 + b12*drvbus + b14*pac45 + b15*pac46
+ b16*pac56 + b18*dbpac11 + b19*dbpac41 + b20*dbpac51 +
b21*dbpac61) *
(b35*sick1 +b36*cough1 +b37*sickcgh1))
+ b178* ((b4*sfcn1 + b5*pac1n1 + b6*pac41 + b7*pac51 + b8*pac61 + b9*pac14
+ b10*pac15 + b11*pac16 + b12*drvbus + b14*pac45 + b15*pac46
+ b16*pac56 + b18*dbpac11 + b19*dbpac41 + b20*dbpac51 +
b21*dbpac61) *
(b40*depres1 +b41*anxious1 +b42*sleep1 +b43*depanx1 +b44*depslp1
+b45*anxsl1))
+ b179* ((b4*sfcn1 + b5*pac1n1 + b6*pac41 + b7*pac51 + b8*pac61 + b9*pac14
+ b10*pac15 + b11*pac16 + b12*drvbus + b14*pac45 + b15*pac46
+ b16*pac56 + b18*dbpac11 + b19*dbpac41 + b20*dbpac51 +
b21*dbpac61) *
tired1)
+ b180* ((b4*sfcn1 + b5*pac1n1 + b6*pac41 + b7*pac51 + b8*pac61 + b9*pac14
+ b10*pac15 + b11*pac16 + b12*drvbus + b14*pac45 + b15*pac46
+ b16*pac56 + b18*dbpac11 + b19*dbpac41 + b20*dbpac51 +
b21*dbpac61) *
(b70*eyepain1 +b71*glasses1 +b72*talk1))
+b181* ((b22*pain1 +b23*rash1 +b24*ubsex1 +b25*painus1 +b26*painrh1)
* (b35*sick1 +b36*cough1 +b37*sickcgh1))
+b182* ((b22*pain1 +b23*rash1 +b24*ubsex1 +b25*painus1 +b26*painrh1)*
(b40*depres1 +b41*anxious1 +b42*sleep1 +b43*depanx1 +b44*depslp1
+b45*anxsl1))
+b183* ((b22*pain1 +b23*rash1 +b24*ubsex1 +b25*painus1 +b26*painrh1)*
tired1)
+b184* ((b22*pain1 +b23*rash1 +b24*ubsex1 +b25*painus1 +b26*painrh1)*
(b70*eyepain1 +b71*glasses1 +b72*talk1))
+b185 *((b35*sick1 +b36*cough1 +b37*sickcgh1)
* (b40*depres1 +b41*anxious1 +b42*sleep1 +b43*depanx1
+b44*depslp1 +b45*anxsl1))
+b186 *((b35*sick1 +b36*cough1 +b37*sickcgh1) * tired1)
+b187 *((b35*sick1 +b36*cough1 +b37*sickcgh1) *
(b70*eyepain1 +b71*glasses1 +b72*talk1))
+b188* ((b40*depres1 +b41*anxious1 +b42*sleep1 +b43*depanx1
+b44*depslp1 +b45*anxsl1) * tired1)
+b189* ((b40*depres1 +b41*anxious1 +b42*sleep1 +b43*depanx1
+b44*depslp1 +b45*anxsl1) *
(b70*eyepain1 +b71*glasses1 +b72*talk1))
+b190* ((b70*eyepain1 +b71*glasses1 +b72*talk1) * tired1));
eta2 = a1 + eta1;
eta3 = a2 + eta2;
eta4 = a3 + eta3;
if EVGFP=1 then do;
p_YLEy = cdf('Normal', eta1);
p_YLEym1 = 0;
end; else
if EVGFP=5 then do;
p_YLEy = 1;
p_YLEym1 = cdf('Normal', eta4);
end;
else do;
p_YLEy = cdf('Normal', eta{EVGFP} );
p_YLEym1 = cdf('Normal', eta{EVGFP-1});
end;
p_YEQy = max(p_YLEy - p_YLEym1, 1E-12);
loglike = log(p_YEQy);
model EVGFP ~ general(loglike);
predict eta1 OUT=Q2 DER; run;
On Mon, 6 Oct 2003, Dale McLerran wrote:
> Susan,
>
> I'm glad that you found the sample code useful. You might
> need to initialize a1, a2, and a3 to positive values. Without
> the initialization, SAS might be computing a zero value for
> the difference of the cdf for adjacent response levels. That
> is, we compute p[i]=cdf(eta[i])-cdf((eta[i-1]). If a[i-1]=0,
> then cdf(eta[i])=cdf(eta[i-1]) and p[i]=0. Since we then
> compute the log-likelihood contribution as log(p[i]), we would
> have computational problems if p[i]=0. Initializing the values
> a1, a2, and a3 so that they are positive would eliminate (at
> the first iteration) the possibility of computing log(0).
>
> We could also check to see if the value of p is zero before
> computing log(p). If the value of p were zero, then we could
> replace p by some (arbitrary) small constant value, say 1E-12.
> Then the logarithm could be obtained without error.
>
> Let me take this opportunity to present revised code which
> should be computationally more robust as well as computationally
> more efficient.
>
>
> proc nlmixed; /* implicit data=_last_ */
> where record ne 11212 and record ne 13294 and
> record ne 13277 and pain1 ne . and evgfp ne . and
> uibow1 ne . and depres1 ne . and anxious1 ne .;
>
> parms a1=1 a2=1 a3=1; /* Could initialize other params */
>
> /* Construct eta1, eta2, ... */
> array eta {4};
> eta1 = b0 + b1*pain1 + b2*uibow1 + b3*depres1 + b4*anxious1 +
> b5*((b1*pain1 + b2*uibow1)*(b3*depres1 + b4*anxious1));
> eta2 = a1 + eta1; /*a's define distance between categories*/
> eta3 = a2 + eta2;
> eta4 = a3 + eta3;
>
>
> /* Compute P(Y=y) = cdf(Y<=y) - cdf(Y<=y-1) */
> /* First construct cdf(Y<=y) and cdf(Y<=y-1) */
> if EVGFP=1 then do; /* Special handling required for Y=1 */
> p_YLEy = cdf('Normal', eta1);
> p_YLEym1 = 0;
> end; else
> if EVGFP=5 then do; /* Special handling required for Y=max */
> p_YLEy = 1;
> p_YLEym1 = cdf('Normal', eta4);
> end;
> else do;
> p_YLEy = cdf('Normal', eta{EVGFP} );
> p_YLEym1 = cdf('Normal', eta{EVGFP-1});
> end;
> /* When constructing p(Y=y), ensure p(Y=y)>0 */
> p_YEQy = max(p_YLEy - p_YLEym1, 1E-12);
>
>
> /* Construct log-likelihood contribution */
> loglike = log(p_YEQy);
>
>
> /* Fit model */
> model EVGFP ~ general(loglike);
> run;
>
>
> I believe that this code should work better. You might again
> want to try the code with your linear probit model just to
> verify that I did not introduce anything wrong. I would note
> that this code is totally untested. But it should work (better).
>
> Best wishes,
>
> Dale
>
>
> --- Susan Stewart <sstewart@nber.org> wrote:
> > Dear Dale,
> >
> > Thank you for your suggestion to use PROC NLMIXED, and the helpful
> > code
> > that you provided. I altered it to accommodate a 5-level response
> > variable, I was able to run a simple linear model and get the same
> > results
> > as I had gotten with the normal ordered probit (using proc logistic).
> > Now
> > I am trying to run the more complicated nonlinear example, and it is
> > getting stuck at certain observations. Do you have an idea why?
> >
> > Here is my SAS code as seen in the .log file. (To troubleshoot, I
> > tried
> > deleting certain observations that were causing execution errors, but
> > it
> > continued to get stuck at a new observation each time. I have not
> > identified anything that these observations have in common in terms
> > of the
> > variables in the equation.) I have tried sorting by the pain variable
> > before proc nlmixed, and the observation numbers that cause errors
> > are
> > higher (i.e. 1028)
> >
> > (My sample size is 1430 and it does not seem to be a memory problem
> > as I
> > ran it on unix with 128mb RAM.)
> >
> > thank you; I appreciate your advice.
> >
> > Susan
> >
> >
> > 463 proc nlmixed; where record ne 11212 and record ne 13294
> > and
> > 464 record ne
> > 465 13277 and pain1 ne . and
> > 466 evgfp ne . and uibow1 ne . and depres1 ne . and anxious1
> > ne .;
> > 467
> > 468 eta1 = b0 + b1*pain1 + b2*uibow1 + b3*depres1 +
> > b4*anxious1 +
> > 469 b5*((b1*pain1 + b2*uibow1)*(b3*depres1 + b4*anxious1));
> > 470
> > 471 eta2 = a1 + eta1; /*a's define distance between
> > categories*/
> > 472 eta3 = a2 + eta2;
> > 473 eta4 = a3 + eta3;
> > 474
> > 475 p_yLE1 = cdf('Normal', eta1);
> > 476 p_yLE2 = cdf('Normal', eta2);
> > 477 p_yLE3 = cdf('Normal', eta3);
> > 478 p_yLE4 = cdf('Normal', eta4);
> > 479 p_yLE5 = 1;
> > 480
> > 481 p_yEQ1 = p_yLE1;
> > 482 p_yEQ2 = p_yLE2 - p_yLE1;
> > 483 p_yEQ3 = p_yLE3 - p_yLE2;
> > 484 p_yEQ4 = p_yLE4 - p_yLE3;
> > 485 p_yEQ5 = p_yLE5 - p_yLE4;
> > 486
> > 487 if EVGFP=1 then loglike=log(p_yEQ1); else
> > 488 if EVGFP=2 then loglike=log(p_yEQ2); else
> > 489 if EVGFP=3 then loglike=log(p_yEQ3); else
> > 490 if EVGFP=4 then loglike=log(p_yEQ4); else
> > 491 if EVGFP=5 then loglike=log(p_yEQ5);
> > 492
> > 493 model EVGFP ~ general(loglike);
> > 494 run;
> >
> > NOTE: A finite difference approximation is used for the derivative of
> > the
> > CDF
> > function at line 475 column 15.
> > NOTE: A finite difference approximation is used for the derivative of
> > the
> > CDF
> > function at line 476 column 15.
> > NOTE: A finite difference approximation is used for the derivative of
> > the
> > CDF
> > function at line 477 column 15.
> > NOTE: A finite difference approximation is used for the derivative of
> > the
> > CDF
> > function at line 478 column 15.
> > NOTE: Execution error for observation 391.
> > NOTE: The PROCEDURE NLMIXED printed page 41.
> >
> >
> -------------------------------------------------------------------------------
> > Output in .lst file:
> >
> > Data Set WORK.D0
> > Dependent Variable EVGFP
> > Distribution for Dependent Variable General
> > Optimization Technique Dual Quasi-Newton
> > Integration Method None
> >
> >
>
>
> =====
> ---------------------------------------
> Dale McLerran
> Fred Hutchinson Cancer Research Center
> mailto: dmclerra@fhcrc.org
> Ph: (206) 667-2926
> Fax: (206) 667-5977
> ---------------------------------------
>
> __________________________________
> Do you Yahoo!?
> The New Yahoo! Shopping - with improved product search
> http://shopping.yahoo.com
>
