橋本・村井（2015）の一要因分散分析＋多重比較のシミュレーション結果をRで再現してみる

Packages

pacman::p_load(
  tidyverse,
  PMCMRplus,
  Rfast,
  future,
  furrr,
  progressr
)

Background

見つけた論文（橋本・村井，2015）の再現です。論文では「分散分析で有意なら多重比較する」という手法において、分散分析と多重比較で判断が異なる場合がどれくらいの確率で生じるかについてシミュレーションしていました。Java言語を使用したと書いてありましたが、私はRしか使えないのでRで再現してみようと思います。シミュレーションの条件は以下の通りです。

• 分析方法：一要因分散分析＋多重比較
• 多重比較法：3手法
  • TukeyのHSD法
  • ペアワイズt検定のBonferroni補正：不偏分散はpoolしたもの。そのままだと長いので、この記事ではBonferroni法と略すことにします。
  • Scheffe法：対比較のみ
• 水準数a：3
  • 3, 4, 5
• 母平均\(\mu_{1}\)のパターン：2
  • 0：すべての水準の母平均が等しい。
  • 0.2：第1群だけ母平均が異なる。
    • 第1群は正規分布\(N(\mu_1, 1)\)から発生。他の群は正規分布\(N(0, 1)\)から。
• 各群のサンプルサイズ：668
  • 対応のないt検定でd = 0.2の差を検出したいときに、有意水準\(\alpha = 0.5\%\)（水準数が最大の5のとき、対比較の組み合わせ数は\({}_5 \mathrm{C}_2=10\)なので、\(5\%\div10=0.5\%\)）で検出力80%になる数。

というわけで、シミュレーションの条件は水準数3パターン×母平均2パターンの全部で6条件です。それぞれの条件で100万回反復して、各検定の有意性の判断の割合を求めていきます。

Preparation

下準備として、使用する関数を決めます。本来はなるべく標準で提供されている既存関数を用いるのが再現性的にもよさそうなのですが、条件ごとに100万回の反復というのはそこそこ時間がかかるので、なるべく時間がかからない関数を選びたいです。もっとも、速さにこだわるのであれば、既存関数からシミュレーションに不要な処理をそぎ落とした関数を自分で書く方がいいのですが……ちょっとめんどくさいので今回は既存関数を使っていきます。（いずれ自力でやってみたいです。）

一番時間がかかりそうな水準数5のベクトルを用意して、ベンチマークを回します。

set.seed(2026)
bench_vec_value5 <- rnorm(668L * 5)
bench_vec_group5 <- factor(1:5) |>
  rep(each = 668L)

one-way ANOVA

一元配置分散分析についてはいくつか種類があります。今回は下記5つを回してみたいと思います。

• oneway.test(): 標準で使える関数。デフォルトではWelchのANOVA（var.equal = FALSE）なので、古典的な一要因分散分析をしたい場合はvar.equal = TRUEにします。
• aov() |> summary(): こちらも標準で使える関数たち。
• aov() |> anova(): これも標準で使える関数たち。分散分析表をdata.frameで出力してくれるので、summary()と比べてp値などの結果が取り出しやすいです。
• Rfast::anova1: 入力がベクトルか行列のときに高速な計算ができるRfastパッケージの関数。
• matrixTests::row_oneway_equalvar: 今回ちょっと探してみたら見つかったパッケージ。row-wiseまたはcol-wiseに計算ができるらしい。列（col）バージョンもあったけど、中身を見たら、入力を転置（t()）してこれに渡していたので、ここでは一旦行方向の関数を採用。

bench::mark(
  "oneway.test" = {
    oneway.test(bench_vec_value5 ~ bench_vec_group5, var.equal = TRUE) |>
      _$p.value |>
      `<`(e1 = _, e2 = .05)
  },
  "aov_summary" = {
    aov(bench_vec_value5 ~ bench_vec_group5) |>
1      summary.aov() |>
      _[[1]] |>
      _[1, "Pr(>F)"] |>
      `<`(e1 = _, e2 = .05)
  },
  "aov_anova" = {
    aov(bench_vec_value5 ~ bench_vec_group5) |>
      anova() |>
      _[1, "Pr(>F)"] |>
      `<`(e1 = _, e2 = .05)
  },
  "matrixTests::row_oneway" = {
     matrixTests::row_oneway_equalvar(
      x = bench_vec_value5,
      g = bench_vec_group5
    ) |>
      _$pvalue |>
      `<`(e1 = _, e2 = .05)
  },
  "Rfast::anova1" = {
    Rfast::anova1(
      x = bench_vec_value5,
      ina = bench_vec_group5
    ) |>
      _["p-value"] |>
      `<`(e1 = _, e2 = .05)
  },
  iterations = 100L,
2  check = \(x, y) all.equal(x, y, check.attributes = FALSE)
)

1

generic関数はクラス名を絞ると、他のクラスを探されないからちょっと速くなる？みたいなのをネットのどこかで見た覚えがあります。

2

Rfast::anova1だけnamed_vectorを返してくるので、check.attributes = FALSEにしてあります。

# A tibble: 5 × 6
  expression                   min   median `itr/sec` mem_alloc `gc/sec`
  <bch:expr>              <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
1 oneway.test              697.5µs  802.6µs     1231.   870.8KB     12.4
2 aov_summary               1.01ms   1.07ms      908.     1.5MB     18.5
3 aov_anova                 1.38ms   1.45ms      670.     941KB     20.7
4 matrixTests::row_oneway  703.4µs  731.2µs     1355.   651.5KB     13.7
5 Rfast::anova1             60.3µs  65.85µs    14660.     102KB      0  

組み込みのoneway.test()もマイクロ秒単位なので十分速いのですが、やはりRfast::anova1が最速ですね。10倍くらい速いです。aov_anovaがaov_summaryよりも若干遅いのには驚きました。（中身を見ると何かわかるんだろうか…？）

なお、matrixTestsもRfastも行列を入力して行（列）ごとに分析を実施する高速な関数を提供している（matrixTest::row_oneway_equalvar()はまさにそう）ので、行列を入力にした場合も試してみます。水準数5で、あまり列数を大きくするとapplyが遅くなるのは目に見えているので、一旦100回分のデータが詰まった行列を出してみます。通常、matrix()は列方向にデータが収納されるので、各列に乱数データが収まっていることになります。

set.seed(2026)
bench_mat_value5 <- matrix(
  rnorm(668L * 5L * 100),
  ncol = 100
)

# check
dim(bench_mat_value5)

[1] 3340  100

入力が行列の場合で列（行）ごとに計算したい場合は、基本はapply系関数を使うと思います。今回はapply()とsapply()を試してみます。中身で使うのはベクトル入力で最速だったRfast::anova1()にします。

