『入門・社会統計学――2ステップで基礎から〔Rで〕学ぶ』

第7章　平均値の差の分散分析(ANOVA)

0 全体の章構成

• 第 1章　1変量の記述統計の基礎
• 第 2章　2変数の関連の記述統計
• 第 3章　推測統計の基礎
• 第 4章　統計的検定の一般型(general form of statistical testing)
• 第 5章　2変数の関連の推定と検定
• 第 6章　2群の母平均の差のt検定
• 第 7章　平均値の差の分散分析(ANOVA)
• 第 8章　単回帰分析(SRA)
• 第9章　重回帰分析（Ⅰ）(MRA)
• 第10章　重回帰分析（Ⅱ）(MRA)
• 第11章　主成分分析と因子分析
• 第12章　ロジスティック回帰分析
• 第13章　階層構造のあるデータの分析

1　1要因分散分析

グループ間の母平均の差を検定する分散分析，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])

g1 <- ggplot(d01, aes(x = group, y = score)) +
         geom_boxplot(aes(color = group))   +
         geom_jitter(alpha = .4, aes(color = group)) + 
         theme_classic() + ylim(min_y, max_y) +
         theme(legend.position = "none") +
         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")

g2 <- ggplot(d02, aes(x = group, y = score)) + 
         geom_boxplot(aes(color = group)) + 
         geom_jitter(alpha = .4, aes(color = group)) + 
         theme_classic() + ylim(min_y, max_y) + 
         theme(legend.position = "none") + 
         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="")

グラフ

実行するたびに結果が変わるので，例としては余り良くない場合もある。
（後から全く同じ結果を再現したい場合には，最初に set.seed( ) を使う。）

gridExtra::grid.arrange(g1, g2, nrow = 1)

対応するF分布を描く

スクリプト

df1 <- ano1$Df[1]
df2 <- ano1$Df[2]

p <- .05; upr <- qf(.99, df1, df2)

curve(df(x, df1, df2), 
      xlim = c(0, upr), 
      bty  = "n", 
      xlab = "", 
      ylab = "", 
      xaxt = "n", 
      col  = "#339900",
      main = paste("対応するF分布(df=", df1, ", ", df2, ")の確率", 
                   round(p*100, 2), "%領域"))
abline(h = 0, col = "#33990080")
segments(x0 <- seq(qf(1-.05, df1, df2), upr, by = .01), 0, 
         x0, df(x0, df1, df2), 
         col = "#ff3333")
axis(side = 1, 
     at = seq(0, upr, by = 0.5), 
     col.axis = "#555555", 
     cex.axis = 0.8, 
     las = 1)
axis(side = 1, 
     at = c(round(qf(1-.05, df1, df2), 2)), 
     col.axis = "#ff0000", 
     cex.axis = 1,
     las = 2)

1-1　偏差平方和(Sum of Squares)の分解

　グループ 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}\]

1-2　ゼロ仮説(Null Hypothesis; H0)とF統計量

偏差平方和の分解の式
\[\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 統計量によってゼロ仮説を検定することが可能になる。

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), xlim = c(0, 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), xlim = c(0, 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), 
      xlim = c(0, 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")

1-3　4群の母平均差のF検定

使用する変数についての基本的確認

　まずは使用する変数の基本的な情報（度数分布や要約統計量，ケース数）を確認しておく必要がある。欠損値NAが含まれる場合はひと手間余計にかかる事が多いので注意しよう。また，本文と異なり，データフレイムの名前はdata01としている。
　また，いちいちデータフレイム名 data01 を書く事を省略したい場合には with( ) を使用している。

data01 <- read.csv("practice.csv")

# まずは各変数の度数分布または基本要約統計量を確認
with(data01, table(job, useNA = "always"))

job
   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( )関数のうち，まずは最初の二つを説明する。

■ 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("正規", "非正規", "自営", "無職"))

with(data01, table(job, JOB, useNA = "ifany")) # 新変数の確認

      JOB
job    正規 非正規 自営 無職 <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)
ANOVA01

