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

[vito muggeo]

A completely (simpler and naive) approach could be to use an algorithm 
which works reasonably well in non-mixed GLM (see the package segmented).

The algorithm works via linear approximations of the kink function, 
therefore extension to mixed models is, at least in principle, 
straightforward. I discuss such an approach in a paper published in the 
Proceedings of the 19th IWSM, and I have (very) raw R code implementing 
the idea.

Let me know off-list if you are interested in the paper and/or code,

best,
vito


estimate the cluster-specific parameters (includ

Ken Beath ha scritto:
> 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
> 
> _______________________________________________
> R-sig-mixed-models at r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> 
> 

-- 
====================================
Vito M.R. Muggeo
Dip.to Sc Statist e Matem `Vianelli'
Università di Palermo
viale delle Scienze, edificio 13
90128 Palermo - ITALY
tel: 091 6626240
fax: 091 485726/485612