# Pothoff & Roy (1964); Jennrich & Schluchter (1986)
# Pinheiro (1994) PhD U WI Madison Statistics
pothoff <- read.table("pothoff2.dat",header=T)
pothoff$person <- factor(pothoff$person)
pothoff$agefac <- factor(pothoff$age)

# split plot
pothoff.fit <- aov(y~agefac*sex+person,pothoff)
print(summary(pothoff.fit))

# split plot using person as random effect
# see Venables & Ripley (1994, sec. 6.7) or Chambers & Hastie (1992, sec. 5.2.1)
pothoff.bfit <- aov(y~agefac*sex+Error(person),pothoff)
print(summary(pothoff.bfit))

# split plot using LME (Splus on ALPHA computers only)
# see Lindstrom & Bates (1988) JASA 83:1014-1022
library(nlme)
# split plot (works, but this is not the main use of lme)
pothoff.lme <- lme(fixed=y~agefac*sex,random=~1,
	cluster=~person,data=pothoff,
	est.method="RML", re.structure="identity")
print(pothoff.lme)
print(anova(pothoff.lme))
print(summary(pothoff.lme))
#print(fitted.values(pothoff.lme))

# random slopes by age*sex (this is in the spirit of lme use)
pothoff.slme <- lme(y~age*sex,random=~age,cluster=~person,data=pothoff,
	est.method="RML", re.structure="identity")
print(pothoff.slme)
print(anova(pothoff.slme))
print(summary(pothoff.slme))
print(anova(pothoff.lme,pothoff.slme))