Call:
   aov(formula = fincome ~ JOB, data = data01)

Terms:
                     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

結果の表示

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)

　aov( )の出力結果を用いて，「グループ平均はすべて等しい」というゼロ仮説の有意性検定を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( )関数

間違い

これも故意に間違い。jobが数量型になっている。

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(lm( ))

因みにsummary( )関数の出力はこうなる。これは回帰分析の時に重要になってくる。

summary(LM01)

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 ***
JOB非正規    -172.55      65.97  -2.616  0.00936 ** 
JOB自営      -173.95      95.31  -1.825  0.06901 .  
JOB無職      -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)のＦ分布の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 平方和の分解の末尾に例示を付け加えておく。

1-4　多重比較(multiple comparison)，事後検定(post-hoc test）

　一つ一つの検定の有意水準をp1，そうした検定を3つ繰り返した場合の最終的な有意水準をp3とすると，p3 = p1 + p1 + p1 - 3*p1*p1 + p1*p1*p1 である。
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回）をかけた結果が出力されている。

等分散を仮定しない→Welch

例えばjob = 1とjob = 2だけでWelchの検定を行うと有意確率は.005345494となる。これに全ての検定のペアの数6をかけると，.03207297となり，上記のボンフェローニ法の結果の1と2の有意確率0.032に一致する。

# job=1とjob=2だけでウェルチ検定
welch01 <- with(data01[is.element(data01$job, c(1, 2)), ], 
                t.test(fincome ~ job, 
                       var.equal = F))
welch01

    Welch Two Sample t-test

data:  fincome by job
t = 2.8242, df = 159.32, p-value = 0.005345
alternative hypothesis: true difference in means between group 1 and group 2 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

non-adjustment

本文で触れた，無調整の検定結果を一応示しておく。

# 共通の標準偏差を推定して検定（等分散性の仮定），有意水準の調整なし
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・不等分散の実行例を示しておく。

Hommel, 等分散

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 

BenjaminiとHochberg, 不等分散

with(data01, 
     pairwise.t.test(fincome, job, 
                     p.adjust.method = "BH",
                     pool.sd = F))

    Pairwise comparisons using t tests with non-pooled SD 

data:  fincome and job 

  1     2     3    
2 0.032 -     -    
3 0.243 0.989 -    
4 0.243 0.886 0.886

P value adjustment method: BH 

テューキー(Tukey)の多重比較

　Tukey法は，例1の2行で実行し信頼区間を図示できるが，より見易くする為にオブジェクトの代入を活用し，グラフに少しオプションを指定したのが例2である。

例1（最低限）

Tukey法の実行と信頼区間の作図

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
非正規-正規 -172.553588 -343.0008  -2.106415 0.0459821
自営-正規   -173.953202 -420.2323  72.325945 0.2636155
無職-正規   -133.902182 -331.0475  63.243119 0.2973623
自営-非正規   -1.399614 -266.1111 263.311891 0.9999991
無職-非正規   38.651407 -181.0871 258.389960 0.9687159
無職-自営     40.051020 -242.5905 322.692510 0.9831950

plot(TukeyHSD(aov(fincome ~ JOB, data01)))

例2（オプション付き）

ANOVA01 <- aov(fincome ~ JOB, data01)
TukeyHSD(ANOVA01)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = fincome ~ JOB, data = data01)

$JOB
                   diff       lwr        upr     p adj
非正規-正規 -172.553588 -343.0008  -2.106415 0.0459821
自営-正規   -173.953202 -420.2323  72.325945 0.2636155
無職-正規   -133.902182 -331.0475  63.243119 0.2973623
自営-非正規   -1.399614 -266.1111 263.311891 0.9999991
無職-非正規   38.651407 -181.0871 258.389960 0.9687159
無職-自営     40.051020 -242.5905 322.692510 0.9831950

par(mar = c(5, 7, 4, 2))
plot(TukeyHSD(ANOVA01), las = 1)

1-5　等分散性の前提とウェルチ検定

one-way( )

