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

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


##(b) 平均 (table 2.1)
tab2.1 <- summaryBy(MATH + HOMEWORK ~ SCHID, data=nels88, FUN=c(mean, length))
tab2.1


##(c) 共分散関数の定義 
scov <- function(dat){
      X    <- t(dat)
      N    <- nrow(dat)
      a    <- rep(1, N)
      Xa   <- rbind(X, a)
      xbar  <- 1/N * X %*% a
      xbar1 <- xbar %*% a
      V     <- X - xbar1
      S     <- 1/N * V %*% t(V)
      return(S)
}


##(d) 分散共分散行列と相関 (table 2.2)
tab2.2 <- data.frame(matrix(nrow=nlevels(nels88$SCHID)*2, ncol=4))
for(i in 1:nlevels(nels88$SCHID)){
  temp <- subset(nels88, SCHID==LETTERS[i], c(MATH, HOMEWORK))
  tab2.2[i*2-1, c(1, 4)] <- c(LETTERS[i], cor(temp)[1,2])
  tab2.2[c(i*2-1, i*2), c(2, 3)] <- scov(temp)
}
tab2.2


##(e) Total or ppoled regression (Table 2.3)
lm0 <- lm(MATH~1, data=nels88)
lm1 <- lm(MATH~HOMEWORK, data=nels88)
summary(lm0); summary(lm1)


##(f) Aggregate regression (Table 2.4)
agglm0 <- lm(MATH.mean ~ 1, data=tab2.1, weights=MATH.length)
agglm1 <- lm(MATH.mean ~ HOMEWORK.mean, data=tab2.1, weights=MATH.length)
summary(agglm0); summary(agglm1)


##(g) Contextual model (Table 2.5)
nels88$home.by.sch <- nels88$SCHID
nels88$home.by.sch <- as.numeric(as.character(factor(nels88$home.by.sch, labels=tab2.1$HOMEWORK.mean)))

contlm0 <- lm(MATH~1, data=nels88)
contlm1 <- lm(MATH~HOMEWORK+home.by.sch, data=nels88)
summary(contlm0); summary(contlm1)


##(h) Cronbach model (Table 2.6) 教科書のinterceptは，おそらく誤植 
nels88$home.bw <- nels88$HOMEWORK - nels88$home.by.sch
nels88$home.bb <- nels88$home.by.sch - mean(nels88$HOMEWORK)

cronlm0 <- lm(MATH~1, data=nels88)
cronlm1 <- lm(MATH~home.bw + home.bb, data=nels88)
summary(cronlm0); summary(cronlm1)


##(g) Analysis of covariance
ancovalm0 <- lm(MATH~SCHID-1, data=nels88)
ancovalm1 <- lm(MATH~SCHID+HOMEWORK-1, data=nels88)
summary(ancovalm0); summary(ancovalm1)