== R code for a mixed model - a worked example ==

Here is an example for comparing Beck depression inventory scores between two treatment groups over three time points (baseline, post and follow-up) which firstly fits the mixed model to obtain regression coefficients corresponding to the group by time interaction and then extracts the design matrix corresponding to the levels of group and time to reproduce the standard errors for each group at each time point (Thanks to Adam Wagner for the code).  

{{{
# Exploring http://stackoverflow.com/questions/14358811/extract-prediction-band-from-lme-fit

# Fit BDI model
print("Fitting BDI model")
bdiModelA   <- lme(bdiScore~Group*time,
                   random=~1|ID,
                   BDIData,
                   method="ML",
                   na.action=na.omit)

summary(bdiModelA)



print("Getting a predicted fit of the fixed terms and SEs")
# Create a dataframe for the predicted fit
names(BDIData)
bdiPredFitData <- expand.grid(Group=levels(BDIData$Group),
                              time=levels(BDIData$time))
# Add in point estimate from model ie add a third column called Pred
bdiPredFitData$Pred <- predict(bdiModelA, newdata=bdiPredFitData, level=0)
bdiPredFitData



# Create a design matrix for these predictions
# NOTE: bdiPredFitData[-ncol(bdiPredFitData)] is simply a way of passing in the data 
#       sans the predictor point estimate (which is in the last column ie column three)
DesignmatBDIPred <- model.matrix(eval(eval(bdiModelA$call$fixed)[-2]), 
                                 bdiPredFitData[-ncol(bdiPredFitData)])


# Compute the standard errors for predictions; the $varFix bit is the variance-covariance matrix of the fixed regression estimates from the mixed model
predvarBDI <- diag(DesignmatBDIPred %*% bdiModelA$varFix %*% t(DesignmatBDIPred))

# Standard error of the fix part of the model adding a fourth column to bdiPredFitData called SE
bdiPredFitData$SE <- sqrt(predvarBDI)

# Standard error of the fix part of the model plus the within-group error SD adding a fifth column to bdiPredFitData called SE2
bdiPredFitData$SE2 <- sqrt(predvarBDI+bdiModelA$sigma^2)
}}}

The predictions and standard errors for both groups at each of the three time point can then be plotted as error bar plot in R using the below syntax:

{{{
# Show the profile
library(plotrix)

# Plot variables
# Labels
grpTitle <- "BDI-II profiles: Mean and 95% CI"
XLabel <- "Time point"
YLabel <- "Average BDI-II score"
subLabel <- "N= observations in each group by time combination"
timeLabels <- c("Baseline", "Post", "Follow-up")
nObsYLine <- 2
groupPch <- c(15, 16)
groupCol <- c(1,2)
# X Locations
xCentrePts <- c(1, 2, 3) 
xLims <- c(0.75, 3.25)
xSplt <- 0.2


# Determine plotting points
xPts <- rep(c(-xSplt, xSplt), 3)
xPts <- rep(xCentrePts, each=2)+xPts
# Calculate the number of observations at each point
bdiObsN <- data.frame(with(subset(BDIData, !is.na(bdiScore)), 
                           table(Group, time)))
# Check that the rows of bdiPredFitData and match up
groupMatchFail <- any((!(bdiObsN$Group==bdiPredFitData$Group)))
timeMatchFail <-  any((!(bdiObsN$time==bdiPredFitData$time))) 
if (groupMatchFail|timeMatchFail){
  stop("The rows of bdiPredFitData and bdiObsN are not in the correct order")
}


# Main plotting details
plotCI(xPts, y=bdiPredFitData$Pred, uiw=1.96*bdiPredFitData$SE,
       xlim=xLims, ylim=c(0,30), xaxt="n", 
       pch=groupPch[as.numeric(bdiPredFitData$Group)],
       col=groupCol[as.numeric(bdiPredFitData$Group)],
       xlab=XLabel, ylab=YLabel, main=grpTitle,
       sub=subLabel)
# Divide up plotting area
abline(v=c(mean(xCentrePts[1:2]), mean(xCentrePts[2:3])), lty=2)
# Label time points
axis(1, xCentrePts, timeLabels)
legend(x="top", c("Control", "Treatment"), pch=groupPch,
       col=groupCol)
# Include observations at each combination of group and time
text(x=xPts, y=nObsYLine, 
     labels=paste("N=", bdiObsN$Freq, sep=""))
}}}