　次のスクリプトの一行目がデフォルトの（等分散性を仮定しない）一要因分散分析，二行目が等分散の仮定を指定した一要因分散分析である。
後者の結果は，先に紹介した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")
library(car)

Loading required package: carData

leveneTest(fincome ~ JOB, data01)

2　2要因分散分析

2-1　平方和の分解

　二要因分散分析でも，（偏差）平方和SSの分解が基本であることには変わらないが，各グループのサイズが異なる場合（非釣り合い型データ，アンバランスデータ，非直交なデータ）には分解の仕方が複数存在し，標準のaov( )関数はその点で問題を生じる。
　まずはデータの準備をする。分割表でも分かる様に，セルのサイズ（グループのサイズ）はまるでばらばらである。
また，投入する変数によって，有効ケース数が変化してしまう（後で対処する）。

## 二元配置の分散分析
data01$SEX <- factor(data01$sex, levels = 1:2,
                     labels = c("男性", "女性")) # 独立変数その1

with(data01, table("性別" = SEX, "従業上の地位" = JOB, useNA="ifany"))

      従業上の地位
性別 正規 非正規 自営 無職 <NA>
  男性  127      9   19    7    3
  女性   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( )

　ここで，“car”パッケイジのAnova( )関数を利用する（パッケイジのインストールには上の1-5を参照）。
model32とmodel33は既に上で作成した分散分析のモデルである。
この，投入順を変えた二つのモデルで，AnovaのタイプⅡ平方和の計算による分散分析を行った結果は，平方和も有意確率も全て一致していて，投入順に左右されない事が確認出来る。
因みに，Anova( )関数でタイプⅢの平方和を用いる事も出来る（type= 3 というオプションを付ける）。

car::Anova(model32)

car::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[!is.na(data01$SEX) & !is.na(data01$JOB) & !is.na(data01$fincome), ], 
     tapply(fincome, SEX, mean, na.rm = T))

    男性     女性 
702.9762 737.0588 

# 従業上の地位での平均
with(data01[!is.na(data01$SEX) & !is.na(data01$JOB) & !is.na(data01$fincome), ], 
     tapply(fincome, JOB, mean, na.rm = T))

    正規   非正規     自営     無職 
804.3103 631.7568 630.3571 670.4082 

# 性別と従業上の地位での平均
with(data01[!is.na(data01$SEX) & !is.na(data01$JOB) & !is.na(data01$fincome), ], 
     tapply(fincome, list(SEX, JOB), mean, na.rm = T))

         正規   非正規     自営     無職
男性 748.5149 387.5000 723.2143 105.0000
女性 932.3864 653.3088 537.5000 734.6591

# 有効ケース数
with(data01[!is.na(data01$SEX) & !is.na(data01$JOB) & !is.na(data01$fincome), ], 
     addmargins(table(SEX, useNA = "ifany")))

SEX
男性 女性  Sum 
 126  170  296 

完備ケース分析 complete.cases( )

　次の様に，有効ケースのみを選び出してから集計・分析する事も出来る。

Full <- with(data01, complete.cases(fincome, SEX, JOB))
d02 <- data01[Full, c("fincome", "SEX", "JOB")]

with(d02, tapply(fincome, SEX, mean))

    男性     女性 
702.9762 737.0588 

with(d02, tapply(fincome, JOB, mean))

    正規   非正規     自営     無職 
804.3103 631.7568 630.3571 670.4082 

with(d02, addmargins(table(SEX, JOB)))

      JOB
SEX    正規 非正規 自営 無職 Sum
  男性  101      6   14    5 126
  女性   44     68   14   44 170
  Sum   145     74   28   49 296

分散説明率(Proportion of variance explained)