bench::mark(
  "apply_anova1" = {
    bench_mat_value5 |>
      apply(
        MARGIN = 2,
        FUN = \(x) {
          Rfast::anova1(x, ina = bench_vec_group5) |>
            _["p-value"] |>
            `<`(e1 = _, e2 = .05)
        }
      )
  },
  "sapply_anova1" = {
    sapply(
      1:ncol(bench_mat_value5),
      FUN = \(x) {
        Rfast::anova1(bench_mat_value5[, x], ina = bench_vec_group5) |>
            _["p-value"] |>
            `<`(e1 = _, e2 = .05)
      }
    )
  },
  "matrixTests::col_oneway_equalvar" = {
    bench_mat_value5 |>
      matrixTests::col_oneway_equalvar(g = bench_vec_group5) |>
      _[, "pvalue"] |>
      `<`(e1 = _, e2 = .05)
  },
  "Rfast::anovas" = {
    bench_mat_value5 |>
      Rfast::anovas(ina = bench_vec_group5) |>
       _[, "p-value"] |>
      `<`(e1 = _, e2 = .05)
  },
  check = \(x, y) all.equal(x, y, check.attributes = FALSE),
  iterations = 100L
)

# A tibble: 4 × 6
  expression                            min   median `itr/sec` mem_alloc `gc/sec`
  <bch:expr>                       <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
1 apply_anova1                       7.12ms   7.41ms      132.   12.78MB    162. 
2 sapply_anova1                      6.34ms   6.66ms      149.   10.26MB     58.0
3 matrixTests::col_oneway_equalvar   4.43ms   4.57ms      216.    7.97MB    111. 
4 Rfast::anovas                      1.22ms   1.26ms      775.    2.71MB     67.4

apply系よりも行列処理系の関数が速いですが、やはりRfast::anovas()が最速です。というわけで、今回は速さを優先してRfastの関数で一要因分散分析をすることにしました。入力をベクトルにするか行列にするかは、いったん後に回します。

Multiple Comparison Procedures

Tukey-HSD

標準でTukeyHSD()が使えるのですが、aov()で作ったオブジェクトしか（？）引数に入りません。そして、aov()はformula形式でモデルを入れるのですが、これが微妙に時間がかかる感じがします。普段の一発の分析なら別に体感できないので標準で使える方を使う方がいいと思うのですが、シミュレーションで何万回も反復しようとするとちょっと気になります。というわけで、ベクトルでそのまま入れられるPMCMRplus::tukeyTest()を試してみます。なおPMCMRplus::tukeyTest()はformula形式でも入れられるので、それも比較してみます。計測はいずれかの対比較ペアのp値が.05を下回るかどうかの判定までです。

bench::mark(
  "TukeyHSD_aov" = {
    aov(bench_vec_value5 ~ bench_vec_group5) |>
      TukeyHSD() |>
      _$bench_vec_group5 |>
      _[, "p adj"] |>
      (\(p) any(p < .05, na.rm = TRUE)) ()
  },
  "PMCMRplus::tukeyTest_default" = {
    PMCMRplus::tukeyTest(
      x = bench_vec_value5,
      g = bench_vec_group5
    ) |>
      _$p.value |>
      (\(p) any(p < .05, na.rm = TRUE)) ()
  },
  "PMCMRplus::tukeyTest_formula" = {
    PMCMRplus::tukeyTest(bench_vec_value5 ~ bench_vec_group5) |>
      _$p.value |>
      (\(p) any(p < .05, na.rm = TRUE)) ()
  },
  iterations = 100L
)

# A tibble: 3 × 6
  expression                        min   median `itr/sec` mem_alloc `gc/sec`
  <bch:expr>                   <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
1 TukeyHSD_aov                   6.06ms   6.14ms     161.     3.47MB    15.9 
2 PMCMRplus::tukeyTest_default   2.74ms   2.83ms     354.   502.85KB     7.22
3 PMCMRplus::tukeyTest_formula  15.52ms  15.77ms      63.1  844.12KB     1.29

ベクトルで入力したPMCMRplus::tukeyTest()が一番早いですね。なお、なぜかformula形式で入れるとTukeyHSDよりも遅くなりますが。というわけで、速さの観点からPMCMRplus::tukeyTest()のベクトル入力を採用しました。

Pairwise t-test

標準で使えるpariwise.t.test()の他にペアワイズでプールした不偏分散を使ってくれるStudentのt検定をしてくれる関数をrstatix::pairwise_t_test以外に知らないです¹。copilotにも聞いて探してみましたが、代替の関数が見つからず、やるとしたら自分で組んでみるかどうかみたいな感じでした。Bioconductorの方は探してないので、もしかしたらあるかもしれない…？ご存知の方がいたら教えてください。今回は、そもそもベクトルを入力できるし、引数でp値調整法も選べることであとからp.adjust()する手間が省けるしというわけでpariwise.t.test(..., p.adjust.method = "bonferroni", pool.sd = TRUE)を使いました。

Scheffe’s test

Scheffe法の対比較については、PMCMRplusとDescToolsに収録されているのを知っていたので、その2パッケージの関数を比較してみました。ベクトル入力とformula入力の両方を試しています。

bench::mark(
  "PMCMRplus_default" = {
    PMCMRplus::scheffeTest(
      x = bench_vec_value5,
      g = bench_vec_group5
    ) |>
      _$p.value |>
      (\(p) any(p < .05, na.rm = TRUE)) ()
  },
  "PMCMRplus_formula" = {
    PMCMRplus::scheffeTest(bench_vec_value5 ~ bench_vec_group5) |>
      _$p.value |>
      (\(p) any(p < .05, na.rm = TRUE)) ()
  },
  "DescTools_default" = {
    DescTools::ScheffeTest(
      x = bench_vec_value5,
      g = bench_vec_group5,
    ) |>
      _$g |>
      _[, "pval"] |>
      (\(p) any(p < .05, na.rm = TRUE)) ()
  },
  "DescTools_default_onlyp" = {
    DescTools::ScheffeTest(
      x = bench_vec_value5,
      g = bench_vec_group5,
      conf.level = NA
    ) |>
      _$g |>
      (\(p) any(p < .05, na.rm = TRUE)) ()
  },
  "DescTools_aov" = {
    DescTools::ScheffeTest(aov(bench_vec_value5 ~ bench_vec_group5)) |>
      _$bench_vec_group |>
      _[, "pval"] |>
      (\(p) any(p < .05, na.rm = TRUE)) ()
  },
  iterations = 100L
)

# A tibble: 5 × 6
  expression                   min   median `itr/sec` mem_alloc `gc/sec`
  <bch:expr>              <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
1 PMCMRplus_default        772.4µs  859.3µs    1119.   409.77KB    22.8 
2 PMCMRplus_formula        13.59ms  13.89ms      71.2  681.85KB     2.97
3 DescTools_default         4.57ms   5.06ms     196.    71.13MB    34.6 
4 DescTools_default_onlyp   4.52ms   4.69ms     212.     3.04MB    21.0 
5 DescTools_aov             4.59ms   4.72ms     210.     3.04MB    18.2 

DescTooks::scheffeTest()は信頼区間を計算しないようにできる（conf.level = NA）ので、もしかしたら速いかも…と思いましたが、一番早いのはPMCMRplus::scheffeTest()のベクトル入力でした。なのでこれを採用しました。

input

使う関数は決まったのですが、乱数データの入力をベクトルにするのか行列にするのかを決めていませんでした。どちらがいいか、これもまたベンチマークを回して決めたいと思います。

vector

使う関数群を一つの関数にまとめて関数にするのですが、この時に乱数データをベクトルで作ってしまいます。

