Hierarchical linear growth

gistfile1.txt

# Multiple subjects -- correlated intercept and slope ---------------------

N <- 50

times <- seq(0,10,1)

interceptmu <- 4 #intercept / starting point

slopemu <- .3 #slope / growth rate

interceptsd <- 2 #sd of individual differences in intercept

slopesd <- .1 #sd of random differences in slope

interceptsloperegression <- -.1 #effect between the two

errorsd <- .2 #random / measurement error

#creating correlated data for the individual intercepts / slopes

intercepts <- rnorm(N,interceptmu,interceptsd) 

slopes <- rnorm(N, slopemu, slopesd) + #proportion of new, slope specific info

  (intercepts-interceptmu)*interceptsloperegression #plus proportion of info from intercepts

#check the subject specific parameters have the right characteristics

#increase N temporarily (eg 10000) to improve precision of these checks

mean(slopes)

mean(intercepts)

sd(slopes)

sd(intercepts)

cor(slopes, intercepts)

#now generate the data

dat <- data.frame( #create data.frame to store our data

  id= rep(1:N,each=length(times)), #make sure each subject has a row in the dataframe for each time of observation

  time= rep(times, N), #repeat the time sequence N subjects times.

  y=NA #fill in a blank observation variable that we fill using a loop later

)

head(dat) #view the first rows of the data.frame

for(i in 1:nrow(dat)){

  id <- dat$id[i] #set a temporary id variable so we can refer to each subjects intercept / slope easily

  dat$y[i] <-  intercepts[id] + slopes[id] * dat$time[i]

}

dat$y <- dat$y + rnorm(length(dat$y),0,errorsd) #include random error

dat$id <- factor(dat$id) #so that R knows that id is a categorical variable

#plot

library(ggplot2)

ggplot(dat,aes(x=time,y=y,colour=id))+

  geom_line() #different plotting styles are possible

#fit

library(lme4)

fit <- lmer(y~time + (time|id),dat) #random slope (with random intercept by default)

summary(fit)