　ここで，効果サイズ(effect size)の指標である分散説明率について紹介しておこう。
分散分析では伝統的に\(\eta^2\)（イータの二乗）と云う単語・統計量が使用され，回帰分析では\(R^2\)（決定係数）が用いられてきた。
\(R^2\)は通常，モデル全体の説明率を示すのに用いられ，\(\eta^2\)は個々の要因の説明率を示すのに用いられる。
（重回帰分析では偏回帰係数のt検定は分散説明率の増加分についてのF検定と同値であるので，実際には回帰分析でも個々の要因の説明率を示すのにも使用可能である。）
　世帯収入を，性別と従業上の地位で説明するモデルを，aov( )とlm( )で実行し，anova( )やAnova( )，heplotsパッケイジ中のetasq( )関数に与えて出力結果を確認する。

anova(aov( ))

## 各種の分散説明率
out01   <- aov(fincome ~ SEX * JOB, data01) # aovの結果
anova01 <- anova(out01)                     # 敢えて type I SS を計算
anova01

sum(anova01$"Sum Sq"[1:3])/sum(anova01$"Sum Sq") # R2を手計算

[1] 0.0840703

summary(lm( ))

out02 <-  lm(fincome ~ SEX * JOB, data01) # lmの結果
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 ***
SEX女性             183.87      81.62   2.253  0.02503 *  
JOB非正規          -361.01     189.88  -1.901  0.05826 .  
JOB自営             -25.30     128.87  -0.196  0.84449    
JOB無職            -643.51     207.03  -3.108  0.00207 ** 
SEX女性:JOB非正規    81.94     209.04   0.392  0.69537    
SEX女性:JOB自営    -369.59     189.30  -1.952  0.05186 .  
SEX女性:JOB無職     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

type II SSと（偏）イータの二乗

car::Anova(out01) # type II SS

heplots::etasq(out01, anova = T) # 偏イータの二乗

heplots::etasq(out01, anova = T, partial = F) # (古典的)イータの二乗

（偏）イータの二乗が何であるかの確認

　以下に，偏イータの二乗とイータの二乗の計算の仕方を実行例とともに示す。
解読してみよう。
定義上必ず偏イータの二乗の方が大きくなる。

# （古典的）イータの二乗の計算方法
eta02 <- heplots::etasq(out01, anova = T, partial = F)
eta02

eta02$"Sum Sq"

[1]  1562196  3428643  1885038 58807357

eta02$"Sum Sq"[c(1:3)]

[1] 1562196 3428643 1885038

sum(eta02$"Sum Sq"[c(1:4)])

[1] 65683234

eta02$"Sum Sq"[c(1:3)]/sum(eta02$"Sum Sq"[c(1:4)])

[1] 0.02378378 0.05219967 0.02869892

※ 2021-01-16の時点ではこの計算結果が一致していたが，2024-11-17時点ではheplots::etasq( )の計算方法が変わったらしく，一致しない。

# 偏イータの二乗の計算方法
eta01 <- heplots::etasq(out01, anova = T)
eta01

eta01$"Sum Sq"

[1]  1562196  3428643  1885038 58807357

eta01$"Sum Sq"[c(1:3)]

[1] 1562196 3428643 1885038

(eta01$"Sum Sq"[c(1:3)] + eta01$"Sum Sq"[4])

[1] 60369553 62236001 60692395

eta01$"Sum Sq"[c(1:3)]/(eta01$"Sum Sq"[c(1:3)] + eta01$"Sum Sq"[4])

[1] 0.02587721 0.05509100 0.03105888

2-2　交互作用項(interaction)

　複数の要因（独立変数）が複合的に従属変数に影響を及ぼす事を，計量分析では「交互作用interaction」と呼ぶが，質的比較分析(QCA, Qualitative Comparative Analysis)でいう「結合因果conjunctural causation」と同種の考えであると言ってよいだろう。
以下，本文のスクリプトとその結果を，より装飾を施した形で紹介してゆく。

# 性別×従業上の地位によって8つの群に分けた時の世帯年収平均
# 本文では小数点以下が表示されているが，不要なので四捨五入する。
with(data01, round(tapply(fincome, list(SEX, JOB), mean, na.rm = T), 0))

     正規 非正規 自営 無職
男性  749    388  723  105
女性  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) と指定したり，幾つか工夫を加えている。

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), 
                     las    = 1,
                     col    = c("blue", "red"))
)

