[R-sig-ME] nlmer: gradient for random change point model

[David Atkins]

Hi all--

I am attempting to implement a nonlinear mixed-effects model in nlmer() 
that I have gotten to work (well, at least run) in nlme().  However, 
nlmer() requires the gradient function, and I've gotten a bit stuck there.

Here's a bit of background, and some data are attached with script below:

In psychotherapy research we often see two phases of improvement: an 
early rapid phase, and a later slower phase.  Thus, a piecewise linear 
model can be a reasonably good fit, except that the breakpoint is 
different for different folk.  The following article describes how it is 
possible to set up a nonlinear mixed-effects model like this with a 
random effect for the breakpoint:

Cudeck, R., & Klebe, K. J. (2002). Multiphase mixed-effects models for 
repeated measures data. Psychological Methods, 7, 41-63.

A biostatistician colleague of mine has used SAS to fit such a model to 
some weekly psychotherapy data of a mutual colleague.  However, NLMIXED 
didn't particularly like our data, and he had to search over a wide grid 
of initial starting values to ever get it to converge.

Given this and Doug's recent efforts with the lme4 library, I thought it 
might be worthwhile to try fitting the model using nlmer().  I've got it 
running in nlme(), though it often bonks and gives error msgs.  At the 
moment, I haven't been successful in getting it to run in nlmer().  This 
  appears to be due to the fact that nlmer() wants a gradient explicitly 
provided (whereas nlme() doesn't), and deriv() doesn't recognize pmax() 
that we use in the nonlinear formula.

The data attached are a small subset of the actual data, which are quite 
a bit larger (241 individuals) and more "messy" (i.e., some mostly 
linear, some rather chaotic growth patterns), which no doubt lead to 
some of the difficulty in model fitting.  Finally, this is my first 
foray into nonlinear mixed-effects models, though I regularly use linear 
and generalized linear mixed-effects models.  Any general pointers or 
obvious "user errors" appreciated.

Thanks in advance.

cheers, Dave
-- 
Dave Atkins, PhD
Associate Professor in Clinical Psychology
Fuller Graduate School of Psychology
Email: datkins at fuller.edu
Phone: 626.584.5554

 > sessionInfo()
R version 2.6.1 (2007-11-26)
i386-pc-mingw32

locale:
LC_COLLATE=English_United States.1252;LC_CTYPE=English_United 
States.1252;LC_MONETARY=English_United 
States.1252;LC_NUMERIC=C;LC_TIME=English_United States.1252

attached base packages:
[1] stats     graphics  grDevices datasets  utils     methods   base

other attached packages:
[1] nlme_3.1-86       Matrix_0.999375-3 lattice_0.17-2

loaded via a namespace (and not attached):
[1] gdata_2.3.1     grid_2.6.1      gtools_2.4.0    lme4_0.999375-0 
tools_2.6.1


load("randomChangePoint.Rdata")

### get a picture of the data
library(lattice)
xyplot(dv ~ weeks | id, data = small.df, type = c("g","p","smooth"),
       pch = 16)

### fit simple linear model with "forced" break at 6 weeks
test.lm <- lm(dv ~ weeks + pmax(0, (weeks - 6)),
                        data = small.df)
summary(test.lm)

### fit nonlinear model
test.nls <- nls(dv ~ b0 + b1*weeks + b2*pmax(0, (weeks - brk)),
                           data = small.df,
                           start = c(b0 = 28, b1 = -4.1, b2 = 4.1, brk = 6))
summary(test.nls)

### use nlme() to fit nonlinear mixed-effects model, though keep it to
### diagonal random-effects for now
library(nlme)

test.nlme <- nlme(dv ~ b0 + b1*weeks + b2*pmax(0, (weeks - brk)),
                       data = small.df,
                       fixed = b0 + b1 + b2 + brk ~ 1,
                       random = pdDiag(b0 + b1 + b2 + brk ~ 1),
                       groups = ~ id,
                       control = list(msVerbose = TRUE, niterEM = 200,
                                     opt = "nlminb", tol = 1e-5),
                       start = coef(test.nls))
summary(test.nlme)

### swap packages
detach("package:nlme")

library(lme4)

test.nlmer <- nlmer(dv ~ b0 + b1*weeks + b2*pmax(0, (weeks - brk))
                         ~  (b0 | id),
                         data = small.df,
                         verbose = TRUE,
                         start = coef(test.nls))
### gives error about gradient, though deriv() doesn't like pmax()