グループ間の母平均の差を検定する分散分析,ANOVA(Analysis of Variance)の初歩について解説する。
mux <- 45
muy <- 50
muz <- 55
n <- 25
sd0 <- 10
multip <- 2
uj <- rnorm(n, mean=0, sd=sd0)
uj <- uj - mean(uj)
x <- data.frame("group"=rep("x", n), "score"=mux + uj)
y <- data.frame("group"=rep("y", n), "score"=muy + uj)
z <- data.frame("group"=rep("z", n), "score"=muz + uj)
x2 <- data.frame("group"=rep("x", n), "score"=mux + multip*uj)
y2 <- data.frame("group"=rep("y", n), "score"=muy + multip*uj)
z2 <- data.frame("group"=rep("z", n), "score"=muz + multip*uj)
d01 <- rbind(x, y, z) # ; str(d01)
d02 <- rbind(x2, y2, z2)
min_y <- min(d01$score, d02$score) -5
max_y <- max(d01$score, d02$score) +5
ano1 <- anova(aov(score ~ group, data=d01))
ano2 <- anova(aov(score ~ group, data=d02))
mean1 <- by(d01$score, d01$group, mean)
mean2 <- by(d02$score, d02$group, mean)
g1.xlab <- sprintf("mean(x, y, z)=(%.1f, %.1f, %.1f)", mean1[1], mean1[2], mean1[3])
g2.xlab <- sprintf("mean(x, y, z)=(%.1f, %.1f, %.1f)", mean2[1], mean2[2], mean2[3])
g <- ggplot(d01, aes(x=group, y=score))
g <- g + geom_boxplot(aes(color=group))
g <- g + geom_jitter(alpha=.4, aes(color=group))
g <- g + theme_classic() + ylim(min_y, max_y)
g <- g + theme(legend.position="none")
g1 <- g + labs(title=sprintf("F(%d, %d)=%.3f, p-value=%.4f", ano1$Df[1], ano1$Df[2], ano1$"F value", ano1$"Pr(>F)"), x=g1.xlab, y="score")
g <- ggplot(d02, aes(x=group, y=score))
g <- g + geom_boxplot(aes(color=group))
g <- g + geom_jitter(alpha=.4, aes(color=group))
g <- g + theme_classic() + ylim(min_y, max_y)
g <- g + theme(legend.position="none")
g2 <- g + labs(title=sprintf("F(%d, %d)=%.3f, p-value=%.4f", ano2$Df[1], ano2$Df[2], ano2$"F value", ano2$"Pr(>F)"), x=g2.xlab, y="")
実行するたびに結果が変わるので,例としては余り良くない場合もある。
gridExtra::grid.arrange(g1, g2, nrow = 1)
分散分析は,級間変動と級内変動の比が核心。
級間変動は,級間分散,級間平方和で考えても同じ。群間,グループ間の(平均の)差を示す。何らかのグループ/カテゴリーに属する事で人々がどの程度違ってくるか,である。
級内変動は,級内分散,級内平方和で考えても同じ。同じグループ/カテゴリーに属する人達の間での個人差を示す。グループの効果を見る観点からすればこれは誤差の扱いになる。
グループ j に属している個人 i の測定値を次のように変形して考えると,第二項の( )内が全体平均からのグループ j の平均の偏差,第三項の( )内がグループ j の平均からの個人 i の測定値の偏差となる。前者がグループ j に属することの効果を表し,後者は誤差として扱われる。
\[y_{ij}=\bar{y}+(\bar{y_j}-\bar{y})+(y_{ij}-\bar{y_j})\]
平均からの偏差の二乗(偏差平方)は次のようになる。
\[(y_{ij}-\bar{y})^2=\{(\bar{y_j}-\bar{y})+(y_{ij}-\bar{y_j})\}^2=(\bar{y_j}-\bar{y})^2+2(\bar{y_j}-\bar{y})(y_{ij}-\bar{y_j})+(y_{ij}-\bar{y_j})^2\]
この偏差平方を,以下のように,ひとまずグループ j の中だけで合計する。\(n_{j}\) はグループ j に属するケース数である。一行目の右辺の 3 つの項の変形の際に,添え字 i を含まない第一項と第二項前半の( )は,i のΣに関しては定数扱いである点に着目する。定数を 1 からまで \(n_{j}\) まで足し合わせるというのは,その定数を\(n_{j}\) 倍するのと同じである。
\[\begin{align*} \sum_{i=1}^{n_{j}}(y_{ij}-\bar{y})^2&=\sum_{i=1}^{n_{j}}(\bar{y_j}-\bar{y})^2+\sum_{i=1}^{n_{j}}2(\bar{y_j}-\bar{y})(y_{ij}-\bar{y_j})+\sum_{i=1}^{n_{j}}(y_{ij}-\bar{y_j})^2\\ &=n_{j}(\bar{y_{j}}-\bar{y})^2+2(\bar{y_{j}}-\bar{y})\sum_{i=1}^{n_{j}}(y_{ij}-\bar{y})+\sum_{i=1}^{n_{j}}(y_{ij}-\bar{y_{j}})^2\\ &=n_{j}(\bar{y_{j}}-\bar{y})^2+\sum_{i=1}^{n_{j}}(y_{ij}-\bar{y_{j}})^2 \end{align*}\]
二行目の第二項のΣの部分は,グループ j の測定値の合計からグループ j の平均の \(n_{j}\) 倍を引くことになるので,以下の通り 0 となる。よって上式の三行目のようになる。 \[\sum_{i=1}^{n_{j}}(y_{ij}-\bar{y_{j}})=\sum_{i=1}^{n_{j}}y_{ij}-\sum_{i=1}^{n_{j}}\bar{y_{j}}=\sum_{i=1}^{n_{j}}y_{ij}-n_{j}\bar{y_{j}}=\sum_{i=1}^{n_{j}}y_{ij}-n_{j}\frac{1}{n_{j}}\sum_{i=1}^{n_{j}}y_{ij}=0\]
これはグループ j の中だけでの合計である。このグループごとの合計を m 個のグループすべてについて合計することによって,全体の偏差平方和は次のようになる。グループの数 m は慣習的に水準数と呼ぶ。
\[\sum_{j=1}^{m}\sum_{i=1}^{n_{j}}(y_{ij}-\bar{y})^2=\sum_{j=1}^{m}n_{j}(\bar{y_{j}}-\bar{y})^2+\sum_{j=1}^{m}\sum_{i=1}^{n_{j}}(y_{ij}-\bar{y_{j}})^2\] 左辺を総平方和\(SS_{T} (SS_{Total})\),右辺第一項を級間平方和(グループ間平方和)\(SS_{B} (SS_{Between})\),右辺第二項を級内平方和(グループ内平方和)\(SS_{W} (SS_{Within})\)と呼ぶ。そして級間平方和が関心の対象であるグループによる影響,級内平方和は誤差としての個人差となる。この平方和の分解が分散分析の根幹である。 \[SS_{Total}=SS_{Between}+SS_{Within}\]
偏差平方和の分解の式
\[\sum_{j=1}^{m}\sum_{i=1}^{n_{j}}(y_{ij}-\bar{y})^2=\sum_{j=1}^{m}n_{j}(\bar{y_{j}}-\bar{y})^2+\sum_{j=1}^{m}\sum_{i=1}^{n_{j}}(y_{ij}-\bar{y_{j}})^2\] 測定値に等分散正規性\(y_{ij}\sim N(\mu,\sigma^2)\)が仮定できるとき,平方和の分解の左辺を母分散\(\sigma^2\)で割っ た式は次のように変形される。
\[\begin{align*} \sum_{j=1}^{m}\sum_{i=i}^{n_{j}}\left(\frac{y_{ij}-\bar{y}}{\sigma}\right)^2&=\sum_{j=1}^{m}\sum_{i=i}^{n_{j}}\left\{\frac{(y_{ij}-\mu)-(\bar{y}-\mu)}{\sigma}\right\}^2\\ &=\sum_{j=1}^{m}\sum_{i=i}^{n_{j}}\left\{\left(\frac{y_{ij}-\mu}{\sigma}\right)^2-\frac{2(y_{ij}-\mu)(\bar{y}-\mu)}{\sigma^2}+\left(\frac{\bar{y}-\mu}{\sigma}\right)^2\right\}\\ &=\sum_{j=1}^{m}\sum_{i=i}^{n_{j}}\left(\frac{y_{ij}-\mu}{\sigma}\right)^2-\frac{2(\bar{y}-\mu)}{\sigma^2}\sum_{j=1}^{m}\sum_{i=i}^{n_{j}}(y_{ij}-\mu)+\sum_{j=1}^{m}\sum_{i=i}^{n_{j}}\left(\frac{\bar{y}-\mu}{\sigma}\right)^2\\ &=\sum_{j=1}^{m}\sum_{i=i}^{n_{j}}\left(\frac{y_{ij}-\mu}{\sigma}\right)^2-\frac{2(\bar{y}-\mu)}{\sigma^2}(n\bar{y}-n\mu)+n\left(\frac{\bar{y}-\mu}{\sigma}\right)^2\\ &=\sum_{j=1}^{m}\sum_{i=i}^{n_{j}}\left(\frac{y_{ij}-\mu}{\sigma}\right)^2-\frac{2n(\bar{y}-\mu)^2}{\sigma^2}+n\left(\frac{\bar{y}-\mu}{\sigma}\right)^2\\ &=\sum_{j=1}^{m}\sum_{i=i}^{n_{j}}\left(\frac{y_{ij}-\mu}{\sigma}\right)^2-n\left(\frac{\bar{y}-\mu}{\sigma}\right)^2\\ &=\sum_{j=1}^{m}\sum_{i=i}^{n_{j}}\left(\frac{y_{ij}-\mu}{\sigma}\right)^2-\left(\frac{\bar{y}-\mu}{\frac{\sigma}{\sqrt{n}}}\right)^2 \sim \chi_{(n-1)}^2 \end{align*}\]
これは,総ケース数と同じ n 個の「標準化された正規変量」の二乗和から,確率変数である標本平均を標準化したもの(第 3 章 1-2 参照),すなわち 1 個の標準正規変量の二乗を引いたものなので,定義上,自由度 n-1 のカイ二乗分布に従う。
自由度ν の\(\chi^2\) 分布とは,互いに独立なν 個の標準正規変量の二乗和として定義される。したがって,カイ二乗 分布で近似できるために等分散正規性が必要とされる。
平方和の分解の右辺の第一項と第二項を母分散で割った\(SS_{Between}/\sigma^2\),\(SS_{Within}/\sigma^2\)もそれぞれ,自由度\(m-1, n-m\)のカイ二乗分布に従う。
\[\begin{align*} \sum_{j=1}^{m}\frac{n_{j}(\bar{y_{j}}-\bar{y})^2}{\sigma^2} &= \sum_{j=1}^{m}\left(\frac{\bar{y_{j}}-\bar{y}}{\frac{\sigma}{\sqrt{n_{j}}}}\right)^2=\sum_{j=1}^{m}\left\{\frac{(\bar{y_{j}}-\mu)-(\bar{y}-\mu)}{\frac{\sigma}{\sqrt{n_{j}}}}\right\}^2\\ &=\sum_{j=1}^{m}\left\{\left(\frac{\bar{y_{j}}-\mu}{\frac{\sigma}{\sqrt{n_{j}}}}\right)^2-\frac{2(\bar{y_{j}}-\mu)(\bar{y}-\mu)}{\left(\frac{\sigma}{\sqrt{n_{j}}}\right)^2}+\left(\frac{\bar{y}-\mu}{\frac{\sigma}{\sqrt{n_{j}}}}\right)^2\right\}\\ &=\sum_{j=1}^{m}\left(\frac{\bar{y_{j}}-\mu}{\frac{\sigma}{\sqrt{n_{j}}}}\right)^2-2(\bar{y}-\mu)\sum_{j=1}^{m}\frac{(\bar{y_{j}}-\mu)}{\frac{\sigma^2}{n_{j}}}+\sum_{j=1}^{m}\frac{(\bar{y}-\mu)^2}{\frac{\sigma^2}{n_{j}}}\\ &=\sum_{j=1}^{m}\left(\frac{\bar{y_{j}}-\mu}{\frac{\sigma}{\sqrt{n_{j}}}}\right)^2-2\frac{(\bar{y}-\mu)}{\sigma^2}\sum_{j=1}^{m}n_{j}(\bar{y_{j}}-\mu)+\frac{1}{\sigma^2}\sum_{j=1}^{m}n_{j}(\bar{y}-\mu)^2\\ &=\sum_{j=1}^{m}\left(\frac{\bar{y_{j}}-\mu}{\frac{\sigma}{\sqrt{n_{j}}}}\right)^2-2\frac{(\bar{y}-\mu)}{\sigma^2}n(\bar{y}-\mu)+\frac{1}{\sigma^2}n(\bar{y}-\mu)^2\\ &=\sum_{j=1}^{m}\left(\frac{\bar{y_{j}}-\mu}{\frac{\sigma}{\sqrt{n_{j}}}}\right)^2-\frac{n(\bar{y}-\mu)^2}{\sigma^2}\\ &=\sum_{j=1}^{m}\left(\frac{\bar{y_{j}}-\mu}{\frac{\sigma}{\sqrt{n_{j}}}}\right)^2-\left(\frac{\bar{y}-\mu}{\frac{\sigma}{\sqrt{n}}}\right)^2 \sim \chi_{(m-1)}^2 \end{align*}\]
最後の辺りで,\(n_{j}\bar{y_{j}}\)がj群の合計,よってそれを j=1 から m まで合計したものは全体の合計に なる事などに気付く事がポイントである。
\[\begin{align*} \sum_{j=1}^{m}\sum_{i=1}^{n_{j}}\left(\frac{y_{ij}-\bar{y_{j}}}{\sigma}\right)^2 &=\sum_{j=1}^{m}\sum_{i=1}^{n_{j}}\left\{\frac{(y_{ij}-\mu)-(\bar{y_{j}}-\mu)}{\sigma}\right\}^2\\ &=\sum_{j=1}^{m}\sum_{i=1}^{n_{j}}\left\{\left(\frac{y_{ij}-\mu}{\sigma}\right)^2-2\frac{(y_{ij}-\mu)(\bar{y_{j}}-\mu)}{\sigma^2}+\left(\frac{\bar{y_{j}}-\mu}{\sigma}\right)^2\right\}\\ &=\sum_{j=1}^{m}\sum_{i=1}^{n_{j}}\left(\frac{y_{ij}-\mu}{\sigma}\right)^2-\frac{2}{\sigma^2}\sum_{j=1}^{m}\sum_{i=1}^{n_{j}}(y_{ij}\bar{y_{j}}-\mu(y_{ij}+\bar{y_{j}})+\mu^2)+\sum_{j=1}^{m}\sum_{i=1}^{n_{j}}\left(\frac{\bar{y_{j}}-\mu}{\sigma}\right)^2\\ &=\sum_{j=1}^{m}\sum_{i=1}^{n_{j}}\left(\frac{y_{ij}-\mu}{\sigma}\right)^2-\frac{2}{\sigma^2}\sum_{j=1}^{m}(n_{j}\bar{y_{j}}^{2}-2\mu n_{j}\bar{y_{j}}+n_{j}\mu^2)+\sum_{j=1}^{m}n_{j}\left(\frac{\bar{y_{j}}-\mu}{\sigma}\right)^2\\ &=\sum_{j=1}^{m}\sum_{i=1}^{n_{j}}\left(\frac{y_{ij}-\mu}{\sigma}\right)^2-\frac{2}{\sigma^2}\sum_{j=1}^{m}n_{j}(\bar{y_{j}}-\mu)^2+\sum_{j=1}^{m}n_{j}\left(\frac{\bar{y_{j}}-\mu}{\sigma}\right)^2\\ &=\sum_{j=1}^{m}\sum_{i=1}^{n_{j}}\left(\frac{y_{ij}-\mu}{\sigma}\right)^2-\sum_{j=1}^{m}n_{j}\left(\frac{\bar{y_{j}}-\mu}{\sigma}\right)^2\\ &=\sum_{j=1}^{m}\sum_{i=1}^{n_{j}}\left(\frac{y_{ij}-\mu}{\sigma}\right)^2-\sum_{j=1}^{m}\left(\frac{\bar{y_{j}}-\mu}{\frac{\sigma}{\sqrt{n_{j}}}}\right)^2 \sim \chi_{(n-m)}^2 \end{align*}\]
これも最後の辺りで,i を含まない項は i についてのΣの中では定数扱いである事に気付くのがポイントである。
ここで,「カイ二乗分布に従う統計量÷その自由度」という統計量が二つあったとき,その比が F 分布として定義されることを利用する。
\[F_{(df_{1},df_{2})}=\frac{\chi_{(df_{1})}^2/df_{1}}{\chi_{(df_{2})}^2/df_{2}}\] 平方和(Sum of Squares, SS)の分解の式の右辺の第一項と第二項をそれぞれの自由度で割った統計量を平均平方(Mean Square, MS)という。これは,自由度一つあたりの平方和という意味である。 \[MS_{Between}=\frac{SS_{Between}}{m-1}\] \[MS_{Within}=\frac{SS_{Within}}{n-m}\] そしてこの二つの平均平方の比は,自由度 m-1,n-m の F 分布に従うことが分かる。 \[\frac{MS_{Between}}{MS_{Within}}=\frac{SS_{Between}/(m-1)}{SS_{Within}/(n-m)}=\frac{(SS_{Between}/\sigma^2)/(m-1)}{(SS_{Within}/\sigma^2)/(n-m)}=F_{(m-1,n-m)}\] 母分散が分母と分子で約分されて消え去り,平均平方の比自体はデータから具体的に計算可能な値になっているので,平均平方の比である F 統計量によってゼロ仮説を検定することが可能になる。
Rでは,自由度\(df_{1}, df_{2}\)のF分布の右端5%領域の限界値は,qf(p=.95, df1, df2)もしくはqf(p=.05, df1, df2, lower.tail=F)
として求められる。グラフは,curve(df(df1, df2), from=0, to=4)
などとすると描ける。以下は,グループの数mが3,総ケース数nが30を想定した数値例である。
m <- 3; n <- 30
qf(p=.95, df1=m-1, df2=n-m)
[1] 3.354131
qf(p=.05, df1=m-1, df2=n-m, lower.tail=F)
[1] 3.354131
curve(df(x, df1=m-1, df2=n-m), from=0, to=4)
グラフにいくつかオプションなどで装飾を施してみよう。
m <- 5; n <- 50; p0 <- .05
F0 <- qf(p=p0, df1=m-1, df2=n-m, lower.tail=F)
END <- qf(min(p0/2, .01), df1=m-1, df2=n-m, lower.tail=F)
curve(df(x, df1=m-1, df2=n-m), from=0, to=END, bty="n",
xlab="F統計量", ylab="確率密度",
main=sprintf("自由度%d,%dのF分布の%.1f%%棄却域", m-1, n-m, p0*100))
axis(side=1, at=F0, label=round(F0, 2), col.axis="red", las=2)
segments(rej <- seq(F0, END, by=.02), 0,
rej, df(rej, df1=m-1, df2=n-m), col="red")
abline(h=0, lty=2, col="grey")
まずは使用する変数の基本的な情報(度数分布や要約統計量,ケース数)を確認しておく必要がある。欠損値NAが含まれる場合はひと手間余計にかかる事が多いので注意しよう。また,本文と異なり,データフレイムの名前はdata01としている。
また,いちいちデータフレイム名 data01 を書く事を省略したい場合には with( ) を使用している。
data01 <- read.csv("practice.csv")
# まずは各変数の度数分布または基本要約統計量を確認
table(data01$job, useNA="always")
1 2 3 4 <NA>
182 92 38 68 4
summary(data01$fincome)
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.0 400.0 600.0 720.1 925.0 2500.0 87
# 次に,従業上の地位別の世帯年収の要約統計
with(data01, tapply(fincome, job, summary))
$`1`
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
75.0 500.0 700.0 804.3 925.0 2500.0 37
$`2`
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.0 325.0 550.0 631.8 925.0 1750.0 18
$`3`
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
75.0 375.0 500.0 630.4 800.0 2500.0 10
$`4`
Min. 1st Qu. Median Mean 3rd Qu. Max. NA's
0.0 200.0 600.0 670.4 925.0 2500.0 19
# 各群のケース数は以下のコマンドで分かる様に思えるが,欠損値を含む
with(data01, tapply(fincome, job, length))
1 2 3 4
182 92 38 68
# 欠損値抜きのケース数を求める( Rコマンダー の出力を参考にした)
with(data01, tapply(fincome, job, function(x) sum(!is.na(x))))
1 2 3 4
145 74 28 49
# tableで各群のケース数を求める事も不可能ではない。
t01 <- with(data01, addmargins(table(job, fincome, useNA="ifany")))
t01
fincome
job 0 25 75 125 200 300 400 500 600 700 800 925 1125 1375 1750 2500 <NA>
1 0 0 2 2 3 6 12 18 18 15 17 20 15 8 5 4 37
2 2 1 2 2 5 7 10 8 7 7 3 7 7 3 3 0 18
3 0 0 1 0 4 2 4 5 3 1 2 2 2 1 0 1 10
4 2 0 2 4 5 4 2 5 4 4 1 4 7 1 3 1 19
<NA> 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3
Sum 5 1 7 8 17 19 28 36 32 27 23 33 31 13 11 6 87
fincome
job Sum
1 182
2 92
3 38
4 68
<NA> 4
Sum 384
t01[, ncol(t01)] - t01[, ncol(t01)-1]
1 2 3 4 <NA> Sum
145 74 28 49 1 297
本文の箱ひげ図は次のスクリプトで描いている。グラフを描く時に with( ) を使うと出力が少しだけすっきりする。
with(data01,
boxplot(fincome ~ job, varwidth=T,
pars=par(bty="n", family="serif"), las=1,
main="従業上の地位別の世帯年収(万円)",
names=c("正規", "非正規", "自営", "無職"))
)
いよいよ(一要因)分散分析を行う。aov( )
関数,lm( )
関数,oneway.test( )
関数のうち,まずは最初の二つを説明する。
ここでも,出力を少しすっきりさせる為に with( ) を使っている。最初は却って面倒に感じるだろうが,慣れると結構便利である。例をよく見て参考にすると良い。
#### aov
aov(fincome ~ job, data01) # 間違い。自由度を見ると分かる。
Call:
aov(formula = fincome ~ job, data = data01)
Terms:
job Residuals
Sum of Squares 1038846 63166253
Deg. of Freedom 1 294
Residual standard error: 463.5204
Estimated effects may be unbalanced
88 observations deleted due to missingness
aov(fincome ~ factor(job), data01) # 正しい。要因型にする必要。
Call:
aov(formula = fincome ~ factor(job), data = data01)
Terms:
factor(job) Residuals
Sum of Squares 1950508 62254591
Deg. of Freedom 3 292
Residual standard error: 461.7366
Estimated effects may be unbalanced
88 observations deleted due to missingness
# 出力が見にくい場合は,文字型の変数を新規作成するのも良い。
data01$JOB <- factor(data01$job, levels=c(1:4), labels=c("1 正規", "2 非正規", "3 自営", "4 無職"))
with(data01, table(job, JOB, useNA="ifany")) # 新変数の確認
JOB
job 1 正規 2 非正規 3 自営 4 無職 <NA>
1 182 0 0 0 0
2 0 92 0 0 0
3 0 0 38 0 0
4 0 0 0 68 0
<NA> 0 0 0 0 4
# aov( )の結果をオブジェクトに格納しておくと便利。名前(ANOVA01)は自由。
ANOVA01 <- aov(fincome ~ JOB, data01)
summary(ANOVA01)
Df Sum Sq Mean Sq F value Pr(>F)
JOB 3 1950508 650169 3.05 0.029 *
Residuals 292 62254591 213201
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
88 observations deleted due to missingness
anova(ANOVA01)
出力結果を用いて,「グループ平均はすべて等しい」というゼロ仮説の有意性検定をF分布を用いて行う事が出来る。この例ではF統計量は約3.05,有意確率は約2.9%となり,5%水準で有意となった。
F0 <- (1950508/3)/(62254591/292) # Mean Squareの比
F0
[1] 3.049565
1 - pf(F0, df1=3, df2=292)
[1] 0.0289779
pf(F0, df1=3, df2=292, lower.tail=F) # これでも結果は同じ
[1] 0.0289779
#### lm
lm(fincome ~ job, data01) # 間違い。数量型になっている。
Call:
lm(formula = fincome ~ job, data = data01)
Coefficients:
(Intercept) job
825.43 -53.15
lm(fincome ~ factor(job), data01) # 要因型にしており正しい。
Call:
lm(formula = fincome ~ factor(job), data = data01)
Coefficients:
(Intercept) factor(job)2 factor(job)3 factor(job)4
804.3 -172.6 -174.0 -133.9
anova(LM01 <- lm(fincome ~ JOB, data01))
summary(LM01) # 因みにsummary( )関数の出力はこうなる
Call:
lm(formula = fincome ~ JOB, data = data01)
Residuals:
Min 1Q Median 3Q Max
-729.3 -304.3 -104.3 190.9 1869.6
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 804.31 38.35 20.976 < 2e-16 ***
JOB2 非正規 -172.55 65.97 -2.616 0.00936 **
JOB3 自営 -173.95 95.31 -1.825 0.06901 .
JOB4 無職 -133.90 76.30 -1.755 0.08031 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 461.7 on 292 degrees of freedom
(88 observations deleted due to missingness)
Multiple R-squared: 0.03038, Adjusted R-squared: 0.02042
F-statistic: 3.05 on 3 and 292 DF, p-value: 0.02898
次のスクリプトを実行すると,このF分布の5%棄却域(赤)とF値(紫)の関係をカラーで図示できる。関数ごとに分けて実行すると関数の意味が分かりやすい。
# F分布
F0 <- (1950508/3)/(62254591/292) # F統計量
DF1 <- 3; DF2 <- 292 # 自由度は関数のオプションの名称と区別する為敢えて大文字にした
curve(df(x, DF1, DF2), xlim=c(0, END <- 4), xaxp=c(0, END, 2),
bty="n", family="serif", xlab="", ylab="",
main=sprintf("自由度(%d, %d)のF分布の5%%棄却領域", DF1, DF2))
segments(x1 <- seq(x2 <- qf(.95, DF1, DF2), END, by=.02), 0,
x1, df(x1, DF1, DF2), col="red")
axis(side=1, at=round(x2, 3), col="red", col.axis="red", las=2, family="serif")
axis(side=1, at=round(F0, 3), col="purple", col.axis="purple", las=2, family="serif")
なお,各種の分散説明率(R-squared, Eta-squared, partial Eta-squared)については,二元配置分散分析の2-1 平方和の分解の末尾に例示を付け加えておく。
一つ一つの検定の有意水準をp1,そうした検定を3つ繰り返した場合の最終的な有意水準をp3とすると,p3 = p1 + p1 + p1 - 3p1p1 + p1p1p1 である。p1を.05(5%)とすると,p3=.142625となる。そこでp1=.05/3とすると,p3は5%以内に収まる。
# 多重比較
p1 <- .05
p3 <- p1*3 - 3*p1^2 + p1^3; p3
[1] 0.142625
p1 <- .05/3
p3 <- p1*3 - 3*p1^2 + p1^3; p3
[1] 0.0491713
世帯年収fincomeと従業上の地位jobで実際にボンフェローニ法を実行してみよう。共通の標準偏差を用いるか否かというオプションpool.sdによって結果が変わる。本文ではデフォルトのpool.sd=Tの結果を述べたが,ここでは両方実行する。(わざと,with( )を使った書き方と使わない書き方で示している。どう異なるか見比べて欲しい。)
# 共通の標準偏差を推定して検定(等分散性の仮定)
pairwise.t.test(data01$fincome, data01$job, p.adjust.method="bonferroni")
Pairwise comparisons using t tests with pooled SD
data: data01$fincome and data01$job
1 2 3
2 0.056 - -
3 0.414 1.000 -
4 0.482 1.000 1.000
P value adjustment method: bonferroni
# 等分散を仮定しないで検定
with(data01, pairwise.t.test(fincome, job, p.adjust.method="bonferroni", pool.sd=F))
Pairwise comparisons using t tests with non-pooled SD
data: fincome and job
1 2 3
2 0.032 - -
3 0.524 1.000 -
4 0.728 1.000 1.000
P value adjustment method: bonferroni
デフォルトはプールした標準偏差を使う,つまり全ての群で等分散を仮定しているが,pool.sd=Fとすると各群での標準偏差を使用する(等分散を仮定しない)。等分散を仮定しないt検定とはつまりWelchの検定であり,上の出力の後半は,各ペアでウェルチの検定を行って計算された有意確率に検定の回数(ここでは6回)をかけた結果が出力されている。例えばjob=1とjob=2だけでWelchの検定を行うと有意確率は.005345494となる。これに全ての検定のペアの数6をかけると,.03207297となり,上記のボンフェローニ法の結果の1と2の有意確率0.032に一致する。
# job=1とjob=2だけでウェルチ検定
welch01 <- with(data01, t.test(fincome[job==1 | job==2] ~ job[job==1 | job==2], var.equal=F))
welch01
Welch Two Sample t-test
data: fincome[job == 1 | job == 2] by job[job == 1 | job == 2]
t = 2.8242, df = 159.32, p-value = 0.005345
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
51.88747 293.21970
sample estimates:
mean in group 1 mean in group 2
804.3103 631.7568
names(welch01)
[1] "statistic" "parameter" "p.value" "conf.int" "estimate"
[6] "null.value" "stderr" "alternative" "method" "data.name"
welch01$"p.value"*6
[1] 0.03207297
本文で触れた,無調整の検定結果を一応示しておく。
# 共通の標準偏差を推定して検定(等分散性の仮定),有意水準の調整なし
with(data01, pairwise.t.test(fincome, job, p.adjust.method="none"))
Pairwise comparisons using t tests with pooled SD
data: fincome and job
1 2 3
2 0.0094 - -
3 0.0690 0.9891 -
4 0.0803 0.6498 0.7145
P value adjustment method: none
指定しうる多重比較の方法は,ボンフェローニ“bonferroni”,ホルム“holm”, ホッホベルク“hochberg”, ホンメル“hommel”, ベンジャミニとホッホベルク“BH”, ベンジャミニとイェクティエリ“BY”である。bonferroni,holm,hochberg,hommelは“family-wise error rate(FWER)”を調整するものであるが,BHとBYは“false discovery rate(FDR)”を統制するものである。以下のサイトの解説が分かり易い。
『大阪大学大学院医学系研究科 老年・腎臓内科学 腎臓内科』の多重比較解説(2019年9月3日閲覧)
hommel・等分散の例と,BH・不等分散の実行例を示しておく。
with(data01, pairwise.t.test(fincome, job, p.adjust.method="hommel"))
Pairwise comparisons using t tests with pooled SD
data: fincome and job
1 2 3
2 0.056 - -
3 0.276 0.989 -
4 0.321 0.989 0.989
P value adjustment method: hommel
with(data01, pairwise.t.test(fincome, job, p.adjust.method="BH"))
Pairwise comparisons using t tests with pooled SD
data: fincome and job
1 2 3
2 0.056 - -
3 0.161 0.989 -
4 0.161 0.857 0.857
P value adjustment method: BH
Tukey法は,例1の2行で実行し信頼区間を図示できるが,より見易くする為にオブジェクトの代入を活用し,グラフに少しオプションを指定したのが例2である。
# Tukey法の実行と信頼区間の作図
TukeyHSD(aov(fincome ~ JOB, data01))
plot(TukeyHSD(aov(fincome ~ JOB, data01)))
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = fincome ~ JOB, data = data01)
$JOB
diff lwr upr p adj
2 非正規-1 正規 -172.553588 -343.0008 -2.106415 0.0459821
3 自営-1 正規 -173.953202 -420.2323 72.325945 0.2636155
4 無職-1 正規 -133.902182 -331.0475 63.243119 0.2973623
3 自営-2 非正規 -1.399614 -266.1111 263.311891 0.9999991
4 無職-2 非正規 38.651407 -181.0871 258.389960 0.9687159
4 無職-3 自営 40.051020 -242.5905 322.692510 0.9831950
ANOVA01 <- aov(fincome ~ JOB, data01)
TukeyHSD(ANOVA01)
par(mar = c(5, 7, 4, 2))
plot(TukeyHSD(ANOVA01), las=1)
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = fincome ~ JOB, data = data01)
$JOB
diff lwr upr p adj
2 非正規-1 正規 -172.553588 -343.0008 -2.106415 0.0459821
3 自営-1 正規 -173.953202 -420.2323 72.325945 0.2636155
4 無職-1 正規 -133.902182 -331.0475 63.243119 0.2973623
3 自営-2 非正規 -1.399614 -266.1111 263.311891 0.9999991
4 無職-2 非正規 38.651407 -181.0871 258.389960 0.9687159
4 無職-3 自営 40.051020 -242.5905 322.692510 0.9831950
次のスクリプトの一行目がデフォルトの(等分散性を仮定しない)一要因分散分析,二行目が等分散の仮定を指定した一要因分散分析である。後者の結果は,先に紹介したaov( )やlm( )と一致している事を見比べよう。
# 一元配置分散分析
oneway.test(fincome ~ job, data01)
One-way analysis of means (not assuming equal variances)
data: fincome and job
F = 3.2034, num df = 3.00, denom df = 89.75, p-value = 0.02699
oneway.test(fincome ~ job, data01, var.equal=T)
One-way analysis of means
data: fincome and job
F = 3.0496, num df = 3, denom df = 292, p-value = 0.02898
本文では簡単に触れているだけの等分散性の検定を実際に行ってみる。標準で存在しているバートレット検定bartlett.test( )
と,追加パッケイジ“car”に入っているルヴィーン検定(レヴィン検定)の二つを行ってみよう。“car”パッケイジは,インストールした事がなければ最初にインストールし,その後も使用する際にはlibrary(car)
として有効化しなければならない。また,bartlettはグループ変数がfactor型になっていなくても自動で要因型だと見做して分析するが,leveneTest( )
はfactor型にしなければならない。
# 等分散性検定
bartlett.test(fincome ~ job, data01)
Bartlett test of homogeneity of variances
data: fincome by job
Bartlett's K-squared = 4.2245, df = 3, p-value = 0.2382
# 初めての場合は,"car"パッケイジをインストール
# install.packages("car", repos="http://cran.ism.ac.jp/")
library(car)
Loading required package: carData
leveneTest(fincome ~ JOB, data01)
二要因分散分析でも,(偏差)平方和SSの分解が基本であることには変わらないが,各グループのサイズが異なる場合(非釣り合い型データ,アンバランスデータ,非直交なデータ)には分解の仕方が複数存在し,標準のaov( )
関数はその点で問題を生じる。まずはデータの準備をする。分割表でも分かる様に,セルのサイズ(グループのサイズ)はまるでばらばらである。
※ スクリプトの書き方を少し変更した(2020.06.21)。
## 二元配置の分散分析
data01$SEX <- factor(data01$sex, levels=1:2,
labels=c("1 男性", "2 女性")) # 独立変数その1
with(data01, table("性別"=SEX, "従業上の地位"=JOB, useNA="ifany"))
従業上の地位
性別 1 正規 2 非正規 3 自営 4 無職 <NA>
1 男性 127 9 19 7 3
2 女性 55 83 19 61 1
次に,それぞれの一要因分散分析と,性別,従業上の地位の順に変数を並べた二要因分散分析の結果を表示する。一要因の場合とは,各独立変数(主効果)の平方和も有意確率も異なる。
anova(model30 <- aov(fincome ~ SEX, data01)) # 性別の一要因分散分析
anova(model31 <- aov(fincome ~ JOB, data01)) # 従業上の地位の一要因分散分析
anova(model32 <- aov(fincome ~ SEX + JOB, data01)) # 二要因分散分析,その1
ここで,独立変数の並べ順を変えてみよう。変数の並べ順を変えただけで,それぞれの独立変数の偏差平方和も有意確率もすっかり変わってしまった。
anova(model33 <- aov(fincome ~ JOB + SEX, data01)) # 二要因分散分析,その2
ここで,“car”パッケイジのAnova( )
関数を利用する。パッケイジのインストールには上の1-5を参照。model32とmodel33は既に上で作成した分散分析のモデルである。この,投入順を変えた二つのモデルで,AnovaのタイプⅡ平方和の計算による分散分析を行った結果は,平方和も有意確率も全て一致していて,投入順に左右されない事が確認出来る。因みに,Anova( )関数でタイプⅢの平方和を用いる事も出来る(type=3というオプションを付ける)。
library("car")
Anova(model32)
Anova(model33)
model32,model33自体は以下の通りであり,Anova( )関数を使わなければtypeⅡの計算は行われない。
model32
Call:
aov(formula = fincome ~ SEX + JOB, data = data01)
Terms:
SEX JOB Residuals
Sum of Squares 84061 3428643 60692395
Deg. of Freedom 1 3 291
Residual standard error: 456.6891
Estimated effects may be unbalanced
88 observations deleted due to missingness
model33
Call:
aov(formula = fincome ~ JOB + SEX, data = data01)
Terms:
JOB SEX Residuals
Sum of Squares 1950508 1562196 60692395
Deg. of Freedom 3 1 291
Residual standard error: 456.6891
Estimated effects may be unbalanced
88 observations deleted due to missingness
各グループのケース数や年収の平均といった基本的な情報を求めるには以下の様にすれば良い。
どの変数にも欠損値が無い様にする為に少し面倒になっているが,この後の「完備ケース分析」にしてしまえばその面倒は無くなる。
with(data01, tapply(fincome[!is.na(JOB)], SEX[!is.na(JOB)], mean, na.rm=T)) # 性別での平均
1 男性 2 女性
702.9762 737.0588
with(data01, tapply(fincome[!is.na(SEX)], JOB[!is.na(SEX)], mean, na.rm=T)) # 従業上の地位での平均
1 正規 2 非正規 3 自営 4 無職
804.3103 631.7568 630.3571 670.4082
with(data01, table(SEX[!is.na(fincome)&!is.na(JOB)], useNA="ifany")) # 有効ケース数
1 男性 2 女性
126 170
次の様に,有効ケースのみを選び出してから集計・分析する事も出来る。(処理をより工夫した。2020.06.21)
Full <- with(data01, complete.cases(fincome, SEX, JOB))
d02 <- data01[Full, c("fincome", "SEX", "JOB")]
with(d02, tapply(fincome, SEX, mean))
1 男性 2 女性
702.9762 737.0588
with(d02, tapply(fincome, JOB, mean))
1 正規 2 非正規 3 自営 4 無職
804.3103 631.7568 630.3571 670.4082
with(d02, addmargins(table(SEX, JOB)))
JOB
SEX 1 正規 2 非正規 3 自営 4 無職 Sum
1 男性 101 6 14 5 126
2 女性 44 68 14 44 170
Sum 145 74 28 49 296
ここで,効果サイズ(effect size)の指標である分散説明率について紹介しておこう。
分散分析では伝統的に\(\eta^2\)(イータの二乗)と云う単語・統計量が使用され,回帰分析では\(R^2\)(決定係数)が用いられてきた。R2は通常,モデル全体の説明率を示すのに用いられ,\(\eta^2\)は個々の要因の説明率を示すのに用いられる。(重回帰分析では偏回帰係数のt検定は分散説明率の増加分についてのF検定と同値であるので,実際には回帰分析でも個々の要因の説明率を示すのにも使用可能である。)
世帯収入を,性別と従業上の地位で説明するモデルを,aov( )
とlm( )
で実行し,anova( )
やAnova( )
,heplotsパッケイジ中のetasq( )
関数に与えて出力結果を確認する。
## 各種の分散説明率
out01 <- aov(fincome ~ SEX * JOB, data01) # aovの結果
out02 <- lm(fincome ~ SEX * JOB, data01) # lmの結果
anova(out01) # type I SS
sum(anova(out01)$"Sum Sq"[1:3])/sum(anova(out01)$"Sum Sq") # R2を手計算
[1] 0.0840703
summary(out02) # lm()の結果からR-squared
Call:
lm(formula = fincome ~ SEX * JOB, data = data01)
Residuals:
Min 1Q Median 3Q Max
-857.39 -253.31 -53.31 190.34 1776.79
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 748.51 44.96 16.647 < 2e-16 ***
SEX2 女性 183.87 81.62 2.253 0.02503 *
JOB2 非正規 -361.01 189.88 -1.901 0.05826 .
JOB3 自営 -25.30 128.87 -0.196 0.84449
JOB4 無職 -643.51 207.03 -3.108 0.00207 **
SEX2 女性:JOB2 非正規 81.94 209.04 0.392 0.69537
SEX2 女性:JOB3 自営 -369.59 189.30 -1.952 0.05186 .
SEX2 女性:JOB4 無職 445.79 228.35 1.952 0.05188 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 451.9 on 288 degrees of freedom
(88 observations deleted due to missingness)
Multiple R-squared: 0.08407, Adjusted R-squared: 0.06181
F-statistic: 3.776 on 7 and 288 DF, p-value: 0.0006164
car::Anova(out01) # type II SS
heplots::etasq(out01, anova=T) # 偏イータの二乗
heplots::etasq(out01, anova=T, partial=F) # (古典的)イータの二乗
以下に,偏イータの二乗とイータの二乗の計算の仕方を実行例とともに示す。解読してみよう。定義上必ず偏イータの二乗の方が大きくなる。
# 偏イータの二乗の計算方法
eta01 <- heplots::etasq(out01, anova=T); eta01
eta01$"Sum Sq"[c(1:3)]/(eta01$"Sum Sq"[c(1:3)] + eta01$"Sum Sq"[4])
[1] 0.02587721 0.05509100 0.03105888
# (古典的)イータの二乗の計算方法
eta02 <- heplots::etasq(out01, anova=T, partial=F); eta02
eta02$"Sum Sq"[c(1:3)]/sum(eta02$"Sum Sq"[1:4])
[1] 0.02378378 0.05219967 0.02869892
複数の要因(独立変数)が複合的に従属変数に影響を及ぼす事を,計量分析では「交互作用intraction」と呼ぶが,質的比較分析(QCA, Qualitative Comparative Analysis)でいう「結合因果conjunctural causation」と同種の考えであると言ってよいだろう。
以下,本文のスクリプトとその結果を,より装飾を施した形で紹介してゆく。
# 性別×従業上の地位によって8つの群に分けた時の世帯年収平均
# 本文では小数点以下が表示されているが,不要なので工夫する。
with(data01, round(tapply(fincome, list(SEX, JOB), mean, na.rm=T), 0))
1 正規 2 非正規 3 自営 4 無職
1 男性 749 388 723 105
2 女性 932 653 538 735
# 交互作用項を含むモデルの二つの書き方
aov(fincome ~ SEX + JOB + SEX:JOB, data01)
Call:
aov(formula = fincome ~ SEX + JOB + SEX:JOB, data = data01)
Terms:
SEX JOB SEX:JOB Residuals
Sum of Squares 84061 3428643 1885038 58807357
Deg. of Freedom 1 3 3 288
Residual standard error: 451.8763
Estimated effects may be unbalanced
88 observations deleted due to missingness
aov(fincome ~ SEX * JOB, data01)
Call:
aov(formula = fincome ~ SEX * JOB, data = data01)
Terms:
SEX JOB SEX:JOB Residuals
Sum of Squares 84061 3428643 1885038 58807357
Deg. of Freedom 1 3 3 288
Residual standard error: 451.8763
Estimated effects may be unbalanced
88 observations deleted due to missingness
# type I平方和による分散分析
anova(aov(fincome ~ SEX * JOB, data01))
anova(aov(fincome ~ JOB * SEX, data01))
# type II平方和による分散分析
Anova(aov(fincome ~ SEX * JOB, data01))
# type III平方和による分散分析
Anova(aov(fincome ~ SEX * JOB, data01), type=3)
平方和の分解方法によって結果が異なる事を確認しておこう。
最後に,交互作用プロットの描き方を紹介する。本文のと類似のグラフが以下のスクリプトで描ける。
with( )でくるんだり,描画する関数を fun=function(x) mean(x, na.rm=T) と指定したり,幾つか工夫を加えた(2020.06.21)。
with(data01,
interaction.plot(JOB, SEX, fincome, bty="l",
fun=function(x) mean(x, na.rm=T),
xlab="従業上の地位", ylab="平均世帯年収(万円)", family="serif",
main="交互作用プロット", ylim=c(0,1000), col=c("blue", "red"))
)
平均は先に計算した通りだが,交互作用プロットの様なグラフではそれぞれのグループのケース数(グループサイズ)が全く分からないので,グループ毎の有効ケース数も計算しておくのが良い。
with(data01, round(tapply(fincome, list(SEX, JOB), mean, na.rm=T), 0))
1 正規 2 非正規 3 自営 4 無職
1 男性 749 388 723 105
2 女性 932 653 538 735
with(data01, tapply(fincome, list(SEX, JOB), function(x) sum(!is.na(x))))
1 正規 2 非正規 3 自営 4 無職
1 男性 101 6 14 5
2 女性 44 68 14 44
性別による幸福度の平均の差の分析結果を再掲しておこう。
x <- data01$sex
y <- c(0:10)[data01$q1700+1]
by(y, x, mean, na.rm=T)
x: 1
[1] 6.36646
------------------------------------------------------------
x: 2
[1] 6.917051
by(y, x, var, na.rm=T)
x: 1
[1] 4.083618
------------------------------------------------------------
x: 2
[1] 3.08568
by(y, x, function(x) sum(!is.na(x)))
x: 1
[1] 161
------------------------------------------------------------
x: 2
[1] 217
var.test(y ~ x)
F test to compare two variances
data: y by x
F = 1.3234, num df = 160, denom df = 216, p-value = 0.05555
alternative hypothesis: true ratio of variances is not equal to 1
95 percent confidence interval:
0.9933775 1.7743088
sample estimates:
ratio of variances
1.323409
t.test(y ~ x, var.equal=T)
Two Sample t-test
data: y by x
t = -2.8252, df = 376, p-value = 0.004977
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.9337912 -0.1673909
sample estimates:
mean in group 1 mean in group 2
6.366460 6.917051
これを,一要因分散分析で行ってみる。この例では等分散を仮定出来る為,関数は3種類あるので全て実行して見せよう。独立変数であるx(元はdata01$sex)はfactor型にすべきだが,二値変数の場合には結果に本質的な違いはない。
oneway.test(y ~ x, var.equal=T)
One-way analysis of means
data: y and x
F = 7.9819, num df = 1, denom df = 376, p-value = 0.004977
anova(aov(y ~ x)) # anova(aov(y ~ factor(x))) と一致する
anova(lm(y ~ x)) # anova(lm(y ~ factor(x)))と一致する
t検定の結果とF検定の結果をそれぞれオブジェクトに格納し,並べて表示させよう。更に,オブジェクトにどんな情報が格納されているのかも一覧を確認しよう。
T01 <- t.test(y ~ x, var.equal=T)
F01 <- anova(aov(y ~ x))
F02 <- oneway.test(y ~ x, var.equal=T)
T01; names(T01)
Two Sample t-test
data: y by x
t = -2.8252, df = 376, p-value = 0.004977
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.9337912 -0.1673909
sample estimates:
mean in group 1 mean in group 2
6.366460 6.917051
[1] "statistic" "parameter" "p.value" "conf.int" "estimate"
[6] "null.value" "stderr" "alternative" "method" "data.name"
F01; names(F01)
[1] "Df" "Sum Sq" "Mean Sq" "F value" "Pr(>F)"
F02; names(F02)
One-way analysis of means
data: y and x
F = 7.9819, num df = 1, denom df = 376, p-value = 0.004977
[1] "statistic" "parameter" "p.value" "method" "data.name"
t分布の自由度とF分布の残差(級内)の自由度,有意確率同士が対応している事が分かる。
t統計量とF統計量の関係は,次の様に確認してみよう。
T01$statistic^2
t
7.981851
F01$"F value"
[1] 7.981851 NA
F02$statistic
F
7.981851
t分布はマイナス領域とプラス領域がある。二乗すれば同じ値になる点が正負の二つある。F分布は全て非負である。どの様に対応しているだろうか。下の様に確認すると,tが負の時と正の時の確率を合わせたものが,(tの二乗に等しい)Fの上側確率に等しい事が分かる。tの両側検定はFの上側検定に対応しているのである。
# t分布の下側確率と上側確率の和
1 - pt(abs(T01$statistic), df=T01$parameter) + pt(-1*abs(T01$statistic), df=T01$parameter)
t
0.004976846
# F分布の上側確率
1 - pf(F01$"F value"[1], df1=F01$Df[1], df2=F01$Df[2])
[1] 0.004976846
等分散を仮定出来ずにウェルチの検定を行った,男女での世帯年収差はどうであろうか。
x <- data01$sex
y <- data01$fincome
T01 <- t.test(y ~ x)
F02 <- oneway.test(y ~ x)
T01
Welch Two Sample t-test
data: y by x
t = -0.756, df = 294.6, p-value = 0.4503
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-142.7534 63.5176
sample estimates:
mean in group 1 mean in group 2
697.4409 737.0588
F02
One-way analysis of means (not assuming equal variances)
data: y and x
F = 0.57153, num df = 1.0, denom df = 294.6, p-value = 0.4503
T01$statistic^2
t
0.5715288
F02$statistic
F
0.5715288
これもやはり,自由度,有意確率が一致し,t統計量の2乗とF統計量が一致している。anova(aov( ))
関数では等分散性が仮定出来ない分散分析を行う事は出来ない。
問題文で要求されていないコマンドも含んだスクリプトを掲載する。
x <- factor(data01$edu2, levels=1:3, labels=c("高校", "短大", "四大"))
y <- data01$income
# まずは変数の分布を確認しておく
table(x, useNA="always")
x
高校 短大 四大 <NA>
128 123 125 8
table(y, useNA="always")
y
0 25 75 125 200 300 400 500 600 700 800 925 1125 1375 1750 2500
45 24 28 29 40 35 31 29 27 18 12 13 6 1 3 1
<NA>
42
# カテゴリー別の平均値,ケース数を求める
by(y, x, mean, na.rm=T)
x: 高校
[1] 283.6777
------------------------------------------------------------
x: 短大
[1] 236.4865
------------------------------------------------------------
x: 四大
[1] 536.7925
by(y, x, function(x) sum(!is.na(x)))
x: 高校
[1] 121
------------------------------------------------------------
x: 短大
[1] 111
------------------------------------------------------------
x: 四大
[1] 106
# 等分散性の検定
bartlett.test(y ~ x)
Bartlett test of homogeneity of variances
data: y by x
Bartlett's K-squared = 33.378, df = 2, p-value = 5.651e-08
# 等分散を仮定しない分散分析と,仮定する分散分析
oneway.test(y ~ x)
One-way analysis of means (not assuming equal variances)
data: y and x
F = 20.878, num df = 2.00, denom df = 210.05, p-value = 5.369e-09
anova(aov(y ~ x))
# 多重比較
TukeyHSD(aov(y ~ x))
Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = y ~ x)
$x
diff lwr upr p adj
短大-高校 -47.1912 -144.5816 50.19921 0.4897002
四大-高校 253.1148 154.5337 351.69587 0.0000000
四大-短大 300.3060 199.6726 400.93929 0.0000000
まず,anova(aov( ))では独立変数の投入順が結果に影響する事を確認しておく(タイプⅠ平方和)。
“car”パッケイジのAnova( )関数(タイプⅡ平方和)ではそうした違いが生じない事を次に確認し,交互作用を含んだモデルも実行する。そして最後に交互作用プロットを描く。
x2 <- factor(data01$sex, levels=1:2, labels=c("男性", "女性"))
# anova(aov( ))では投入順で結果が変わる
anova(aov(y ~ x + x2))
anova(aov(y ~ x2 + x))
# install.packages("car", repos="https://cran.ism.ac.jp")
library("car")
Anova(aov(y ~ x + x2))
Anova(aov(y ~ x2 + x)) # "car"のAnova( )なら順番に左右されない
Anova(aov(y ~ x * x2)) # 交互作用項を含める
MAX <- max(tapply(y, list(x, x2), mean, na.rm=T)) # y軸の最大値を自動で設定させる
interaction.plot(x, x2, y, bty="l",
fun=function(x) mean(x, na.rm=T),
xlab="最終学歴", ylab="平均個人年収(万円)", family="serif",
main="交互作用プロット", ylim=c(0, MAX), col=c("blue", "red"))
交互作用を含むモデルを取り上げ,aov( ),anova( ),Anova( )の結果も比較してみる。
model1 <- aov(y ~ x * x2)
model1 # aov( )のオブジェクトの主な内容
Call:
aov(formula = y ~ x * x2)
Terms:
x x2 x:x2 Residuals
Sum of Squares 5658984 7807360 93424 25290094
Deg. of Freedom 2 1 2 332
Residual standard error: 275.9982
Estimated effects may be unbalanced
46 observations deleted due to missingness
names(model1)
[1] "coefficients" "residuals" "effects" "rank"
[5] "fitted.values" "assign" "qr" "df.residual"
[9] "na.action" "contrasts" "xlevels" "call"
[13] "terms" "model"
anova(model1) # anova( )の主な内容
names(anova(model1))
[1] "Df" "Sum Sq" "Mean Sq" "F value" "Pr(>F)"
Anova(model1) # Anova( )の主な内容
names(Anova(model1))
[1] "Sum Sq" "Df" "F value" "Pr(>F)"
model1の“coefficients”を表示させると以下の通りである。
model1$coefficients
(Intercept) x短大 x四大 x2女性 x短大:x2女性 x四大:x2女性
462.7551020 61.3189720 195.9750567 -300.9495465 -79.0769085 0.3589226
カテゴリごとの平均は以下の通りである。
tapply(y, list(x, x2), mean)
男性 女性
高校 NA 161.8056
短大 NA NA
四大 NA NA
これをよく見比べると,model1$coefficients
の(Intercept)は,高校―男性カテゴリの平均個人年収に一致している。また,「x2女性」の-300.9495465は,高校―男性と高校―女性の差にぴったり一致している。「x短大」の61.3189720は高校―男性と短大―男性の差にぴったり一致する。最終学歴においては「高校」が,性別においては「男性」が「基準カテゴリ」となっており,そこからの差が“coefficients”として計算されているのである。これは「ダミー変数を用いて,交互作用項も投入した重回帰分析」と同じである。
ちょっと難しいかも知れないが,この coefficients を次の様に合計すると,各カテゴリの平均が求められる。いずれも最後の列の Sum が,上で調べた各カテゴリの平均年収に一致している事を確認しよう。
b0 <- model1$coefficients # 係数をオブジェクトに格納する
b <- matrix(c(rep(b0[1], 3), 0, b0[2:3], rep(b0[4], 3), 0, b0[5:6]), ncol=4)
rownames(b) <- c("高校", "短大", "四大")
colnames(b) <- c("基準", "学歴", "女性", "女性×学歴")
b # 係数行列から,基準値(高校―男性),学歴の主効果,性別の主効果,学歴×性別の交互作用効果を構成
基準 学歴 女性 女性×学歴
高校 462.7551 0.00000 -300.9495 0.0000000
短大 462.7551 61.31897 -300.9495 -79.0769085
四大 462.7551 195.97506 -300.9495 0.3589226
addmargins(b[, 1:2])[1:3,] # 男性の各カテゴリはこの様に計算する
基準 学歴 Sum
高校 462.7551 0.00000 462.7551
短大 462.7551 61.31897 524.0741
四大 462.7551 195.97506 658.7302
addmargins(b)[1:3,] # 女性の各カテゴリはこの様に計算する
基準 学歴 女性 女性×学歴 Sum
高校 462.7551 0.00000 -300.9495 0.0000000 161.8056
短大 462.7551 61.31897 -300.9495 -79.0769085 144.0476
四大 462.7551 195.97506 -300.9495 0.3589226 358.1395
以下の様にも計算出来る。この方がダミー変数を投入した重回帰分析の考え方をよく表現している。
sum(b0*c(1, 0, 0, 0, 0, 0)) # 男性―高校
[1] 462.7551
sum(b0*c(1, 1, 0, 0, 0, 0)) # 男性―短大・高専
[1] 524.0741
sum(b0*c(1, 0, 1, 0, 0, 0)) # 男性―四大
[1] 658.7302
sum(b0*c(1, 0, 0, 1, 0, 0)) # 女性―高校
[1] 161.8056
sum(b0*c(1, 1, 0, 1, 1, 0)) # 女性―短大・高専
[1] 144.0476
sum(b0*c(1, 0, 1, 1, 0, 1)) # 女性―四大
[1] 358.1395