ケィオスの時系列解析メモランダム

時系列解析、生体情報解析などをやわらかく語ります

【Rで統計】分散分析ANOVAの考え方:データのばらつきの要因の分解

実験や観測データの解析では「複数の群(グループ)の平均値に差があるか」 を知りたいことが頻繁にあります。たとえば、

  • 3 種類の薬 A、B、C の効果の違い
  • 複数条件下での生理指標の差
  • 学年・群・処理条件ごとの測定値の比較

などです。
 そのような多群比較では、すべての条件どうしを直接比べようとすると、統計的に困った問題が生じます。たとえば、A、B、C の3 群の平均値を比較する場合、「A と B」「B と C」「A と C」という すべての組み合わせについて t 検定を繰り返すことが考えられます。しかし、このような方法を取ると、本当は差がないにもかかわらず、偶然のばらつきを「有意な差」と判断してしまう確率(第 I 種の誤り)が高くなります。
 これは、帰無仮説を棄却する p 値の有意水準を 0.05 に設定していたとしても、実際には差がまったく存在しない群どうしの比較を独立に 100 回繰り返せば、期待値としては約 5 回程度、誤って帰無仮説を棄却してしまうことを意味しています。つまり、各検定では誤りの確率が 5% に抑えられていても、検定の回数を重ねることで、見かけ上の有意差が現れる可能性が累積的に高まるのです。
 そのような理由から、複数の群を比較する際に、個々の組み合わせごとに検定を繰り返すことは望ましくないと考えられています。そこで、まず全体として群の間に差があるかどうかを 一度の検定で評価する方法として、ANOVA(Analysis of Variance:分散分析)が用いられます。

1. ANOVA の本質的な考え方

 ANOVA(分散分析)は名前に「分散」とありますが、研究者が本当に知りたいのは、出発点から一貫して 「群によって平均が変わるか」です。にもかかわらず「分散」という言葉と枠組みが中心になるのは、統計学が 観測誤差を扱う学問として成立してきた歴史と、平均差を「誤差と比較できる形」に落とし込む必要があるためです。
 歴史的には、観測誤差を二乗して足し合わせる 最小二乗法が早くから標準化されました。ルジャンドル(Adrien-Marie Legendre)は 1805 年に最小二乗法(method of least squares)を公表し、ガウス(Carl Friedrich Gauss)は 1809 年の著書 Theoria motus corporum coelestium において、その理論的基礎を体系的に論じました。この「二乗して足す」発想は、平均との差を扱いやすい量(平方和)に変換するもので、後の分散分析にも直結します。
 その後、R.A. Fisher がこの考え方を 実験計画と結びつけ、分散を成分に分けて検定する体系を整えました。フィッシャー(Ronald Aylmer Fisher)は 1918 年に「分散(variance)」という概念を明確に導入し(集団遺伝学の文脈)、1921 年には作物データの解析で分散分析の考え方を実データに適用し、1923 年に続編を発表しています。さらに 1925 年の著書で分散分析が広く知られるようになり、1935 年の『The Design of Experiments』で実験計画と検定の思想が決定的に普及しました。
 この流れを踏まえると、ANOVA が立てている問いは「分散」そのものではなく、次のように言い換えるのが正確です。

観測された全体のばらつきのうち、 群平均の違い(=私たちが見たい“平均差”)で説明できる部分はどれくらいか。 それは、個体差や測定誤差といった 偶然のばらつきに比べて十分に大きいか。

要するに ANOVA は、平均差を見たいという目的を、誤差(偶然のばらつき)と比較可能な形に翻訳するために、「分散(平方和の分解)」という言葉と手続きを採用している、ということです。

■ Rでやってみる

 RでANOVAを実行するには以下のコマンドを実行するだけです。

# 一元配置 ANOVA
fit_aov <- aov(y ~ group, data = dat)

# 結果の表示
summary(fit_aov)

サンプルデータは、この記事の最後「【付録】Rスクリプト」で与えた関数 generate_and_plot_3group_data() を使い、

dat <- generate_and_plot_3group_data(
  n_vec     = c(A = 18, B = 18, C = 18),
  mu_vec    = c(A = 0.0, B = 0.5, C = -1.0),
  var_total = 1.0
)

