###########################################
#Kreft ,& de Leeuw (1998).
#Varying and random coefficient models.
#In Introducing Multilevel Modeling (pp.36--56).
#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)


#-----------------------------------
##Slope-as-outcome analysis
#-----------------------------------

##(a) Table 3.2
tab3.2 <- data.frame(matrix(ncol=7, nrow=11))
colnames(tab3.2) <- c("intercept", "intercept.se", "slope", "slope.se", "r", "N", "Pu.Pr")
rownames(tab3.2) <- c(levels(nels88$SCHID), "Total")
for(i in 1:10){
  lm.temp <- lm(MATH~HOMEWORK, data=nels88, subset=SCHID==levels(SCHID)[i])
  lm.temp.coef <- summary(lm.temp)$coef
  cor.temp <- cor.test(~ MATH + HOMEWORK, data=nels88, subset=SCHID==levels(SCHID)[i])
  pu.temp <- nels88$PUBLIC[is.element(nels88$SCHID, levels(nels88$SCHID)[i])][1]
  tab3.2[i, 1:7] <- c(lm.temp.coef[1, 1:2], lm.temp.coef[2, 1:2], cor.temp$estimate, cor.temp$parameter+2, pu.temp)
}

lm.temp <- lm(MATH~HOMEWORK, data=nels88)
lm.temp.coef <- summary(lm.temp)$coef
cor.temp <- cor.test(~ MATH + HOMEWORK, data=nels88)
tab3.2[11, 1:6] <- c(lm.temp.coef[1, 1:2], lm.temp.coef[2, 1:2], cor.temp$estimate, cor.temp$parameter+2)
tab3.2


##(b) Figure 3.7
y <- seq(20, 100, length=9)
x <- seq(0,  7, length=9)
plot(y~x, type="n", main="OLS regression lines for 10 schools")
for(i in 1:11){
  if(i != 11){
    abline(coef=tab3.2[i, c(1, 3)])
  }else{
    abline(coef=tab3.2[i, c(1, 3)], lty=2)
  }
}
graphics.off()


##(c) two macro equation
sl.as.outlm1 <- lm(intercept~Pu.Pr, data=tab3.2[-11,])
sl.as.outlm2 <- lm(slope~Pu.Pr, data=tab3.2[-11,])

predict(sl.as.outlm1, newdata=data.frame(Pu.Pr=c(0, 1)))  #the effects of private and public sector
predict(sl.as.outlm2, newdata=data.frame(Pu.Pr=c(0, 1)))  #the effects of private and public sector



#-----------------------------------
##Random coefficient results
#-----------------------------------
library(lme4)
library(arm)
library(nlme)

##Table 3.3
random.lmer1 <- lmer(MATH~HOMEWORK + (HOMEWORK|SCHID), data=nels88)
random.lmer1

random.lme2  <- lme(MATH~HOMEWORK, random=~HOMEWORK | SCHID, data=nels88)
intervals(random.lme2)


##Table 3.4
random.lmer2 <- lmer(MATH~HOMEWORK + PUBLIC + (HOMEWORK|SCHID), data=nels88)
random.lmer2
anova(random.lmer1, random.lmer2)

random.lme2  <- lme(MATH~HOMEWORK + PUBLIC, random=~HOMEWORK | SCHID, data=nels88)
intervals(random.lme2)


##Table 3.5
random.lmer3 <- lmer(MATH~HOMEWORK * PUBLIC + (HOMEWORK|SCHID), data=nels88)
random.lmer3

random.lme3  <- lme(MATH~HOMEWORK * PUBLIC, random=~HOMEWORK | SCHID, data=nels88)
intervals(random.lme3)