平均は先に計算した通りだが，交互作用プロットの様なグラフではそれぞれのグループのケース数（グループサイズ）が全く分からないので，グループ毎の有効ケース数も計算しておくのが良い。

# 各グループ平均
with(data01, round(tapply(fincome, list(SEX, JOB), mean, na.rm=T), 0))

     正規 非正規 自営 無職
男性  749    388  723  105
女性  932    653  538  735

# 各グループの有効ケース数
with(data01, tapply(fincome, list(SEX, JOB), function(x) sum(!is.na(x))))

     正規 非正規 自営 無職
男性  101      6   14    5
女性   44     68   14   44

発展1　2群の母平均の差のt検定と1要因分散分析

　性別による幸福度の平均の差の分析結果を再掲しておこう。

data01$happiness <- c(0:10)[data01$q1700 + 1]

with(data01, 
     by(happiness, sex, mean, na.rm = T))

sex: 1
[1] 6.36646
------------------------------------------------------------ 
sex: 2
[1] 6.917051

with(data01, 
     by(happiness, sex, var,  na.rm = T))

sex: 1
[1] 4.083618
------------------------------------------------------------ 
sex: 2
[1] 3.08568

with(data01, 
     by(happiness, sex, function(x) sum(!is.na(x))))

sex: 1
[1] 161
------------------------------------------------------------ 
sex: 2
[1] 217

var.test(happiness ~ sex, data01)

    F test to compare two variances

data:  happiness by sex
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 

T01 <- t.test(happiness ~ sex, var.equal = T, data01)
T01

    Two Sample t-test

data:  happiness by sex
t = -2.8252, df = 376, p-value = 0.004977
alternative hypothesis: true difference in means between group 1 and group 2 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型にすべきだが，二値変数の場合には結果に本質的な違いはない。

F02 <- oneway.test(happiness ~ sex, var.equal = T, data01)
F01 <- anova(aov(happiness ~ sex, data01)) # anova(aov( )) と一致する
L01 <- anova(lm(happiness ~ sex, data01))  # anova(lm(  ))と一致する
F02

    One-way analysis of means

data:  happiness and sex
F = 7.9819, num df = 1, denom df = 376, p-value = 0.004977

F01

L01

　t検定の結果とF検定の結果を並べて表示させよう。更に，オブジェクトにどんな情報が格納されているのかも一覧を確認しよう。

T01; names(T01)

    Two Sample t-test

data:  happiness by sex
t = -2.8252, df = 376, p-value = 0.004977
alternative hypothesis: true difference in means between group 1 and group 2 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:  happiness and sex
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分布の下側確率と上側確率の和
pt(   abs(T01$statistic), df = T01$parameter, lower.tail = F) +    pt(-1*abs(T01$statistic), df = T01$parameter)

          t 
0.004976846 

# F分布の上側確率
pf(F01$"F value"[1], df1 = F01$Df[1], df2 = F01$Df[2], lower.tail = F)

[1] 0.004976846

　等分散を仮定出来ずにウェルチの検定を行った男女での世帯年収差はどうであろうか。

T02 <- t.test(fincome ~ sex,      data01)
F03 <- oneway.test(fincome ~ sex, data01)
T02

    Welch Two Sample t-test

data:  fincome by sex
t = -0.756, df = 294.6, p-value = 0.4503
alternative hypothesis: true difference in means between group 1 and group 2 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 

F03

    One-way analysis of means (not assuming equal variances)

data:  fincome and sex
F = 0.57153, num df = 1.0, denom df = 294.6, p-value = 0.4503

T02$statistic^2

        t 
0.5715288 

F03$statistic

        F 
0.5715288 

これもやはり，自由度，有意確率が一致し，t統計量の2乗とF統計量が一致している。anova(aov( ))関数では等分散性が仮定出来ない分散分析を行う事は出来ない。

第7章の【練習問題】の解答

1) の答え

問題文で要求されていないコマンドも含んだスクリプトを掲載する。

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

2) の答え

まず，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"))

3) の答え

交互作用を含むモデルを取り上げ，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