マルチレベルモデル

http://www.ats.ucla.edu/stat/r/examples/mlm_imm/immch4.htm を参考にした。

演習用に作ったデータなのでこの結果を用いて実践の改善を図ろうとしないこと

データ

Kreft & de Leeuw (1998) のサンプルデータ
- http://www.ats.ucla.edu/stat/mplus/examples/imm/
もとはNELS:88（米国の大規模縦断調査）のデータ
このデータを演習用に改変した
この演習用データは以下から取得可能 https://dl.dropboxusercontent.com/u/7515903/dat.csv

コードブック

変数	内容
schid	学校番号
stuid	生徒番号
ses	生徒の家庭のSES
meanses	学校の平均SES
homework	1週間あたり宿題時間
white	人種：白人を1，それ以外を0
parented	親の学歴
public	公立私立：公立を1
ratio	PT比
percmin	（不明）
math	数学得点
sex	性別
race	民族
sctype	（学校のタイプ？）
cstr	（不明）
scsize	学校規模
urban	都市
region	地域

一般的な回帰分析の問題点

たとえば，3校を対象に調査を行い 宿題に取り組む時間 と 数学の得点 との関係を検討したとする。学校別に2変数の関係をプロットすると以下の通りとなり，まとめて回帰直線を描くと太線のようになる。

# データ読み込み
dat <- read.csv("dat.csv")
# 学校別のサブセット
s26537 <- subset(dat, schid == 26537)
s47583 <- subset(dat, schid == 47583)
s54344 <- subset(dat, schid == 54344)
# 学校別に色分けした散布図
plot(s26537$homework, s26537$math, pch = 16, col = 2, xlab = "", ylab = "", 
    xlim = c(0, 7), ylim = c(0, 75))
par(new = T)
plot(s47583$homework, s47583$math, pch = 16, col = 3, xlab = "", ylab = "", 
    xlim = c(0, 7), ylim = c(0, 75))
par(new = T)
plot(s54344$homework, s54344$math, pch = 16, col = 4, xlab = "", ylab = "", 
    xlim = c(0, 7), ylim = c(0, 75))
# 回帰直線
lmres <- lm(math ~ homework, data = dat)  #回帰係数を求めて格納
abline(lmres, lwd = 2, lty = 3)

plot of chunk HLM2 plot1

この図では3校とも同様の散布状況ではなく，学校ごとに違いがあるように見える。
学校ごとに回帰直線を描くと以下の色つきの直線の通りとなる。

# 学校別に色分けした散布図
plot(s26537$homework, s26537$math, pch = 16, col = 2, xlab = "", ylab = "", 
    xlim = c(0, 7), ylim = c(0, 75))
par(new = T)
plot(s47583$homework, s47583$math, pch = 16, col = 3, xlab = "", ylab = "", 
    xlim = c(0, 7), ylim = c(0, 75))
par(new = T)
plot(s54344$homework, s54344$math, pch = 16, col = 4, xlab = "", ylab = "", 
    xlim = c(0, 7), ylim = c(0, 75))
# 回帰直線
lm26537 <- lm(math ~ homework, data = s26537)
lm47583 <- lm(math ~ homework, data = s47583)
lm54344 <- lm(math ~ homework, data = s54344)

abline(lm26537, lwd = 2, col = 2)
abline(lm47583, lwd = 2, col = 3)
abline(lm54344, lwd = 2, col = 4)
abline(lmres, lwd = 2, lty = 3)

plot of chunk HLM2 plot2

このように，学校別で見ると回帰直線が平行であり，かつ傾きの勾配が学校別で見た場合とまとめて見た場合とで異なることがある。
この場合の学校をまとめた分析結果は宿題に取り組む時間の効果が過大評価されているといえる。たとえばこの結果を鵜呑みにして「宿題に取り組む時間を増やすと学力向上」と指導しても，実際にはその向上度合いが小さいといったことが起こると，子どもはがっかりする。
このデータの場合，そもそも数学の得点が学校間で差があることに加え，宿題に取り組む時間も学校ごとに異なる。
このようなデータに対しては，データの階層構造を仮定した分析を行う方がよい。

