分散分析、共分散分析、（重）回帰分析を一般線形モデルで統一的に把握する。
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つの分散分析。平方和のタイプなどが問題になってくる。ここでは、性別と学部によって達成度を比較。
まずは、性別のみの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つの場合には結果は同じ）。
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 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.2a, 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.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 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.2b, 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では先に投入するよりも、後に投入した方が平方和の値が小さくなることが確認できる。書籍で示されている偏η2乗は、タイプ2平方和を使って計算したもののようだ。
Rでは、factor型の変数は自動的にダミー変数に変換してくれるが、その変換方法にはいくつかパターンがある。その関係で、単純にcarパッケージのAnova関数でtype3平方和を計算すると、おかしな値となる。そこで、以下のオプションをあらかじめ指定しておく。
options(contrasts = c("contr.sum", "contr.sum"));
optionsは、Rの基本的な挙動を決定する値を管理するもののようである。その中のcontrastsが、カテゴリカル変数をダミー変数に変換する際の方式を決定している。デフォルトの値は、以下の通り。
> options()$contrasts
unordered ordered
"contr.treatment" "contr.poly"
この値に対して、上記の設定では「contr.sum」という方式に変更している。2つ設定しているのは、デフォルト値からもわかるように前者が名義変数で、後者は順序カテゴリ変数の処理方法を指定している。
名義変数のデフォルト値「contr.treatment」は、一番シンプルな変換で、factorのlevelesの最初の水準を基準（0）として、それ以外のカテゴリは対応するダミー変数にそれぞれ1を割り当てる方式である。
順序カテゴリ変数のデフォルト値「contr.poly」は、順序カテゴリデータに用いられるもので、カテゴリの順序に応じて数値が割り当てられる。contr.poly()関数にカテゴリ数を与えると、どのようにデザイン行列がつくられるのかがわかる。1列目から順に、1次、2次、3次の関数になっていて、列方向の合計は0になる。どのような効果があるのかは、今度よく考えてみることにして、とりあえず順序カテゴリ変数向けの方式だということだけ把握しておこう。
タイプ3平方和を計算する際に指定する「contr.sum」は、「contr.treatment」と似ているが、ダミー変数に変換したときの各変数の列方向の合計が0になるようになっている。タイプ3平方和を計算する際には、この条件が必要なようだ。
> contr.sum(5); [,1] [,2] [,3] [,4] 1 1 0 0 0 2 0 1 0 0 3 0 0 1 0 4 0 0 0 1 5 -1 -1 -1 -1
カテゴリカルデータからダミー変数を生成する方式に関して、optionsで指定しておくと以降はすべての分析に影響する。値の推定においては、どの方式でも同じ結果を得ることができるが、切片（Intercept）と係数を組み合わせて推定値を計算する際にややこしくなるので、タイプ3平方和を計算するとき以外は、デフォルトのcontr.treatmentに戻しておくのが無難かもしれない。なお、lm関数のオプションでも個別に方式を指定することができる。
lm(式, data=データフレーム, contrasts=list(カテゴリ変数名="contr.treatment"))
contrastsオプションに、方式を指定するカテゴリ変数名と方式名をセットにしたリストを渡してやればよい。
データフレームの属性に、あらかじめデザイン行列を埋め込むこともできるようだ。この場合、contr.treatment関数、contr.sumなど目的の方式のデザイン行列を作成し、contrasts関数で目的の変数に埋め込む。この方法だと、baseオプションで基準とするカテゴリを指定できる。
contrasts(df1$VAL1)<-contr.treatment(6, base=5); > contrasts(df1$VAL1) 1 2 3 4 6 a 1 0 0 0 0 b 0 1 0 0 0 c 0 0 1 0 0 d 0 0 0 1 0 e 0 0 0 0 0 f 0 0 0 0 1
本線に戻って、書籍の分析結果を追いかける。分散分析の平方和のタイプ以外はシンプルな回帰分析なので、Rでの分析でも結果は一致する。Rでちょっと引っかかるのは、標準化した分析とセンタリングした分析の方法か（使える人にとっては特に問題ないのだと思うが）。
標準化（p.25, モデル1.6.2Z）した回帰分析。自主学習＋興味のモデルについて、変数を標準化して分析する。
> m1.6.2Z <- lm(scale(ACHIEVE) ~ scale(STUDY) + scale(INTEREST), data=df1);
> summary(m1.6.2Z);
Call:
lm(formula = scale(ACHIEVE) ~ scale(STUDY) + scale(INTEREST),
data = df1)
Residuals:
Min 1Q Median 3Q Max
-0.44360 -0.17371 -0.04769 0.19323 0.70300
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.049e-15 4.553e-02 0.00 1
scale(STUDY) 6.416e-01 4.646e-02 13.81 2.86e-15 ***
scale(INTEREST) 6.529e-01 4.646e-02 14.05 1.74e-15 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.2732 on 33 degrees of freedom
Multiple R-squared: 0.9296, Adjusted R-squared: 0.9254
F-statistic: 218 on 2 and 33 DF, p-value: < 2.2e-16
> Anova(m1.6.2Z, type=3);
Anova Table (Type III tests)
Response: scale(ACHIEVE)
Sum Sq Df F value Pr(>F)
(Intercept) 0.0000 1 0.00 1
scale(STUDY) 14.2335 1 190.73 2.864e-15 ***
scale(INTEREST) 14.7405 1 197.52 1.744e-15 ***
Residuals 2.4627 33
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
scale関数を用いて、変数を標準化してから分析する。書籍では標準化した変数が用意されているが、関数で標準化してみた。得られた結果は、書籍と一致（指数表記のためちょっと見にくいが）。分散分析は、optionsを指定のうえで、書籍の注記にしたがってtype3をあえて指定。type2でも結果は同じ。
ちなみに、この方法で変数を標準化すると、欠損値の扱いの関係でSPSSなどと結果が一致しないことがある。その場合は、まず変数そのままでlm関数に投入し、そこから得られたオブジェクトの「オブジェクト名$model」に分析に利用されたデータセットが入っているので、これをscaleで標準化して再分析すると、SPSSの標準化係数（ベータ）に一致する値が得られる。ちょっと面倒だが。
続いて、モデル1.7.1および1.7.1C。自主学習と興味の交互作用を投入したモデル。
> m1.7.1 <- lm(ACHIEVE ~ STUDY + INTEREST + STUDY:INTEREST, data=df1);
> summary(m1.7.1);
Call:
lm(formula = ACHIEVE ~ STUDY + INTEREST + STUDY:INTEREST, data = df1)
Residuals:
Min 1Q Median 3Q Max
-2.4284 -0.8941 -0.1191 1.0809 2.8352
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 67.1227 7.4203 9.046 2.48e-10 ***
STUDY -0.7472 0.9497 -0.787 0.43720
INTEREST -0.2873 1.0352 -0.278 0.78315
STUDY:INTEREST 0.3871 0.1320 2.932 0.00617 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.287 on 32 degrees of freedom
Multiple R-squared: 0.9445, Adjusted R-squared: 0.9393
F-statistic: 181.7 on 3 and 32 DF, p-value: < 2.2e-16
> Anova(m1.7.1, type=3);
Anova Table (Type III tests)
Response: ACHIEVE
Sum Sq Df F value Pr(>F)
(Intercept) 135.435 1 81.8265 2.485e-10 ***
STUDY 1.025 1 0.6190 0.437195
INTEREST 0.127 1 0.0770 0.783150
STUDY:INTEREST 14.229 1 8.5969 0.006174 **
Residuals 52.965 32
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
交互作用項が入っているので、タイプ2平方和とは値が異なる。optionsを指定すると、書籍で示されるタイプ3平方和と一致。
続いて、モデル1.7.1C。興味と自主学習の両変数を平均で中心化して分析する。中心化は、scale関数で可能であるが、オプションを指定しないと標準化されるので、「scale(変数, scale = F)」として、中心化のみ行うようにすればよい。上のように、モデル式中にscale関数を書いても良いが、交互作用項のかき方が煩雑になるのでセンタリング変数を作ってから分析した方が、スマートかも。ここでは、とりあえず無理やりモデル式中でセンタリングする。
> m1.7.1C <- lm(ACHIEVE ~ scale(STUDY, scale=F) + scale(INTEREST, scale=F) + scale(STUDY, scale=F):scale(INTEREST, scale=F), data=df1);
> summary(m1.7.1C);
Call:
lm(formula = ACHIEVE ~ scale(STUDY, scale = F) + scale(INTEREST,
scale = F) + scale(STUDY, scale = F):scale(INTEREST, scale = F),
data = df1)
Residuals:
Min 1Q Median 3Q Max
-2.4284 -0.8941 -0.1191 1.0809 2.8352
Coefficients:
Estimate Std. Error t value
(Intercept) 82.4412 0.2164 380.885
scale(STUDY, scale = F) 2.2741 0.1591 14.298
scale(INTEREST, scale = F) 2.7018 0.1735 15.570
scale(STUDY, scale = F):scale(INTEREST, scale = F) 0.3871 0.1320 2.932
Pr(>|t|)
(Intercept) < 2e-16 ***
scale(STUDY, scale = F) 1.9e-15 ***
scale(INTEREST, scale = F) < 2e-16 ***
scale(STUDY, scale = F):scale(INTEREST, scale = F) 0.00617 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.287 on 32 degrees of freedom
Multiple R-squared: 0.9445, Adjusted R-squared: 0.9393
F-statistic: 181.7 on 3 and 32 DF, p-value: < 2.2e-16
> Anova(m1.7.1C, type=3);
Anova Table (Type III tests)
Response: ACHIEVE
Sum Sq Df F value
(Intercept) 240118 1 1.4507e+05
scale(STUDY, scale = F) 338 1 2.0442e+02
scale(INTEREST, scale = F) 401 1 2.4243e+02
scale(STUDY, scale = F):scale(INTEREST, scale = F) 14 1 8.5969e+00
Residuals 53 32
Pr(>F)
(Intercept) < 2.2e-16 ***
scale(STUDY, scale = F) 1.901e-15 ***
scale(INTEREST, scale = F) < 2.2e-16 ***
scale(STUDY, scale = F):scale(INTEREST, scale = F) 0.006174 **
Residuals
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
無理にセンタリングを式中の関数で行っているので、出力の各項のラベルが長くなってしまい、見づらい。しかし、結果は書籍に一致しているので、良しとする。
2次以上の項を含む回帰分析。
モデル1.8.1の分析例。
> m1.8.1 <- lm(ACHIEVE ~ STUDY + I(STUDY^2), df1);
> summary(m1.8.1);
Call:
lm(formula = ACHIEVE ~ STUDY + I(STUDY^2), data = df1)
Residuals:
Min 1Q Median 3Q Max
-6.2599 -2.2677 -0.2599 2.3280 7.7479
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 42.0018 11.1533 3.766 0.000651 ***
STUDY 9.0956 3.2013 2.841 0.007645 **
I(STUDY^2) -0.4767 0.2210 -2.157 0.038417 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3.531 on 33 degrees of freedom
Multiple R-squared: 0.5692, Adjusted R-squared: 0.5431
F-statistic: 21.8 on 2 and 33 DF, p-value: 9.238e-07
> Anova(m1.8.1, type=2);
Anova Table (Type II tests)
Response: ACHIEVE
Sum Sq Df F value Pr(>F)
STUDY 100.64 1 8.0724 0.007645 **
I(STUDY^2) 57.98 1 4.6510 0.038417 *
Residuals 411.40 33
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
二次以上の項を投入する場合は、I()関数を使う必要がある。ここでは、STUDYの2次項をモデルに入れるためにI(STUDY^2)としている。
1.8.1のデータをセンタリングして分析。データの加工方法がスマートではないかもしれないが（いずれ調査しよう）、センタリングした変数を作成したうえで分析する。
> C_STUDY <- scale(df1$STUDY, scale=F);
> m1.8.1C <- lm(df1$ACHIEVE ~ C_STUDY + I(C_STUDY^2));
> summary(m1.8.1C);
Call:
lm(formula = df1$ACHIEVE ~ C_STUDY + I(C_STUDY^2))
Residuals:
Min 1Q Median 3Q Max
-6.2599 -2.2677 -0.2599 2.3280 7.7479
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 83.8151 0.8382 99.993 < 2e-16 ***
C_STUDY 1.7337 0.4269 4.061 0.000283 ***
I(C_STUDY^2) -0.4767 0.2210 -2.157 0.038417 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3.531 on 33 degrees of freedom
Multiple R-squared: 0.5692, Adjusted R-squared: 0.5431
F-statistic: 21.8 on 2 and 33 DF, p-value: 9.238e-07
> Anova(m1.8.1C, type=2);
Anova Table (Type II tests)
Response: df1$ACHIEVE
Sum Sq Df F value Pr(>F)
C_STUDY 205.64 1 16.495 0.0002825 ***
I(C_STUDY^2) 57.98 1 4.651 0.0384173 *
Residuals 411.40 33
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
まず、scale関数でSTUDYをセンタリングして変数に代入する。lmのモデル式は、ベクトルそのものを指定しているので、データフレームを指定しなくても計算してくれる。
変数を対数変換したうえで回帰する方法。こちらも変数変換の手順だけ確認。
> m1.9.1 <- lm(log(ACHIEVE) ~ log(STUDY) + log(INTEREST), df1);
> summary(m1.9.1);
Call:
lm(formula = log(ACHIEVE) ~ log(STUDY) + log(INTEREST), data = df1)
Residuals:
Min 1Q Median 3Q Max
-0.028015 -0.015089 -0.003859 0.012266 0.046677
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 3.57837 0.04553 78.59 < 2e-16 ***
log(STUDY) 0.16991 0.01345 12.63 3.45e-14 ***
log(INTEREST) 0.23991 0.02026 11.84 2.01e-13 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.019 on 33 degrees of freedom
Multiple R-squared: 0.9164, Adjusted R-squared: 0.9113
F-statistic: 180.9 on 2 and 33 DF, p-value: < 2.2e-16
> Anova(m1.9.1);
Anova Table (Type II tests)
Response: log(ACHIEVE)
Sum Sq Df F value Pr(>F)
log(STUDY) 0.057605 1 159.58 3.445e-14 ***
log(INTEREST) 0.050605 1 140.19 2.007e-13 ***
Residuals 0.011913 33
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
log関数を使って、各変数を対数変換するだけ。log関数に単純に変数を渡すと、自然対数になる。
名義変数と連続変数を同時に扱うモデル。共分散分析と呼ばれているものであるが、一般線形モデルの枠組みではカテゴリカル変数をダミーデータとして投入した重回帰分析と同じ。分散分析の検定力を上げるために連続尺度の共変量を考慮すると考えるのか、カテゴリカル変数と連続変数を合わせて、目的変数の値を予測するのか、という分析視点の違いと考えられる。本質的に同じ分析なので、目的に応じて分析の見方を変えればよいと思う。
ここでは、学部と自主学習で達成度を予想するモデル。
> m1.10.1 <- lm(ACHIEVE ~ FACULTY + STUDY, df1);
> summary(m1.10.1);
Call:
lm(formula = ACHIEVE ~ FACULTY + STUDY, data = df1)
Residuals:
Min 1Q Median 3Q Max
-6.0532 -2.5351 -0.9895 2.2294 7.7917
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 64.3910 3.2036 20.100 < 2e-16 ***
FACULTYS -1.1399 1.4814 -0.769 0.447
STUDY 2.4225 0.4507 5.374 6.11e-06 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3.738 on 33 degrees of freedom
Multiple R-squared: 0.5171, Adjusted R-squared: 0.4879
F-statistic: 17.67 on 2 and 33 DF, p-value: 6.067e-06
> anova(m1.10.1);
Analysis of Variance Table
Response: ACHIEVE
Df Sum Sq Mean Sq F value Pr(>F)
FACULTY 1 90.25 90.25 6.4588 0.01592 *
STUDY 1 403.61 403.61 28.8845 6.109e-06 ***
Residuals 33 461.11 13.97
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
> Anova(m1.10.1, type=3);
Anova Table (Type III tests)
Response: ACHIEVE
Sum Sq Df F value Pr(>F)
(Intercept) 5645.1 1 403.9986 < 2.2e-16 ***
FACULTY 8.3 1 0.5921 0.4471
STUDY 403.6 1 28.8845 6.109e-06 ***
Residuals 461.1 33
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
分散分析の結果については、Rで出力した結果と書籍に示される結果は一致する。しかし、いくつかわからない点もある。
まず、p.35の表1-21について。学部（FACULTY)の平均の推定値については、summaryに示される（Intercept）がリファレンスであるL学部で、この値にFACULTYSの係数を加えたものS学部の平均の推定値となる。ここまではOK。よくわからないのは、標準誤差の値。本書のp.9では、誤差の平均平方を観測数で除した値の平方根とされているが、ここではそのように計算した値とは異なる。ただし、Rの出力の（Intercept）の標準誤差の値と同一の値が示されている。これはSPSSの出力でも同様。また、L学部とS学部の標準誤差も、分散分析の場合とは異なり両者で違う値が示されている。これは、カテゴリカルな変数だけではなく連続変数を同時に扱うため、回帰直線の傾きがあるためだと思う（もう少し数学的理解が必要。回帰式の切片の標準誤差の算出方法を調べるとなんとなくわかる気がするが、完全には理解できていない。それぞれの回帰直線のX,Y変数の平均が異なるため、計算上は標準誤差が異なることになるのだとおもうのだが、、）。他方の学部の切片の標準誤差を求めるには、リファレンスとなる学部を入れ替えて再計算すればよい。具体的には、factore(df$FACULTY, levels=c(“S”,“L”))などとして、S学部をリファレンスにする。基準となるカテゴリをスマートに指定するオプションがあるような気もする。
書籍で示される数値を求める方法に終始したが、切片の誤差は、summaryの出力（およびp.35表1-21）では切片が0という帰無仮説を検定するだけで、特に意味があるようには思えない。問題は、L学部とS学部の間で切片に有意な差が認められるか、という点。これはp.37の表1-24で検定されている。これは、lmのsummaryの（Intercept）の一行下のFACULTYSの推定値、標準誤差、t値、p値に一致する。したがって、学部間の切片に差があるかを検定している。この検定は、タイプ3平方和による分散分析のFACULTYのp値と完全に一致しており、表1-24のt値を2乗した値が分散分析のF値となっている。つまり、両者の検定は等価であり、分散分析で学部の効果として検討されているのは、切片の違いということである。
共分散分析と一般線形モデルとの違いは、平方和のタイプの違いともいえる。タイプ1平方和は、共分散分析に対応する。すなわち、目的のカテゴリカル変数の検定力を上げるために、連続変数を考慮する。なので、最初に検定の目的であるカテゴリカル変数を投入し、同時に連続変数を投入することでモデル全体としての誤差の平方平均を小さくして検定力を向上させる。ここで、目的は最初に投入したカテゴリカル変数の検定なので、用いるのはタイプ1の平方和である。投入順序も、当然ながら目的のカテゴリカル変数を第1番目にしなければ意味がないだろう。タイプ3平方和を用いると、他の投入変数の影響を全て取り除いたうえで計算されるため、他の変数との兼ね合いでその変数の効果が検定される。
本書の例でいえば、共分散分析（タイプ1平方和）では自主学習の効果を考慮すると学部の効果は有意だが、一般線形モデル（タイプ3平方和）ならば有意ではない。これは、前者の共分散分析では、自主学習の多さ（平均的な水準）もひっくるめて学部の達成度の平均を検定しており、後者は学部の違いのみの効果を検定しているために有意にならないのだと思う。学部のみの違いによる達成度の差＋自主学習との相関部分（学部の違いによって自主学習が増える分）を含めた達成度の違いの検定がタイプ1平方和、タイプ3平方和では他の変数で説明される部分はすべて取り除いたうえでの平方和なので、純粋な学部の違いのみの効果を検定している。
前節のように、連続変数とカテゴリカル変数をモデルに同時に投入すると、水準ごとに切片が異なるモデルで分析できるが、水準ごとに傾きが常に同じとは限らない。そこで、水準ごとに傾きも変化するモデルを考えるために、両者の交互作用項を導入する。
p.39の表1-25、表1-26、表1-27の分析例。まず表1-25の数値を拾う。
> options(contrasts = c("contr.treatment", "contr.poly"));
> m1.11.1 <- lm(ACHIEVE ~ FACULTY + STUDY + FACULTY:STUDY, df1);
> summary(m1.11.1);
Call:
lm(formula = ACHIEVE ~ FACULTY + STUDY + FACULTY:STUDY, data = df1)
Residuals:
Min 1Q Median 3Q Max
-5.537 -2.632 -0.537 2.392 7.463
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 58.7830 3.7323 15.750 < 2e-16 ***
FACULTYS 15.8978 6.9754 2.279 0.0295 *
STUDY 3.2431 0.5329 6.086 8.46e-07 ***
FACULTYS:STUDY -2.1480 0.8621 -2.492 0.0181 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3.474 on 32 degrees of freedom
Multiple R-squared: 0.5956, Adjusted R-squared: 0.5577
F-statistic: 15.71 on 3 and 32 DF, p-value: 1.86e-06
平均の推定値の計算で混乱しないように、念のためcontrastsをデフォルト（contr.treatment）に戻してから分析。表1-25の1行目、2行目の切片については、（Intercept）と（Intercept）にFACULTYSを加えた値（書籍の切片[L]は誤植があるようだ。SPSSでもRと同じ値になった）。標準誤差は、（Intercept）の方はRの値で良いが、FACULTYS（S学部）については、学部のカテゴリのlevelesを入れ替えて再分析する必要がある。しかし、切片の標準誤差はあまり使い道もないので、特に計算する必要はないかもしれない。
次にSTUDYの係数。こちらは連続変数なので、係数は傾きとなる。こでの分析では、L学部を基準とするので、STUDYの係数はL学部の傾きとなり、S学部の傾きはSTUDYの係数に交互作用項（FACULTYS:STUDY）の値を加えたもの。ただし、L学部の標準誤差はSTUDY行の値であるが、S学部の標準誤差はS学部を基準カテゴリとして分析をやり直さなければ得られない。
次に、表1-26の値を拾う。これはタイプ3平方和を用いて分散分析を行えばよいので、簡単。contrastsをcontr.sumに変更してから分析すればよい。
> options(contrasts = c("contr.sum", "contr.sum"));
> m1.11.1s <- lm(ACHIEVE ~ FACULTY + STUDY + FACULTY:STUDY, df1);
> Anova(m1.11.1s, type=3);
Anova Table (Type III tests)
Response: ACHIEVE
Sum Sq Df F value Pr(>F)
(Intercept) 4418.2 1 366.0924 < 2.2e-16 ***
FACULTY 62.7 1 5.1944 0.02947 *
STUDY 305.6 1 25.3226 1.81e-05 ***
FACULTY:STUDY 74.9 1 6.2079 0.01809 *
Residuals 386.2 32
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
こちらは特に問題なし。
続いて表1-27の値を拾う。これは学部による切片の違い、自主学習の効果（傾き≠0といえるかどうか）、学部によって自主学習の効果（傾き）が異なるかどうかを検定したもので、検定方法は異なるが表1-26と等価（p値は完全に一致）。表1-26はF値による検定で、表1-27はt値による検定であるが、後者のt値を2乗すれば前者のF値になる。
値を拾っていくと、「α1-α2」は切片の違いを検定しているので、最初の回帰係数表（summaryの出力）のFACULTYS行の値と一致する。符号は逆転しているが、特に問題はない。次に「β」は自主学習単独の効果であるが、最初の回帰係数表のSTUDY行は、基準のL学部の傾きとなってしまうので、この検定の値とは一致しない。これは、タイプ3平方和を求める時と同様にcontrastsをcontr.sumにして分析した結果の回帰係数表（下記）のSTUDY行に一致する。ダミー変数に変換する際に、列ごとの合計が0になるようにコーディングされているので、STUDY行の出力がL学部でもなくS学部でもなく、自主学習単独の効果を検定することになっているのだと思う。
最後の「ω1-ω2」は、最初の回帰係数表の交互作用項の値と一致する。contrastsをcontr.sumに変更した場合の交互作用項とは、推定値と標準誤差が異なるが、t値とp値は一致しているので検定内容としては等価だろう。基準の取り方が違うだけである。
こんな面倒なことをせずとも分散分析表の結果を用いれば事足りるが、一応値を拾うことができたので、良しとしよう。分散分析の検定結果と、回帰係数の検定結果の対応関係を考えてみることで、分散分析と回帰分析を一般線形モデルで統一的に捉えるという意味が、（数学的にはともかくとしても）枠組みとして理解が進む気がする。また、重回帰分析単独で考えれば、カテゴリカル変数の各水準の効果（係数の値）だけに注目しがちだが、分散分析表も合わせてみることで、各水準を通して見た変数全体としての効果量（偏η2乗）も意味ある知見だということに気が付く。分散分析は、カテゴリ間のどこかに差があることを見出すものであるため、カテゴリ数が3つ以上になるとどこに差があるのかは、多重比較などで別途検討する必要がある。ただ、比較対象の基準が特定のカテゴリならば、回帰係数表の係数の検定だけでも十分かもしれない（と思う）。
> summary(m1.11.1s);
Call:
lm(formula = ACHIEVE ~ FACULTY + STUDY + FACULTY:STUDY, data = df1)
Residuals:
Min 1Q Median 3Q Max
-5.537 -2.632 -0.537 2.392 7.463
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 66.7319 3.4877 19.134 < 2e-16 ***
FACULTY1 -7.9489 3.4877 -2.279 0.0295 *
STUDY 2.1691 0.4311 5.032 1.81e-05 ***
FACULTY1:STUDY 1.0740 0.4311 2.492 0.0181 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 3.474 on 32 degrees of freedom
Multiple R-squared: 0.5956, Adjusted R-squared: 0.5577
F-statistic: 15.71 on 3 and 32 DF, p-value: 1.86e-06