Prev: Missing values
Next: Exporting the results
From: Eli Y. Kling on 23 Dec 2009 07:25
The (non-)convergence of NLMIXED (& NLP) is an issue that haunted me
in the last two projects I worked on.

I am reading a very interesting paper by Serroyen et al (Nonlinear
Models for Longitudinal Data, The American Statistician, Nov 09, Vol
63 #4, pp378-388.

They use NLMIXED to module Orange tree circumference. They said on
page 380: “To improve numerical stability, and hence convergence,
while fitting the model, both the response values, trunk
circumference, and the covariate values, age, were divided by 100.”

I tired it out:

data trees;
input day tree c;
ct=c/100;
dayt=day/100;
datalines;
118 1 30
484 1 58
664 1 87
1004 1 115
1231 1 120
1372 1 142
1582 1 145
118 2 33
484 2 39
664 2 111
1004 2 156
1231 2 172
1372 2 203
1582 2 203
118 3 30
484 3 51
664 3 75
1004 3 108
1231 3 115
1372 3 139
1582 3 140
118 4 32
484 4 62
664 4 112
1004 4 167
1231 4 179
1372 4 209
1582 4 214
118 5 30
484 5 49
664 5 81
1004 5 125
1231 5 142
1372 5 174
1582 5 177
;
run;

*The variable names here are different to the ones used in the
appendix of the paper as I wrote the nlmixed before I looked at their
implementation;

*Using original values C & Day;
proc nlmixed data=trees qpoints=20;
e=exp(-(Day-b2)/b3);
m=(b1+U)/(1+e);
model C ~ normal(m,s2);
random U ~ normal(0,s1) subject=tree out=nlm_U;
run;quit;

*Using values/100 Ct & Dayt;
proc nlmixed data=trees qpoints=20;
e=exp(-(Dayt-b2)/b3);
m=(b1+U)/(1+e);
model Ct ~ normal(m,s2);
random U ~ normal(0,s1) subject=tree out=nlm_U;
run;quit;


When using the original values I get :
WARNING: The final Hessian matrix is not positive definite, and
therefore the estimated covariance matrix is not full rank and may be
unreliable. The variance of some parameter estimates is zero or some
parameters are linearly related to other parameters.


When using the values divided by 100 the NLMIXED converged and
reported values similar (but not spot on) to the ones listed in table
2 of the paper.

Interesting behaviour. I do not think the division by 100 made any
difference to the positiveness of the Hessian. Could anybody shed
light on this?

With regard,

Eli
From: Dale McLerran on 28 Dec 2009 18:18
--- On Wed, 12/23/09, Eli Y. Kling <eli.kling(a)GMAIL.COM> wrote:

> From: Eli Y. Kling <eli.kling(a)GMAIL.COM>
> Subject: NLMIXED numerical stability
> To: SAS-L(a)LISTSERV.UGA.EDU
> Date: Wednesday, December 23, 2009, 4:25 AM
> The (non-)convergence of NLMIXED
> (& NLP) is an issue that haunted me
> in the last two projects I worked on.
>
> I am reading a very interesting paper by Serroyen et al
> (Nonlinear
> Models for Longitudinal Data, The American Statistician,
> Nov 09, Vol
> 63 #4, pp378-388.
>
> They use NLMIXED to module Orange tree circumference. They
> said on
> page 380: “To improve numerical stability, and hence
> convergence,
> while fitting the model, both the response values, trunk
> circumference, and the covariate values, age, were divided
> by 100.”
>
> I tired it out:
>
> data trees;
> input day tree c;
> ct=c/100;
> dayt=day/100;
> datalines;
> � 118 1� 30
> � 484 1� 58
> � 664 1� 87
> 1004 1 115
> 1231 1 120
> 1372 1 142
> 1582 1 145
> � 118 2� 33
> � 484 2� 39
> � 664 2 111
> 1004 2 156
> 1231 2 172
> 1372 2 203
> 1582 2 203
> � 118 3� 30
> � 484 3� 51
> � 664 3� 75
> 1004 3 108
> 1231 3 115
> 1372 3 139
> 1582 3 140
> � 118 4� 32
> � 484 4� 62
> � 664 4 112
> 1004 4 167
> 1231 4 179
> 1372 4 209
> 1582 4 214
> � 118 5� 30
> � 484 5� 49
> � 664 5� 81
> 1004 5 125
> 1231 5 142
> 1372 5 174
> 1582 5 177
> ;
> run;
>
> *The variable names here are different to the ones used in
> the
> appendix of the paper as I wrote the nlmixed before I
> looked at their
> implementation;
>
> *Using original values C & Day;
> proc nlmixed data=trees qpoints=20;
> e=exp(-(Day-b2)/b3);
> m=(b1+U)/(1+e);
> model C ~ normal(m,s2);
> random U ~ normal(0,s1) subject=tree out=nlm_U;
> run;quit;
>
> *Using values/100 Ct & Dayt;
> proc nlmixed data=trees qpoints=20;
> e=exp(-(Dayt-b2)/b3);
> m=(b1+U)/(1+e);
> model Ct ~ normal(m,s2);
> random U ~ normal(0,s1) subject=tree out=nlm_U;
> run;quit;
>
>
> When using the original values I get :
> WARNING: The final Hessian matrix is not positive definite,
> and
> therefore the estimated covariance matrix is not full rank
> and may be
> unreliable.� The variance of some parameter estimates
> is zero or some
> parameters are linearly related to other parameters.
>
>
> When using the values divided by 100 the NLMIXED converged
> and
> reported values similar (but not spot on) to the ones
> listed in table
> 2 of the paper.
>
> Interesting behaviour. I do not think the division by 100
> made any
> difference to the positiveness of the Hessian. Could
> anybody shed
> light on this?
>
> With regard,
>
> Eli
>

Eli,

The result can be readily explained. Note that the initial
parameter estimates for b1, b2, and b3 are all 1. Now, if
we evaluate exp(-(Day-b2)/b3) at b2=b3=1 for the observed Day
values, we obtain

Day Day-1 exp(-(Day-1))
118 117 1.5400882849875E-51
484 483 1.7209380649407E-210
664 663 1.155469531661E-288
1004 1003 0.0
1231 1230 0.0
1372 1371 0.0
1582 1581 0.0


I should state that these are the values which are obtained for
SAS running on a 32-bit OS (Win-XP). Results would certainly
be different on a 64-bit OS. Notice that the value of
exp(-(Day-1)) is exactly 0.0 for the four values with Day>1000.
Thus, you have lost the ability to differentiate between the
response for those four days. It should also be noted that for
Day=664, the value is very near the machine limits for the
smallest positive number. Again, running on 32-bit Win-XP,
I find that the smallest positive number which SAS can
represent is approximately 1E-321. While you can represent
the value of the function at Day=664 (at the initial
parameter values), you may not be able to properly represent
values of derivatives of the function with respect to the
parameters of the model.

By rescaling the data, the problems with exponents of extremely
large negative values is alleviated. We obtain different
values of exp(-(Dayt-1)) for all seven values of Dayt. This
allows the model to compute all functionals without loss of
information. This results in numerical stability.

Dale

---------------------------------------
Dale McLerran
Fred Hutchinson Cancer Research Center
mailto: dmclerra(a)NO_SPAMfhcrc.org
Ph: (206) 667-2926
Fax: (206) 667-5977
---------------------------------------
From: Eli Y. Kling on 29 Dec 2009 06:09
Very interesting.

Thanks for that.

Happy New Year,

Eli
 | 
Pages: 1
Prev: Missing values
Next: Exporting the results