とすれば、生成できます。データフレーム dat に含まれる群(group)によって平均値 y が異なるかを、一元配置 ANOVA によって検定できます。aov() が検定モデルを作成し、summary() が ANOVA 表を表示します。
 今回得られた結果は次の通りです。

             Df Sum Sq Mean Sq F value   Pr(>F)    
group        2  19.33   9.666   22.53 9.74e-08 ***
Residuals   51  21.88   0.429                     

ここでは、group の効果が大きいこと(p 値 Pr(>F) がきわめて小さいこと)が分かります。つまり、このデータでは「3つの群の平均はすべて同じ」という仮説は支持されず、少なくともどこかの群で平均が異なると判断されます。
 本記事ではRの使い方を詳しく述べるのではなく、ANOVAの考え方に焦点を当てて説明していきます。

【重要】ANOVA が仮定している前提条件

 ただし、この比較が統計的に正当化されるためには、いくつかの前提となる仮定が置かれています。代表的なものは次の3点です。

  1. 独立性(independence) 各観測値は互いに独立である。

  2. 正規性(normality) 各群における誤差(群平均からのずれ)が、正規分布に従うと仮定します。実務的には、サンプルサイズが十分大きい場合、多少の非正規性は許容できます。

  3. 等分散性(homoscedasticity) すべての群で分散が等しいと仮定します。ANOVA では、群内分散をまとめて「誤差分散」として扱うため、この仮定が特に重要になります。

 これらの仮定は、ANOVA が「分散の比」を用いて平均差を評価するために求められる前提です。実データでは厳密に満たされないことも多く、その場合には分散の不等性に配慮した手法(Welch の ANOVA など)や、非パラメトリック手法を検討する必要があります。そのような方法についても、今後の記事で説明していきます。

2 2種類の「ばらつき」

 フィッシャーがANOVAで行った最も重要な発想は、観測されたデータ全体のばらつきを、性質の異なる二つの成分に分けて考えることでした。ANOVA では、データのばらつきを次の二種類に分解します。

3群のデータ比較の例。水平破線:全体平均、水平実線:グループ平均。

(1) 群間分散(between-group variance)

 これは、群ごとの平均値の違いによって生じるばらつきです。これは上図の各群の水平なカラー実線のばらつきで、全体の水平な破線(全体平均)からのずれです。もし条件 A、B、C がそれぞれ異なる効果を持つなら、各群の平均値は互いにずれ、そのずれがデータ全体のばらつきの一部を作ります。

(2) 群内分散(within-group variance)

 一方で、同じ条件の下に集められたデータであっても、測定値は完全に一致することはありません。個体差、測定誤差、環境の微妙な違いなどにより、同じ群の中にも必ずばらつきが生じます。上図の例では、各群(各色)における水平実線からの観測データのばらつきです。
 このような、条件とは無関係に生じるばらつきが群内分散です。これは「ノイズ」や「誤差」に相当する成分であり、平均値の違いを判断するときの基準になります。

 これら 2 種類のばらつきは、ANOVA の結果表(Rでの分析結果)の中では次のように数値化されています。

             Df Sum Sq Mean Sq F value   Pr(>F)    
group        2  19.33   9.666   22.53 9.74e-08 ***
Residuals   51  21.88   0.429                     

ここで、group の行が 群間分散(between-group variance) に対応します。Sum Sq は「群平均どうしの違いによって生じた平方和」、Df はその自由度(群数 − 1)であり、それらを割った Mean Sq群間分散です。この例では、

  • 群間平方和:19.33
  • 自由度:2
  • 群間分散(Mean Sq):9.666

となっています。
 一方、Residuals の行が 群内分散(within-group variance) に対応します。これは、同じ群の中で生じた個体差や測定誤差によるばらつきを集めたもので、Sum Sq が群内平方和、Df が群内の自由度(全サンプル数 − 群数)です。その比である Mean Sq群内分散に相当します。この例では、

  • 群内平方和:21.88
  • 自由度:51
  • 群内分散(Mean Sq):0.429

です。
 ANOVAでは、群間分散と群内分散を使い F 値

\displaystyle{
F
= \frac{\text{群間分散}}{\text{群内分散}}
}

を計算します。次節では、このF 値について説明していきます。

3. ANOVA における F 値の正体

 すでに導入したように、ANOVA の検定統計量は F 値で、形式的には次の比として定義されます。