f_analyse_vec <- function(n_levels, mu_g1, n_per_group, vec_group, alpha = .05) {
  # data generation
1  vec_value <- c(
    rnorm(n_per_group, mean = mu_g1, sd = 1),
    rnorm(n_per_group * (n_levels - 1), mean = 0, sd = 1)
  )
  # analysis
  p_oneway <- Rfast::anova1(
    x = vec_value,
    ina = vec_group
  ) |>
    _["p-value"]
  p_tukey <- PMCMRplus::tukeyTest(
    x = vec_value,
    g = vec_group
  ) |>
    _$p.value
  p_bonferroni <- pairwise.t.test(
    x = vec_value,
    g = vec_group,
    p.adjust.method = "bonf",
    pool.sd = TRUE
  ) |>
    _$p.value
  p_scheffe <- PMCMRplus::scheffeTest(
    x = vec_value,
    g = vec_group
  ) |>
    _$p.value
  # result
  c(
    oneway = p_oneway < alpha,
    tukey = any(p_tukey < alpha, na.rm = TRUE),
    bonferroni = any(p_bonferroni < alpha, na.rm = TRUE),
    scheffe = any(p_scheffe < alpha, na.rm = TRUE)
  )
}

1

ここ。

チェックしてみましょう。

f_analyse_vec(
  n_levels = 5,
  mu_g1 = 0,
  n_per_group = 668L,
  vec_group = bench_vec_group5,
  alpha = .05
)

oneway.p-value          tukey     bonferroni        scheffe 
         FALSE          FALSE          FALSE          FALSE 

OKですね。それで、この関数をreplicate()で反復して²、集計までする関数を作ります。

f_analyse_batch <- function(n_sim, n_levels, mu_g1, n_per_group, vec_group, alpha = .05) {
  replicate(
    n = n_sim,
    expr = f_analyse_vec(
      n_levels = n_levels,
      mu_g1 = mu_g1,
      n_per_group = n_per_group,
      vec_group = vec_group,
      alpha = alpha
    )
  ) |>
1    t() |>
    data.frame() |>
2    rename(oneway = 1) |>
3    pivot_longer(
      cols = -oneway,
      names_to = "method",
      values_to = "anysig"
    ) |>
    count(oneway, method, anysig)
}

1

f_analyse_vecは長さ4のベクトルを戻すので、replicate()の戻り値は、4×n行の行列になります。あとあとでdplyr::count()を使いたいので、転置してdf化します。

2

Rfast::anova1()の戻り値は名前付きベクトルで、p値はp-valueという名前で返ってきています。それが引き継がれてoneway.p-valueという列名になっているので（1個前の結果参照）、短くします。

3

tidyr::pivot_longer()の引数values_trans_formに無名関数を入れて値をfactor化して、dplyr::count(..., .drop = FALSE)で0件のケースを集計しようか迷いました。でも最後にtidyr::complete()を1回挟んで拾ってあげればいいということに気づいたので、そうします。

試しに100回で確認。

f_analyse_batch(
  n_sim = 100,
  n_levels = 5,
  mu_g1 = 0,
  n_per_group = 668L,
  vec_group = bench_vec_group5,
  alpha = .05
)

# A tibble: 9 × 4
  oneway method     anysig     n
  <lgl>  <chr>      <lgl>  <int>
1 FALSE  bonferroni FALSE     92
2 FALSE  bonferroni TRUE       2
3 FALSE  scheffe    FALSE     94
4 FALSE  tukey      FALSE     92
5 FALSE  tukey      TRUE       2
6 TRUE   bonferroni TRUE       6
7 TRUE   scheffe    FALSE      4
8 TRUE   scheffe    TRUE       2
9 TRUE   tukey      TRUE       6

100回分の結果の集計ができました。

matrix

まず、多重比較の関数群だけ1つにまとめます。

f_analyse_mcp <- function(x, g) {
  p_tukey <- PMCMRplus::tukeyTest(
    x = x,
    g = g
  ) |>
    _$p.value
  p_bonf <- pairwise.t.test(
    x = x,
    g = g,
    p.adjust.method = "bonf",
    pool.sd = TRUE
  ) |>
    _$p.value
  p_scheffe <- PMCMRplus::scheffeTest(
    x = x,
    g = g
  ) |>
    _$p.value
  c(
    tukey = any(p_tukey < .05, na.rm = TRUE),
    bonferroni = any(p_bonf < .05, na.rm = TRUE),
    scheffe = any(p_scheffe < .05, na.rm = TRUE) 
  )
}

チェックします。

f_analyse_mcp(x = bench_vec_value5, g = bench_vec_group5)

     tukey bonferroni    scheffe 
      TRUE       TRUE      FALSE 

OKですね。これを関数の中に入れます。関数の中でデータ行列を作り、分散分析はRfast::anovas()で行列のまま処理、多重比較はsapply()（apply()でも可）でデータ行列の行で処理する感じです。

f_analyse_mat <- function(n_sim, n_levels, mu_g1, n_per_group, vec_group, alpha = .05) {
  # data generation
  mat_value <- replicate(
    n = n_sim,
    expr = c(
      rnorm(n_per_group, mean = mu_g1, sd = 1),
      rnorm(n_per_group * (n_levels - 1), mean = 0, sd = 1)
    )
  )
  # analysis
  oneway <- Rfast::anovas(
    mat_value,
    ina = vec_group
  ) |>
    _[, "p-value"] |>
    `<`(e1 = _, e2 = .05)
  mcp <- sapply(
    X = 1:n_sim,
    FUN = \(x) f_analyse_mcp(mat_value[, x], g = vec_group)
  )
  # result
  cbind(oneway, t(mcp)) |>
    data.frame() |>
    pivot_longer(
      cols = -oneway,
      names_to = "method",
      values_to = "anysig"
    ) |>
    count(oneway, method, anysig)
}

チェックします。

f_analyse_mat(
  n_sim = 100,
  n_levels = 5,
  mu_g1 = 0,
  n_per_group = 668L,
  vec_group = bench_vec_group5,
  alpha = .05
)

# A tibble: 10 × 4
   oneway method     anysig     n
   <lgl>  <chr>      <lgl>  <int>
 1 FALSE  bonferroni FALSE     88
 2 FALSE  scheffe    FALSE     88
 3 FALSE  tukey      FALSE     87
 4 FALSE  tukey      TRUE       1
 5 TRUE   bonferroni FALSE      6
 6 TRUE   bonferroni TRUE       6
 7 TRUE   scheffe    FALSE     10
 8 TRUE   scheffe    TRUE       2
 9 TRUE   tukey      FALSE      4
10 TRUE   tukey      TRUE       8

OKですね。100回分の結果が出ています。

Comparison

どっちの方が高速なのでしょうか。正直なところ、使っている関数はほぼ同じなので変わらないような気がしますが、行列を入力する方が多重比較の関数をsapplyしているところが気になります（実はapply系はfor loop に比べてそんなに速くない説ですね。replicate()の中身はsapply()だろ！という声は聞かなかったことにします）。100回で試してみます。時間がかかりそうなのでiterations = 10程度にしておきます。

