#================================================
# Multilevel data
# Example: Smoking prevention and cessation data
#===============================================

# Read smoking data from the web:
data.name="http://www.stat.ufl.edu/~aa/glm/data/Smoking.dat"
smoking=read.table(data.name,header=T)
head(smoking)

# Use the nlme library 
# (alternatively one may use the lme4 library, cf page 298 in the text book)
library(nlme)

# Fit model with random effect for school and class:
fit.smoke=lme(y~PTHK+SC+TV,random=list(~1|school,~1|class),data=smoking) 
summary(fit.smoke)

# A linear model fit that disregards the hierarchical structure of the data, 
#will give too small standard errors:
summary(lm(y~PTHK+SC+TV,data=smoking))