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