bench::mark(
  "input_vector" = f_analyse_batch(
    n_sim = 100,
    n_levels = 5,
    mu_g1 = 0,
    n_per_group = 668L,
    vec_group = bench_vec_group5
  ),
  "input_matrix" = f_analyse_mat(
    n_sim = 100,
    n_levels = 5,
    mu_g1 = 0,
    n_per_group = 668L,
    vec_group = bench_vec_group5
  ),
  check = FALSE,
  iterations = 10L
)

Warning: Some expressions had a GC in every iteration; so filtering is disabled.

# A tibble: 2 × 6
  expression        min   median `itr/sec` mem_alloc `gc/sec`
  <bch:expr>   <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
1 input_vector    460ms    476ms      1.89     114MB     4.16
2 input_matrix    449ms    461ms      1.98     119MB     4.56

あまり変わらないですね…。microbenchmarkでも試しに計測してみましょう。

microbenchmark::microbenchmark(
  "input_vector" = f_analyse_batch(
    n_sim = 100,
    n_levels = 5,
    mu_g1 = 0,
    n_per_group = 668L,
    vec_group = bench_vec_group5
  ),
  "input_matrix" = f_analyse_mat(
    n_sim = 100,
    n_levels = 5,
    mu_g1 = 0,
    n_per_group = 668L,
    vec_group = bench_vec_group5
  ),
  times = 10L,
  check = "equivalent",
  setup = set.seed(1)
)

Unit: milliseconds
         expr      min       lq     mean   median       uq      max neval cld
 input_vector 460.9722 462.2968 495.6739 464.7887 496.7971 607.5561    10   a
 input_matrix 448.9897 452.1344 482.6125 455.4639 460.6882 600.4662    10   a

ベンチを回すたびに値は若干上下しますが、まあ傾向は変わらないですね。結局のところ、Rfastの計算が速すぎること、相対的に遅い多重比較の関数の呼び出し回数は結局同じなので、どちらでやっても変わらないってところでしょうか。おそらくもっと速くしたいなら、多重比較の関数を軽量化した自作関数にしないとダメですね。ベンチの結果を見ると行列を入力する方が若干RAMを食いがちなので、今回はベクトルを入力する方を使いたいと思います。

Simulation

シミュレーションの事前準備はできたので、条件を設定して実際に回してみましょう。Packagesにもあるように、今回はfutureとfurrrパッケージを使って、かつ時間がかかるのでプログレスバーを出しながらやりたいと思います。

Set up

シミュレーションの条件は以下のdfの通りです。

df_design <- expand_grid(
  n_levels = 3:5,
  mu_g1 = c(0, .2)
)

df_design

# A tibble: 6 × 2
  n_levels mu_g1
     <int> <dbl>
1        3   0  
2        3   0.2
3        4   0  
4        4   0.2
5        5   0  
6        5   0.2

定数を設定していきます。有意水準は5%、1群あたりのサンプルサイズは668です。各群のインデックスベクトルはlist化して1つにまとめておきます。

num_alpha <- .05
int_n_per_group <- 668L

list_vec_group <- unique(df_design$n_levels) |>
  lapply(
    FUN = \(x) {
      seq_len(x) |>
        rep(each = int_n_per_group) |>
        as.factor()
    }
  ) |>
  setNames(nm = 3:5)

# check
list_vec_group

$`3`
  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [48] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [95] 1 1 1 1 1 1
 [ reached 'max' / getOption("max.print") -- omitted 1904 entries ]
Levels: 1 2 3

$`4`
  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [48] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [95] 1 1 1 1 1 1
 [ reached 'max' / getOption("max.print") -- omitted 2572 entries ]
Levels: 1 2 3 4

