library(nlme)
library(lattice) 

data(MathAchieve)
head(MathAchieve)

data(MathAchSchool)
head(MathAchSchool)

attach(MathAchieve)
mses <- tapply(SES, School, mean) 
mses[as.character(MathAchSchool$School[1:10])]
detach(MathAchieve)

Bryk <- as.data.frame(MathAchieve[, c("School", "SES", "MathAch")])
names(Bryk) <- c("school", "ses", "mathach")
sample20 <- sort(sample(7185, 20)) # 20 randomly sampled students
Bryk[sample20, ]

Bryk$meanses <- mses[as.character(Bryk$school)]
Bryk$cses <- Bryk$ses - Bryk$meanses
sector <- MathAchSchool$Sector
names(sector) <- row.names(MathAchSchool)
Bryk$sector <- sector[as.character(Bryk$school)]
Bryk[sample20,]



attach(Bryk)
Cat <- sample(unique(school[sector=='Catholic']), 20)
Cat.20 <- groupedData(mathach ~ ses | school,
                      data=Bryk[is.element(school, Cat),])
Pub <- sample(unique(school[sector=='Public']), 20)
Pub.20 <- groupedData(mathach ~ ses | school,
                      data=Bryk[is.element(school, Pub),])

#################################################################

trellis.device(color=F)
xyplot(mathach ~ ses | school, data=Cat.20, main="Catholic",
	panel=function(x, y){
		panel.xyplot(x, y)
		panel.loess(x, y, span=1)
		panel.lmline(x, y, lty=2)
	}
)

xyplot(mathach ~ ses | school, data=Pub.20, main="Public",
	panel=function(x, y){
		panel.xyplot(x, y)
		panel.loess(x, y, span=1)
		panel.lmline(x, y, lty=2)
	}
)

####################################################################
cat.list <- lmList(mathach ~ ses | school, subset = sector=='Catholic',
			   	data=Bryk)
pub.list <- lmList(mathach ~ ses | school, subset = sector=='Public',
				data=Bryk)

plot(intervals(cat.list), main='Catholic')
plot(intervals(pub.list), main='Public')

cat.coef <- coef(cat.list); cat.coef[1:10,]
pub.coef <- coef(pub.list); pub.coef[1:10,]

###################################################################

old <- par(mfrow=c(1,2))
boxplot(cat.coef[,1], pub.coef[,1], main='Intercepts',
		names=c('Catholic', 'Public'))
boxplot(cat.coef[,2], pub.coef[,2], main='Slopes',
		names=c('Catholic', 'Public'))
par(old)

###################################################################

Bryk$sector <- factor(Bryk$sector, levels=c('Public', 'Catholic'))
contrasts(Bryk$sector)


bryk.lme.1 <- lme(mathach ~ meanses*cses + sector*cses,
			random = ~ cses | school,
			data=Bryk)
summary(bryk.lme.1)

bryk.lme.2 <- update(bryk.lme.1,
	random = ~ 1 | school) # omitting random effect of cses
anova(bryk.lme.1, bryk.lme.2)

bryk.lme.3 <- update(bryk.lme.1,
	random = ~ cses - 1 | school) # omitting random intercept
anova(bryk.lme.1, bryk.lme.3)

summary(bryk.lme.2)
summary(bryk.lme.2)