###########################################
#Kreft ,& de Leeuw (1998).
#Analyses.
#In Introducing Multilevel Modeling (pp.57--104).
#Thousand Oaks: Sage.
#
#データの所在: Textbook Examples Introduction to Multilevel Modeling
#http://www.ats.ucla.edu/stat/examples/imm/default.htm
###########################################


##データの読み込み
library(doBy)
library(lme4)
library(arm)
nels88b <- read.table("http://blue.zero.jp/yokumura/R/multilevel/nels88b.csv", sep=",", header=T)


##Table 4.1
stu.ex.var <- c("SES", "HOMEWORK", "WHITE", "PUBLIC", "MEANSES", "PERCMIN", "RATIO")
nels88b[,stu.ex.var]
cor(nels88b[,stu.ex.var])



##---------------------------------#
# Session 1 (pp.63--71)
##---------------------------------#

##Table 4.2
ses1.model0 <- lmer(MATH~(1|SCHID), data=nels88b)
ses1.model0
residual <- unlist(sigma.hat(ses1.model0)$sigma) ^ 2
ses1.model0.corr <- as.numeric(residual[2] / sum(residual))
ses1.model0.corr


##Table 4.3
ses1.model1 <- lmer(MATH~HOMEWORK+(1|SCHID), data=nels88b)
ses1.model1
residual <- unlist(sigma.hat(ses1.model1)$sigma) ^ 2
ses1.model1.corr <- as.numeric(residual[2] / sum(residual))
ses1.model1.corr
anova(ses1.model0, ses1.model1)


##Table 4.4
ses1.model2 <- lmer(MATH~HOMEWORK+(HOMEWORK|SCHID), data=nels88b)
ses1.model2
anova(ses1.model1, ses1.model2)


##Table 4.5
ses1.model3 <- lmer(MATH~HOMEWORK+ PARENTED + (HOMEWORK|SCHID), data=nels88b)
ses1.model3


##Table 4.6
ses1.model4 <- lm(MATH~HOMEWORK + PARENTED, data=nels88b)
summary(ses1.model4)



##---------------------------------#
# Session 2 (pp.71--90)
##---------------------------------#

##Table 4.7
ses2.model2 <- lmer(MATH~HOMEWORK + SCSIZE + (HOMEWORK|SCHID), data=nels88b)
ses2.model2
anova(ses1.model2, ses2.model2)


##Table 4.8
ses2.model3 <- lmer(MATH~HOMEWORK + PUBLIC + (HOMEWORK|SCHID), data=nels88b)
ses2.model3
anova(ses1.model2, ses2.model3)


##Table 4.9
attach(nels88b)
  nels88b$HOMEWORK.PUBLIC <- HOMEWORK * PUBLIC
  nels88b$HOMEWORK.MEANSES <- HOMEWORK * MEANSES
  nels88b$HOMEWORK.RATIO <- HOMEWORK * RATIO
detach(nels88b)

varlist <- c("PUBLIC", "MEANSES", "RATIO", "HOMEWORK", "HOMEWORK.PUBLIC", "HOMEWORK.MEANSES", "HOMEWORK.RATIO")
cor(nels88b[, varlist])


##Table 4.10
ses2.model4 <- lmer(MATH~HOMEWORK * PUBLIC + (HOMEWORK|SCHID), data=nels88b)
ses2.model4
anova(ses1.model0, ses2.model4)
anova(ses1.model2, ses2.model4)


##Table 4.12
ses2.model5 <- lmer(MATH~HOMEWORK + WHITE + PUBLIC + (HOMEWORK|SCHID), data=nels88b)
anova(ses1.model2, ses2.model5)


##Table 4.13
ses2.model6 <- lmer(MATH~HOMEWORK + WHITE + PUBLIC + (HOMEWORK+WHITE|SCHID), data=nels88b)
ses2.model6
anova(ses2.model5, ses2.model6)


##Table 4.14
ses2.model7 <- lmer(MATH~HOMEWORK + WHITE + PUBLIC + MEANSES+ (HOMEWORK|SCHID), data=nels88b)
ses2.model7
anova(ses2.model5, ses2.model7)


##Table 4.15
ses2.model8 <- lmer(MATH~HOMEWORK + WHITE + MEANSES+ (HOMEWORK|SCHID), data=nels88b)
ses2.model8
anova(ses2.model7, ses2.model8)
deviance(ses2.model5); deviance(ses2.model8)


##Table 4.16
ses2.model9 <- lmer(MATH~HOMEWORK * MEANSES + WHITE + (HOMEWORK|SCHID), data=nels88b)
ses2.model9
anova(ses2.model8, ses2.model9)


##Table 4.17
ses2.model10 <- lmer(MATH~HOMEWORK + WHITE + SES+ MEANSES + (HOMEWORK|SCHID), data=nels88b)
ses2.model10
anova(ses2.model8, ses2.model10)



##---------------------------------#
# Session 3 (pp.90--98)
##---------------------------------#

##Table 4.19
ses3.model1 <- lmer(MATH~ SES + (1|SCHID), data=nels88b)
ses3.model1
anova(ses1.model0, ses3.model1)


##Table 4.20 (Error)
ses3.model2 <- lmer(MATH~ SES + (SES|SCHID), data=nels88b)
ses3.model2


##Table 4.21
ses3.model3 <- lmer(MATH~ SES + PERCMIN + (1|SCHID), data=nels88b)
ses3.model3
anova(ses3.model1, ses3.model3)


##Table 4.22
ses3.model4 <- lmer(MATH~ SES + PERCMIN + MEANSES+ (1|SCHID), data=nels88b)
ses3.model4
anova(ses3.model1, ses3.model4)
anova(ses3.model3, ses3.model4)



##---------------------------------#
# Session 4 (pp.98--102)
##---------------------------------#

##Table 4.25
ses4.model1 <- lmer(MATH~ HOMEWORK + RATIO + (HOMEWORK|SCHID), data=nels88b)
ses4.model1
anova(ses1.model2, ses4.model1)


##Table 4.26
ses4.model2 <- lmer(MATH~ HOMEWORK + HOMEWORK:RATIO + (HOMEWORK|SCHID), data=nels88b)
ses4.model2
anova(ses1.model2, ses4.model2)