5 Rによる分析

フリーの統計解析ソフトであるR(R Development Core Team, 2005)には，混合モデルを分析するためのパッケージnmleが用意されている。このパッケージにあるlme関数を用いることで，階層的線形モデルと同様の分析を行うことができる ¹ 。

nlmeパッケージをインストールすると，サンプルデータとしてHSBデータがついてくるので，それを用いて今回行ったのと同じ分析を行うことができる。 ²

5.1 データ

SASの場合と同様，データファイルは１つにまとめる必要がある。以下に，HSBのデータを整形したものの最初の10ケース目までを示した。

   school    ses mathach   meanses        cses sector
1    1224 -1.528   5.876 -0.434383 -1.09361702 Public
2    1224 -0.588  19.708 -0.434383 -0.15361702 Public
3    1224 -0.528  20.349 -0.434383 -0.09361702 Public
4    1224 -0.668   8.781 -0.434383 -0.23361702 Public
5    1224 -0.158  17.898 -0.434383  0.27638298 Public
6    1224  0.022   4.583 -0.434383  0.45638298 Public
7    1224 -0.618  -2.832 -0.434383 -0.18361702 Public
8    1224 -0.998   0.523 -0.434383 -0.56361702 Public
9    1224 -0.888   1.527 -0.434383 -0.45361702 Public
10   1224 -0.458  21.521 -0.434383 -0.02361702 Public

5.2 Rプログラム

#----------------------------------------------------------------------#
#                          RによるHSBデータの分析
#                            By Taichi OKUMURA
#----------------------------------------------------------------------#

# nlmeパッケージの読み込み
library(nlme)

# レベル1のデータの読み込み
data(MathAchieve)
MathAchieve[1:10,]

# レベル2のデータの読み込み
data(MathAchSchool)
MathAchSchool[1:10,]

# 学校ごとのSESの平均をmsesという新しい変数にする。
attach(MathAchieve)
mses <- tapply(SES, School, mean)
detach(MathAchieve)

# HSBというデータセット
HSB <- as.data.frame(MathAchieve[, c("School", "SES", "MathAch")])
names(HSB) <- c("school", "ses", "mathach")

# 学校ごとのSESの平均をmsesという新しい変数にする。
HSB$meanses <- mses[as.character(HSB$school)]
HSB$cses <- HSB$ses - HSB$meanses

# sector変数を取り出してHSBに投入
sector <- MathAchSchool$Sector
names(sector) <- row.names(MathAchSchool)
HSB$sector <- sector[as.character(HSB$school)]
HSB$sector <- factor(HSB$sector, levels=c('Public', 'Catholic'))

# result.hsbに分析結果を格納
result.hsb <- lme(mathach ~ sector + meanses + cses + 
                            sector*cses + meanses*cses,
                  random =  ~ cses | school, 
                  data = HSB)
# 結果の表示
summary(result.hsb)

5.3 Rのアウトプット

上記のプログラムを実行したときのアウトプットを以下に示す。 ³

Linear mixed-effects model fit by REML /* REML推定（デフォルト） */
 Data: HSB 
       AIC      BIC    logLik 
  46523.66 46592.45 -23251.83 /* LogLik（対数尤度）の値に-2をかけたものがdeviance */

Random effects:
 Formula: ~cses | school
 Structure: General positive-definite, Log-Cholesky parametrization
            StdDev    Corr  
(Intercept) 1.5426150 (Intr) /* レベル2の分散の平方根（切片） */
cses        0.3182015 0.391  /* レベル2の分散の平方根（傾き） */
Residual    6.0597955        /* レベル1の分散の平方根 */

/* レベル2の偏回帰係数の推定結果 */
Fixed effects: mathach ~ sector + meanses + cses + sector * cses + meanses * cses 
                        Value Std.Error   DF  t-value p-value
(Intercept)         12.127931 0.1992920 7022 60.85510   0e+00  /* G00 */
sectorCatholic       1.226579 0.3062733  157  4.00485   1e-04  /* G01 */
meanses              5.332875 0.3691684  157 14.44564   0e+00  /* G02 */
cses                 2.945041 0.1556005 7022 18.92694   0e+00  /* G10 */
sectorCatholic:cses -1.642674 0.2397800 7022 -6.85076   0e+00  /* G11 */
meanses:cses         1.039230 0.2988971 7022  3.47688   5e-04  /* G12 */
 Correlation: 
                    (Intr) sctrCt meanss cses   sctrC:
sectorCatholic      -0.699                            
meanses              0.256 -0.356                     
cses                 0.075 -0.053  0.019              
sectorCatholic:cses -0.052  0.077 -0.027 -0.696       
meanses:cses         0.019 -0.026  0.074  0.293 -0.351

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-3.1592608 -0.7231893  0.0170471  0.7544510  2.9582205 

Number of Observations: 7185 /* データに含まれる生徒の数 */
Number of Groups: 160        /* 学校の数 */

Footnotes

http://socserv.socsci.mcmaster.ca/jfox/Books/Companion/appendix-mixed-models.pdfに詳しい解説がある。
ただし，以下のプログラムで示したようにデータを整形する必要がある。
" $/*$ " 〜 " $*/$ "の部分はコメントである。