$`5`
  [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [48] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
 [95] 1 1 1 1 1 1
 [ reached 'max' / getOption("max.print") -- omitted 3240 entries ]
Levels: 1 2 3 4 5

繰り返し回数は100万回です。
次にブロックサイズを考えます。最初自分でコードを書いたときは、作った関数をそのまま100万回回せばいいやーと思って試しに1000回くらいで回したのですが、プログレスバーのフリッカーがひどくて進捗状況が逆に分からなくなる感じでした。Copilotに聞いたら、1回の処理でプログレスバー1回更新だと処理速度が速すぎて表示がおかしくなるかもということと、furrrで並列処理するなら100万を各workerに分散させるよりも、ある程度の処理数をまとめてブロックにして分散させる数を多くしない方がいいかもという点を教えてもらったので、今回はその方針をとってみた感じです。先ほどベンチマークを回したときの100回だと約500ミリ秒くらいで結果が返ってきているのですが、それだと小さいような気もします。1000回くらいだと処理速度とメモリ使用量的にも悪くなさそうな気がするので、ブロックサイズは1000にしました。

int_n_iter <- 1e6L
int_block_size <- 1000L

# check
c(int_n_iter = int_n_iter, int_block_size = int_block_size) |>
  format(trim = TRUE, scientific = FALSE, big.mark = ",")

    int_n_iter int_block_size 
   "1,000,000"        "1,000" 

並列処理用のインデックスdfを用意します。

df_jobs <- df_design |>
  crossing(
    block_id = seq_len(int_n_iter / int_block_size)
  ) |>
  mutate(
    n_sim = int_block_size
  )

df_jobs

# A tibble: 6,000 × 4
   n_levels mu_g1 block_id n_sim
      <int> <dbl>    <int> <int>
 1        3     0        1  1000
 2        3     0        2  1000
 3        3     0        3  1000
 4        3     0        4  1000
 5        3     0        5  1000
 6        3     0        6  1000
 7        3     0        7  1000
 8        3     0        8  1000
 9        3     0        9  1000
10        3     0       10  1000
# ℹ 5,990 more rows

プログレスバーを準備をします。furrr::future_*()を利用し並列処理をプログレスバーありで使うときは、引数.progress = TRUEにすれば簡易的なものが出てきます。ただし、Helpによるとこの引数は廃止予定らしく、よりロバストに使えるprogressrパッケージを使うことが推奨されています。progressrは様々なパッケージが提供しているプログレスバー表示をサポートしています。その中でも今回はcliパッケージの表示を使うことにします。cliはコンソールに色付きでいろいろ表示してくれるパッケージで、purrrのmap系関数で引数.progress = TRUEとすると表示されるプログレスバーはこれを利用しています。

1progressr::handlers(
2  progressr::handler_cli(
      format = paste(
        "{cli::pb_bar}",
        "{cli::pb_percent}",
        "[{cli::pb_elapsed}]",
        "{cli::pb_eta_str}"
      ),
      clear = FALSE
    ),
  global = TRUE
)

1

progressr::handlers()にprogressr::handler_*()を設定して、プログレスバーを表示させるようにします。

2

cliパッケージのプログレスバー表示を使用できるようにします。引数formatに文字列で入力することで、好きな表示にすることができます。今回は上から順にプログレスバー、完了した割合（％）、経過時間、完了までの時間です。clear = FALSEで終わった後に消えないようにしています。

シードを固定して、並列化します。future::multisesson()は特に設定しないと利用可能なCPUスレッド数をすべて使ってくれますが、今回は全コア数－1にします。このPCは8C16Tなので15になるはずですね。

# seed
set.seed(2026)

# multi-session
1n_workers <- parallelly::availableCores() - 1
2future::plan(future::multisession(workers = n_workers))

# check
3future::nbrOfWorkers()

1

parallellyパッケージはfutureパッケージをインストールすると一緒に入るはず。

2

futureパッケージを読み込んでいるなら、plan(mutltisession, workers = ...)で通ります。

3

workerが増えたかの確認は、future::plan()かfuture::nbrOfWorkers()で立ち上がっているworkerの数を見ることで確認できます。ちなみにマルチセッション化した後のfuture::plan()の戻り値は長いです。

[1] 15

Run Simulation

シミュレーションを回します。処理用のインデックスdfをfurrr::future_pmap_dfr()を使って処理します。purrrやfurrrのpmap系は引数.xにdfを持ってくると、dfの各行の要素を引数.fで作られた関数の引数に入れていきます。例えばdf_jobsの一行目は、n_levels = 3, mu_g1 = 0, block_id = 1, n_sim = 500なので、それが.fの無名関数の引数にそのまま入るという仕組みです。また、今回は処理中にプログレスバーを出したいので、progressr::with_progress()の第1引数exprの中に処理内容を書いていきます。こうすることで、並列化していてもうまくプログレスバーが表示されます。

progressr::with_progress({
1  pb <- progressr::progressor(steps = nrow(df_jobs))
  
  temp_res <- future_pmap_dfr(
    .l = df_jobs,
    .f = \(n_levels, mu_g1, block_id, n_sim) {
      # group vector
      vec_group <- list_vec_group[[as.character(n_levels)]]
      # processing
      res <- f_analyse_batch(
        n_sim = n_sim,
        n_levels = n_levels,
        mu_g1 = mu_g1,
        n_per_group = int_n_per_group,
        vec_group = vec_group,
        alpha = num_alpha
      ) |>
        mutate(
          n_levels = n_levels,
          mu_g1 = mu_g1,
          block_id = block_id,
          .before = 1
        )
      # progress bar update
      pb()
      # return
      res
    },
    .options = furrr_options(
2      seed = TRUE
    )
  )
})

1

プログレスバーをオブジェクトとして登録します。引数stepsはプログレスバーの100%となる数を入れます。今回はdf_jobsの各行が並列処理で割り振られるので、行数を入れます。

2

並列化しているときに乱数を固定したい場合は、必ずfurrr::future_*()の引数.optionsでfurrr_options(seed = ...)を設定します。今回のように事前にシードを設定してればTRUEで、もしくはシード値を1つまたは.xの要素数だけ入れることもできます。

処理が終わったら、忘れないうちに必ずfuture::plan(future::sequential())で並列化を終わらせます。

1future::plan(future::sequential()); future::nbrOfWorkers()

1

futureパッケージを読み込んでいるなら、plan(sequential)で通ります。

[1] 1

かかった時間は44分27秒でした。シングルコアで回すよりは速いはずですが、思ったほど速くならなかった印象です。

結果はdfで返ってきているのですが、各ブロックごとの集計がまとまっているだけなので、全体で合算して最終結果にします。

res_df_table <- 
  temp_res |>
  reframe(
    n = sum(n),
    .by = c(n_levels, mu_g1, oneway, method, anysig)
  ) |>
1  complete(
    n_levels, mu_g1, oneway, method, anysig,
    fill = list(n = 0)
  ) |>
  mutate(
    across(
      .cols = c(oneway, anysig),
      .fns = \(x) factor(x, levels = c(FALSE, TRUE), labels = c("ns", "sig"))
    ),
    method = factor(method, levels = c("tukey", "bonferroni", "scheffe")),
  ) |>
  mutate(
    pct = n / int_n_iter,
    mcse = sqrt(pct * (1 - pct) / int_n_iter),
    .by = c(n_levels, mu_g1)
  ) |>
  pivot_wider(
    names_from = n_levels,
    names_glue = "lvs{n_levels}_{.value}",
    names_vary = "slowest",
    values_from = c(n, pct, mcse),
    values_fill = 0
  ) |>
  arrange(mu_g1, oneway, method, anysig)

1

分散分析が非有意でScheffe法で有意になるケースは観測されないはずなので、dplyr::count()を使うとケースがその組み合わせがそもそも存在しないことになります。tidyr::complete()で選択した変数であり得る組み合わせを列挙してもらって、観測されなかったケースについては引数fillに設定した値を入れてもらいます。欠損を埋める値は、named list形式で入れる必要があります。

Result

Quartoで作っている都合上、qmdからhtmlへレンダリングするときにRのチャンクはすべて再計算されます。ということは、この投稿を作るのに再度45分以上かかることが確定です。さすがにそれはイヤなので、あらかじめ回しておいたデータを読み込んでおくことにします。シード値を固定してあるので、おそらく問題はないでしょう。（この記事のフォルダにrdsファイルもついてくるはずなので、需要はないと思いますが欲しい人がいたらgithubの/posts/の中のこの記事のフォルダを探してください。）

res_df_table <- readr::read_rds("anova_mcp_sig_res_df.rds")

one-way ANOVA

まず初めに、分散分析の有意性判断を見てみましょう。

res_df_table |>
  filter(oneway == "sig") |>
  reframe(
    across(
      .cols = ends_with("_n"),
1      .fns = \(x) sum(x) / 3
    ),
    .by = mu_g1
  ) |>
  mutate(
    across(
      .cols = -mu_g1,
      .fns = list(
        pr = \(x) x / int_n_iter,
        pr_mcse = \(x) {
          pr <- x / int_n_iter
          sqrt(pr * (1 - pr) / int_n_iter)
        }
      )
    )
  ) |>
2  relocate(mu_g1, sort(tidyselect::peek_vars())) |>
  gt::gt() |>
  gt::fmt_percent(
    columns = contains("pr"),
    decimals = 4
  ) |>
  gt::tab_spanner_delim(
    delim = "_",
    columns = starts_with("lvs"),
    split = "first",
    limit = 1
  ) |>
  gt::cols_label(
    mu_g1 ~ gt::md("$\\mu_{1}$"),
    ends_with("n_pr") ~ "%",
    ends_with("_mcse") ~ "mcse(%)"
  ) |>
  gt::text_transform(
    fn = \(x) str_replace_all(x, pattern = "lvs", replacement = "a = "),
    locations = gt::cells_column_spanners()
  )

1

各手法ごとに分散分析が有意判定だった数は同じなので、合計して÷3すれば分散分析が有意になった数が出てきますね。

2

tidyselect::peek_vars()はselectの文脈で現状で利用可能な変数名を返してくれるらしい。一個上の処理で_nの列、_pctの列、…の順番になっていたので、lvs...でソートするのに使ってみました。
（参考：https://stackoverflow.com/questions/29873293/dplyr-order-columns-alphabetically-in-r ）

\(\mu_{1}\)	a = 3			a = 4			a = 5
	n	%	mcse(%)	n	%	mcse(%)	n	%	mcse(%)
0.0	50009	5.0009%	0.0218%	50443	5.0443%	0.0219%	49950	4.9950%	0.0218%
0.2	973365	97.3365%	0.0161%	975313	97.5313%	0.0155%	974066	97.4066%	0.0159%

母平均に差がない場合はどの水準数でも5%程度になりました。元論文もおおむね5%程度だったのでOKですね。母平均に差があるときは検定力ということになりますが、どの水準数でも97%程度でした。

MCP

次に元論文とほぼ同じ表を作ってみます。gtパッケージで表を作りますが、共通している部分は関数化しときます。

f_make_gt <- \(x) {
  x |>
    gt::gt() |>
    gt::fmt_percent(
    columns = matches("(pct$|mcse$)"),
    decimals = 4
    ) |>
    gt::tab_spanner_delim(
      delim = "_",
      columns = starts_with("lvs")
    ) |>
    gt::cols_label(
      method = "Procedures"
    ) |>
    gt::text_transform(
      fn = \(x) {
        str_replace_all(
          x,
          pattern = c(
            "pct" = "%",
            "mcse" = "mcse(%)"
          )
        )
      },
      locations = gt::cells_column_labels(
        columns = matches("(pct$|mcse$)")
      )
    ) |>
    gt::text_transform(
      fn = \(x) str_replace(x, pattern = "lvs", replacement = "a = "),
      locations = gt::cells_column_spanners()
    ) |>
    gt::tab_style(
      style = gt::cell_text(align = "left"),
      locations = gt::cells_body(columns = method)
    ) |>
    gt::text_transform(
      fn = \(x) str_to_title(x),
      locations = gt::cells_body(columns = method)
    ) |>
    gt::cols_hide(
      columns = c(mu_g1, oneway, anysig)
    )
}

Table 1

まずは元論文のTable 1、すべての水準の母平均が等しいときに、一要因分散分析では有意だけど多重比較では有意差ありな対比較ペアがない場合です。

res_df_table |>
  filter(mu_g1 == 0, oneway == "sig", anysig == "ns") |>
  f_make_gt() |>
  gt::tab_header(
    title = gt::md("$\\mu_{1}$: **0** / ANOVA: **sig.** / MCP: **n.s.**")
  )

\(\mu_{1}\): 0 / ANOVA: sig. / MCP: n.s.
Procedures	a = 3			a = 4			a = 5
	n	%	mcse(%)	n	%	mcse(%)	n	%	mcse(%)
Tukey	5265	0.5265%	0.0072%	8191	0.8191%	0.0090%	10863	1.0863%	0.0104%
Bonferroni	8026	0.8026%	0.0089%	12260	1.2260%	0.0110%	15334	1.5334%	0.0123%
Scheffe	11988	1.1988%	0.0109%	23370	2.3370%	0.0151%	32278	3.2278%	0.0177%

水準数3のBonferroni法が0.02%ほど元論文より高いくらいで、他は0.1%の桁までは同じになりました。再現成功とっていいでしょう。

この表は分散分析で有意、すなわちType I error を犯している状態で、多重比較をすることによって正しい判断ができている割合になります。水準数が増えるほど割合が上がっているので、元論文でも指摘されているように、水準数が多いほどType I errorの割合を抑えられているかもしれないといえるでしょう。水準数をもう少し増やして検証してみるのもいいかもしれません。

参考までに、すべての水準の母平均が等しいときの各多重比較法のType I errorの割合は以下の通りです。

res_df_table |>
  filter(mu_g1 == 0, anysig == "sig") |>
  reframe(
    across(
      .cols = ends_with("_n"),
      .fns = \(x) sum(x)
    ),
    .by = method
  ) |>
  mutate(
    across(
      .cols = -method,
      .fns = list(
        pct = \(x) {
          (x / int_n_iter) |>
            num(digits = 4, label = "%", scale = 100)
        }
      )
    )
  ) |>
  relocate(method, sort(tidyselect::peek_vars()))

# A tibble: 3 × 7
  method     lvs3_n lvs3_n_pct lvs4_n lvs4_n_pct lvs5_n lvs5_n_pct
  <fct>       <int>          %  <int>          %  <int>          %
1 tukey       50084     5.0084  50448     5.0448  50178     5.0178
2 bonferroni  43821     4.3821  41872     4.1872  40094     4.0094
3 scheffe     38021     3.8021  27073     2.7073  17672     1.7672

Tukey法が5%程度を維持、Bonferroni法はそれより少し低く、Scheffe法が一番低い結果になりました。つまり、多重比較を単体で行ってもType I error（多重比較ならfamilywise errorと言う方が正確でしょうか）は抑えられているわけです。

Table 2

次に、元論文のTable 2、第1群の母平均だけ異なるときに、一要因分散分析では有意だけど多重比較では有意な対比較ペアがなしという判断になる場合です。

res_df_table |>
  filter(mu_g1 == 0.2, oneway == "sig", anysig == "ns") |>
  f_make_gt() |>
  gt::tab_header(
    title = gt::md("$\\mu_{1}$: **0.2** / ANOVA: **sig.** / MCP: **n.s.**")
  )

\(\mu_{1}\): 0.2 / ANOVA: sig. / MCP: n.s.
Procedures	a = 3			a = 4			a = 5
	n	%	mcse(%)	n	%	mcse(%)	n	%	mcse(%)
Tukey	4940	0.4940%	0.0070%	4902	0.4902%	0.0070%	4091	0.4091%	0.0064%
Bonferroni	7786	0.7786%	0.0088%	8999	0.8999%	0.0094%	8441	0.8441%	0.0091%
Scheffe	11983	1.1983%	0.0109%	23345	2.3345%	0.0151%	38428	3.8428%	0.0192%

水準数3のScheffe法が若干低めに出てますが（と言っても四捨五入すれば同じです）、これも0.1%の桁までほぼ同じですので再現成功と言っていいでしょう。

中身に関していえば、実践上、これが一番困るパターンです。分散分析の結果を見て「じゃあどこに差があるんだろう」と多重比較の結果を見たら、どの対比較ペアでも有意差ありと判断できず、どう結果を書いたらいいんだろう…となってしまうパターンです。割合が低い方がありがたいのですが、学生のデータを見ていると年に1回は必ず見かける事象です。しかし、今回の条件設定（正規分布、等分散）だとそんなに多く発生する訳ではなさそうですね。
水準数との関係性はどうなのでしょう？Scheffe法は水準数が増えるにつれて割合が高くなっているように見えますが、Tukey法はその逆のような気もします。一方でBonferroni法は水準数が増えるにつれ検出力が落ちていくので、この割合が高くなっていく傾向がみられてもいいのですが、そういうわけでもなさそう（？）。巷では「水準数が5以下ならBonferroni、それ以上ならHolmを使え！」というどこ出典なのそれ基準があったりしますが、これについては水準数をもっと増やしてみないとわかりませんね。

Table 3

次に、元論文のTable 3、第1群の母平均だけ異なるときに、一要因分散分析では有意ではないけど多重比較だと有意差ありの対比較ペアがある場合です。

res_df_table |>
  filter(mu_g1 == 0.2, oneway == "ns", anysig == "sig") |>
  f_make_gt() |>
  gt::tab_header(
    title = gt::md("$\\mu_{1}$: **0.2** / ANOVA: **n.s.** / MCP: **sig.**")
  )

\(\mu_{1}\): 0.2 / ANOVA: n.s. / MCP: sig.
Procedures	a = 3			a = 4			a = 5
	n	%	mcse(%)	n	%	mcse(%)	n	%	mcse(%)
Tukey	1704	0.1704%	0.0041%	2569	0.2569%	0.0051%	3991	0.3991%	0.0063%
Bonferroni	627	0.0627%	0.0025%	1173	0.1173%	0.0034%	2030	0.2030%	0.0045%
Scheffe	0	0.0000%	0.0000%	0	0.0000%	0.0000%	0	0.0000%	0.0000%

0.1%の桁までは同じになりました。これも再現できたといっていいでしょう。

Scheffe法は分散分析の結果と整合するため0件になるのはいいとして、元論文でも指摘されているように、Tukey法とBonferroni法はそれ単体で実施していれば検出できていた差を、先に分散分析を行うことによって無視してしまったことになります。つまり、もともと各群を対比較するつもりがあっても、先に分散分析を行うという慣習に則り実行することで、見つけられたであろう差を逃しているわけです。割合としてはかなり低い事象ですが、Table 2の条件と同様にこれはこれで困る事象で少ないに越したことはありません。一方で、水準数が増えるほどその割合も増えているように見えます。

参考までに、多重比較法単体での検定力を出しておきます。

res_df_table |>
  filter(mu_g1 == 0.2, anysig == "sig") |>
  reframe(
    across(
      .cols = ends_with("_n"),
      .fns = \(x) sum(x)
    ),
    .by = method
  ) |>
  mutate(
    across(
      .cols = -method,
      .fns = list(
        pct = \(x) {
          (x / int_n_iter) |>
            num(digits = 4, label = "%", scale = 100)
        }
      )
    )
  ) |>
  relocate(method, sort(tidyselect::peek_vars()))

# A tibble: 3 × 7
  method     lvs3_n lvs3_n_pct lvs4_n lvs4_n_pct lvs5_n lvs5_n_pct
  <fct>       <int>          %  <int>          %  <int>          %
1 tukey      970129    97.0129 972980    97.2980 973966    97.3966
2 bonferroni 966206    96.6206 967487    96.7487 967655    96.7655
3 scheffe    961382    96.1382 951968    95.1968 935638    93.5638

Table 1のところでも見たように、Type I errorの生じにくさと検定力は裏表の関係です。Tukey法が一番検定力が高く、次いでBonferroni法、Scheffe法の順になりました。

All results

最後に今回の結果をすべて載せておきます。まずはコンソール表示版。割合の列名はpctですが％の計算はしてないです（列名はprとかにしておけばよかったかも。癖であとで表を作るときにどうせgtでパーセント表記にしちゃうからいいだろって思っていつもpctにしてしまうんですよね）。

（個人的メモ）：チャンク単体でR.optionsを変えるにはこういう書き方をするっぽい⇩。

```{r}
#| attr-output: 'style="max-height: 350px;"'
1#| R.options: !expr list(width = 140)
res_df_table |>
  print(n = 99)
```

1

ここ。微妙に横に長いので横の表示を広くしたい。

# A tibble: 24 × 13
   mu_g1 oneway method     anysig lvs3_n lvs3_pct lvs3_mcse lvs4_n lvs4_pct lvs4_mcse lvs5_n lvs5_pct lvs5_mcse
   <dbl> <fct>  <fct>      <fct>   <int>    <dbl>     <dbl>  <int>    <dbl>     <dbl>  <int>    <dbl>     <dbl>
 1   0   ns     tukey      ns     944651 0.945    0.000229  941361  0.941   0.000235  938959  0.939   0.000239 
 2   0   ns     tukey      sig      5340 0.00534  0.0000729   8196  0.00820 0.0000902  11091  0.0111  0.000105 
 3   0   ns     bonferroni ns     948153 0.948    0.000222  945868  0.946   0.000226  944572  0.945   0.000229 
 4   0   ns     bonferroni sig      1838 0.00184  0.0000428   3689  0.00369 0.0000606   5478  0.00548 0.0000738
 5   0   ns     scheffe    ns     949991 0.950    0.000218  949557  0.950   0.000219  950050  0.950   0.000218 
 6   0   ns     scheffe    sig         0 0        0              0  0       0              0  0       0        
 7   0   sig    tukey      ns       5265 0.00526  0.0000724   8191  0.00819 0.0000901  10863  0.0109  0.000104 
 8   0   sig    tukey      sig     44744 0.0447   0.000207   42252  0.0423  0.000201   39087  0.0391  0.000194 
 9   0   sig    bonferroni ns       8026 0.00803  0.0000892  12260  0.0123  0.000110   15334  0.0153  0.000123 
10   0   sig    bonferroni sig     41983 0.0420   0.000201   38183  0.0382  0.000192   34616  0.0346  0.000183 
11   0   sig    scheffe    ns      11988 0.0120   0.000109   23370  0.0234  0.000151   32278  0.0323  0.000177 
12   0   sig    scheffe    sig     38021 0.0380   0.000191   27073  0.0271  0.000162   17672  0.0177  0.000132 
13   0.2 ns     tukey      ns      24931 0.0249   0.000156   22118  0.0221  0.000147   21943  0.0219  0.000146 
14   0.2 ns     tukey      sig      1704 0.00170  0.0000412   2569  0.00257 0.0000506   3991  0.00399 0.0000630
15   0.2 ns     bonferroni ns      26008 0.0260   0.000159   23514  0.0235  0.000152   23904  0.0239  0.000153 
16   0.2 ns     bonferroni sig       627 0.000627 0.0000250   1173  0.00117 0.0000342   2030  0.00203 0.0000450
17   0.2 ns     scheffe    ns      26635 0.0266   0.000161   24687  0.0247  0.000155   25934  0.0259  0.000159 
18   0.2 ns     scheffe    sig         0 0        0              0  0       0              0  0       0        
19   0.2 sig    tukey      ns       4940 0.00494  0.0000701   4902  0.00490 0.0000698   4091  0.00409 0.0000638
20   0.2 sig    tukey      sig    968425 0.968    0.000175  970411  0.970   0.000169  969975  0.970   0.000171 
21   0.2 sig    bonferroni ns       7786 0.00779  0.0000879   8999  0.00900 0.0000944   8441  0.00844 0.0000915
22   0.2 sig    bonferroni sig    965579 0.966    0.000182  966314  0.966   0.000180  965625  0.966   0.000182 
23   0.2 sig    scheffe    ns      11983 0.0120   0.000109   23345  0.0233  0.000151   38428  0.0384  0.000192 
24   0.2 sig    scheffe    sig    961382 0.961    0.000193  951968  0.952   0.000214  935638  0.936   0.000245 

gtのinteractive table版。そのままだと横に長くなるので、水準数を列にしてあります。use_filters = TRUEで絞り込みできるようにしたのですが、母平均の列は0とだけ打つと0も0.2もどっちも引っかかるので、0は0.0に変更しました。

res_df_table |>
  pivot_longer(
    cols = starts_with("lvs"),
    names_to = c("n_levels", ".value"),
    names_pattern = "lvs(\\d)_(.+)"
  ) |>
  relocate(n_levels, .after = mu_g1) |>
  arrange(mu_g1, n_levels) |>
  mutate(
    mu_g1 = factor(mu_g1, levels = c(0, 0.2), labels = c("0.0", "0.2"))
  ) |>
  gt::gt() |>
  gt::fmt_number(
    columns = c(pct, mcse),
    decimals = 6
  ) |>
  gt::opt_interactive(
    use_filters = TRUE,
    use_page_size_select = TRUE,
    page_size_default = 12,
    page_size_values = seq(12, 12 * 6, by = 12)
  )

Conclusion

橋本・村井 （2015） による一要因分散分析＋多重比較のシミュレーションをRで再現してみました。使用言語は違いましたが、再現は成功したといえるでしょう。言語が異なっても再現できたというのはいいですね。時間をかけて100万回回した甲斐がありました。乱数発生に関してはRのデフォルトは元論文と同じメルセンヌ・ツイスタ法なので、シード値さえわかれば元論文をピッタリ再現できるのでしょうか…？ちょっと気になります。

一要因分散分析＋多重比較という分析手法についてですが、個人的には多重比較を最初からやればよくねと最近思うようになってきました。心理学だと（心理学に限らず？）、最終的にはどの群と群の間に差があるか知りたい場合の方がほとんどなはずです。となると最初から対比較をすればよくね？となるわけです。２つの検定を行う分、判断の相違みたいな面倒なことだって生じるわけですし、適切な多重比較法を用いればそれ単体でFamilywise errorも抑えられます。

多重比較する前に一要因分散分析をするという慣行については、 永田 （1998, p. 98） はFisherのPLSD法³の名残ではないかと指摘しています。妥当な多重比較法が開発されていなった時代のやり方が今も受け継がれているようです。（なお、 永田 （1998） は、Scheffe法のような手順に分散分析が含まれている手法を使うのでないなら、分散分析は不要どころか原則行わないと述べています。） 橋本・村井 （2015） が主張するようにこの手法は考え直されてもいいのではないかと思います。

今後また時間があるときに、条件や設定を増やして「一要因分散分析＋多重比較」手法のシミュレーションに挑戦してみようと思います。

Session Information

Notesessioninfo

R version 4.5.3 (2026-03-11 ucrt)
Platform: x86_64-w64-mingw32/x64
Running under: Windows 11 x64 (build 26100)

Matrix products: default
  LAPACK version 3.12.1

locale:
[1] LC_COLLATE=Japanese_Japan.utf8  LC_CTYPE=Japanese_Japan.utf8    LC_MONETARY=Japanese_Japan.utf8
[4] LC_NUMERIC=C                    LC_TIME=Japanese_Japan.utf8    

time zone: Asia/Tokyo
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] progressr_1.0.0    furrr_0.4.0        future_1.75.0      Rfast_2.1.5.2      RcppParallel_6.2.0
 [6] zigg_0.0.2         Rcpp_1.1.2         PMCMRplus_1.9.12   lubridate_1.9.5    forcats_1.0.1     