\displaystyle{
F
= \frac{\text{群間分散}}{\text{群内分散}}
=
\frac{\text{群間平均平方(MS}_{\text{between}}\text{)}}
{\text{群内平均平方(MS}_{\text{within}}\text{)}}
}

ここで重要なのは、分散とは単なるばらつきではなく、平方和を自由度で割った量だという点です。

■ データと記号の定義

  • 群の数: K
  •  k に含まれるデータ数:[texn _ k]
  • 観測値: y _ {ki}(群 ki 番目)
  • 群平均:
\displaystyle{
  \bar{y}_k = \frac{1}{n_k}\sum_{i=1}^{n_k} y_{ki}
}
  • 全体平均:
\displaystyle{
  \bar{y} = \frac{1}{N}\sum_{k=1}^{K}\sum_{i=1}^{n_k} y_{ki},
  \quad
  N=\sum_{k=1}^{K} n_k
}

■ 平方和の分解(ANOVA の核心)

 ANOVA の出発点は、全体のばらつき(平方和)を二つに分解できるという事実です。  

 全平方和(Total Sum of Squares)を

\displaystyle{
SS_{\text{total}}
=
\sum_{k=1}^{K}\sum_{i=1}^{n_k}
\left(y_{ki}-\bar{y}\right)^2

}

と定義します。これは「すべてのデータが全体平均からどれだけずれているか」を表します。

 また、群間平方和(Between-group SS)を

\displaystyle{
SS_{\text{between}}
=
\sum_{k=1}^{K}
n_k\left(\bar{y}_k-\bar{y}\right)^2

}

と定義します。これは群平均が全体平均からどれだけずれているかを測った量で、平均差に対応する成分です。

 さらに、群内平方和(Within-group SS)を

\displaystyle{
SS_{\text{within}}
=
\sum_{k=1}^{K}\sum_{i=1}^{n_k}
\left(y_{ki}-\bar{y}_k\right)^2
}

と定義します。これは、同じ群の中でどれくらいばらつくのが普通かを表し、誤差や個体差に対応します。

 以上の分解の関係は、

\displaystyle{
SS_{\text{total}}
=
SS_{\text{between}}
+
SS_{\text{within}}
}

となります。 この「足し算としての分解」が、ANOVA を可能にしている数学的基盤です。

ANOVAにおける、観測されたばらつきの要因の分解。要因は「群間のベースラインの変動」と「群内の変動」に分けられる。

■ 分散(平均平方)への変換

平方和そのものは、データ数に依存して大きくなるため、自由度で割って平均平方(Mean Square)に変換します。

  • 群間平均平方:
\displaystyle{
  MS_{\text{between}}
  =
  \frac{SS_{\text{between}}}{K-1}
}
  • 群内平均平方:
\displaystyle{
  MS_{\text{within}}
  =  \frac{SS_{\text{within}}}{N-K}
}

■ F 値の定義と意味

 以上を用いて、F 値は

\displaystyle{
F=\frac{MS_{\text{between}}}{MS_{\text{within}}}
}

と定義されます。

 F 値の解釈は、

  •  F \approx 1 群平均のずれは、群内のばらつきと同程度 → 平均差は誤差の範囲内

  •  F \gg 1 群平均のずれが、誤差に比べて明らかに大きい → 条件(群)の効果が疑われる

帰無仮説

\displaystyle{
H_0:\ \mu_1=\mu_2=\cdots=\mu_K
}

が正しければ、理論的に

\displaystyle{
\mathbb{E}[F] = 1
}

となることが、F 分布の性質から保証されています。

 F 値が評価しているのは、あくまで

\displaystyle{
SS_{\text{between}}
=
\sum_{k=1}^{K}
n_k(\bar{y}_k-\bar{y})^2
}

という 合計量です。この中には、

  • どの群が高いのか
  • どの群どうしが違うのか

という情報は個別には含まれていません。そのため、ANOVA が教えてくれるのは

「少なくともどこかの群平均は異なるらしい」

という事実だけです。

4. ANOVA 後の事後検定(概要と次回予告)

 ANOVA で帰無仮説(「すべての群平均は等しい」)が棄却された場合、次に知りたくなるのは「どの群とどの群の間に差があるのか」という点です。この問いに答えるために用いられるのが 事後検定(post-hoc test) です。
 ただし、ここで注意しなければならないのは、

  • 単純に t 検定を繰り返してはいけない
  • 比較回数が増えることで 多重比較問題 が生じる

