メモ:マルチレベルモデル
差分
このページの2つのバージョン間の差分を表示します。
両方とも前のリビジョン前のリビジョン次のリビジョン | 前のリビジョン | ||
メモ:マルチレベルモデル [2015/12/14 12:04] – [重(単)回帰分析] Wiki Editor | メモ:マルチレベルモデル [2017/02/11 09:57] (現在) – Wiki Editor | ||
---|---|---|---|
行 1: | 行 1: | ||
====== マルチレベルモデルに関するメモ ====== | ====== マルチレベルモデルに関するメモ ====== | ||
+ | |||
+ | いろいろ調べたりして、ちょっとまとめた資料を作成する機会があったので、なくさないように添付しておく。検証用に、サンプルデータもおいておく。サンプルデータは、意図的に生成したものであり、PDF中の設定も架空のもの。 | ||
+ | |||
+ | * {{: | ||
+ | * {{: | ||
+ | |||
[安藤正人(2011)『マルチレベルモデル入門 実習:継時データ分析』ナカニシヤ出版]に基づいて、マルチレベルモデルを勉強したメモ。 | [安藤正人(2011)『マルチレベルモデル入門 実習:継時データ分析』ナカニシヤ出版]に基づいて、マルチレベルモデルを勉強したメモ。 | ||
- | 本書は、SPSS/ | + | 本書は、SPSS/ |
===== R環境の準備 ===== | ===== R環境の準備 ===== | ||
- | 必要なライブラリをあらかじめインストールしておく。線形モデルを扱う関数(lm)は標準で入っているが、分散分析でType2/ | + | 必要なライブラリをあらかじめインストールしておく。線形モデルを扱う関数(lm)は標準で入っているが、分散分析でType2/ |
| | ||
install.packages(" | install.packages(" | ||
install.packages(" | install.packages(" | ||
| | ||
+ | # | ||
+ | install.packages(" | ||
+ | # | ||
+ | install.packages(" | ||
+ | | ||
library(car); | library(car); | ||
library(lme4); | library(lme4); | ||
+ | library(lmerTest); | ||
+ | library(pbkrtest); | ||
+ | |||
+ | |||
サンプルデータも読み込んでおく。 | サンプルデータも読み込んでおく。 | ||
行 17: | 行 32: | ||
df1 <- read.csv(" | df1 <- read.csv(" | ||
- | ===== 一般線形モデル | + | - [[メモ: |
- | 分散分析、共分散分析、(重)回帰分析を一般線形モデルで統一的に把握する。 | + | - [[メモ:マルチレベルモデル:混合モデル|混合モデル(マルチレベルモデル)]] |
- | + | ||
- | ==== 1元配置分散分析 ==== | + | |
- | lm関数で学科間の平均得点について分散分析する。 | + | |
- | + | ||
- | >m1.2.1 <- lm(ACHIEVE ~ DEPART, data=df1); | + | |
- | > summary(m1.2.1); | + | |
- | + | ||
- | Call: | + | |
- | lm(formula = ACHIEVE ~ DEPART, data = df1) | + | |
- | + | ||
- | Residuals: | + | |
- | | + | |
- | -7.000 -3.542 | + | |
- | + | ||
- | Coefficients: | + | |
- | Estimate Std. Error t value Pr(>|t|) | + | |
- | (Intercept) | + | |
- | DEPARTb | + | |
- | DEPARTc | + | |
- | DEPARTd | + | |
- | DEPARTe | + | |
- | DEPARTf | + | |
- | --- | + | |
- | Signif. codes: | + | |
- | + | ||
- | Residual standard error: 4.478 on 30 degrees of freedom | + | |
- | Multiple R-squared: | + | |
- | F-statistic: | + | |
- | + | ||
- | > anova(m1.2.1); | + | |
- | + | ||
- | Analysis of Variance Table | + | |
- | + | ||
- | Response: ACHIEVE | + | |
- | Df Sum Sq Mean Sq F value Pr(>F) | + | |
- | DEPART | + | |
- | Residuals 30 601.50 | + | |
- | --- | + | |
- | Signif. codes: | + | |
- | + | ||
- | 学科(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の各出力と書籍の内容は一致する。 | + | |
- | + | ||
- | ==== ネストした分散分析 ==== | + | |
- | 下位カテゴリが上位カテゴリにもれなく含まれるような場合の分散分析。上位カテゴリで説明される分散と下位カテゴリで説明される分散を分けてみる。両者の合計は、上位カテゴリだけで分析した場合に説明される分散の大きさと一致する。 | + | |
- | + | ||
- | ここでは、学科は学部にネストされているデータを用いる。まずは、学部だけの分散分析の結果。 | + | |
- | + | ||
- | | + | |
- | > summary(m1.3.1); | + | |
- | + | ||
- | Call: | + | |
- | lm(formula = ACHIEVE ~ FACULTY, data = df1) | + | |
- | + | ||
- | Residuals: | + | |
- | Min 1Q Median | + | |
- | -7.9444 -4.3194 -0.1111 | + | |
- | + | ||
- | Coefficients: | + | |
- | Estimate Std. Error t value Pr(> | + | |
- | (Intercept) | + | |
- | FACULTYS | + | |
- | --- | + | |
- | Signif. codes: | + | |
- | + | ||
- | Residual standard error: 5.043 on 34 degrees of freedom | + | |
- | Multiple R-squared: | + | |
- | F-statistic: | + | |
- | + | ||
- | > anova(m1.3.1); | + | |
- | Analysis of Variance Table | + | |
- | + | ||
- | Response: ACHIEVE | + | |
- | Df Sum Sq Mean Sq F value Pr(>F) | + | |
- | FACULTY | + | |
- | Residuals 34 864.72 | + | |
- | --- | + | |
- | Signif. codes: | + | |
- | + | ||
- | 書籍と一致。 | + | |
- | + | ||
- | 続いて、ネストした分散分析。 | + | |
- | > m1.3.2 <- lm(ACHIEVE ~ FACULTY/ | + | |
- | > summary(m1.3.2); | + | |
- | + | ||
- | Call: | + | |
- | lm(formula = ACHIEVE ~ FACULTY/ | + | |
- | + | ||
- | Residuals: | + | |
- | | + | |
- | -7.000 -3.542 | + | |
- | + | ||
- | Coefficients: | + | |
- | | + | |
- | (Intercept) | + | |
- | FACULTYS | + | |
- | FACULTYL: | + | |
- | FACULTYS: | + | |
- | FACULTYL: | + | |
- | FACULTYS: | + | |
- | FACULTYL: | + | |
- | FACULTYS: | + | |
- | FACULTYL: | + | |
- | FACULTYS: | + | |
- | FACULTYL: | + | |
- | FACULTYS: | + | |
- | --- | + | |
- | Signif. codes: | + | |
- | + | ||
- | Residual standard error: 4.478 on 30 degrees of freedom | + | |
- | Multiple R-squared: | + | |
- | F-statistic: | + | |
- | + | ||
- | > anova(m1.3.2); | + | |
- | Analysis of Variance Table | + | |
- | + | ||
- | Response: ACHIEVE | + | |
- | Df Sum Sq Mean Sq F value Pr(>F) | + | |
- | FACULTY | + | |
- | FACULTY: | + | |
- | Residuals | + | |
- | --- | + | |
- | Signif. codes: | + | |
- | + | ||
- | 各学科の平均の推定値の計算は、やや面倒。学部はLとSがあり、それぞれa, | + | |
- | + | ||
- | 各学部の平均の推定値に関しては、Rでは直接出力はされない。切片の値は学部の平均の推定値ではなく、この場合はL学部a学科の推定値。学部の推定値を求めるには、学部だけで1元配置分散分析をして平均の推定値を求めて、標準誤差はネストしたモデルに基づいて計算する必要がある。 | + | |
- | + | ||
- | その他、各数値は整理されていないので値を拾うか計算が必要。平方和、平方平均は、分散分析表(anovaの結果)から拾える。モデル全体のR2乗、調整R2乗はlmの結果のsummary下部に表示。F値およびp値の検定に関する数値も、表示されたものを拾う。書籍にあるもので計算しなければならないのは、学部および学科個別の偏η2乗値で、分散分析表から誤差の平方和と学部および学科変数の平方和をもってきて、[学部の平方和/ | + | |
- | + | ||
- | 平均の推定値の標準誤差は、一般線形モデルでは先に書いたように誤差の平均平方を観測数で除した値となる。したがって、学部レベルの平均の推定値の標準誤差は、このデータ(学部、学科ともに全ての同水準カテゴリで観測数が同じバランスデータ)については、誤差の平均平方を各学部の観測数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 | + | |
- | -8.8333 -4.4722 | + | |
- | + | ||
- | Coefficients: | + | |
- | Estimate Std. Error t value Pr(> | + | |
- | (Intercept) | + | |
- | GENDERm | + | |
- | --- | + | |
- | Signif. codes: | + | |
- | + | ||
- | Residual standard error: 5.251 on 34 degrees of freedom | + | |
- | Multiple R-squared: | + | |
- | F-statistic: | + | |
- | + | ||
- | > anova(m1.4.1); | + | |
- | Analysis of Variance Table | + | |
- | + | ||
- | Response: ACHIEVE | + | |
- | Df Sum Sq Mean Sq F value Pr(>F) | + | |
- | GENDER | + | |
- | Residuals 34 937.61 | + | |
- | + | ||
- | 単純な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 | + | |
- | -7.0000 -2.8056 -0.2222 | + | |
- | + | ||
- | Coefficients: | + | |
- | Estimate Std. Error t value Pr(> | + | |
- | (Intercept) | + | |
- | FAC_GENLm | + | |
- | FAC_GENSf | + | |
- | FAC_GENSm | + | |
- | --- | + | |
- | Signif. codes: | + | |
- | + | ||
- | Residual standard error: 3.711 on 32 degrees of freedom | + | |
- | Multiple R-squared: | + | |
- | F-statistic: | + | |
- | + | ||
- | > anova(m1.4.2); | + | |
- | Analysis of Variance Table | + | |
- | + | ||
- | Response: ACHIEVE | + | |
- | Df Sum Sq Mean Sq F value Pr(>F) | + | |
- | FAC_GEN | + | |
- | Residuals 32 440.67 | + | |
- | --- | + | |
- | Signif. codes: | + | |
- | + | ||
- | 続いて、性別、学部および両者の交互作用を含めた2元配置分散分析。モデルの式は、「*」で直接効果と全ての組み合わせでの交互作用を検討する。省略しなければ[ACHIEVE ~ FACULTY+GENDER+FACULTY: | + | |
- | + | ||
- | > m1.4.3 <- lm(ACHIEVE ~ FACULTY*GENDER, | + | |
- | > summary(m1.4.3); | + | |
- | + | ||
- | Call: | + | |
- | lm(formula = ACHIEVE ~ FACULTY * GENDER, data = df1) | + | |
- | + | ||
- | Residuals: | + | |
- | Min 1Q Median | + | |
- | -7.0000 -2.8056 -0.2222 | + | |
- | + | ||
- | Coefficients: | + | |
- | | + | |
- | (Intercept) | + | |
- | FACULTYS | + | |
- | GENDERm | + | |
- | FACULTYS: | + | |
- | --- | + | |
- | Signif. codes: | + | |
- | + | ||
- | Residual standard error: 3.711 on 32 degrees of freedom | + | |
- | Multiple R-squared: | + | |
- | F-statistic: | + | |
- | + | ||
- | > anova(m1.4.3); | + | |
- | Analysis of Variance Table | + | |
- | + | ||
- | Response: ACHIEVE | + | |
- | Df Sum Sq Mean Sq F value Pr(>F) | + | |
- | FACULTY | + | |
- | GENDER | + | |
- | FACULTY: | + | |
- | Residuals | + | |
- | --- | + | |
- | Signif. codes: | + | |
- | + | ||
- | 各変数の平均の推定値の求め方は、学部、学科についてはそれぞれで1元配置分散分析をした値で、標準誤差はこれまでに見たように誤差の平均平方を観測数で除した値とする。学部×学科については、L学部女性が基準になっているので、L学部男性(基準に対して性別のみ違う)ならGENDERの係数を加え、S学部女性(基準に対して学部のみ違う)ならFACULTYの係数を加え、S学部男性(基準に対して、学部・性別とも異なる)なら、GENDER+FACULTY+GENDER: | + | |
- | + | ||
- | 上記の性別×学部変数と、性別+学部+交互作用での2つの分析は、モデルとして説明される分散の総和は同じ(R2乗や平方和が一致)。つまり、性別×学部変数で説明された情報量を、性別の効果、学部の効果、および交互作用に分割して分析したのが、後者のモデル。要素に分割してみると、交互作用による効果が大きいことがわかる。 | + | |
- | + | ||
- | 分散分析の基本的な考え方は、説明変数によって目的変数の分散をどの程度説明できるか、という点に尽きる。したがって、説明変数の投入方法を変更してみても、モデル全体として見れば説明される分散(情報量)は等価。この点は、マルチレベルモデルを理解する上で、大事だと思われる。すなわち、書籍で後に解説されていくことになるが、説明変数を固定効果として投入しても変量効果として投入しても、モデルとしてその説明変数が目的変数の分散を説明する割合は変わらない。なので、変量効果そのものに関心がなく、その変数の影響を統制したいということであれば、固定効果としてモデルに投入しても、概ね同じ分析結果が得られる(と思う)。 | + | |
- | + | ||
- | + | ||
- | ==== 重(単)回帰分析 ==== | + | |
- | 単回帰分析のパラメータの推定結果については、SPSSの出力(書籍の数値)との対応は容易に判断できるので、省略。 | + | |
- | + | ||
- | === 平方和のタイプ === | + | |
- | 重回帰分析(ないしは、2元以上の分散分析=説明変数が複数の一般線形モデル)のモデルについて回帰分析を行う際には、平方和のタイプに注意をする必要がある。平方和のタイプには1, | + | |
- | + | ||
- | * タイプ1平方和 | + | |
- | * モデル式で与えられた順に変数を投入し、その都度増加する平方和。変数の投入順序が変わると、変数間に相関がある場合には平方和の値も変化する。これは、変数間で情報が共有されている部分(相関)があると、先に投入された変数の平方和に、相関(情報の重なり)部分が含まれるため、後に投入された変数の平方和が小さくなる。 | + | |
- | * タイプ2平方和 | + | |
- | * モデル式での変数の投入順序に関わらず、他の変数が全て投入された状態で最後に目的の変数が投入されたときの平方和を、それぞれの変数について計算する。したがって、投入順序に関係なく同じ値が得られる。また、計算方法から考えれば、タイプ2平方和は、他の変数が既に投入された状態で計算されているため、他の変数の影響が統制された値(他の変数との相関部分を除く)が得られているはず。 | + | |
- | * タイプ3平方和 | + | |
- | * タイプ2の計算と同様だが、交互作用も投入したうえで、目的の変数の平方和を計算する。 | + | |
- | * RのcarパッケージAnova関数でタイプ3平方和を計算する際には、オプションを指定しないとSPSSなどの計算結果と一致しない(後述)。 | + | |
- | + | ||
- | p.24の表1-10、表1-11の確認。達成度を興味と自主学習で説明するモデルで、変数の投入順序を入れ替えたとき、平方和のタイプがどう変化するか。 | + | |
- | + | ||
- | 自主学習を先に投入 | + | |
- | > m1.6.2a <- lm(ACHIEVE ~ STUDY + INTEREST, data=df1); | + | |
- | > anova(m1.6.2a); | + | |
- | Analysis of Variance Table | + | |
- | + | ||
- | Response: ACHIEVE | + | |
- | Df Sum Sq Mean Sq F value Pr(>F) | + | |
- | STUDY 1 485.58 | + | |
- | INTEREST | + | |
- | Residuals 33 67.19 2.04 | + | |
- | --- | + | |
- | Signif. codes: | + | |
- | > Anova(m1.6.2a, | + | |
- | Anova Table (Type II tests) | + | |
- | + | ||
- | Response: ACHIEVE | + | |
- | Sum Sq Df F value Pr(>F) | + | |
- | STUDY | + | |
- | INTEREST | + | |
- | Residuals | + | |
- | --- | + | |
- | Signif. codes: | + | |
- | + | ||
- | 興味を先に投入 | + | |
- | > m1.6.2b <- lm(ACHIEVE ~ INTEREST + STUDY, data=df1); | + | |
- | > anova(m1.6.2b); | + | |
- | Analysis of Variance Table | + | |
- | + | ||
- | Response: ACHIEVE | + | |
- | Df Sum Sq Mean Sq F value Pr(>F) | + | |
- | INTEREST | + | |
- | STUDY 1 388.36 | + | |
- | Residuals 33 67.19 2.04 | + | |
- | --- | + | |
- | Signif. codes: | + | |
- | > Anova(m1.6.2b, | + | |
- | Anova Table (Type II tests) | + | |
- | + | ||
- | Response: ACHIEVE | + | |
- | Sum Sq Df F value Pr(>F) | + | |
- | INTEREST | + | |
- | STUDY | + | |
- | Residuals | + | |
- | --- | + | |
- | Signif. codes: | + | |
- | + | ||
- | タイプ2の計算結果は、順序に関わらず同じ値だが、タイプ1では先に投入するよりも、後に投入した方が平方和の値が小さくなることが確認できる。書籍で示されている偏η2乗は、タイプ2平方和を使って計算したもののようだ。 | + | |
- | + | ||
- | + | ||
- | == RのcarパッケージAnova()関数でのタイプ3平方和 == | + | |
- | Rでは、factor型の変数は自動的にダミー変数に変換してくれるが、その変換方法にはいくつかパターンがある。その関係で、単純にAnova関数でtype3平方和を計算すると、おかしな値となる。そこで、以下のオプションをあらかじめ指定しておく。 | + | |
- | + | ||
- | options(contrasts = c(" | + | |
- | + | ||
メモ/マルチレベルモデル.txt · 最終更新: 2017/02/11 09:57 by Wiki Editor