[R] nlme
Jim Lindsey
jlindsey at alpha.luc.ac.be
Thu Jan 13 13:36:19 CET 2000
After a lot of struggle and considerable help from my Chinese
colleague (who is on the R list), I have managed to fit a few models
with nlme.
I have discovered that nlme is extremely sensitive to starting
values. I have never experienced this with any other nonlinear
optimization involving nlm.
One of the reasons for some of the weird error messages was that the
default value of stepmax for nlm is far too large.
I still have not been able to fit more than one covariate in fixed:
fixed=list(p1+p2~1,p3~sex)
works but
fixed=list(p1+p2~1,p3~sex+bilirubin)
does not, as pointed out in my previous email, giving the error:
Error in fixed[[nm]][[3]] != "1" : comparison (2) is possible only for vector types
How can one easily tell from the summary() output exactly how many
parameters have been estimated by nlme (including all variances)
except by dividing the AIC by 2 and comparing to the -logLik?
I do not understand the following:
Without a random effect,
> print(summary(resmls <- nls(logconcm~p1-p3+log(dose/(exp(p1)-exp(p2))*
+ (exp(-exp(p2)*tm)-exp(-exp(p1)*tm))),
+ control=nls.control(minFactor=0.00001),
+ data=datam,start=list(p1=0.5,p2=-4,p3=4),trace=F)))
Formula: logconcm ~ p1 - p3 + log(dose/(exp(p1) - exp(p2)) * (exp(-exp(p2) *
tm) - exp(-exp(p1) * tm)))
Parameters:
Estimate Std. Error t value
p1 0.52120 0.05981 8.715
p2 -3.88560 0.03004 -129.342
p3 4.00112 0.02319 172.507
Residual standard error: 0.4604 on 856 degrees of freedom
Correlation of Parameter Estimates:
p1 p2
p2 -0.3756
p3 0.5583 -0.6538
the volume is estimated to be
> exp(4.00112)
[1] 54.65933
a reasonable value. However, when volume is random,
> print(summary(resm2bf <- nlme(logconcm~p2-p3+log(dose/(exp(p1)-exp(p2))*
+ (exp(-exp(p2)*tm)-exp(-exp(p1)*tm))),start=list(fixed=c(-0.3,-5,0)),
+ fixed=list(p1+p2+p3~1),random=pdDiag(p3~1),weights=varExp(2),
+ groups=~subj,data=datam,method="ML")))
Nonlinear mixed-effects model fit by maximum likelihood
Model: logconcm ~ p2 - p3 + log(dose/(exp(p1) - exp(p2)) * (exp(-exp(p2) * tm) - exp(-exp(p1) * tm)))
Data: datam
AIC BIC logLik
511.0899 539.6245 -249.5449
Random effects:
Formula: p3 ~ 1 | subj
p3 Residual
StdDev: 0.1998116 0.2759416
Variance function:
Structure: Exponential of variance covariate
Formula: ~fitted(.)
Parameter estimates:
expon
-0.42927
Fixed effects: list(p1 + p2 + p3 ~ 1)
Value Std.Error DF t-value p-value
p1 -0.395378 0.03740172 839 -10.57110 <.0001
p2 -3.772868 0.02390463 839 -157.83002 <.0001
p3 0.488964 0.06350886 839 7.69914 <.0001
Correlation:
p1 p2
p2 -0.473
p3 -0.608 0.499
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-8.3650734 -0.3426087 0.1233143 0.6027582 3.2001694
Number of Observations: 859
Number of Groups: 18
Warning message:
NA/Inf replaced by maximum positive value
the volume is estimated to be
> exp(0.488964)
[1] 1.630626
which is medically meaningless. (Unfortunately, from summary() there
is no way to know that the second model does fit far better than the
first. These are not the same data as in previous emails: here
metabolite, previously parent drug.) This estimate does not change if
all three parameters are random. For comparison, non-normal
autoregression models fit better than these normal random coefficient
models and do not distort the estimate of the volume.
Jim
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !) To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
More information about the R-help
mailing list