#R program for the Singer Article on Hierchical Models

#download the appropriate package

library(nlme)

#Examples using the hsb12 data set.

#The following section of code is not in the article but it is 
#always a good idea to get a look at the data if possible. 
#Let us look at the data and create a groupedData object.
#Let's look at only a few of the many schools in this vast dataset

#inter.only <- groupedData( mathach ~ 1 | school, data=hsb12small) 

#plot(inter.only)
#gsummary(inter.only, FUN= mean)

#Table 1, Random intercept

hsb12.int <- groupedData( mathach ~ 1 | school, data=hsb12) 
basic <- lme(mathach ~ 1, data=hsb12.int, random = ~ 1 | school) 
print(summary(basic))

#Table 2, random intercept, predictor = meanses

hsb12.meanses <- groupedData(mathach ~ 1 | school, data=hsb12, outer= ~ meanses)
fit1 <- lme(mathach ~ meanses, data=hsb12.meanses, random = ~ 1 | school)
print(summary(fit1))

# compute cses

hsb12$cses <- hsb12$ses - hsb12$meanses
hsb12.cses <- groupedData( mathach ~ cses | school, data=hsb12)  

# use unstructured correlation (default in lme function, corr = NULL)

fit2 <- lme(mathach ~ cses, data=hsb12.cses, random = ~ cses | school )
print(summary(fit2))

#table 4

hsb12.big <- groupedData( mathach ~ meanses + sector + cses | school, data=hsb12)  
fit3 <- lme(mathach ~ meanses + sector + cses + meanses:cses + sector:cses, 
             data=hsb12.big, random = ~ cses | school)
print(summary(fit3))

fit4 <- lme(mathach ~ meanses + sector + cses + meanses:cses + sector:cses, 
             data=hsb12.big, random = ~ 1 | school)
print(summary(fit4))

#Example 2:  Willett data

#grouping the data so that we can graph it 
willettg <- groupedData( y ~ time | id, data=willett)

#Looking at the data:
trellis.device()
print(plot(willettg))

#table on p. 342
wfit1 <- lme(y ~ time, data=willettg, random = ~ time | id)
print(summary(wfit1))

trellis.device()
print(plot( wfit1, id ~ resid(.), abline = 0))

trellis.device()
print(plot( augPred(wfit1), layout=c(4,3), col=8 ))

#next model with both level data

#centering the covar variable
willett$ccovar <- willett$covar - mean(willett$covar)

#grouping the data  
willettg1 <- groupedData( y ~ time | id, data=willett, outer=~ccovar)

#table on p. 344
wfit2 <- lme(y ~ time*ccovar, data=willettg1, random = ~ time | id)
print(summary(wfit2))

#windows(6,6)
trellis.device()
print(plot( wfit2, id ~ resid(.), abline = 0))

trellis.device()
print(plot( augPred(wfit2), layout=c(4,3), col=8 ))

#Strictly repeated measure, table on p. 344

#compound symmetry
fit.gls.cs <- gls(y~time, data=willett, corr=corCompSymm(, form=~time | id))
print(summary(fit.gls.cs))

#autoregressive
fit.gls.ar1 <- gls(y~time, data=willett, corr=corAR1(, form=~time | id))

#comparing the models with different covariance structures, can't test the dif since same df
print(anova(fit.gls.cs, fit.gls.ar1))

#NOTE: the unstructured is not available!!!

#table on p. 248
wfit2 <- lme( y ~ time*ccovar, data=willettg1, random = ~ time | id, 
              corr=corAR1(, form =~ time|id))
print(summary(wfit2))

trellis.device()
print(plot( augPred(wfit2), layout=c(4,3), col=8))