文書の過去の版を表示しています。
マルチレベルモデルに関するメモ
[安藤正人(2011)『マルチレベルモデル入門 実習：継時データ分析』ナカニシヤ出版]に基づいて、マルチレベルモデルを勉強したメモ。 本書は、SPSS/SASをベースに分析が進められているが、この内容をRでできるだけ再現してみる。Rについても、併せて勉強する。
R環境の準備
必要なライブラリをあらかじめインストールしておく。線形モデルを扱う関数(lm)は標準で入っているが、分散分析でType2/3平方和を計算する関数や混合モデルのためのライブラリが無いので、carおよびlme4をインストールしておく。
install.packages("car");
install.packages("lme4");
library(car);
library(lme4);
サンプルデータも読み込んでおく。
df1 <- read.csv("sample.csv");
一般線形モデル
分散分析、共分散分析、（重）回帰分析を一般線形モデルで統一的に把握する。
1元配置分散分析
lm関数で学科間の平均得点について分散分析する。
>m1.2.1 <- lm(ACHIEVE ~ DEPART, data=df1);
> summary(m1.2.1);
Call:
lm(formula = ACHIEVE ~ DEPART, data = df1)
Residuals:
Min 1Q Median 3Q Max
-7.000 -3.542 0.000 2.625 7.500
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 85.500 1.828 46.772 <2e-16 ***
DEPARTb -6.833 2.585 -2.643 0.0129 *
DEPARTc -6.833 2.585 -2.643 0.0129 *
DEPARTd -3.167 2.585 -1.225 0.2301
DEPARTe 1.500 2.585 0.580 0.5661
DEPARTf -2.500 2.585 -0.967 0.3413
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.478 on 30 degrees of freedom
Multiple R-squared: 0.3701, Adjusted R-squared: 0.2652
F-statistic: 3.526 on 5 and 30 DF, p-value: 0.01261
> anova(m1.2.1);
Analysis of Variance Table
Response: ACHIEVE
Df Sum Sq Mean Sq F value Pr(>F)
DEPART 5 353.47 70.694 3.5259 0.01261 *
Residuals 30 601.50 20.050
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
学科（DEPART）変数はカテゴリカルなので、Rでは自動的に第1番目のカテゴリを基準にダミー変数になって分析される。この場合は学科aが基準になり、Coefficientsの（Intercept）の値が学科aの平均の推定値となる。その他の学科については、Interceptにそれぞれの係数を加えて計算すると、書籍に示される推定値に一致する。
標本平均の標準誤差（Std. Error）については、従来の分散分析では学科の不偏分散を観測数で除した値の平方根となるため、学科ごとに異なる値になるが、一般線形モデルでは分析モデルの誤差の平方平均に基づいて決定されるため、学科間で値は共通となる。この例では、各学科の平均の推定値の標準誤差は、誤差の平均平方（20.050）を各学科における観測数6で除した値の平方根になる。Rでの分析では、この値は（Intercept）のStd. Errorの値（1.828）となる。Rで表示される各学科の係数の標準誤差の値については、書籍では言及なし。
平方和、R2乗、調整済みR2乗、F値、p値などについては、Rの各出力と書籍の内容は一致する。
ネストした分散分析
下位カテゴリが上位カテゴリにもれなく含まれるような場合の分散分析。上位カテゴリで説明される分散と下位カテゴリで説明される分散を分けてみる。両者の合計は、上位カテゴリだけで分析した場合に説明される分散の大きさと一致する。
ここでは、学科は学部にネストされているデータを用いる。まずは、学部だけの分散分析の結果。
> m1.3.1 <- lm(ACHIEVE ~ FACULTY, data=df1);
> summary(m1.3.1);
Call:
lm(formula = ACHIEVE ~ FACULTY, data = df1)
Residuals:
Min 1Q Median 3Q Max
-7.9444 -4.3194 -0.1111 3.8889 12.0556
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 80.944 1.189 68.096 <2e-16 ***
FACULTYS 3.167 1.681 1.884 0.0682 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 5.043 on 34 degrees of freedom
Multiple R-squared: 0.09451, Adjusted R-squared: 0.06787
F-statistic: 3.549 on 1 and 34 DF, p-value: 0.06817
> anova(m1.3.1);
Analysis of Variance Table
Response: ACHIEVE
Df Sum Sq Mean Sq F value Pr(>F)
FACULTY 1 90.25 90.250 3.5485 0.06817 .
Residuals 34 864.72 25.433
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
書籍と一致。
続いて、ネストした分散分析。
> m1.3.2 <- lm(ACHIEVE ~ FACULTY/DEPART, data=df1);
> summary(m1.3.2);
Call:
lm(formula = ACHIEVE ~ FACULTY/DEPART, data = df1)
Residuals:
Min 1Q Median 3Q Max
-7.000 -3.542 0.000 2.625 7.500
Coefficients: (6 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 85.5000 1.8280 46.772 <2e-16 ***
FACULTYS -2.5000 2.5852 -0.967 0.3413
FACULTYL:DEPARTb -6.8333 2.5852 -2.643 0.0129 *
FACULTYS:DEPARTb NA NA NA NA
FACULTYL:DEPARTc -6.8333 2.5852 -2.643 0.0129 *
FACULTYS:DEPARTc NA NA NA NA
FACULTYL:DEPARTd NA NA NA NA
FACULTYS:DEPARTd -0.6667 2.5852 -0.258 0.7983
FACULTYL:DEPARTe NA NA NA NA
FACULTYS:DEPARTe 4.0000 2.5852 1.547 0.1323
FACULTYL:DEPARTf NA NA NA NA
FACULTYS:DEPARTf NA NA NA NA
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 4.478 on 30 degrees of freedom
Multiple R-squared: 0.3701, Adjusted R-squared: 0.2652
F-statistic: 3.526 on 5 and 30 DF, p-value: 0.01261
> anova(m1.3.2);
Analysis of Variance Table
Response: ACHIEVE
Df Sum Sq Mean Sq F value Pr(>F)
FACULTY 1 90.25 90.250 4.5012 0.04224 *
FACULTY:DEPART 4 263.22 65.806 3.2821 0.02404 *
Residuals 30 601.50 20.050
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
各学科の平均の推定値の計算は、やや面倒。学部はLとSがあり、それぞれa,b,c学科とd,e,f学科が含まれる。学部は、L学科がリファレンス、学科はa学科がリファレンスなので、L学部a学科の推定値は（Intercept）の値。標準誤差は、（Intercept）値で共通（学科の平均の推定値の標準誤差。学部については以下に示す方法で計算の必要あり）。b,c学科は、（Intercept）の値に、それぞれの係数の値を加えればよい。d,e,f学科は、S学部なのでFACULTYSの係数と各学科の係数の両方を加えた数値が、推定値となる。ただし、最後のカテゴリのf学科については、S学部に所属する学科のリファレンスとなるようで、FACULTYSの値を加えるだけでよい。
各学部の平均の推定値に関しては、Rでは直接出力はされない。切片の値は学部の平均の推定値ではなく、この場合はL学部a学科の推定値。学部の推定値を求めるには、学部だけで1元配置分散分析をして平均の推定値を求めて、標準誤差はネストしたモデルに基づいて計算する必要がある。
その他、各数値は整理されていないので値を拾うか計算が必要。平方和、平方平均は、分散分析表（anovaの結果）から拾える。モデル全体のR2乗、調整R2乗はlmの結果のsummary下部に表示。F値およびp値の検定に関する数値も、表示されたものを拾う。書籍にあるもので計算しなければならないのは、学部および学科個別の偏η2乗値で、分散分析表から誤差の平方和と学部および学科変数の平方和をもってきて、[学部の平方和/（学科の平方和＋誤差の平方和）]で、各変数の効果量を求めることができる。具体的には、学部の場合[90.25/(90.25+601.50)=0.130]となる。全体の平方和で除するのではなく、対象となる変数の平方和＋誤差の平方和で除するのがポイント。SPSSでは、自動的に計算されて出力される。
平均の推定値の標準誤差は、一般線形モデルでは先に書いたように誤差の平均平方を観測数で除した値となる。したがって、学部レベルの平均の推定値の標準誤差は、このデータ（学部、学科ともに全ての同水準カテゴリで観測数が同じバランスデータ）については、誤差の平均平方を各学部の観測数18で除した値の平方根となる。学科についても、計算手続きは同じ。カテゴリ間で観測数が異なるアンバランスデータに関しては、当然ながら標準誤差はカテゴリごとに異なることになる。
以上で、一通り書籍に示された値をRで得ることができた。なお、1節前のネスト無しの1変数の場合は、変数が1つなので偏η2乗はR2乗と一致するので計算の手間はない。
2元配置分散分析
説明変数が2つの分散分析。平方和のタイプなどが問題になってくる。ここでは、性別と学部によって達成度を比較。
まずは、性別のみの1元配置分散分析。
> m1.4.1 <- lm(ACHIEVE ~ GENDER, data=df1);
> summary(m1.4.1);
Call:
lm(formula = ACHIEVE ~ GENDER, data = df1)
Residuals:
Min 1Q Median 3Q Max
-8.8333 -4.4722 0.9722 3.0278 10.1667
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 83.222 1.238 67.236 <2e-16 ***
GENDERm -1.389 1.750 -0.793 0.433
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 5.251 on 34 degrees of freedom
Multiple R-squared: 0.01818, Adjusted R-squared: -0.0107
F-statistic: 0.6296 on 1 and 34 DF, p-value: 0.433
> anova(m1.4.1);
Analysis of Variance Table
Response: ACHIEVE
Df Sum Sq Mean Sq F value Pr(>F)
GENDER 1 17.36 17.361 0.6296 0.433
Residuals 34 937.61 27.577
単純な1元配置分散分析なので、先述の要領で値を読み取る。
次に、性別と学部を組み合わせて作成した性別×学部変数による1元配置分散分析。
> m1.4.2 <- lm(ACHIEVE ~ FAC_GEN, data=df1);
> summary(m1.4.2);
Call:
lm(formula = ACHIEVE ~ FAC_GEN, data = df1)
Residuals:
Min 1Q Median 3Q Max
-7.0000 -2.8056 -0.2222 1.7222 8.0000
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 85.000 1.237 68.716 < 2e-16 ***
FAC_GENLm -8.111 1.749 -4.637 5.69e-05 ***
FAC_GENSf -3.556 1.749 -2.033 0.0505 .
FAC_GENSm 1.778 1.749 1.016 0.3171
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3.711 on 32 degrees of freedom
Multiple R-squared: 0.5386, Adjusted R-squared: 0.4953
F-statistic: 12.45 on 3 and 32 DF, p-value: 1.468e-05
> anova(m1.4.2);
Analysis of Variance Table
Response: ACHIEVE
Df Sum Sq Mean Sq F value Pr(>F)
FAC_GEN 3 514.31 171.435 12.449 1.468e-05 ***
Residuals 32 440.67 13.771
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
続いて、性別、学部および両者の交互作用を含めた2元配置分散分析。モデルの式は、「*」で直接効果と全ての組み合わせでの交互作用を検討する。省略しなければ[ACHIEVE ~ FACULTY+GENDER+FACULTY:GENDER]になる。
> m1.4.3 <- lm(ACHIEVE ~ FACULTY*GENDER, data=df1);
> summary(m1.4.3);
Call:
lm(formula = ACHIEVE ~ FACULTY * GENDER, data = df1)
Residuals:
Min 1Q Median 3Q Max
-7.0000 -2.8056 -0.2222 1.7222 8.0000
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 85.000 1.237 68.716 < 2e-16 ***
FACULTYS -3.556 1.749 -2.033 0.0505 .
GENDERm -8.111 1.749 -4.637 5.69e-05 ***
FACULTYS:GENDERm 13.444 2.474 5.434 5.61e-06 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3.711 on 32 degrees of freedom
Multiple R-squared: 0.5386, Adjusted R-squared: 0.4953
F-statistic: 12.45 on 3 and 32 DF, p-value: 1.468e-05
> anova(m1.4.3);
Analysis of Variance Table
Response: ACHIEVE
Df Sum Sq Mean Sq F value Pr(>F)
FACULTY 1 90.25 90.25 6.5537 0.01539 *
GENDER 1 17.36 17.36 1.2607 0.26987
FACULTY:GENDER 1 406.69 406.69 29.5330 5.612e-06 ***
Residuals 32 440.67 13.77
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
各変数の平均の推定値の求め方は、学部、学科についてはそれぞれで1元配置分散分析をした値で、標準誤差はこれまでに見たように誤差の平均平方を観測数で除した値とする。学部×学科については、L学部女性が基準になっているので、L学部男性（基準に対して性別のみ違う）ならGENDERの係数を加え、S学部女性（基準に対して学部のみ違う）ならFACULTYの係数を加え、S学部男性（基準に対して、学部・性別とも異なる）なら、GENDER+FACULTY+GENDER:FACULTY（交互作用）を全て加える。
上記の性別×学部変数と、性別＋学部＋交互作用での2つの分析は、モデルとして説明される分散の総和は同じ（R2乗や平方和が一致）。つまり、性別×学部変数で説明された情報量を、性別の効果、学部の効果、および交互作用に分割して分析したのが、後者のモデル。要素に分割してみると、交互作用による効果が大きいことがわかる。
分散分析の基本的な考え方は、説明変数によって目的変数の分散をどの程度説明できるか、という点に尽きる。したがって、説明変数の投入方法を変更してみても、モデル全体として見れば説明される分散（情報量）は等価。この点は、マルチレベルモデルを理解する上で、大事だと思われる。すなわち、書籍で後に解説されていくことになるが、説明変数を固定効果として投入しても変量効果として投入しても、モデルとしてその説明変数が目的変数の分散を説明する割合は変わらない。なので、変量効果そのものに関心がなく、その変数の影響を統制したいということであれば、固定効果としてモデルに投入しても、概ね同じ分析結果が得られる（と思う）。
重（単）回帰分析
単回帰分析のパラメータの推定結果については、SPSSの出力（書籍の数値）との対応は容易に判断できるので、省略。
平方和のタイプ
重回帰分析（ないしは、2元以上の分散分析＝説明変数が複数の一般線形モデル）のモデルについて回帰分析を行う際には、平方和のタイプに注意をする必要がある。平方和のタイプには1,2,3があり、それぞれ計算方法が異なる（ただし説明変数が1つの場合には結果は同じ）。
- タイプ1平方和
- モデル式で与えられた順に変数を投入し、その都度増加する平方和。変数の投入順序が変わると、変数間に相関がある場合には平方和の値も変化する。これは、変数間で情報が共有されている部分（相関）があると、先に投入された変数の平方和に、相関（情報の重なり）部分が含まれるため、後に投入された変数の平方和が小さくなる。
- タイプ2平方和
- モデル式での変数の投入順序に関わらず、他の変数が全て投入された状態で最後に目的の変数が投入されたときの平方和を、それぞれの変数について計算する。したがって、投入順序に関係なく同じ値が得られる。また、計算方法から考えれば、タイプ2平方和は、他の変数が既に投入された状態で計算されているため、他の変数の影響が統制された値（他の変数との相関部分を除く）が得られているはず。
- タイプ3平方和
- タイプ2の計算と同様だが、交互作用も投入したうえで、目的の変数の平方和を計算する。
p.24の表1-10、表1-11の確認。達成度を興味と自主学習で説明するモデルで、変数の投入順序を入れ替えたとき、平方和のタイプがどう変化するか。
自主学習を先に投入
> m1.6.1a <- lm(ACHIEVE ~ STUDY + INTEREST, data=df1);
> anova(m1.6.1a);
Analysis of Variance Table
Response: ACHIEVE
Df Sum Sq Mean Sq F value Pr(>F)
STUDY 1 485.58 485.58 238.48 < 2.2e-16 ***
INTEREST 1 402.19 402.19 197.52 1.744e-15 ***
Residuals 33 67.19 2.04
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> Anova(m1.6.1a, type=2);
Anova Table (Type II tests)
Response: ACHIEVE
Sum Sq Df F value Pr(>F)
STUDY 388.36 1 190.73 2.864e-15 ***
INTEREST 402.19 1 197.52 1.744e-15 ***
Residuals 67.19 33
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
興味を先に投入
> m1.6.1b <- lm(ACHIEVE ~ INTEREST + STUDY, data=df1);
> anova(m1.6.1b);
Analysis of Variance Table
Response: ACHIEVE
Df Sum Sq Mean Sq F value Pr(>F)
INTEREST 1 499.42 499.42 245.27 < 2.2e-16 ***
STUDY 1 388.36 388.36 190.73 2.864e-15 ***
Residuals 33 67.19 2.04
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> Anova(m1.6.1b, type=2);
Anova Table (Type II tests)
Response: ACHIEVE
Sum Sq Df F value Pr(>F)
INTEREST 402.19 1 197.52 1.744e-15 ***
STUDY 388.36 1 190.73 2.864e-15 ***
Residuals 67.19 33
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
タイプ2の計算結果は、順序に関わらず同じ値だが、タイプ1では先に投入するよりも、後に投入した方が平方和の値が小さくなることが確認できる。