[R] nlme

[Jim Lindsey]

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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._