最終更新日 : 2008.10.1
1要因の分散分析モデル
一般にt検定というと、それは2群(独立2標本)の平均値の差の検定のことを指すことが多いです。そして2群以上の平均値の差の検定は分散分析を用いるのだと教わったかもしれません。実際、多くの統計学の参考書にはそのように書いてあることでしょう。ところが、一般線形モデルの傘下ではt検定も分散分析も同じ理屈で考えることができます。ここでは1要因の分散分析モデルとして、t検定と分散分析(一元配置分散分析)で扱うデータを解析してみましょう。
■ t検定を一般線形モデルで解析する
仮説 : 中学生の1日の平均摂取カロリーと高校生の1日の平均摂取カロリーには差がある。
データ形式 :
1日の摂取カロリー
(単位:kcal)中学生 高校生 No. 1 1387 1044 No. 2 1398 1410 No. 3 1380 1599 No. 4 1257 1207 No. 5 1562 1308
さて、一般線形モデルとして扱うためにはこの表にまとめられたデータをデータフレームの形にまとめ直す必要があります。
学校 摂取カロリー 1
2
3
4
5
6
7
8
9
101
1
1
1
1
2
2
2
2
21387
1298
1380
1257
1562
1044
1410
1599
1207
1308
ここで用いられている変数は2つです。1つは学校という変数で、これはカテゴリカル型で中学生(=1)と高校生(=2)という2つの水準をもっています。もう1つは摂取カロリーという変数で、これは連続型です。この2つの変数をモデル式に書くと次式のように書くことができます。
摂取カロリー = 学校
これは「摂取カロリーは学校によって説明される」という意味です。従属変数が連続型、独立変数がカテゴリカル型(2水準)なので、これは1要因の分散分析モデルです。これを実際にRで解析してみましょう。
> school <- rep(c(1, 2), c(5, 5)) > school [1] 1 1 1 1 1 2 2 2 2 2 > school <- as.factor(school) > school [1] 1 1 1 1 1 2 2 2 2 2 Levels: 1 2 > calorie <- c(1387, 1298, 1380, 1257, 1562, 1044, 1410, 1599, 1207, 1308) > calorie [1] 1387 1298 1380 1257 1562 1044 1410 1599 1207 1308 > result <- lm(calorie ~ school) > summary(aov(result)) Df Sum Sq Mean Sq F value Pr(>F) school 1 9986 9986 0.3476 0.5717 Residuals 8 229800 28725
得られた分散分析表において、Pr(<F)の部分を見ると0.5717となっているのが確認できます。これはp値とよばれるもので、この値が有意水準αの値よりも小さければ帰無仮説を棄却することになるわけです。
ちなみに、最初に述べた仮説は実験者がこうなるであろうと考えた仮説であって、帰無仮説ではありません。統計的仮説検定では帰無仮説はH0で表され、対立仮説をH1で表します。
H0 : 中学生と高校生で摂取カロリーの平均値に差はない
H1 : 中学生と高校生で摂取カロリーの平均値に差がある
今回の結果によれば、仮に有意水準をα=0.05と設定しておいたのであればp=0.5717>0.05となるので、帰無仮説を採択するということになります。ところで有意水準αは「100回検定を繰り返し行ったら、そのうちn回は誤った結論を導く危険がある」ということを意味しており、これを特に危険率といったりもします。α=0.05なら100回に5回は失敗する可能性があるというわけです。また、有意水準は検定を行う前に設定しておくべきものです。
■ 分散分析を一般線形モデルで解析する
繰り返しになりますが、2群の平均値の差の検定はt検定を用い、2群以上の平均値の差の検定を行う場合は分散分析を用いると教わったかもしれません。しかし、一般線形モデルで解析しようとした場合、これらはどちらも同じ1要因の分散分析モデルとして表されます。ただ要因がもつ水準数(これはカテゴリ、あるいは属性ともいうのでした)が異なるというだけです。
仮説 : 期末テストの総合点が良い生徒は家庭での勉強時間が長い。
データ形式 :
この仮説に基づいてデータを解析する方法は2つあります。1つは生徒各々に家庭での勉強時間を質問紙などで聞き、勉強時間そのものを連続型のデータであるとして扱うことです。もう1つは勉強時間をいくつかの水準に区別してカテゴリカル型のデータとして扱うという方法です。今回は分散分析モデルとして解析するので後者の方法を採用することにします。
今回用いる変数は期末テストの総合点(連続型)と家庭での勉強時間(カテゴリカル型:3水準)の2つです。なお、勉強時間の水準は全くしない(=1)、1時間から2時間程度(=2)、2時間以上(=3)とします。
勉強時間 総合得点 1
2
3
4
5
6
7
8
9
10
11
12
13
14
151
1
1
1
1
2
2
2
2
2
3
3
3
3
3324
250
207
213
278
381
370
319
334
416
424
383
426
424
361
解析するモデル式は次のようになります。このモデル式が意味するのは「総合得点は勉強時間によって説明される」ということです。もう少し難しい言葉でいうならば「総合得点の変動は勉強時間の水準ごとの平均値を当てはめることによって説明される」といえます。詳しくは後章の「分散分析の基本原理」を読めば分かりますので、ここでは詳しい説明は省きます(この時点ではあまり難しいことを考えなくて良いです)。
総合得点 = 勉強時間
このモデル式をRのlm()を用いて解析してみましょう。
> learn.time <- rep(c(1, 2, 3), c(5, 5, 5)) > learn.time <- as.factor(learn.time) > total.score <- c( + 324, 250, 207, 213, 278, + 381, 370, 319, 334, 416, + 424, 383, 426, 424, 361) > result2 <- lm(total.score ~ learn.time) > summary(aov(result2)) Df Sum Sq Mean Sq F value Pr(>F) learn.time 2 59735 29867 18.955 0.0001932 *** Residuals 12 18908 1576
さて、出力された分散分析表よりp=0.0001932という値が確認できたので、有意水準を0.05とするならば帰無仮説を棄却することになります。したがって、どうやら勉強時間によって成績に差がでるということが認められたわけです。ただし、この分散分析表からだけでは、どの水準がどれほど総合成績に影響しているかを把握することはできません。
H0 : 3つの水準ごとの平均値に差はない
H1 : 3つの水準ごとの平均値に差がある
■ 係数表の見方
一般線形モデルとして解析することの利点として、例えばRではlm()という関数を用いることになりますが、モデル解析の結果から係数表とよばれるものを得ることができるということがあげられます。それをどのようにして出力させるかというと、次のようにします(実際にはもう少し多くの情報が出力されますが、ここには係数表(Coefficients:という部分)のみを掲載してあります)。
> summary(result2) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 254.40 17.75 14.331 6.54e-09 *** learn.time2 109.60 25.11 4.366 0.00092 *** learn.time3 149.20 25.11 5.943 6.79e-05 ***
このままでは見にくいという人もいるでしょうから、最初は次のように数式で表してあげると分かりやすいでしょう。
実はこれ、learn.timeの最初の水準を基準として他の水準2、3がどれほど従属変数の増減に寄与しているかということを意味しています。だからこの場合だとlearn.time1は0となっていますから、これを基準にlearn.time2は109.60だけlearn.time1より従属変数の増減に寄与しているといえます(符号がプラスなので、当然、増加的に寄与するということ)。
言い方を変えれば、水準2(勉強時間が1時間から2時間の生徒)の平均値は水準1(勉強時間が0時間の生徒)の平均値より109.60だけ高いということです。同じように水準3(勉強時間が2時間以上の生徒)の平均値は水準1の平均値より149.20だけ高いということです。
もう少し数学的な言い方をするなら、定数項の254.40はlearn.time1を基準としたときの平均値である、ということです。実際に水準1に属するデータの平均値を計算してみましょう。
> mean(c(324, 250, 207, 213, 278)) [1] 254.4
すると、この値は定数項の254.40と一致するのが確認できます。同じように水準2の平均値は
> mean(c(381, 370, 319, 334, 416)) [1] 364
となり、これは254.40+109.60=360と一致します。最後に水準3についても、
> mean(c(424, 383, 426, 424, 361)) [1] 403.6
と、これは254.40+149.20=403.6と一致します。
以上をまとめると、係数表からは各水準の平均値が推定することができるということになります。また、dummy.coer()という関数を利用すると、上で提示した数式により忠実な形で係数が表示されます。
> dummy.coef(result2) Full coefficients are (Intercept): 254.4 learn.time: 1 2 3 0.0 109.6 149.2
■ 標準誤差の意味
先ほどの係数表を再掲します。
Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 254.40 17.75 14.331 6.54e-09 *** learn.time2 109.60 25.11 4.366 0.00092 *** learn.time3 149.20 25.11 5.943 6.79e-05 ***
ここにStd. Errorという部分がありますね。これは標準誤差とよばれているもので、簡単にいえば推定値の推定精度を表す指標です。そしてこの値が小さいほどよい推定をしているといえます。なかなか難しいことですが、Estimate / Std. Error = t valueとなり、例えばlearn.time2については109.60 / 25.11 = 4.364795(≒4.366)になることが確認できます。t値が大きいほどそれに対するp値は小さくなりますから、標準誤差が小さいほど(0に近いほど)t値も大きくなることがわかるでしょう。だから標準誤差の値が小さいほどよい推定であるといえるわけです。
ところでこのt値は一体何のための統計量なのでしょうか。それは「係数は0である」という帰無仮説を検定するための統計量です。だから「learn.time2の係数は0である」という帰無仮説についての統計量は先ほどやったように、
109.60 / 25.11 = 4.364795
となり、これに対するp値は0.00092であり、この値は0.05よりも小さい。したがって、「learn.time2の係数は0ではない」ということがいえるわけです。ちなみに、定数項を含めたいずれの項のt値も自由度12のt分布に従います。自由度は全データ数から推定したパラメータの数を引いた値になります。今回の例では全データ数が15で、推定したパラメータは水準1から3の平均値3つですから15-3=12となるわけです(Interceptはlearn.time1の平均値であることに注意)。