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

[Dieter Menne]

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)
}