[11] stringr_1.6.0      dplyr_1.2.1        purrr_1.2.2        readr_2.2.0        tidyr_1.3.2       
[16] tibble_3.3.1       ggplot2_4.0.3      tidyverse_2.0.0   

loaded via a namespace (and not attached):
 [1] gld_2.6.8            sandwich_3.1-3       readxl_1.5.0         rlang_1.3.0         
 [5] magrittr_2.0.5       multcomp_1.4-31      otel_0.2.0           matrixStats_1.5.0   
 [9] e1071_1.7-17         compiler_4.5.3       BWStest_0.2.3        systemfonts_1.3.2   
[13] vctrs_0.7.3          kSamples_1.2-12      pkgconfig_2.0.3      fastmap_1.2.0       
[17] labeling_0.4.3       utf8_1.2.6           rmarkdown_2.31       markdown_2.0        
[21] tzdb_0.5.0           haven_2.5.5          ragg_1.5.2           xfun_0.60           
[25] cachem_1.1.0         Rmpfr_1.1-2          litedown_0.10        jsonlite_2.0.0      
[29] SuppDists_1.1-9.9    matrixTests_0.2.3.1  gmp_0.7-5.1          parallel_4.5.3      
[33] DescTools_0.99.60    R6_2.6.1             stringi_1.8.9        RColorBrewer_1.1-3  
[37] parallelly_1.48.0    boot_1.3-32          cellranger_1.1.0     knitr_1.51          
[41] zoo_1.9-0            base64enc_0.1-6      pacman_0.5.1         Matrix_1.7-6        
[45] splines_4.5.3        timechange_0.4.0     tidyselect_1.2.1     rstudioapi_0.19.0   
[49] yaml_2.3.12          codetools_0.2-20     listenv_1.0.0        lattice_0.22-9      
[53] withr_3.0.3          S7_0.2.2             evaluate_1.0.5       survival_3.8-9      
[57] proxy_0.4-29         bench_1.1.4          xml2_1.6.0           pillar_1.11.1       
[61] generics_0.1.4       hms_1.1.4            commonmark_2.0.0     scales_1.4.0        
[65] rootSolve_1.8.2.4    globals_0.19.1       class_7.3-24         glue_1.8.1          
[69] lmom_3.3             tools_4.5.3          data.table_1.18.4    reactable_0.4.5     
[73] Exact_3.3            fs_2.1.0             mvtnorm_1.4-2        grid_4.5.3          
[77] crosstalk_1.2.2      cli_3.6.6            textshaping_1.0.5    profmem_0.7.0       
[81] expm_1.0-0           gt_1.3.0             gtable_0.3.6         reactR_0.6.1        
[85] sass_0.4.10          digest_0.6.39        TH.data_1.1-5        htmlwidgets_1.6.4   
[89] farver_2.1.2         memoise_2.0.1        htmltools_0.5.9      lifecycle_1.0.5     
[93] httr_1.4.8           multcompView_0.1-12  microbenchmark_1.5.0 MASS_7.3-66         

References

橋本 貴充・村井 潤一郎（2015）．分散分析で有意になり多重比較で有意差がない確率　日本心理学会大会発表論文集，79, 1EV–063.

永田 靖（1998）．多重比較法の実際　応用統計学，27(2), 93–108. https://doi.org/10.5023/jappstat.27.93

Footnotes

1. rstatix::pairwise_t_test()はformula形式の入力なのとパッケージがそもそもdfのパイプライン運用を前提している感があるので、今回は採用していないです。普段の分析だと使いやすいパッケージです。Rfast::ttest()はどうやらWelchしかできないみたいなので、今回の目的だと残念ながら使えないです。

2. replicate()を使うかforを使うかは悩むところです。replicate()の中身はsapply()なのですが、入力がデカすぎるときはapply系よりもforの方が速いという話があります（例えばこの記事とか https://qiita.com/Atsushi776/items/c31f2345b9c698354c81 ）。

3. Rでいうとoneway.test(..., var.equal = TRUE)で有意なら、pairwise.t.test(..., p.adjustment.method = "none", pool.sd = TRUE)をするという手法です。