myscratchpad

「基礎から学ぶマルチレベルモデル」をRで実行
# Kreft, I., & de Leeuw, J. (1998). Introducing multilevel modeling. Newbury Park, CA: Sage. (クレフト, I., デレウ, J. 小野寺孝義 (編訳) (2006). 基礎から学ぶマルチレベルモデル―入り組んだ文脈から新たな理論を創出するための統計手法　ナカニシヤ出版
# 第3章

uclaや奥村先生のサイトを参考にした。というかほとんどコピペ。感謝。
上記サイトではlme4パッケージのlmer関数が主だが、lme (nlmeパッケージ) を使えばいいじゃないか、的な記事もあったのでとりあえずlmeで書き直してみた。
RjpWikiのオブジェクト一覧も有用。感謝

library(foreign)
dat0 <- data.frame(read.spss("http://www.ats.ucla.edu/stat/spss/examples/mlm_imm/imm10.sav"))
names(dat0) <- tolower(names(dat0))
data.frame(names(dat0))
dat <- dat0[c(1,2,11,5,8,19)]
head(dat)
dat.agmean <- aggregate(dat[3:5], list(dat[,1]), mean)
dat.agmean

# 傾きを結果変数とするアプローチ
## 表3.2　10校の宿題に対する数学の成績のOLS回帰　p41
by(dat[3:4], dat[1], cor) # 相関係数
by(dat, dat$schid, function(x) coef(summary(lm(math~homework, x)))) # 学校ごとにデータを分割し、math~homeworkの回帰分析を行ったもの
summary(lm(math~homework, dat)) #p24 表2.3と同じ

## マクロ回帰の結果
bylmres <- by(dat, dat$schid, function(x) lm(math~homework, x))
a <- sapply(bylmres, coef)
intc <- a[1,]
slp <- a[2,]
pbl <- factor(dat.agmean[,4])
summary(lm(intc~pbl))
summary(lm(slp~pbl))

## 表3.3　ランダム係数モデル p42
# レベル1 Y=math, X = homework, i=個人, j=学校, ε=誤差
## β_0j, β_1j はそれぞれa_j, b_jとも表記される (教科書中) 。学校の切片よび傾きを表す。
## 最初の添え字 (i) はミクロ水準、後の添え字は (j) はマクロ水準を示す。つまりiはネストされる変数、jはネストする変数。この例では個人がi, 学校がj。また、例えば日ごとの反復測定ならiが日、jが個人となる
## 表3.3、水準1の分散はε_ijであり、重回帰分析でいうところの誤差項にあたる。VarCorr関数のResidualsにあたる。
# Y_ij = β_0j + β_1jX_ij + ε_ij

# レベル2 # γ₀₀=切片の全体の平均, γ₁₀=全体の傾き、u=誤差
## VarCorr関数で出力される分散成分、および教科書表3.3の「水準2」の「推定値」とはここのu_0j, u_1jの誤差分散を示す。
# β_0j = γ₀₀ + u_0j
# β_1j = γ₁₀ + u_1j
library(nlme)
lme.res <- lme(math ~ homework, random=~homework|schid, data=dat, method="ML")
summary(lme.res)
VarCorr(lme.res)
getVarCov(lme.res) # 分散
logLik(lme.res)*-2 # 逸脱度
intervals(lme.res) # 信頼区間

## 表3.4 p44
lme.res2 <- lme(math ~ homework+public, random=~homework|schid, data=dat, method="ML")
summary(lme.res2)
VarCorr(lme.res2)
getVarCov(lme.res2) # 分散
logLik(lme.res2)*-2 # 逸脱度
intervals(lme.res2) # 信頼区間


## 表3.5 p45
lme.res3 <- lme(math~homework+public+homework:public, random=~homework|schid,  data=dat, method="ML")
summary(lme.res3)
VarCorr(lme.res3)
getVarCov(lme.res3) # 分散
logLik(lme.res3)*-2 # 逸脱度
intervals(lme.res3) # 信頼区間

## 表3.6 p46
lme.res4 <- lme(math~homework, random=~homework|schid, data=dat, method="ML")
summary(lme.res4)
VarCorr(lme.res4)
getVarCov(lme.res4) # 分散
logLik(lme.res4)*-2 # 逸脱度
intervals(lme.res4) # 信頼区間
cf.v <- coef(lme.res4)
cf.v


## 図3.8
n<-nrow(cf.v)
x<-cbind(0, 7)
y<-cbind(seq(1:n), seq(1:n))
y[,1]<-cf.v[,1]
y[,2]<-cf.v[,1]+cf.v[,2]*7
plot(x, y[1,], type="l", ylab="predicte", ylim=c(20,100))
for (i in 2:n){points(x, y[i,], type="l", lty=i)}