# Data in Table 9.3, page 177 timber <- c(0,2.38,4.34,6.64,8.05,9.78,10.97,12.05,12.98,13.94,14.74,16.13,17.98,19.52,19.97,0,2.69,4.75,7.04,9.2,10.94,12.23,13.19,14.08,14.66,15.37,16.89,17.78,18.41,18.97,0,2.85,4.89,6.61,8.09,9.72,11.03,12.14,13.18,14.12,15.09,16.68,17.94,18.22,19.40,0,2.46,4.28,5.88,7.43,8.32,9.92,11.1,12.23,13.24,14.19,16.07,17.43,18.36,18.93,0,2.97,4.68,6.66,8.11,9.64,11.06,12.25,13.35,14.54,15.53,17.38,18.76,19.81,20.62,0,3.96,6.46,8.14,9.35,10.72,11.84,12.85,13.83,14.85,15.79,17.39,18.44,19.46,20.05,0,3.17,5.33,7.14,8.29,9.86,11.07,12.13,13.15,14.09,15.11,16.69,17.69,18.71,19.54,0,3.36,5.45,7.08,8.32,9.91,11.06,12.21,13.16,14.05,14.96,16.24,17.34,18.23,18.87) timber <- matrix(timber, byrow=T, ncol=15) x <- c(0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0,1.2,1.4,1.6,1.8) slippage <- rep(x,8) loads <- as.vector(t(timber)) specimen <- rep(1:8,rep(15,8)) timber.dat <- data.frame(specimen, slippage, loads) # loading the nlme package # You may have to install this package from the CRAN web site first: # If so, type: install.packages("nlme") library(nlme) attach(timber.dat) # Random intercept model: timber.lme <- lme(loads~slippage, random=~1|specimen, data=timber.dat, method = "ML") # Random intercept and slope model: timber.lme1 <- lme(loads~slippage, random=~slippage|specimen, data=timber.dat, method = "ML") # Comparing the two models: anova(timber.lme, timber.lme1) # Getting the estimates in Table 9.4 of the book: summary(timber.lme1) # Plotting the observed data alongside the predicted values from the model: predictions <- matrix(predict(timber.lme1),ncol=15,byrow=T) par(mfrow=c(1,2)) matplot(x,t(timber), type='l', col=1, xlab="Slippage", ylab="Load", lty=1, ylim=c(0,25)) title("(a)") matplot(x, t(predictions), type='l', col=1, xlab="Slippage", ylab="Load", lty=1, ylim=c(0,25)) title("(b)") # Random intercept and slope model, with a quadratic term: slippage.squared <- slippage^2 timber.lme2 <- lme(loads~slippage+slippage.squared, random=~slippage|specimen, data=timber.dat,method="ML") # Comparing the model without the quadratic term to the model with the quadratic term: anova(timber.lme1, timber.lme2) # Getting the estimates in Table 9.5 of the book: summary(timber.lme2) # Plotting the observed data alongside the predicted values from the model: predictions <- matrix(predict(timber.lme2),ncol=15,byrow=T) par(mfrow=c(1,2)) matplot(x,t(timber), type='l', col=1, xlab="Slippage", ylab="Load", lty=1, ylim=c(0,25)) title("(a)") matplot(x, t(predictions), type='l', col=1, xlab="Slippage", ylab="Load", lty=1, ylim=c(0,25)) title("(b)") # Random intercept and slope model, with a quadratic term, constrained to go through origin (see Problem 9.1 in book): timber.lme2.noint <- lme(loads~0+slippage+slippage.squared, random=~0+slippage|specimen, data=timber.dat,method="ML") # Comparing the quadratic model without the intercepts to the quadratic model with the intercepts: anova(timber.lme2.noint, timber.lme2) # Getting the estimates: summary(timber.lme2.noint) # Plotting the observed data alongside the predicted values from the model: predictions <- matrix(predict(timber.lme2.noint),ncol=15,byrow=T) par(mfrow=c(1,2)) matplot(x,t(timber), type='l', col=1, xlab="Slippage", ylab="Load", lty=1, ylim=c(0,25)) title("(a)") matplot(x, t(predictions), type='l', col=1, xlab="Slippage", ylab="Load", lty=1, ylim=c(0,25)) title("(b)") # Is this a better result???