[R-sig-ME] Gastric emptying data: nlmer vs. nlme

[Douglas Bates]

Dieter,

Thanks for the very thorough description of the problem and for
including a reproducible example.  I enclose a modified version of
your script and the output using the development version of the lme4
package (https://svn.R-project.org/R-packages/branches/gappy-lmer/)

You will see that the development version does better on the first
model but still ends up giving a "false convergence" message.  I then
fit ge0A.nlmer and added another model ge0B.nlmer that removes the
random effect for kappa.  The (conservative) p-value for a test of H0:
ge0B.nlmer versus Ha: ge0A.nlmer is

> pchisq(979.76761 - 973.73274, df = 2)
[1] 0.9510734

so incorporating random effects for kappa is, at best, marginally significant.

The reason that the first model is so difficult to fit is because
there are 6 variance-covariance parameters and only 8 levels of
subj:treat.  You are trying to estimate too many variance-covariance
parameters from too few groups.  The likelihood surface will be very
flat and the parameter estimates will be ill-defined.

The reason that nlmer from the 0.99875-8 release of lme4 gave up has
to do with the calculation of the conditional modes of the random
effects to evaluate the Laplace approximation to the deviance.  In
that version I retained the values of the conditional modes of the
random effects (these are the values the maximum the conditional
density of the random effects given the data and the current values of
the model parameters - the Laplace approximation is evaluated at these
values and the adaptive Gauss-Hermite approximation is centered around
these values) between evaluations of the deviance.  That's a good idea
because the vector of the conditional modes for a different set of
parameters is going to be similar to the current values so you have
good starting estimates.  However, it is a bad idea in that the
evaluation of the approximation to the deviance will not only depend
on the values of the parameters and the data but also on where the
last value was taken.  This means that the function value being
optimized is not reproducible and that causes a lot of problems in a
derivative-free optimization.

To avoid this I now start each evaluation of the conditional modes at
the same point (all random effects start at zero) so the evaluation is
reproducible.

On 10/1/07, Dieter Menne <dieter.menne at menne-biomed.de> wrote:
> Dear Group,
>
> I have tried to redo published nlme-fits of gastric emptying data recorded by
> MRI with nlmer. nlme needs time, but give reasonable results, informing me about
> the known correlation between two nonlinear parameters, while nlmer gave up
> after two iterations without telling me why.
>
> I hope I got the translation to nlmer syntax correctly. What's wrong?
>
> Dieter Menne
>
>
> ----------------------------------------------------
> # Gastric Emptying data from
> # Goetze et. al (University Hospital of Zürich)
> # http://ajpgi.physiology.org/cgi/content/abstract/00498.2005v1
> # Demo data set from http://www.menne-biomed.de/gastempt
> # R-Version: 2.5.1, i386, mingw32
> # lme4_0.99875-8.zip
>
> useNlme = TRUE
>
> # Get data
> ge = read.table("http://www.menne-biomed.de/gastempt/gastempt.csv",
>   sep=",",header=TRUE)
> #ge = read.table("gastempt.csv",sep=",",header=TRUE)
> ge$subj = as.factor(ge$subj)
> ge$subjtreat = as.factor(paste(ge$subj,ge$treat,sep="."))
> ge$t = as.double(ge$t)
>
>
> EmptInit= function(mCall,LHS,data){ # dummy, not explicitely used
>   stop("Should not be called")
> }
>
> # Standard LinExp model for gastric emptying
>
> SSEmptLinExp=selfStart(~v0*(1+kappa*t/tempt)*exp(-t/tempt),
>   initial=EmptInit, parameters= c("v0","kappa","tempt"))
>
> # nlme final is 643,1.4,74
> start = list(fixed=c(v0=600,kappa=1.0,tempt=70))
>
> if (useNlme) {
>   library(nlme)
>   contr = nlmeControl(pnlsTol=0.3,pnlsMaxIter=200,msMaxIter=200)
>   ge0.nlme=nlme(v~SSEmptLinExp(t,v0,kappa,tempt),
>     fixed =list(v0+kappa+tempt~1),
>     groups = ~subjtreat,
>     control=contr, start=start,  data=ge, verbose=F
>     )
>   summary(ge0.nlme)
>   # Note correlation between kappa and tempt
>   ge0A.nlme = update(ge0.nlme,random=v0+kappa~1)
>   summary(ge0A.nlme)
>   anova(ge0A.nlme,ge0.nlme)
> }
>
> if (!useNlme) {
>   library(lme4) # lme4_0.99875-8.zip
>   # This stops after two iterations without telling me details
>   ge0.nlmer = nlmer(
>     v~SSEmptLinExp(t,v0,kappa,tempt)~(v0+kappa+tempt|subjtreat),
>     data=ge,  start=start,verbose=T)
>   show(ge0.nlmer)
>   # This is close to the result of ge0A.nlme)
>   ge0A.nlmer = nlmer(
>     v~SSEmptLinExp(t,v0,kappa,tempt)~(v0+kappa|subjtreat),
>     data=ge,  start=start,verbose=T)
>   show(ge0.nlmer)
> }
>
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
>