という点です。事後検定は、この問題を考慮した手法として設計されています。

 代表的な事後検定法には以下のようなものがあります。

Tukey の HSD 検定

  • すべての群ペアを同時に比較
  • 多重比較に対してバランスがよく、標準的に用いられる方法

Bonferroni 法

  • 有意水準を比較回数で割って調整
  • 偽陽性を強く抑える一方で、保守的になりやすい

Holm 法

  • Bonferroni 法を改良した逐次的な方法
  • 実務でよく用いられる現実的な選択肢

Dunnett 検定

  • 対照群(control)と各群のみを比較
  • 薬効試験など、明確な基準群がある場合に有効

 ここでは、それぞれの手法の考え方と位置づけのみを整理しました。具体的な数式の意味、検定統計量の分布、R による実装方法、使い分けの判断基準については、次回以降の記事で順を追って詳しく説明していきます。

5. まとめ

 次回以降の記事では、今回扱った理論を踏まえつつ、数値実験を交えて、ANOVA を「計算手順」ではなく、「考え方」として理解していきたいと思います。

【付録】Rスクリプト

 以下の関数 generate_and_plot_3group_data() は、3群(A, B, C)のサンプルデータを 指定した例数・平均・全体分散にもとづいて乱数生成し、その結果を 群内ばらつきと群間差が直感的に分かる図として描画します。
 入力は、各群の例数ベクトル n_vec、平均ベクトル mu_vec、目標とする全体分散 var_total です。まず、指定した群平均から生じる「群間のばらつき(群平均のズレ)」を計算し、全体分散が var_total になるように、残りを「群内ばらつき(同じ群の中の散らばり)」として割り当てて標準偏差を決めます。その標準偏差を用いて、各群の観測値 y を正規乱数で生成し、dat = data.frame(group, y) を作ります。
 図では、点が各観測値、太い横線が群平均、破線が全体平均を表します。各点から群平均への薄い縦線が「群内のズレ」を可視化し、群平均から全体平均へ向かう矢印が「群間のズレ」を示します。横方向の位置はランダムに揺らさず、群内で等間隔に配置して重なりを避けています。最後に、生成したデータ dat を(表示せずに)返すので、この dat を別スクリプトで ANOVA などの解析にそのまま利用できます。

