二元配置要因分散分析2
大学生の自習時間に関する学部別学年別調査データ(架空データ)を用いる。
ある大学で学生の一日あたりの平均自習時間を調査した。学生を学部ごと、学年ごとに無作為で4人抽出し、1ヶ月の平均時間を調査した。
このデータから学部・学年の違いが自習時間に影響するかどうかついて分析せよ。
| 時間 | 学部 | 学年 |
|---|---|---|
| 94 | 経済 | 2回生 |
| 112 | 文 | 1回生 |
| 132 | 理工 | 4回生 |
| 131 | 社会 | 2回生 |
| … | ||
| 経済 | 法 | 文 | 理工 | 社会 | 全学部 | |
|---|---|---|---|---|---|---|
| 1回 | 121.50 | 116.50 | 113.75 | 112.25 | 115.75 | 115.95 |
| 2回 | 103.25 | 104.25 | 118.75 | 124.50 | 122.50 | 114.65 |
| 3回 | 113.75 | 112.25 | 111.00 | 112.75 | 106.25 | 111.20 |
| 4回 | 90.00 | 88.50 | 106.00 | 129.50 | 96.50 | 102.10 |
| 全学年 | 107.13 | 105.38 | 112.38 | 119.75 | 110.25 |
縦方向、横方向の組み合わせにおいて、平均に差があるのか。あるとすればどの組み合わせか。
分析対象
アンケート調査における元データをそのまま用いる。
その上で今回用いるデータは各学部＊各学年の人数がすべて同じ（4人）である事を前提としている(「釣り合い型データ」)。
data <- read.csv("http://kyoto-edu.sakura.ne.jp/weblesson/statistics/data/self-study1.csv", fileEncoding = "utf-8")
summary(data)
value <- data[[1]] #自習時間を検定変数に
factor1 <- factor(data[[3]]) #学年を要因1
factor2 <- factor(data[[2]]) #学部を要因2
#平均値
mean.table <- tapply(value,list(factor1,factor2),mean) #セルごとの平均(クロス集計)
factor1.mean <- tapply(value,factor1,mean) #要因1:学年単位の平均
factor2.mean <- tapply(value,factor2,mean) #要因2:学部単位の平均
all.mean <- mean(value) #全平均
#データの大きさ
factor1.N <- tapply(value,factor1,length) #要因1のサンプルサイズ
factor2.N <- tapply(value,factor2,length) #要因2のサンプルサイズ
N <- length(value) #サンプルサイズ
#項目数
factor1.num <- length(factor1.N) #要因1数
factor2.num <- length(factor2.N) #要因2数
all.num <- factor1.num*factor2.num #表のセル数
#自由度
factor1.df <- factor1.num - 1 #要因1自由度
factor2.df <- factor2.num -1 #要因2自由度
mix.df <- factor2.df*factor1.df #交互作用自由度
residual.df <- length(value) -1 - (factor2.df + factor1.df + mix.df) #残差自由度
二つの要因偏差
二元配置分散分析には新たに「交互作用」が加わる。
- 個々の値 ← 全体の平均 + 要因1ごとの平均差 + 要因2ごとの平均差 + 要因1要因2の交互作用 + 誤差(残差)
交互作用というのは複数の要因が組み合わさって出てくる作用である。例えば4学年は自習時間が減る傾向があるとして、特定の学部だけは逆に増えるというようなケースでは学部と学年の交互作用がみられるということである。ただし各カテゴリ内にデータが一つしかない時は交互作用は検出できない(理論上残差がなくなる)。
#学年ごとの平均差(要因1偏差) factor1.effects <- sapply(factor1,function(x)mean(data$time[which(factor1==x)])) - all.mean #学部ごとの平均差(要因2偏差) factor2.effects <- sapply(factor2,function(x)mean(data$time[which(factor2==x)])) - all.mean #要因1*要因2を合わせて交互要因列を作成 mix <- paste(data$grade,data$faculty) data2 <- transform(data,mix=mix) #交互作用-クロス表各セルの平均値から各要因の影響を除いた値と全体平均との偏差 mix.effects <- sapply(data2$mix,function(x)mean(data$time[which(data2$mix==x)])) - factor1.effects - factor2.effects - all.mean #残差 ← 個々の値 - 全体の平均 - 学年平均差 - 学部平均差 - 交互作用 residual.value <- value - all.mean - factor1.effects - factor2.effects - mix.effects
ここの計算で要因1*要因2にはいるデータの数がすべて等しいことを前提としている。
二つの要因分散と残差分散
偏差から変動、分散を求めていく。
#変動(偏差平方和) total.ss <- sum((value - all.mean)^2) factor1.ss <- sum(factor1.effects^2) factor2.ss <- sum(factor2.effects^2) mix.ss <- sum(mix.effects^2) residual.ss <- sum(residual.value^2) #分散 factor1.ms <- factor1.ss / factor1.df factor2.ms <- factor2.ss / factor2.df mix.ms <- mix.ss / mix.df residual.ms <- residual.ss / residual.df
F検定
#分散比(F値) factor1.F <- factor1.ms / residual.ms factor2.F <- factor2.ms / residual.ms mix.F <- mix.ms / residual.ms #有意性(p値) factor1.p <- 1 - pf(factor1.F,factor1.df,residual.df) factor2.p <- 1 - pf(factor2.F,factor2.df,residual.df) mix.p <- 1 - pf(mix.F,mix.df,residual.df)
| 変動要因 | 変動 | 自由度 | 分散 | F | p |
|---|---|---|---|---|---|
| 学年 | 2341.45 | 3 | 780.4833 | 4.68173 | 0.005276914 |
| 学部 | 2010.7 | 4 | 502.675 | 3.015296 | 0.02472752 |
| 交互作用 | 4495.3 | 12 | 374.6083 | 2.247088 | 0.02012015 |
| 残差 | 10002.5 | 60 | 166.7083 |
※交互作用に注目。
事後検定
- 単純主効果検定（交互作用が有意の時）
-
- 学部を固定して各学部に対して学年と時間の一元配置分散分析
- 学年を固定して各学年に対して学部と時間の一元配置分散分析
- 多重比較検定
-
- 最初の二元配置分散分析で有意と出た主効果に対して多重比較検定
grade N df mean 1st 20 19 115.950 2nd 20 19 114.650 3rd 20 19 111.200 4th 20 19 102.100 faculty N df mean economy 16 15 107.125 law 16 15 105.375 literature 16 15 112.375 science 16 15 119.750 sociology 16 15 110.250 - 単純主効果検定（一元配置分散分析）で有意と出たすべての主効果で多重比較検定
economy law literature science sociology 1st 121.50 116.50 113.75 112.25 115.75 2nd 103.25 104.25 118.75 124.50 122.50 3rd 113.75 112.25 111.00 112.75 106.25 4th 90.00 88.50 106.00 129.50 96.50
- 最初の二元配置分散分析で有意と出た主効果に対して多重比較検定
クロス表中の組み合わせを一つずつ精査して、有意差がある組み合わせをピックアップしていく。
非釣り合い型データの分析
偏差平方和を計算する過程で、個々の要素に入るデータ数が等しい（釣り合い型）ことが前提とされていた。心理学や理系分野においてはカテゴリーごとに計画的にデータを採取することは可能だが、とりわけアンケート調査においてそれをコントロールするのは難しい。となるとこのままではこの分析は使えないということになる。
ということで平方和計算に一工夫する必要がある。この工夫のやり方（流儀）によって複数のタイプの平方和が提唱されている。そのなかでも特にType1,Type2,Type3が使われている。
Type1はどちらを要因1とするか、その投入順によって結果が変わってしまう。これは少し使いづらい。というのでType2またはType3が実用的。SPSSはType3を使っている。
Rで普及しているaov関数は残念ながらType1を採用しているために、社会学の文脈ではあまり嬉しくない。
twoway.factorial.anova関数(socialStatisticsBasic.Rで読み込まれる自作関数)
- twoway.factorial.anova(formula,data,type,interaction)
-
- formula=「従属変数（検定変数）~因子1*因子2」
- data=モデル中の変数を含むデータフレーム(formulaが実体を持つ時は省略可)
- type=平方和のタイプ(1または2)
- interaction=グラフを出力するか(TRUE/FALSE 初期値はFALSE)
アンケート調査データの分析に特化している。2元配置までで反復測定分散分析は出来ない。平方和のタイプは1か2のみ。初期はType2。
実験
非釣り合い型のデータの分析をしてみる。Type1,Type2双方でfactor1とfactor2を入れ替えて再度計算を行い、結果を比較せよ。
data <- read.csv("http://kyoto-edu.sakura.ne.jp/weblesson/statistics/data/self-study2.csv",na.strings=99999, fileEncoding = "utf-8")
value <- data[[1]] #自習時間を検定変数に
factor1 <- factor(data[[3]]) #学年を要因1
factor2 <- factor(data[[2]]) #学部を要因2
#一元配置分散分析
source("http://kyoto-edu.sakura.ne.jp/weblesson/statistics/socialStatisticsBasic.R", encoding="UTF-8")
output.twoway.anova1.1 <- twoway.factorial.anova(value~factor1*factor2,type=1)
output.twoway.anova1.2 <- twoway.factorial.anova(value~factor2*factor1,type=1)
output.twoway.anova2.1 <- twoway.factorial.anova(value~factor1*factor2,type=2)
output.twoway.anova2.2 <- twoway.factorial.anova(value~factor2*factor1,type=2)
#結果出力
output.twoway.anova1.1
output.twoway.anova1.2
output.twoway.anova2.1
output.twoway.anova2.2
Rでの二元配置分散分析
data <- read.csv("http://kyoto-edu.sakura.ne.jp/weblesson/statistics/data/self-study1", fileEncoding = "utf-8")
summary(data)
value <- data[[1]] #自習時間を検定変数に
factor1 <- factor(data[[3]]) #学年を要因1
factor2 <- factor(data[[2]]) #学部を要因2
#分散分析表
summary(aov(value~factor1*factor2))
R標準のaov関数は平方和タイプ1にしか対応していないので「非釣り合い型」データの分析には不向きである。