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

[David Atkins]

looks like the attachment didn't come through, so data can be found:

http://www.fuller.edu/sop/faculty/datkins/cp_content/randomChangePoint.Rdata


-------- Original Message --------
Subject: nlmer: gradient for random change point model
Date: Tue, 11 Dec 2007 15:26:17 -0800
From: David Atkins <datkins at fuller.edu>
To: r-sig-mixed-models at r-project.org


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



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