Fitting piecewise linear curves with lmer

Question

(cross-posted from the R-sig ME list - maybe this is too basic for it)

Hello I would like to fit piecewise linear growth curves to a number of subjects (each with a number of observations at different abscissa), where all the slopes for a common interval (e.g, 2 <= x < 3) presumable arise from the same normal distribution, and therefore should be "shrunk" towards each other. What is the best way to about doing this in lme4? I have a toy, cleaned up (although I am not very good with R) example of the type of data I am looking at. Thank you for your time

set.seed(1)
# our growth functions are piecewise linear with these breaks
breaks <- seq(0,5)
mean.slopes <- seq(5,2,-.6)

# generate the true slope parameters for 20 subjects (variance isn't
    really constant or known, just simplifying)...
true.slopes <- sapply(mean.slopes,function(x) rnorm(n=20,sd=3,mean=x))

# ...and the intercepts - they really don't come from any particular
    distribution
true.intercepts <- runif(n=20,min=30,max=50)

# generate sample data - we don't necessarily have the same number of
observations per subject
mat <- matrix(nr=20*10,nc=3)
for (i in 1:20)
{
   xs <- runif(n=10,min=0,max=6)
   slopes <- true.slopes[i,]
   ys <- true.intercepts[i] + rnorm(n=10,mean=0,sd=0.9) +
   sapply(xs, function(x) sum(slopes[1:floor(x)]) + (x - floor(x)) *
   slopes[ceiling(x)])
   id <- rep(i,10)
   mat[(1+10*(i-1)):(10*i),] <- cbind(xs,ys,id)
}
mat <- as.data.frame(mat)
names(mat) <- c("x","y","id")
# qplot(x=x,y=y,colour=as.factor(id),data=mat)

Eventually I might like to shrink the slopes for each individual towards each other, in the sense that if for a given subject the slope in the interval [2,3) is very small compared to the slopes of other subjects in the same interval, then the same would hold for the slope in [1,2), but I would like to get the basics first.

Answer 1

Apologies for only posting code with no comments. I wrote this but don't have time now to comment and decided it was best to save this here and comment later.

# packages needed
library(lattice)
library(latticeExtra)
library(lme4)

# make a picture of the data
mat <- mat[order(mat$id, mat$x),]
mat$id <- factor(mat$id)
plot1 <- xyplot(y~x, group=id, data=mat, type="b")
plot1

# function to set up knots
knot <- function(x, knot) {(x-knot)*(x>knot)}
knots <- function(x, knots) {
  out <- sapply(knots, function(k) knot(x, k))
  colnames(out) <- knots
  out
}

# add knots to data frame
mat$knot <- knots(mat$x, 0:5)

# piecewise curves with no random effects
m1 <- lm(y~ knot, data=mat)
summary(m1)

# get predicted values
matX <- data.frame(x=0:6)
matX$knot <- knots(matX$x, 0:5)
matX$predict <- predict(m1, newdata=matX)

# plot of data and predicted values
plot1 <- xyplot(y~x, group=id, data=mat, type="b")
plot2 <- xyplot(predict~x, data=matX, type="b", col="black", lwd=3)
plot1+plot2

# piecewise curves with random effects
m2 <- lmer(y~ knot + (knot|id), data=mat)
summary(m2)
ranef(m2)

# get predicted values
ids <- unique(mat$id)
matXid <- expand.grid(id=unique(mat$id), x=0:6)
matXid$knot <- knots(matXid$x, 0:5)
matXid$predict <- rowSums(as.matrix(coef(m2)$id[matXid$id,]) * cbind(1,matXid$knot))

# plot of data and predicted values
plot1 <- xyplot(y~x|id, data=mat, type="b", as.table=TRUE)
plot2 <- xyplot(predict~x|id, data=matXid, type="b", col="black", lwd=1)