RE: posthoc step
From: jerry.nedelman@pharma.novartis.com
Subject: RE: [NMusers] posthoc step
Date: Sun, December 12, 2004 11:50 pm
Tom: SAS doesn't get past it either. That was with the Residual SD fixed.
When the Residual SD was estimated (starting always from 0.2), the estimate
of k was 10 for initial guesses from 10 to 3. But the estimated Residual SD
was weird then, essentially plus or minus infinity. When the initial guess
of k was 2 or 1, SAS failed to converge, although it stopped with k near 0.94.
Jerry
Results:
Initial k Final k SSE Res. SD
10 9.78 448.1 Fixed
7.25 9.78 448.1 Fixed
7.2 0.94 21.6 Fixed
Initial k Final k SSE Res. SD Initial Res. SD Final Res. SD
10 10 8.49E-14 Estimated 0.2 14526239
7 10 2.86E-19 Estimated 0.2 -7.92E+09
3 10 3.61E-20 Estimated 0.2 2.23E+10
2 0.9847 21.0218 Estimated 0.2 0.3671 Failure to converge
1 0.9578 21.7932 Estimated 0.2 0.1957 Failure to converge
Code:
data one;
******************************************;
* dat: an indicator for whether it's "real data" or an extra obs for the parameter.
* t: the t of y = exp(-kt)
* obs: observation (real data or theta)
* "Real data" generated by S-Plus: 10*exp(-c(1,2,3)) + 0.2*rnorm(3,0,1)
******************************************;
input dat t obs;
cards;
1 1 3.87
1 2 1.66
1 3 0.44
0 999 10
;
run;
**********************************;
title 'Weighted least squares, MAP, initial=10, fix residual sd';
**********************************;
proc nlin data=one;
y = dat*obs/0.2 + (1-dat)*obs/2;
model y = dat*( 10*exp(-k*t)/0.2 ) +
(1-dat)*( k/2 );
parm k=10;
run;
**********************************;
title 'Weighted least squares, MAP, initial=7.25, fix residual sd';
**********************************;
proc nlin data=one;
y = dat*obs/0.2 + (1-dat)*obs/2;
model y = dat*( 10*exp(-k*t)/0.2 ) +
(1-dat)*( k/2 );
parm k=7.25;
run;
**********************************;
title 'Weighted least squares, MAP, initial=7.2, fix residual sd';
**********************************;
proc nlin data=one;
y = dat*obs/0.2 + (1-dat)*obs/2;
model y = dat*( 10*exp(-k*t)/0.2 ) +
(1-dat)*( k/2 );
parm k=7.2;
run;
**********************************;
title 'Weighted least squares, MAP, initial=10, estimate residual sd';
**********************************;
proc nlin data=one;
y = dat*obs/res_sd + (1-dat)*obs/2;
model y = dat*( 10*exp(-k*t)/res_sd ) +
(1-dat)*( k/2 );
parm k=10 res_sd=0.2;
run;
**********************************;
title 'Weighted least squares, MAP, initial=7, estimate residual sd';
**********************************;
proc nlin data=one;
y = dat*obs/res_sd + (1-dat)*obs/2;
model y = dat*( 10*exp(-k*t)/res_sd ) +
(1-dat)*( k/2 );
parm k=7.25 res_sd=0.2;
run;
**********************************;
title 'Weighted least squares, MAP, initial=3, estimate residual variance';
**********************************;
proc nlin data=one;
y = dat*obs/res_sd + (1-dat)*obs/2;
model y = dat*( 10*exp(-k*t)/res_sd ) +
(1-dat)*( k/2 );
parm k=3 res_sd=0.2;
run;
**********************************;
title 'Weighted least squares, MAP, initial=2, estimate residual sd';
**********************************;
proc nlin data=one;
y = dat*obs/res_sd + (1-dat)*obs/2;
model y = dat*( 10*exp(-k*t)/res_sd ) +
(1-dat)*( k/2 );
parm k=2 res_sd=0.2;
run;
**********************************;
title 'Weighted least squares, MAP, initial=1, estimate residual sd';
**********************************;
proc nlin data=one;
y = dat*obs/res_sd + (1-dat)*obs/2;
model y = dat*( 10*exp(-k*t)/res_sd ) +
(1-dat)*( k/2 );
parm k=1 res_sd=0.2;
run;