generate_and_plot_3group_data <- function(
  n_vec,
  mu_vec,
  var_total,
  seed  = 1,
  width = 0.35,
  main  = "Sample Data from Three Groups"
) {

  # -------------------------
  # 0) 前処理・チェック
  # -------------------------
  if (!is.null(seed)) set.seed(seed)

  grp_levels <- names(n_vec)
  if (is.null(grp_levels)) {
    stop("n_vec には group 名(names)を付けてください。")
  }

  n_vec  <- n_vec[grp_levels]
  mu_vec <- mu_vec[grp_levels]

  if (length(n_vec) != length(mu_vec)) {
    stop("n_vec と mu_vec の長さが一致していません。")
  }

  N <- sum(n_vec)

  # -------------------------
  # 1) 群内分散の決定
  #    var_total = var_between + var_within
  # -------------------------
  mu_bar <- sum(n_vec * mu_vec) / N
  var_between <- sum(n_vec * (mu_vec - mu_bar)^2) / N

  var_within <- var_total - var_between
  if (var_within <= 0) {
    stop(sprintf(
      "var_between = %.4f が var_total = %.4f を上回っています。\nmu_vec の差を小さくするか var_total を大きくしてください。",
      var_between, var_total
    ))
  }
  sd_within <- sqrt(var_within)

  # -------------------------
  # 2) データ生成
  # -------------------------
  group <- factor(rep(grp_levels, times = n_vec), levels = grp_levels)

  y <- numeric(N)
  idx_start <- 1
  for (g in grp_levels) {
    idx_end <- idx_start + n_vec[g] - 1
    y[idx_start:idx_end] <- rnorm(n_vec[g], mean = mu_vec[g], sd = sd_within)
    idx_start <- idx_end + 1
  }

  dat <- data.frame(group = group, y = y)

  # -------------------------
  # 3) 可視化用の量
  # -------------------------
  mean_all <- mean(dat$y)
  mean_by_group <- tapply(dat$y, dat$group, mean)

  # x 座標(等間隔配置)
  x_jit <- numeric(nrow(dat))
  for (g in levels(dat$group)) {
    idx <- which(dat$group == g)
    n_g <- length(idx)
    offs <- seq(-width, width, length.out = n_g)
    x_center <- which(levels(dat$group) == g)
    x_jit[idx] <- x_center + offs
  }

  # -------------------------
  # 4) 描画
  # -------------------------
  ylim <- range(dat$y, mean_all, mean_by_group) + c(-0.8, 0.8)

  par(mar = c(5, 5, 3, 1), las = 1)
  plot(
    NA, NA,
    xlim = c(0.5, length(grp_levels) + 0.5),
    ylim = ylim,
    xlab = "Group", ylab = "Value",
    xaxt = "n",
    main = main
  )
  axis(1, at = seq_along(grp_levels), labels = grp_levels)

  # 全体平均
  abline(h = mean_all, lty = 2, lwd = 2)

  # 群平均
  cols_mean <- c("red", "navy", "forestgreen")
  for (g in levels(dat$group)) {
    xg <- which(levels(dat$group) == g)
    segments(
      x0 = xg - 0.45, y0 = mean_by_group[g],
      x1 = xg + 0.42, y1 = mean_by_group[g],
      lwd = 3, col = cols_mean[xg]
    )
  }

  # 群内ばらつき(縦線)
  cols_within_by_group <- c(
    A = rgb(1,   0.5, 0.5, alpha = 0.4),
    B = rgb(0,   0.4, 1,   alpha = 0.4),
    C = rgb(0.3, 1,   0.4, alpha = 0.4)
  )
  cols_within <- cols_within_by_group[as.character(dat$group)]

  for (i in seq_len(nrow(dat))) {
    g <- as.character(dat$group[i])
    segments(
      x0 = x_jit[i], y0 = dat$y[i],
      x1 = x_jit[i], y1 = mean_by_group[g],
      col = cols_within[i], lwd = 1
    )
  }

  # 点
  cols_points <- c(A = 2, B = 4, C = 3)
  points(
    x_jit, dat$y,
    pch = 16, cex = 1.1,
    col = cols_points[as.character(dat$group)]
  )

  # 群平均 → 全体平均
  for (g in levels(dat$group)) {
    xg <- which(levels(dat$group) == g)
    if (abs(mean_by_group[g] - mean_all) > 0.15) {
      arrows(
        x0 = xg - 0.4, y0 = mean_by_group[g],
        x1 = xg - 0.4, y1 = mean_all,
        lwd = 2, length = 0.12, angle = 20, code = 1
      )
    }
  }

  # -------------------------
  # 5) データを返す
  # -------------------------
  invisible(dat)
}

 以下のRスクリプトは、先に定義した generate_and_plot_3group_data() 関数を実際に呼び出し、3グループのサンプルデータを生成・可視化する例です。
 ここでは、各グループ A・B・C の例数をすべて 18 に設定し、群平均をそれぞれ A = 0.0、B = 0.5、C = −1.0 と与えています。また、var_total = 1.0 によって、全観測値をまとめたときの分散が概ね 1 になるようにデータが生成されます(標本誤差のため厳密には一致しません)。
 この1行を実行すると、指定した条件に基づいて乱数データが生成され、同時に図が描画されます。図には、各観測点、群平均、全体平均が重ねて表示され、同じ群の中のばらつき(群内ばらつき)と、群平均どうしの違い(群間差)を視覚的に確認できます。
 戻り値として得られる dat は、group(因子)と y(数値)からなるデータフレームです。この dat をそのまま用いて、後続のスクリプトで ANOVA や回帰分析などを行うことができます。図による直感的理解と、数理的な解析とを切り分けて扱える点が、この使い方のポイントです。

dat <- generate_and_plot_3group_data(
  n_vec     = c(A = 18, B = 18, C = 18),
  mu_vec    = c(A = 0.0, B = 0.5, C = -1.0),
  var_total = 1.0
)

※ もし記事の中で「ここ違うよ」という点や気になるところがあれば、気軽に指摘していただけると助かります。質問や「このテーマも取り上げてほしい」といったリクエストも大歓迎です。必ず対応するとは約束できませんが、できるだけ今後の記事で扱いたいと思います。
 それと、下の はてなブログランキング はあまり信用できる指標ではなさそうですが(笑)、押してもらえるとシンプルに励みになります。気が向いたときにポチッとしていただけたら嬉しいです。