ランダム切片モデル

考え方

上図を見ると，傾きは3校とも同じだが切片が異なると考えられる。
そこで，切片は学校ごとに異なる （ランダム切片） が，傾きは各校同じであると仮定したモデルを組み， どの学校でも共通と思われる宿題に取り組む時間が数学の得点に与える影響の大きさ を推定したい。
一般的な回帰分析は以下の式で表現される。 \[ Y = \beta_0 + \beta_1 X_1 + r \]
ランダム切片モデルは切片\( \beta_0 \)が集団によって異なるため，以下の式で表現する。 \[ Y = \beta_0 + \beta_1 X_1 + r \] \[ \beta_0 = \gamma_{00} + u_0 \]
この例では\( Y \)が数学得点(math)，\( X \)が宿題時間(homework)，\( r \)が誤差となる。

Rで推定

lme4パッケージを用いる。無い場合はインストール(install.packages(“lme4”))する。
以下はmathをhomeworkで予測させるという点ではlmと同じ。(1|schid)が切片にあたり，schidでランダムであることを示している。

library(lme4)

## Loading required package: lattice
## Loading required package: Matrix

res <- lmer(math ~ homework + (1 | schid), data = dat)
summary(res)

## Linear mixed model fit by REML ['lmerMod']
## Formula: math ~ homework + (1 | schid) 
##    Data: dat 
## 
## REML criterion at convergence: 370.2 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  schid    (Intercept) 32.0     5.66    
##  Residual             48.6     6.97    
## Number of obs: 55, groups: schid, 3
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)   40.899      3.705   11.04
## homework       3.653      0.629    5.81
## 
## Correlation of Fixed Effects:
##          (Intr)
## homework -0.396

なお，一般的な回帰分析でmathをhomeworkで予測させた場合の結果は以下の通りとなる。ランダム切片モデルの方が回帰係数が小さく，上図のようなデータとのあてはまりが良さそうだということがわかる。

summary(lmres)

## 
## Call:
## lm(formula = math ~ homework, data = dat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -11.676  -6.934  -0.297   4.445  21.324 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   38.918      1.990   19.55  < 2e-16 ***
## homework       4.379      0.709    6.18  9.3e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Residual standard error: 8.27 on 53 degrees of freedom
## Multiple R-squared: 0.419,   Adjusted R-squared: 0.408 
## F-statistic: 38.2 on 1 and 53 DF,  p-value: 9.33e-08

補遺

このデータでは以下のような特徴が見られる。

mathとses，mathとmeansesとの相関が高い。

cor(dat$math, dat$ses)

## [1] 0.4814

cor(dat$math, dat$meanses)

## [1] 0.4006

sesとhomeworkとの相関もそれなりにある。

cor(dat$ses, dat$homework)

## [1] 0.1965

学校別で見ると以下のような相関が見られる。

homeworkとmath

cor(s26537$homework, s26537$math)

## [1] 0.7235

cor(s47583$homework, s47583$math)

## [1] 0.5631

cor(s54344$homework, s54344$math)

## [1] 0.5871

sesとmath

cor(s26537$ses, s26537$math)

## [1] 0.1195

cor(s47583$ses, s47583$math)

## [1] 0.5125

cor(s54344$ses, s54344$math)

## [1] 0.3738

sesとhomework

cor(s26537$ses, s26537$homework)

## [1] 0.188

cor(s47583$ses, s47583$homework)

## [1] 0.4574

cor(s54344$ses, s54344$homework)

## [1] -0.2506

ランダム傾きモデル

\[ Y = \beta_0 + \beta_1 X_1 + r \] \[ \beta_1 = \gamma_{00} + u_0 \]

lmer(math~homework+(0+homework|schid), data=dat)