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

[Ken Beath]

On 13/12/2007, at 2:41 AM, David Atkins wrote:

>
> Ken--
>
> Many thanks for the reply and the pointer to Mary Lindstrom's post.  
> However, I can't seem to get her function to work with nls().  I  
> assume that the underscores in her function should be replaced with  
> assignments (ie, <-)?
>
> My slightly modified version of her function is below, and data  
> attached.  Any thoughts?
>
> cheers, Dave
>
> # http://www.biostat.wustl.edu/archives/html/s-news/2000-04/msg00209.html
> #
> # It is also possible to fit a linear break point model where the  
> break
> # point is a true parameter in the model The model function below will
> # work in nls.  It includes a very short quadratic piece where the  
> break
> # point is to keep the derivative w.r.t. the break point  
> continuous.  -
> # Mary Lindstrom
>
> hockey <- function(x, alpha1, beta1, beta2, brk, eps = range(x)/100) {

should be eps = diff(range(x)/100)

>
>       ## alpha1 is the intercept of the left line segment
>       ## beta1 is the slope of the left line segment
>       ## beta2 is the slope of the right line segment
>       ## brk is location of the break point
>       ## 2*eps is the length of the connecting quadratic piece
>
>       ## reference: Bacon & Watts "Estimating the Transition Between
>       ## Two Intersecting Straight Lines", Biometrika, 1971
>
>        x1 <- brk - eps
>        x2 <- brk + eps
>        b <- (x2*beta1 - x1*beta2)/(x2 - x1)
>        cc <- (beta2 - b)/(2*x2)
>        a <- alpha1 + beta1*x1 - b*x1 - cc*x1^2
>        alpha2 <- - beta2*x2 + (a + b*x2 + cc*x2^2)
>
>        lebrk <- (x <= brk - eps)
>        gebrk <- (x >= brk + eps)
>        eqbrk <- (x > brk - eps & x < brk + eps)
>
>        result <- rep(0,length(x))
>        result[lebrk] <- alpha1 + beta1*x[lebrk]
>        result[eqbrk] <- a + b*x[eqbrk] + cc*x[eqbrk]^2
>        result[gebrk] <- alpha2 + beta2*x[gebrk]
>        result
> }
>
> test.nls <- nls(dv ~ hockey(weeks, alpha1, beta1, beta2, brk),
>                          data = small.df)
>
> Dave Atkins, PhD
> Associate Professor in Clinical Psychology
> Fuller Graduate School of Psychology
> Email: datkins at fuller.edu
> Phone: 626.584.5554
>

The zero derivatives are due to the interval around the breakpoint not  
containing any data, I suggest placing the initial changepoint in the  
centre of the data and then increasing eps until it works.

I expect things will get even worse with a mixed model so it would be  
worthwhile plotting observed versus predicted for each subject to make  
sure the breakpoints are in sensible places. I think this may be a  
difficult model to fit.

Ken