# turnip.R
library(nlme)
library(car)
turnipdata <- read.table(file="turnip.dat",header=T) # read the data in
head(turnipdata)
turnipdata$Leaf <- factor(turnipdata$leaf) # Leaf is a factor
turnip2 <- groupedData(calcium~1|Leaf, data=turnipdata)
plot(turnip2)
m1 <- lme(calcium~1,data=turnip2, random= ~1|Leaf)
summary(m1)
# notice that the summary of model m1 gives estimates of sigma_a and
# sigma, not sigma^2_a and sigma^2. These estimates are given in the
# following portion of the summary:
# Random effects:
# Formula: ~1 | Leaf
# (Intercept) Residual
#StdDev: 0.2690357 0.08125321
# These values agree with the estimates given by SAS (see turnip.sas) since
# 0.2690357^2 = .07238, and 0.08125321^2 = .006602
# Confidence intervals for sigma_a and sigma are given by the following
# command. Note these are for sigma_a and sigma, not sigma^2_a and sigma^2
# Furthermore, they are obtained with a different method than the one used in
# SAS and presented in the lecture notes, so these intervals don't agree with
# those given by PROC MIXED. E.g., you can't square the endpoints of these
# intervals to get the intervals given by SAS. I recommend using the intervals
# that are given by SAS and in my lecture notes rather than those given by
# the intervals() function as illustrated below.
intervals(m1,which="var-cov")