[R] nlme

[Douglas Bates]

Jim Lindsey <jlindsey at alpha.luc.ac.be> writes:

> Among others, datam contains the columns: logconc, tm, dose, subj, bilirubin.
> None of these are factor variables.

You didn't list `age' in there.  Did you mean to have it listed?

I would have converted `subj' to a factor.  Actually, it will be
converted to a factor anyway before being used in groups = ~ subj.

> The following compartment models work (the first still has not
> converged after 100 interations):
> 
> res1 <- nlme(logconc~p2+p3+log(dose/(exp(p1)-exp(p2))*
> 	(exp(-exp(p2)*tm)-exp(-exp(p1)*tm))),start=list(fixed=c(5,-2,-0.1)),
> 	fixed=list(p1+p2+p3~1),control=list(maxIter=100),
> 	groups=~subj,data=datam,verbose=T,method="ML")
> 
> res2 <- nlme(logconc~p2+p3+log(dose/(exp(p1)-exp(p2))*
> 	(exp(-exp(p2)*tm)-exp(-exp(p1)*tm))),start=list(fixed=c(5,-2,-0.1)),
> 	fixed=list(p1+p2+p3~1),random=pdDiag(p1+p2+p3~1),
> 	groups=~subj,data=datam,verbose=T,method="ML")
> 
> res3 <- nlme(logconc~p2+p3+log(dose/(exp(p1)-exp(p2))*
> 	(exp(-exp(p2)*tm)-exp(-exp(p1)*tm))),start=list(fixed=c(5,-2,-0.1,0)),
> 	fixed=list(p1+p2~1,p3~bilirubin),random=pdDiag(p1+p2+p3~1),
> 	groups=~subj,data=datam,verbose=T,method="ML")
> 
> However, when I add a second covariate,
> 
> res4 <- nlme(logconc~p2+p3+log(dose/(exp(p1)-exp(p2))*
> 	(exp(-exp(p2)*tm)-exp(-exp(p1)*tm))),
> 	start=list(fixed=c(5,-2,-0.1,0,0)),
> 	fixed=list(p1+p2~1,p3~bilirubin+age),random=pdDiag(p1+p2+p3~1),
> 	groups=~subj,data=datam,verbose=T,method="ML")
> 
> I get the error:
> 
> Error in fixed[[nm]][[3]] != "1" : comparison (2) is possible only for vector types

Determining the order of the starting values is a tricky business.
I'll have to try it to figure out what is going on.

I think I know of a way of avoiding a lot of these problems but
haven't had the opportunity to code it up and try it out yet.

> If I try instead
> 
> res4a <- nlme(logconc~p2+p3+p4*bilirubin+p5*age+log(dose/(exp(p1)-exp(p2))*
> 	(exp(-exp(p2)*tm)-exp(-exp(p1)*tm))),
> 	start=list(fixed=c(5,-2,-0.1,0,0)),
> 	fixed=list(p1+p2+p3+p4+p5~1),random=pdDiag(p1+p2+p3~1),
> 	groups=~subj,data=datam,verbose=T,method="ML")
> 
> it runs.
> 
> Can anyone explain to me what I am doing wrong?
> 
> Sorry but I cannot supply the data as I am under nondisclosure to a
> pharmaceutical company.

I appreciate that you can't disclose the data you have but could you
at least send the output of
 str(datum)
so we can see exactly what it looks like?  Otherwise it will be hard
to debug the problem.  Also please report the system type, the version
of R, and the version of the nlme package.

If you wish you can blank out those few individual data values that
are printed by str(datum).  I just want to see exactly what the
structure is.

>   How does one make the variance depend on a function of the mean as
> is usual in PK modelling? I can only find how to make it depend on
> covariates in the docs?

Which docs?  According to the help page for `varPower', for example,
the default formula is `~ fitted(.)', which is what you are looking
for.

   Usage:

        varPower(value, form, fixed)

   Arguments:
    ... 

       form: an optional one-sided formula of the form `~ v',
             or `~ v | g', specifying a variance covariate `v'
             and, optionally, a grouping factor `g' for the
             coefficients. The variance covariate must evaluate
             to a numeric vector and may involve expressions
             using `"."', representing  a fitted model object
             from which fitted values (`fitted(.)') and residu-
             als (`resid(.)') can be extracted (this allows the
             variance covariate to be updated during the opti-
             mization of an object function). When a grouping
             factor is present in `form', a different coeffi-
             cient value is used for each of its levels. Sev-
             eral grouping variables may be simultaneously
             specified, separated by the `*' operator, like in
             `~ v | g1 * g2 * g3'. In this case, the levels of
             each grouping variable are pasted together and the
             resulting factor is used to group the observa-
             tions. Defaults to `~ fitted(.)'  representing a
             variance covariate given by the fitted values of a
             fitted model object and no grouping factor.


>   It would be nice if nls printed out the log likelihood so
> that the fits could be compared with those from nlme.

You can obtain the log-likelihood for the fitted nls model `fm' with
 logLik(fm)

> I also find it very annoying that nlme stops with an error when the
> iteration limit is exceeded (as in res1 above), returning nothing so
> that one cannot even inspect the partially converged results.

I'm afraid we disagree on the approach.  When you have not converged
you don't have estimates and you should not be doing analysis as if
you had estimates.  I have seen too many situations where people have
used PROC NLIN in SAS to fit a nonlinear regression model, had it fail
to converge giving one brief message to that effect and then carry on
to print several pages of results exactly as if it had converged.  The
user happily walks off with the results never noticing that they are
nonsense.
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._