Delacre et al. (2017) の等分散性の検定のシミュレーション (Fig. 1) を再現してみた＋α

R memo
statistics
analyse
simulation
Published

2026-08-26

Abstract

本論に進む前にさらっと等分散性の検定のシミュレーションをしてたのでそれの再現をしてみたのと、シミュレーション条件を拡張してみて回してみた。

Packages

pacman::p_load(
  tidyverse,
  Rfast,
  matrixTests,
  future,
  furrr
)

Background

見つけた論文のシミュレーションを再現してみようシリーズ。今回は Delacre et al. （2017） 1です。論文の内容としては、各種t検定（Student, Welch, Yuen）をシミュレーションで比較して、心理学においてはWelchのt検定を基本的には使うべきだと主張するものです。

その問題提起部分として、2段階t検定、すなわち等分散性の検定を行ってその結果次第でStudentとWelchのどちらかを使うという手法はよくないことの説明がされていました。その根拠の一つとして、等分散性の検定（Levene検定とBrown-Forsythe検定）の検出力についてシミュレーションで示したグラフ（Fig. 1）が載っていたので、それを再現してみようと思います。なお、各種t検定のシミュレーションはRでやったっぽいですが、等分散性の検定のシミュレーションは何でやったかは書いてありません。（グラフのデザインを見る限り、Rのデフォルトのプロットで作った感じがするので、多分Rでやってみたんじゃないでしょうか。）

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

  • 等分散性の検定：2種
    • Levene検定：SPSSの出力でよくみかけるやつ。平均値からの偏差の絶対値を使う。
    • Brown-Forsythe検定：Levene検定とほぼ同じ。Levene検定が平均値からの偏差の絶対値を使うところを中央値からの偏差の絶対値を使う。
  • サンプルサイズn：15パターン。
    • 1群当たり10から80で増分は5。
      • Rでいうとseq(10, 80, 5)
  • 標準偏差の比SDR：3パターン
    • 1.1, 1.5, 2

以上の条件から、論文では\(15\times3=45\)パターンを各100万回反復し、検出力を経験的に求めていました。せっかくなので、ここでは以下のようにもう少し条件を増やして回してみたいと思います。

  • 等分散性の検定：2種⇒3種
    • 【追加】F検定：2群のSDの比をとる。Rでいうとvar.test。2群以上に一般化されたのがBartlett検定（bartrett.test）。これらの手法は、検出力が高いが正規性からの逸脱に脆弱なのが知られているということで、論文では門前払いされています。
  • サンプルサイズn：15⇒40
    • 1群当たり5から100で増分5。：なんかキリのいい数字で。
  • 標準偏差の比SDR：11パターン
    • 1から2まで増分は.1：せっかくなので比が同じである1を含めてみます。
      • 1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2
      • データ発生は第1群は正規分布\(N~(0,1)\)で固定、第2群を正規分布\(N~(0, sd)\)からにしたいと思います。

グラフの形を再現するなら多分10万回あれば十分だとは思いますが、あとあとで計算のパワーに余裕があることが分かったので各条件100万回の反復にしました。断りがない限り有意水準\(\alpha = .05\)とします。

Preparation

まず、使用する関数を決めます。ベンチマーク用にデータを作っておきます。要素数が一番多い1群当たり100人の条件でデータを作ります。

set.seed(2026)
test_vec_val <- rnorm(100 * 2)
test_vec_group <- factor(1:2) |>
  rep(each = 100L)
test_mat <- matrix(
  rnorm(200 * 10),
  ncol = 10
)

F-test

等分散性のためのF検定（2群）は標準でvar.test()が用意されています。これを使えばいいのですが、Rfastにより速いvar2test()が用意されています。速度を比較してみましょう。

bench::mark(
  "var.test" = {
    var.test(
      x = test_vec_val[test_vec_group == "1"],
      y = test_vec_val[test_vec_group == "2"]
    ) |>
      _$p.value
  },
  "Rfast::var2test" = {
    Rfast::var2test(
      x = test_vec_val[test_vec_group == "1"],
      y = test_vec_val[test_vec_group == "2"]
    ) |>
      _[["p-value"]]
  }
)
# A tibble: 2 × 6
  expression           min   median `itr/sec` mem_alloc `gc/sec`
  <bch:expr>      <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
1 var.test          99.8µs  106.7µs     9133.   138.3KB     10.4
2 Rfast::var2test   23.7µs   25.9µs    37129.    31.4KB     11.1

標準の関数と比べてRfast::var2test()が速くて省メモリですね。なお、Rfast::var2tests()という行列データに対応したバージョンもあります。

bench::mark(
  "apply_vartest" = {
    apply(
      test_mat,
      MARGIN = 2,
      FUN = \(x) {
        var.test(
          x = x[test_vec_group == "1"],
          y = x[test_vec_group == "2"]
        ) |>
          _$p.value
      }
    )
  },
  "Rfast::var2tests" = {
    Rfast::var2tests(
      x = test_mat,
      ina = test_vec_group
    ) |>
      _[, "pvalue"]
  }
)
# A tibble: 2 × 6
  expression            min   median `itr/sec` mem_alloc `gc/sec`
  <bch:expr>       <bch:tm> <bch:tm>     <dbl> <bch:byt>    <dbl>
1 apply_vartest      1.02ms   1.06ms      925.   156.9KB     12.8
2 Rfast::var2tests   82.8µs   88.5µs    10936.    72.7KB     10.5

これがエグいほど速いです。後述するLeveneとBF検定の高速版も行列データの入力に対応しているので、今回は行列対応版の方を使います。

Levene Test & Brown-Forsythe test

両検定とも標準では収録されていないです（多分）。自分が入れているパッケージの中では、分散分析で使うcarや何でも入っているDescTools、なんで入れたか覚えていないlawstatなどに収録されていました。引数をmeanにするかmedianにするかで、両検定を使い分けることができます。 その中でもmatrixTestsの関数が行列データの入力に対応していてかつ高速ということで、これを採用したいと思います。2 3

# Levene test
matrixTests::col_levene(
  x = test_mat,
  g = test_vec_group
)
   obs.tot obs.groups df.between df.within   statistic      pvalue
1      200          2          1       198 4.509681583 0.034944874
2      200          2          1       198 8.574689346 0.003808122
3      200          2          1       198 0.023500009 0.878320086
4      200          2          1       198 0.604011691 0.437979007
5      200          2          1       198 0.760892566 0.384105885
6      200          2          1       198 0.025220600 0.873980474
7      200          2          1       198 1.552297246 0.214268830
8      200          2          1       198 0.009257321 0.923447153
9      200          2          1       198 2.019766555 0.156836044
10     200          2          1       198 2.112022002 0.147729170
# BF test
matrixTests::col_brownforsythe(
  x = test_mat,
  g = test_vec_group
)
   obs.tot obs.groups df.between df.within  statistic      pvalue
1      200          2          1       198 4.22455048 0.041155763
2      200          2          1       198 7.53182883 0.006619116
3      200          2          1       198 0.02435853 0.876135127
4      200          2          1       198 0.58566542 0.445011416
5      200          2          1       198 0.73956431 0.390840541
6      200          2          1       198 0.00722156 0.932363342
7      200          2          1       198 1.56804198 0.211968702
8      200          2          1       198 0.00751811 0.930992007
9      200          2          1       198 2.02811747 0.155985832
10     200          2          1       198 2.08520115 0.150313082

Analyse function

というわけで、どの手法に対しても乱数データを行列で入れて高速に計算してもらいたいと思います。

f_analyse <- function(mat, g) {
  p_vartest <- Rfast::var2tests(
    x = mat,
    ina = g
  ) |>
    _[, "pvalue"]
  p_levene <- matrixTests::col_levene(
    x = mat,
    g = g
  ) |>
    _[, "pvalue"]
  p_bf <-  matrixTests::col_brownforsythe(
    x = mat,
    g = g
  ) |>
    _[, "pvalue"]
  cbind(vartest = p_vartest, levene = p_levene, bf = p_bf)
}

# check
f_analyse(test_mat, g = test_vec_group)
          vartest      levene          bf
 [1,] 0.077580928 0.034944874 0.041155763
 [2,] 0.007063534 0.003808122 0.006619116
 [3,] 0.726876959 0.878320086 0.876135127
 [4,] 0.222244350 0.437979007 0.445011416
 [5,] 0.245936901 0.384105885 0.390840541
 [6,] 0.909507055 0.873980474 0.932363342
 [7,] 0.366067277 0.214268830 0.211968702
 [8,] 0.679857459 0.923447153 0.930992007
 [9,] 0.171250194 0.156836044 0.155985832
[10,] 0.300277072 0.147729170 0.150313082

OKですね。バッチ処理用の関数も用意します。

f_analyse_batch <- function(n_sim, n_each, sd_g2, vec_group , sig_level = .05) {
  # data generation
  temp_mat <- replicate(
    n_sim,
    expr = c(
      rnorm(n = n_each, mean = 0, sd = 1),
      rnorm(n = n_each, mean = 0, sd = sd_g2)
    )
  )
  # analyse
  f_analyse(temp_mat, g = vec_group) |>
    `<`(e1 = _, e2 = sig_level) |>
    colSums() |>
1    tibble::enframe(name = "procedures", value = "n_rej")
}
1
名前付きベクトルをn行2列df（1列目名前、2列目値）にするのに便利。ここだけベンチ回したら、なんか(\(x) data.frame(procedures = names(x), n_rej = x)) ()より5倍くらい速かった。

チェック。

system.time(
  f_analyse_batch(
    n_sim = 1000,
    n_each = 100,
    sd_g2 = 1,
    vec_group = test_vec_group
  ) |>
    print()
)
# A tibble: 3 × 2
  procedures n_rej
  <chr>      <dbl>
1 vartest       58
2 levene        50
3 bf            48
   user  system elapsed 
   0.05    0.00    0.04

OKですね。十分速い。

Simulation

下準備はできたので、シミュレーションの条件を設定します。

Set up

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

df_design <- expand_grid(
  n_each = seq(5, 100, 5),
  sd_g2 = seq(1, 2, .1)
)

# check
df_design
# A tibble: 220 × 2
   n_each sd_g2
    <dbl> <dbl>
 1      5   1  
 2      5   1.1
 3      5   1.2
 4      5   1.3
 5      5   1.4
 6      5   1.5
 7      5   1.6
 8      5   1.7
 9      5   1.8
10      5   1.9
# ℹ 210 more rows

群分けベクトルのリストも作っておきます

list_vec_group <- unique(df_design$n_each) |>
  setNames(nm = _) |>
  lapply(
    FUN = \(x) {
      factor(1:2) |>
        rep(each = x)
    }
  )

# check
head(list_vec_group)
$`5`
 [1] 1 1 1 1 1 2 2 2 2 2
Levels: 1 2

$`10`
 [1] 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2
Levels: 1 2

$`15`
 [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
Levels: 1 2

$`20`
 [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
Levels: 1 2

$`25`
 [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 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[49] 2 2
Levels: 1 2

$`30`
 [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 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
[49] 2 2 2 2 2 2 2 2 2 2 2 2
Levels: 1 2

反復数は1条件当たり100万です。ブロックサイズですが、1%（= 10000）にしておきます。4

int_n_iter <- 1e6
int_n_block_size <- int_n_iter * .01

# check
c(n_iter = int_n_iter, block_size = int_n_block_size) |>
  format(scientific = FALSE, big.mark = ",") 
     n_iter  block_size 
"1,000,000" "   10,000"

割り振り用のdfを作ります5。

df_jobs <- df_design |>
  crossing(batch_id = seq_len(int_n_iter / int_n_block_size)) |>
  mutate(n_sim = int_n_block_size)

# check
df_jobs
# A tibble: 22,000 × 4
   n_each sd_g2 batch_id n_sim
    <dbl> <dbl>    <int> <dbl>
 1      5     1        1 10000
 2      5     1        2 10000
 3      5     1        3 10000
 4      5     1        4 10000
 5      5     1        5 10000
 6      5     1        6 10000
 7      5     1        7 10000
 8      5     1        8 10000
 9      5     1        9 10000
10      5     1       10 10000
# ℹ 21,990 more rows

プログレスバーを用意します。

progressr::handlers(
1  list(
    progressr::handler_cli(
      show_after = 0,
      format = paste(
        "{cli::pb_bar}",
        "{cli::pb_percent}",
        "[in {cli::pb_elapsed}]",
        "| {cli::pb_eta_str}"
      ),
      clear = FALSE
    ),
    progressr::handler_beepr(
      update = NA_integer_,
      finish = 8L
    )
  ),
  global = TRUE
)
1
progressr::handlers()にlist形式で複数のhandlerを入れると、入れたhandlerを全部実行してくれます。今回はbeeprを使って終了時に音を鳴らしてもらいます。

Run simulation

では、実行しましょう…というか、毎度同じく、ポストを作る前に既に実行してあるので、コードだけ。 今回もfutureとfurrrで並列化して実行します。

# seed
set.seed(2026)

# multisession
future::plan(
  future::multisession,
  workers = parallelly::availableCores() - 1
)

# check
future::nbrOfWorkers()

実行部分。

progressr::with_progress({
  pb <- progressr::progressor(steps = nrow(df_jobs))
  
  res_df <- furrr::future_pmap_dfr(
    .l = df_jobs,
    .f = \(n_each, sd_g2, batch_id, n_sim) {
      temp_vec_group <- list_vec_group[[as.character(n_each)]]
      res <- f_analyse_batch(
        n_sim = n_sim,
        n_each = n_each,
        sd_g2 = sd_g2,
        vec_group = temp_vec_group,
        sig_level = .05
      ) |>
        mutate(
          n_each = n_each,
          sd_g2 = sd_g2,
          batch_id = batch_id,
          .before = 1
        )
      # update progress bar
      pb()
      # return
      res
    },
    .options = furrr::furrr_options(seed = TRUE)
  ) |>
    reframe( 
      err = sum(n_rej) / int_n_iter,
      .by = c(n_each, sd_g2, procedures)
    )
})

# end parallelization and check
future::plan(future::sequential); future::nbrOfWorkers()

かかった時間は12分50秒でした。条件数・反復数の割には思ったよりも速かったです。行列入力対応の高速な関数はシミュレーションのときにありがたいですね。

Result

結果を読み込んでおきます。（需要があるかわかりませんが、結果はgithubのこのポストのディレクトリの中に入っているので、欲しい人がいたらpost/から探してください。容量節約のためgzipにしてありますが、readr::write_rds()を使って圧縮したので、readr::read_rds()で読み込めます。）

res_df <- readr::read_rds("rep_delacre_etal_2017_fig1.gz")

Replication of Fig. 1

元論文のグラフを再現してみましょう。元論文と同じ条件のデータだけ抽出して作ります。今回のシミュレーションでは群2のSDがそのまま群間のSD比（SDR）になります。

Code 
res_df |>
  filter(
    between(n_each, 10, 80),
    sd_g2 %in% c(1.1, 1.5, 2)
  ) |>
  filter_out(procedures == "vartest") |>
  mutate(
    n_total = n_each * 2,
    across(
      .cols = c(sd_g2, procedures),
      .fns = \(x) as_factor(x)
    )
  ) |>
  ggplot(aes(x = n_total, y = err, group = interaction(sd_g2, procedures))) +
  geom_hline(
    yintercept = c(.8, .95),
    linetype = "dashed",
    color = "red"
  ) +
  geom_line(aes(linetype = procedures)) +
  geom_point(
    aes(shape = sd_g2),
    size = 2.5
  ) +
  scale_x_continuous(
    name = "total sample size N",
    breaks = seq(20, 160, 20),
    expand = expansion(add = 5)
  ) +
  scale_y_continuous(
    name = "Empirical Rejection Rate (%)",
    limits = c(0, 1),
    breaks = seq(0, 1, .2),
    labels = scales::label_number(scale = 100),
    expand = expansion(add = c(0, .05))
  ) +
  scale_linetype_manual(
    name = "Method",
    values = c("solid", "22"),
    labels = c(
      "levene" = "**MEAN** centered (Levene)",
      "bf" = "**MEDIAN** centered (Brown-Forsythe)"
    ),
    guide = guide_legend(
      nrow = 2,
1      order = 1
    )
  ) +
  scale_shape_manual(
    name = "SDR",
    values = c(19, 15, 12),
    guide = guide_legend(
      nrow = 2,
      order = 0
    )
  ) +
  theme_classic() +
  labs(
    title = "Replication of Fig. 1 of Delacre et al. (2017)",
    caption = paste0(
      "1,000,000 iterations per condition. ",
      "<i>n</i><sub>1</sub> = <i>n</i><sub>2</sub> = <i>N</i>/2. ",
      "Group 1: rnorm(<i>n</i><sub>1</sub>, mean = 0, sd = 1); ",
      "Group 2: rnorm(<i>n</i><sub>2</sub>, mean = 0, sd = SDR)."
    )
  ) +
  theme(
    plot.title.position = "plot",
    panel.border = element_rect(color = "black"),
    legend.position = "top",
    legend.text = ggtext::element_markdown(),
    plot.caption = ggtext::element_markdown(size = 7)
  )
1
グラフ上の複数の凡例の順番を変えるには、ggplot2::guide_legend(order = ...)で正の整数を渡せばいいとのこと6。

元論文のデザインに似せてみました。元論文では数値は示されていませんが、グラフを見比べる限り再現成功と言っていいでしょう。

All Results

全ての結果を確認してみましょう。見づらいけど一枚絵にした版。

Code 
res_df |>
  mutate(
    n_total = n_each * 2,
    across(
      .cols = c(sd_g2, procedures),
      .fns = \(x) as_factor(x)
    )
  ) |>
  ggplot(aes(
    x = n_total,
    y = err,
    group = interaction(sd_g2, procedures),
    colour = sd_g2,
    shape = procedures,
    linetype = procedures
  )) +
  geom_hline(
    yintercept = c(.05, .8, .95),
    linetype = "dashed"
  ) +
  geom_line() +
  geom_point() +
  scale_x_continuous(
    breaks = c(10, seq(20, 200, 20)),
    minor_breaks = seq(20, 200, 10),
    expand = expansion(add = 2)
  ) +
  scale_y_continuous(
    limits = c(0, 1),
    breaks = c(.05, .95, seq(0, 1, .1)),
    minor_breaks = seq(0, 1, .05),
    labels = scales::label_number(scale = 100),
    expand = expansion(add = c(0, .01))
  ) +
  scale_color_discrete(
    guide = guide_legend(
      nrow = 2,
      byrow = TRUE
    )
  ) +
  scale_shape_discrete(
    labels = c(
      "vartest" = "F test",
      "levene" = "Levene",
      "bf" = "Brown-Forsythe"
    ),
    guide = guide_legend(
      nrow = 2,
      byrow = TRUE
    )
  ) +
  scale_linetype_discrete(
    labels = c(
      "vartest" = "F test",
      "levene" = "Levene",
      "bf" = "Brown-Forsythe"
    ),
    guide = guide_legend(
      nrow = 2,
      byrow = TRUE
    )
  ) +
  labs(
    title = "Simulation Result: all",
    x = "total sample size *N*",
    y = "Empirical Rejection Rate (%)",
    shape = "Test",
    linetype = "Test",
    color = "SD of Group 2",
    caption = paste0(
      "1,000,000 iterations per condition. ",
      "<i>n</i><sub>1</sub> = <i>n</i><sub>2</sub> = <i>N</i>/2. ",
      "Group 1: rnorm(<i>n</i><sub>1</sub>, mean = 0, sd = 1); ",
      "Group 2: rnorm(<i>n</i><sub>2</sub>, mean = 0, sd = SD)."
    )
  ) +
  theme_bw() +
  theme(
    plot.title.position = "plot",
    axis.title.x = ggtext::element_markdown(),
    legend.position = "top",
    legend.title = element_text(size = 8),
    legend.text = element_text(size = 7),
    plot.caption = ggtext::element_markdown(size = 6)
  )

各手法とも大体似たような結果になりました。検出力はF検定＞Levene検定＞Brown-Forsythe検定の順ですね。

手法ごとに分けた版。

Code 
res_df |>
  mutate(
    n_total = n_each * 2,
    across(
      .cols = c(sd_g2, procedures),
      .fns = \(x) as_factor(x)
    )
  ) |>
  ggplot(aes(
    x = n_total,
    y = err,
    group = interaction(sd_g2, procedures),
    colour = sd_g2
  )) +
  geom_hline(
    yintercept = c(.05, .8, .95),
    linetype = "dashed"
  ) +
  geom_line() +
  geom_point() +
  scale_x_continuous(
    breaks = c(10, seq(20, 200, 20)),
    minor_breaks = seq(20, 200, 10),
    expand = expansion(add = 2)
  ) +
  scale_y_continuous(
    limits = c(0, 1),
    breaks = c(.05, .95, seq(0, 1, .1)),
    minor_breaks = seq(0, 1, .05),
    labels = scales::label_number(scale = 100),
    expand = expansion(add = c(0, .01))
  ) +
  scale_color_discrete(
    guide = guide_legend(
      nrow = 1,
      byrow = TRUE
    )
  ) +
  labs(
    title = "Simulation Result: per procedure",
    x = "total sample size *N*",
    y = "Empirical Rejection Rate (%)",
    color = "SD of Group 2",
    caption = paste0(
      "1,000,000 iterations per condition. ",
      "<i>n</i><sub>1</sub> = <i>n</i><sub>2</sub> = <i>N</i>/2. ",
      "Group 1: rnorm(<i>n</i><sub>1</sub>, mean = 0, sd = 1); ",
      "Group 2: rnorm(<i>n</i><sub>2</sub>, mean = 0, sd = SD)."
    )
  ) +
  facet_wrap(
    facet = vars(procedures),
    labeller = labeller(
      procedures = c(
        "vartest" = "F test",
        "levene" = "Levene test",
        "bf" = "Brown-Forsythe test"
      )
    )
  ) +
  theme_bw() +
  theme(
    plot.title.position = "plot",
    axis.title.x = ggtext::element_markdown(),
    legend.position = "top",
    legend.title = element_text(size = 8),
    legend.text = element_text(size = 7),
    plot.caption = ggtext::element_markdown(size = 6)
  )

サンプルサイズ\(N=10\)のとき挙動が興味深いですね。F検定とLevene検定はSD比が大きくなるにつれて検出率も高くなっていますが、BF検定はそれほど変わっていないように見受けられます。またLevene検定は\(SD=1.1\)のとき、サンプルサイズ\(N=10\)から\(N=20\)で単純に検出力が高まっているわけではないみたいです。不思議。 あと、SD比が大きくないときに、BF検定はなんか\(N=30\)のあたりがへこんでいるというか、なんかグラフが滑らかじゃない気がします。謎。反復1万回だとまああるあるだよねと思えるんですが、100万回してこれだと何か理由があるのか気になります。

標準偏差の条件ごとにプロットを作って、もう少し見やすくします。

Code 
res_list_plot <- res_df |>
  mutate(
    n_total = n_each * 2,
    across(
      .cols = c(sd_g2, procedures),
      .fns = \(x) as_factor(x)
    )
  ) |>
  group_by(sd_g2) |>
  group_map(
    .f = \(x, idx) {
      temp_gg <- x |>
        ggplot(
          aes(
            x = n_total,
            y = err,
            group = procedures,
            color = procedures,
            shape = procedures
          )
        )
      if(idx$sd_g2 == 1) {
        temp_gg <- temp_gg +
          geom_hline(
            yintercept = .05,
            linetype = "dashed"
          ) +
          scale_y_continuous(
            limits = c(0, 1),
            breaks = c(.05, seq(0, 1, .1)),
            minor_breaks = seq(0, 1, .05),
            labels = scales::label_number(scale = 100),
            expand = expansion(add = c(0, .01))
          )
      } else {
        temp_gg <- temp_gg +
          geom_hline(
            yintercept = c(.8, .95),
            linetype = "dashed"
          ) +
          scale_y_continuous(
            limits = c(0, 1),
            breaks = c(seq(0, 1, .1), .95),
            minor_breaks = seq(0, 1, .05),
            labels = scales::label_number(scale = 100),
            expand = expansion(add = c(0, .01))
          )
      }
      temp_gg +
        geom_line() +
        geom_point(size = 2.5) +
        scale_x_continuous(
          breaks = c(10, seq(20, 200, 20)),
          minor_breaks = seq(30, 190, 20),
          expand = expansion(add = 2)
        ) +
        scale_color_discrete(
          name = "Method",
          labels = c(
            "vartest" = "F test",
            "levene" = "Levene Test",
            "bf" = "Brown-Forsythe Test"
          ),
          guide = guide_legend(nrow = 1)
        ) +
        scale_shape_discrete(
          name = "Method",
          labels = c(
            "vartest" = "F test",
            "levene" = "Levene Test",
            "bf" = "Brown-Forsythe Test"
          ),
          guide = guide_legend(nrow = 1)
        ) +
        labs(
          title = paste0(
            "Simulation Result: SD of Group 2 = ",
            idx$sd_g2
          ),
          x = "total sample size *N*",
          y = "Empirical Rejection Rate (%)",
          color = "SD Ratio",
          caption = paste0(
            "1,000,000 iterations per condition. ",
            "<i>n</i><sub>1</sub> = <i>n</i><sub>2</sub> = <i>N</i>/2. ",
            "Group 1: rnorm(<i>n</i><sub>1</sub>, mean = 0, sd = 1); ",
            "Group 2: rnorm(<i>n</i><sub>2</sub>, mean = 0, sd = SD)."
          )
        ) +
        theme_bw() +
        theme(
          plot.title.position = "plot",
          axis.title.x = ggtext::element_markdown(),
          legend.position = "top",
          legend.title = element_text(size = 8),
          legend.text = element_text(size = 7),
          plot.caption = ggtext::element_markdown(size = 6)
        )
    }
  )

まずは\(SD=1\)の条件。この条件では群1のSDと群2のSDは同じなので、ERRはType I error rateを意味しています。

Code 
res_list_plot[1]

サンプルサイズが少ないときの挙動が各手法で異なりました。F検定は5％程度を維持できていたのに対して、Levene検定は若干5％を上回り気味、Brown-Forsythe検定は5％下回り気味なように見受けられます。実際のERRは以下の通りです。

Code 
res_df |>
  filter(sd_g2 == 1) |>
  mutate(
    mcse = sqrt(err * (1 - err) / int_n_iter),
    .by = c(n_each, procedures)
  ) |>
  pivot_wider(
    names_from = procedures,
    names_vary = "slowest",
    names_glue = "{procedures}_{.value}",
    values_from = c(err, mcse),
    values_fn = \(x) num(x, scale = 100, label = "%")
  ) |>
  relocate(sd_g2, .before = 1)
# A tibble: 20 × 8
   sd_g2 n_each vartest_err vartest_mcse levene_err levene_mcse bf_err bf_mcse
   <dbl>  <dbl>           %            %          %           %      %       %
 1     1      5        5.01       0.0218       7.37      0.0261  0.644 0.00800
 2     1     10        4.99       0.0218       5.98      0.0237  3.57  0.0186 
 3     1     15        4.97       0.0217       5.66      0.0231  3.31  0.0179 
 4     1     20        5.03       0.0219       5.51      0.0228  4.03  0.0197 
 5     1     25        5.00       0.0218       5.39      0.0226  4.00  0.0196 
 6     1     30        4.98       0.0218       5.31      0.0224  4.28  0.0202 
 7     1     35        4.97       0.0217       5.26      0.0223  4.27  0.0202 
 8     1     40        4.98       0.0218       5.21      0.0222  4.41  0.0205 
 9     1     45        5.03       0.0219       5.23      0.0223  4.45  0.0206 
10     1     50        5.00       0.0218       5.19      0.0222  4.55  0.0208 
11     1     55        5.00       0.0218       5.20      0.0222  4.58  0.0209 
12     1     60        4.99       0.0218       5.14      0.0221  4.59  0.0209 
13     1     65        4.99       0.0218       5.11      0.0220  4.58  0.0209 
14     1     70        5.00       0.0218       5.16      0.0221  4.68  0.0211 
15     1     75        4.98       0.0218       5.11      0.0220  4.66  0.0211 
16     1     80        5.02       0.0218       5.12      0.0220  4.71  0.0212 
17     1     85        4.98       0.0217       5.08      0.0220  4.69  0.0211 
18     1     90        5.00       0.0218       5.10      0.0220  4.73  0.0212 
19     1     95        5.03       0.0218       5.12      0.0220  4.74  0.0213 
20     1    100        4.99       0.0218       5.08      0.0220  4.74  0.0213

個人的な感覚として±1%未満に収まっていればOKなんじゃないかなというところではありますが、こうしてみると、Levene検定は少し高めに出がちでBF検定は低めに出がちですね。生のp値を保存しておいてECDFのグラフを描いてみたらもっと面白いことがわかるかもしれませんが、今回はやっていないのでこれ以上は何とも言えません。

群間でSDが異なる条件を見てみます。この場合のERRは検出力を意味しています。

Code 
res_list_plot[-1]

当然ですが、サンプルサイズが少ないときは検出力は低く、サンプルサイズが大きくなるにしたがって検出力も高くなりました。また、SD比が大きければ大きいほど、ある程度の検出力（例えば80%だったり95%だったり）に達するのに必要なサンプルサイズは少ないです。 逆にいえば、SD比が小さめのときにはより多くのサンプルサイズがないと正しく判断できないわけです。 Delacre et al. （2017） は、SD比が1.1のときでさえWelchのt検定を使った方がいいのに、Levene検定は総サンプルサイズ160でも十数パーセントの検出力しかないことを指摘しています。検出力がその程度なら、2段階t検定を行ったときにはたいていStudentのt検定を実施することになりますが、等分散ではないのでType I error rateはαからズレてしまうでしょう。（実際には今回みたいに正規分布データかつ各群等サンプルサイズならStudentのt検定はまだ耐えなのですが、等サンプルサイズでない場合にはStudentは異分散に脆弱です。）

Extra (Simulation 2)

等サンプルサイズじゃなくするくらいなら簡単にできないか？と思ったのと、サンプルサイズを等しくしないならSDも大小の条件をつくらないとなぁと思いやってみました。条件ごとの反復数は100万だと確実に時間がかかるので、10万回にしました。MCSEを考えると10万でも十分精度は高いと思います。

Set up

群1のサンプルサイズ\(n_1\)は変わらずで、群2のサンプルサイズ\(n_2\)は群1の半分と2倍を追加します。半分の場合、下一桁が5だと割り切れなくなってしまうのと、matrixTests::col_brownforsythe()が\(n_i\geq3\)を要求するので、小数点以下は切り上げします。例えば\(n_1=5\)のとき、\(n_2=3\)です。

また群2のSDに1.1から2の逆数を用意して、\(1:0.5\)から\(1:2\)まで条件を増やします。

df_design2 <- expand_grid(
  n_g1 = seq(5, 100, 5),
  sd_g2 = c(
    1 / seq(2, 1, -.1),
    seq(1.1, 2, .1)
  )
) |>
  mutate(
    n_g2 = map(
      .x = n_g1,
      .f = \(x) ceiling(x * c(1/2, 1, 2))
    ),
    .after = n_g1
  ) |>
  unnest(cols = n_g2) |>
  mutate(
    name_cond = str_c(n_g1, n_g2, sep = "_")
  ) |>
  arrange(n_g1, n_g2, sd_g2)

# check
df_design2
# A tibble: 1,260 × 4
    n_g1  n_g2 sd_g2 name_cond
   <dbl> <dbl> <dbl> <chr>    
 1     5     3 0.5   5_3      
 2     5     3 0.526 5_3      
 3     5     3 0.556 5_3      
 4     5     3 0.588 5_3      
 5     5     3 0.625 5_3      
 6     5     3 0.667 5_3      
 7     5     3 0.714 5_3      
 8     5     3 0.769 5_3      
 9     5     3 0.833 5_3      
10     5     3 0.909 5_3      
# ℹ 1,250 more rows

グループ分け用のベクトルのリストも更新です。

list_vec_group2 <- df_design2 |>
  distinct(n_g1, n_g2, name_cond) |>
  mutate(
    res = map2(
      .x = n_g1,
      .y = n_g2,
      .f = \(x, y) {
        rep(1:2, times = c(x, y)) |>
          as.factor()
      }
    ) |>
      set_names(name_cond)
  ) |>
  pull(res)

# check partially
head(list_vec_group2)
$`5_3`
[1] 1 1 1 1 1 2 2 2
Levels: 1 2

$`5_5`
 [1] 1 1 1 1 1 2 2 2 2 2
Levels: 1 2

$`5_10`
 [1] 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2
Levels: 1 2

$`10_5`
 [1] 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2
Levels: 1 2

$`10_10`
 [1] 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2
Levels: 1 2

$`10_20`
 [1] 1 1 1 1 1 1 1 1 1 1 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
Levels: 1 2

分析に使う関数はそのまま流用で、バッチ処理の方だけ更新します。

f_analyse_batch2 <- function(n_sim, n_g1, n_g2, sd_g2, vec_group, sig_level = .05) {
  # data generation
  temp_mat <- replicate(
    n_sim,
    expr = c(
      rnorm(n_g1, mean = 0, sd = 1),
      rnorm(n_g2, mean = 0, sd = sd_g2)
    )
  )
  # analyse
  f_analyse(temp_mat, g = vec_group) |>
    (\(x) colSums(x < sig_level, na.rm = TRUE)) () |>
    tibble::enframe(name = "procedures", value = "n_rej")
}

反復数は10万でブロックサイズは25000にします。5万だと並列化した時にRAMが足りなくなりそうでした。

int_n_iter2 <- 100000
int_n_block_size2 <- 25000

df_jobs2 <- df_design2 |>
  crossing(batch_id = seq_len(int_n_iter2 / int_n_block_size2)) |>
  mutate(n_sim = int_n_block_size2)

# check
df_jobs2
# A tibble: 5,040 × 6
    n_g1  n_g2 sd_g2 name_cond batch_id n_sim
   <dbl> <dbl> <dbl> <chr>        <int> <dbl>
 1     5     3 0.5   5_3              1 25000
 2     5     3 0.5   5_3              2 25000
 3     5     3 0.5   5_3              3 25000
 4     5     3 0.5   5_3              4 25000
 5     5     3 0.526 5_3              1 25000
 6     5     3 0.526 5_3              2 25000
 7     5     3 0.526 5_3              3 25000
 8     5     3 0.526 5_3              4 25000
 9     5     3 0.556 5_3              1 25000
10     5     3 0.556 5_3              2 25000
# ℹ 5,030 more rows

プログレスバーの構成はそのままで。

Run Simulation

ということで回してきました。

# seed
set.seed(202608)

# parallelization and check
plan(multisession, workers = parallelly::availableCores() - 1); nbrOfWorkers()

# run
progressr::with_progress({
  pb <- progressr::progressor(steps = nrow(df_jobs2))
  
  res_df2 <- future_pmap_dfr(
    .l = df_jobs2,
    .f = \(n_g1, n_g2, sd_g2, name_cond, batch_id, n_sim) {
      # group vector
      temp_vec_group <- list_vec_group2[[name_cond]]
      # analyse
      res <- f_analyse_batch2(
        n_sim = n_sim,
        n_g1 = n_g1,
        n_g2 = n_g2,
        sd_g2 = sd_g2,
        vec_group = temp_vec_group,
        sig_level = .05
      ) |>
        mutate(
          n_g1 = n_g1,
          n_g2 = n_g2,
          sd_g2 = sd_g2,
          name_cond = name_cond,
          batch_id = batch_id,
          .before = 1
        )
      # update pb
      pb()
      # return
      res
    },
    .options = furrr::furrr_options(seed = TRUE)
  ) |>
    reframe(
      err = sum(n_rej) / int_n_iter2,
      .by = c(n_g1, n_g2, sd_g2, procedures)
    )
})

# end parallelization and check
plan(sequential); nbrOfWorkers()

8分半ちょっとかかりました。条件が多いとやっぱり時間かかりますね。

Result (Sim. 2)

この結果も先ほどの結果のファイルと同じところにあるので、欲しい人がいたらどうぞ。（rdsファイルそのままだと166KBもあったのでgzipに圧縮してあります。先ほどと同様にreadr::write_rds()で解凍して読み込めます。）

res_df2 <- readr::read_rds("rep_delacre_etal_2017_fig1_sim2.gz")

出来ればなるべく一枚絵でグラフを出したいんですが、条件が多すぎて難しいです。とりあえずサンプルサイズの関係性ごとにまずは出すことにします。

Code 
list_plot_sim2_per_n <- res_df2 |>
  mutate(
    condition = case_when(
      n_g1 > n_g2 ~ "<i>n</i><sub>1</sub> > <i>n</i><sub>2</sub>",
      n_g1 == n_g2 ~ "<i>n</i><sub>1</sub> = <i>n</i><sub>2</sub>",
      n_g1 < n_g2 ~ "<i>n</i><sub>1</sub> < <i>n</i><sub>2</sub>"
    ),
    n_ratio = str_c(n_g1, n_g2, sep = ":") |>
      as_factor(),
    across(
      .cols = c(condition, sd_g2, procedures),
      .fns = \(x) as_factor(x)
    )
  ) |>
  group_by(condition) |>
  group_map(
    .f = \(x, idx) {
      x |>
        ggplot(
          aes(
            x = n_ratio,
            y = err,
            group = interaction(sd_g2, procedures),
            color = sd_g2,
            shape = procedures
          )
        ) +
        geom_hline(
          yintercept = c(.05, .8, .95),
          linetype = "dashed"
        ) +
        geom_line() +
        geom_point() +
        scale_x_discrete(
          expand = expansion(add = .2)
        ) +
        scale_y_continuous(
          limits = c(0, 1),
          breaks = c(.05, .95, seq(0, 1, .1)),
          minor_breaks = seq(0, 1, .05),
          labels = scales::label_number(scale = 100),
          expand = expansion(add = c(0, .01))
        ) +
        scale_color_discrete(
          guide = guide_legend(
            nrow = 2,
            byrow = TRUE
          ),
          labels = c(
            paste0("1/", seq(2, 1.1, -.1)),
            seq(1, 2, .1)
          )
        ) +
        scale_shape_discrete(
          guide = "none"
        ) +
        facet_wrap(
          vars(procedures),
          labeller = labeller(
            procedures = c(
              "vartest" = "F test",
              "levene" = "Levene",
              "bf" = "Brown-Forsythe"
            )
          )
        ) +
        labs(
          title = paste(
            "Result of Simulation 2:",
            idx$condition
          ),
          x = "Sample size <i>n</i><sub>1</sub>:<i>n</i><sub>2</sub>",
          y = "Empirical Rejection Rate (%)",
          color = "SD of Group 2",
          caption = paste0(
            "100,000 iterations per condition. ",
            "Group 1: rnorm(<i>n</i><sub>1</sub>, mean = 0, sd = 1); ",
            "Group 2: rnorm(<i>n</i><sub>2</sub>, mean = 0, sd = SD)."
          )
        ) +
        theme_bw() +
        theme(
          plot.title = ggtext::element_markdown(),
          plot.title.position = "plot",
          axis.title.x = ggtext::element_markdown(),
          axis.text.x = element_text(
            hjust = 1,
            angle = 60,
            size = 7
          ),
          legend.position = "top",
          legend.title = element_text(size = 8),
          legend.text = element_text(size = 7),
          legend.key.spacing.y = unit(0, "pt"),
          legend.box.spacing = unit(1, "pt"),
          plot.caption = ggtext::element_markdown(size = 6)
        )
    }
  )

list_plot_sim2_per_n

（普通に見づらくてよくない。）反復数は1条件当たり10万回に減っていますが、\(n_1=n_2\)のグラフを見る限り問題なさそうですね。しかもこの条件は、SDの大きさでほぼきれいに対称関係になっていますね。

さて、各グラフで横軸の\(n_1\)を固定した時、\(n_1 > n_2\)の条件ではグラフの左上のスペースが広めです。当然ですが、単純にサンプルサイズが減っているので、等サンプルサイズに比べれば検出力も低めに出がちです。また、どの検定も群2のSDが大きい場合、すなわちサンプルサイズが小さい群のSDが大きい場合の方に検出力が高くなっています。そして、その差の大きさはF検定＞Levene検定＞BF検定のように見受けられます。

\(n_1 < n_2\)の条件ではサンプルサイズが大きくなっているので、等サンプルサイズの条件よりも検出力も高めに出がちです。また、どの検定も群2のSDが小さい場合の方が検出力が高くなっています。その差の大きさの関係も同じくF検定＞Levene検定＞BF検定のように見受けられますが、こちらの条件の方が差が小さいように見えます。

SDごとに分けてもう少し見やすくします。まずは\(SD = 1\)のときで、これは実質的にはType I errorの確率。

Code 
list_plot_sim2_per_sd <- res_df2 |>
  mutate(
    condition = case_when(
      n_g1 > n_g2 ~ "<i>n</i><sub>1</sub> > <i>n</i><sub>2</sub>",
      n_g1 == n_g2 ~ "<i>n</i><sub>1</sub> = <i>n</i><sub>2</sub>",
      n_g1 < n_g2 ~ "<i>n</i><sub>1</sub> < <i>n</i><sub>2</sub>"
    ),
    n_ratio = str_c(
      str_pad(n_g1, width = 3, side = "left"),
      str_pad(n_g2, width = 3, side = "right"),
      sep = ":"
    ) |>
      as_factor(),
    across(
      .cols = c(condition, sd_g2, procedures),
      .fns = \(x) as_factor(x)
    )
  ) |>
  nest(.by = sd_g2) |>
  mutate(
    SD = c(paste0("1/", seq(2, 1.1, -.1)), seq(1, 2, .1)),
    .after = 1
  ) |>
  reframe(
    res = map2(
      .x = data,
      .y = SD,
      .f = \(x, y) {
        temp_gg <- x |>
          ggplot(aes(x = n_ratio, y = err, group = procedures, colour = procedures))
        if(y == "1") {
          temp_gg <- temp_gg +
            geom_hline(
              yintercept = .05,
              linetype = "dashed"
            ) +
            scale_y_continuous(
              limits = c(0, 1),
              breaks = c(.05, seq(0, 1, .1)),
              minor_breaks = seq(0, 1, .05),
              labels = scales::label_number(scale = 100),
              expand = expansion(add = c(0, .01))
            )
        } else {
          temp_gg <- temp_gg +
            geom_hline(
              yintercept = c(.8, .95),
              linetype = "dashed"
            ) +
            scale_y_continuous(
              limits = c(0, 1),
              breaks = c(seq(0, 1, .1), .95),
              minor_breaks = seq(0, 1, .05),
              labels = scales::label_number(scale = 100),
              expand = expansion(add = c(0, .01))
            )
        }
        temp_gg +
          geom_line() +
          geom_point() +
          scale_color_discrete(
            labels = c(
              "vartest" = "F test",
              "levene" = "Levene",
              "bf" = "Brown-Forsythe"
            )
          ) +
          facet_wrap(
            facet = vars(condition),
            scales = "free_x"
          ) +
          labs(
            title = paste(
              "Result of Simulation 2: SD of Group 2 =",
              y
            ),
            x = "Sample size <i>n</i><sub>1</sub>:<i>n</i><sub>2</sub>",
            y = "Empirical Rejection Rate (%)",
            color = "Test",
            caption = paste0(
              "100,000 iterations per condition. ",
              "Group 1: rnorm(<i>n</i><sub>1</sub>, mean = 0, sd = 1); ",
              "Group 2: rnorm(<i>n</i><sub>2</sub>, mean = 0, sd = SD)."
            )
          ) +
          theme_bw() +
          theme(
            plot.title.position = "plot",
            axis.title.x = ggtext::element_markdown(),
            axis.text.x = element_text(
              angle = 60,
              hjust = 1,
              size = 7
            ),
            strip.text = ggtext::element_markdown(),
            legend.position = "top",
            legend.title = element_text(size = 8),
            legend.text = ggtext::element_markdown(size = 7),
            legend.key.spacing.y = unit(0, "pt"),
            legend.box.spacing = unit(1, "pt"),
            plot.caption = ggtext::element_markdown(size = 6)
          )
      }
    ) |>
      setNames(nm = SD)
  ) |>
  pull(res)

list_plot_sim2_per_sd[["1"]]

F検定はサンプルサイズの等しさに関わらず5%程度を維持していました。一方で、Levene検定はサンプルサイズが少ないとちょっと高めに出がち、BF検定はサンプルサイズが少ないと低めに出がちという結果になりました。Levene検定よりもBF検定の方が頑健と評されるのはこういったところも理由の一つでしょうか（外れ値とかのこともあるだろうけど）。\(n_1\)が少なめのときのERRを見てみます。

Code 
res_df2 |>
  filter(n_g1 <= 40, sd_g2 == 1) |>
  mutate(
    mcse = sqrt(err * (1 - err) / int_n_iter2)
  ) |>
  pivot_wider(
    names_from = procedures,
    names_vary = "slowest",
    names_glue = "{procedures}_{.value}",
    values_from = c(err, mcse),
    values_fn = \(x) num(x, scale = 100, label = "%")
  ) |>
  select(-sd_g2) |>
  print(n = 99)
# A tibble: 24 × 8
    n_g1  n_g2 vartest_err vartest_mcse levene_err levene_mcse bf_err bf_mcse
   <dbl> <dbl>           %            %          %           %      %       %
 1     5     3        5.05       0.0692       7.33      0.0824  0.094 0.00969
 2     5     5        4.91       0.0683       7.40      0.0828  0.617 0.0248 
 3     5    10        5.07       0.0694       6.35      0.0771  2.46  0.0490 
 4    10     5        5.01       0.0690       6.28      0.0767  2.43  0.0487 
 5    10    10        4.99       0.0688       6.02      0.0752  3.64  0.0592 
 6    10    20        4.90       0.0683       5.53      0.0723  3.74  0.0600 
 7    15     8        4.98       0.0688       5.89      0.0744  3.46  0.0578 
 8    15    15        5.09       0.0695       5.66      0.0730  3.33  0.0568 
 9    15    30        4.94       0.0686       5.44      0.0717  3.95  0.0616 
10    20    10        5.03       0.0691       5.60      0.0727  3.81  0.0605 
11    20    20        5.04       0.0692       5.53      0.0723  4.07  0.0625 
12    20    40        5.05       0.0693       5.35      0.0712  4.34  0.0644 
13    25    13        5.00       0.0689       5.54      0.0724  3.73  0.0599 
14    25    25        4.98       0.0688       5.37      0.0713  3.97  0.0617 
15    25    50        4.94       0.0685       5.15      0.0699  4.27  0.0639 
16    30    15        4.98       0.0688       5.44      0.0717  3.97  0.0618 
17    30    30        5.00       0.0689       5.19      0.0702  4.13  0.0629 
18    30    60        4.95       0.0686       5.09      0.0695  4.42  0.0650 
19    35    18        5.00       0.0689       5.34      0.0711  4.12  0.0629 
20    35    35        5.05       0.0693       5.24      0.0704  4.28  0.0640 
21    35    70        5.07       0.0694       5.34      0.0711  4.64  0.0665 
22    40    20        4.94       0.0685       5.28      0.0707  4.26  0.0639 
23    40    40        5.05       0.0693       5.32      0.0710  4.47  0.0653 
24    40    80        4.92       0.0684       5.12      0.0697  4.60  0.0663

\(n_1=5,10\)のときのLevene検定は6~7%台、一方BF検定は~3%でした。F検定が安定して5%付近にいるのと比べてみると、興味深いです。

次は検出力の方です。

Code 
list_plot_sim2_per_sd[!names(list_plot_sim2_per_sd) == "1"]

当然ながら、全体のサンプルサイズが大きい条件の方が、グラフの傾きが大きいですね。

では、総サンプルサイズが同じ場合で比較するとどうでしょうか。今回の場合、総サンプルサイズが30の倍数のときに総サンプルサイズが揃えられる（例：\(N=30\)のとき、\((n_1,n_2)=(20,10), (15,15), (10,20)\)）ので見てみましょう。左右対称にするために、グラフの横軸は群2のSDをlog2でとっています。

Code 
res_df2 |>
  filter(
    (n_g1 + n_g2) %in% seq(30, 150, 30)
  ) |>
  mutate(
    n_total = n_g1 + n_g2,
    greater = case_when(
      n_g1 > n_g2 ~ "g1",
      n_g1 == n_g2 ~ "eq",
      n_g1 < n_g2 ~ "g2"
    ) |>
      factor(levels = c("g1", "eq", "g2")),
    procedures = factor(
      procedures,
      levels = c("vartest", "levene", "bf")
    )
  ) |>
  group_by(n_total) |> 
  group_map(
    .f = \(x, idx) {
      x |>
        ggplot(
          aes(
            x = sd_g2,
            y = err,
            group = greater,
            color = greater
          )
        ) +
        geom_vline(
          xintercept = 1,
          linetype = "dashed"
        ) +
        geom_line() +
        geom_point() +
        scale_x_continuous(
          transform = "log2",
          breaks = c(
            1 / seq(2, 1.1, -.1),
            seq(1, 2, .1)
          ),
          minor_breaks = NULL,
          labels = c(
            paste0("1/", seq(2, 1.1, -.1)),
            seq(1, 2, .1)
          ),
          expand = expansion(add = .02)
        ) +
        scale_y_continuous(
          limits = c(0, 1),
          breaks = seq(0, 1, .1),
          labels = scales::label_number(scale = 100),
          expand = expansion(add = c(0, .01))
        ) +
        scale_color_discrete(
          labels = c(
            "g1" = "<i>n</i><sub>1</sub> > <i>n</i><sub>2</sub>",
            "eq" = "<i>n</i><sub>1</sub> = <i>n</i><sub>2</sub>",
            "g2" = "<i>n</i><sub>1</sub> < <i>n</i><sub>2</sub>"
          )
        ) +
        facet_wrap(
          vars(procedures),
          labeller = labeller(
            procedures = c(
              "vartest" = "F test",
              "levene" = "Levene",
              "bf" = "Brown-Forsythe"
            )
          )
        ) +
        labs(
          title = paste(
            "Result of Sim. 2: <i>N</i> =",
            idx$n_total
          ),
          x = "SD of Group 2",
          y = "Empirical Rejection Rate (%)",
          color = "Sample Size",
          caption = paste0(
            "100,000 iterations per condition. ",
            "Group 1: rnorm(<i>n</i><sub>1</sub>, mean = 0, sd = 1); ",
            "Group 2: rnorm(<i>n</i><sub>2</sub>, mean = 0, sd = SD).",
            "<br>The x-axis is plotted on a log2 scale."
          )
        ) +
        theme_bw() +
        theme(
          plot.title = ggtext::element_markdown(),
          plot.title.position = "plot",
          axis.text.x = element_text(
            angle = 60,
            hjust = 1,
            size = 7
          ),
          legend.position = "top",
          legend.title = element_text(size = 8),
          legend.text = ggtext::element_markdown(size = 7),
          legend.key.spacing.y = unit(0, "pt"),
          legend.box.spacing = unit(1, "pt"),
          plot.caption = ggtext::element_markdown(size = 6)
        )
    }
  )

点線は等分散のときです。群2のSDが相対的に小さいとき（グラフの左半分）、どの手法も群2のサンプルサイズが小さい場合の検出力が一番低いです。群2のSDが相対的に大きいとき（グラフの右半分）は、サンプルサイズの大小について入れ替わった関係になってます。\(N=30\)で\(1/1.6\leq SD\leq1.6\)のときを除いて（というかここは何でこうなったんだろう？）、基本的には等サンプルサイズのときが一番検出力が高くなっているように見受けられます。 検出力の最大最小差については、F検定＞Levene検定＞BF検定になっているように見受けられます。

Conclusion

今回は Delacre et al. （2017） のFig. 1の、等分散性の検定のシミュレーションをRで再現＋αをしてみました。グラフの再現についてはうまくできたので特に言うことなしです。高速に処理してくれる関数があったので、条件数と反復数の割には速くできたと思います。

＋αの方ですが、やはりシミュレーションの条件が多くなると、グラフでコンパクトにあらわすのが難しいです。可視化の仕方をもっと考えないといけませんね。それかあきらめて表にするか。あと生のp値を保存していないので、ECDFが書けないのもちょっと損な気がします。これについてはSimDesignでシミュレーションを組んだ方が絶対に効率がいいんですが、HTMLの記事にするときに圧倒的量の結果のrdsファイルをどうするかが悩ましいところ。

中身に関していえば、小サンプルサイズor異分散の程度が大きくないときは等分散性の検定の検出力もそれほど高まらないのがよくわかりました。それゆえに不等分散に脆弱なStudentのt検定を選択しやすくなるため、2段階t検定で使うのはよくないとする主張もわかります。これについては左右非対称なゆがんだ分布も持ってきて確かめる必要がありますね。元論文で引用されていた Nordstokke & Zumbo （2007） は\(\chi^2\)分布を使って歪度のパターンを作っていたり、 Brown-Forsythe検定の Brown & Forsythe （1974） では他の等分散性の検定やデータ発生の分布を使っていたりしていました。今までやったシミュレーションは基本的に正規分布以外使ったことがないので、ここら辺はもっと勉強して取り組みたいです。

なお、Delacre et al. （2017） の本論はt検定のシミュレーションだったので、それらもやってみたいですね。（だいぶ前に正規分布データで、不等サンプルサイズ、異分散におけるStudentとWelchの比較のシミュレーションしたことはあります。）ただ、Appendixファイルを見たら条件がかなり多かったので時間がかかりそうです。これもSimDesignでやる方が多分楽。ただ、2段階t検定に関しては Zimmerman （1996） や Zimmerman （2004） のシミュレーションがあるので、それの再現すれば事足りる気もします。読み直してそのうち再現をやりたいです。

Session Information

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] furrr_0.4.0         future_1.75.0       matrixTests_0.2.3.1 Rfast_2.1.5.2      
 [5] RcppParallel_6.2.0  zigg_0.0.2          Rcpp_1.1.2          lubridate_1.9.5    
 [9] forcats_1.0.1       stringr_1.6.0       dplyr_1.2.1         purrr_1.2.2        
[13] readr_2.2.0         tidyr_1.3.2         tibble_3.3.1        ggplot2_4.0.3      
[17] tidyverse_2.0.0    

loaded via a namespace (and not attached):
 [1] gtable_0.3.6       xfun_0.60          htmlwidgets_1.6.4  bench_1.1.4        tzdb_0.5.0        
 [6] vctrs_0.7.3        tools_4.5.3        generics_0.1.4     parallel_4.5.3     pacman_0.5.1      
[11] pkgconfig_2.0.3    RColorBrewer_1.1-3 S7_0.2.2           lifecycle_1.0.5    compiler_4.5.3    
[16] farver_2.1.2       textshaping_1.0.5  codetools_0.2-20   litedown_0.10      htmltools_0.5.9   
[21] yaml_2.3.12        profmem_0.7.0      pillar_1.11.1      parallelly_1.48.0  commonmark_2.0.0  
[26] tidyselect_1.2.1   digest_0.6.39      stringi_1.8.9      listenv_1.0.0      labeling_0.4.3    
[31] fastmap_1.2.0      grid_4.5.3         cli_3.6.6          magrittr_2.0.5     utf8_1.2.6        
[36] withr_3.0.3        scales_1.4.0       timechange_0.4.0   rmarkdown_2.31     matrixStats_1.5.0 
[41] globals_0.19.1     otel_0.2.0         ggtext_0.1.2       ragg_1.5.2         hms_1.1.4         
[46] evaluate_1.0.5     knitr_1.51         markdown_2.0       rlang_1.3.0        gridtext_0.1.6    
[51] glue_1.8.1         xml2_1.6.0         rstudioapi_0.19.0  jsonlite_2.0.0     R6_2.6.1          
[56] systemfonts_1.3.2

References

Brown, M. B., & Forsythe, A. B. (1974). Robust Tests for the Equality of Variances. Journal of the American Statistical Association, 69(346), 364–367. https://doi.org/10.2307/2285659
Delacre, M., Lakens, D., & Leys, C. (2017). Why Psychologists Should by Default Use Welch’s t-test Instead of Student’s t-test. International Review of Social Psychology, 30(1). https://doi.org/10.5334/irsp.82
Delacre, M., Lakens, D., & Leys, C. (2022). Correction: Why Psychologists Should by Default Use Welch’s t-test Instead of Student’s t-test. International Review of Social Psychology, 35(1). https://doi.org/10.5334/irsp.661
Nordstokke, D. W., & Zumbo, B. D. (2007). A Cautionary Tale about Levene’s Tests for Equal Variances. Journal of Educational Research & Policy Studies, 7(1), 1–14.
Zimmerman, D. W. (1996). Some Properties of Preliminary Tests of Equality of Variances in the Two-Sample Location Problem. The Journal of General Psychology, 123(3), 217–231. https://doi.org/10.1080/00221309.1996.9921274
Zimmerman, D. W. (2004). A note on preliminary tests of equality of variances. British Journal of Mathematical and Statistical Psychology, 57(1), 173–181. https://doi.org/10.1348/000711004849222

Footnotes

  1. 論文の訂正箇所については、 Delacre et al. （2022） 参照のこと。↩︎

  2. ちなみに*_f_var()という2群の等分散性のF検定をやってくれる関数もあります。ただし、入力を群1のデータ行列と群2のデータ行列に分けなきゃいけないのと、Rfastの方が約9倍速かったので、採用しませんでした。matrixTestsは戻り値でdfを構成している分、速くなりきらないかもしれません↩︎

  3. Rfastにもvartests()をというほとんど同じ機能の関数があるのですが、この記事を書いている時点でのバージョン（2.1.5.2）だと、正しい結果が返ってきません。Issueで修正をお願いしました。↩︎

  4. ベンチをいろいろ回した結果、実行速度とメモリ消費量の観点から1000～5万くらいがよさそうでした。10万とかにしても十分速いのですが、メモリ消費量がGB単位で大きくなっていって並列化した時にRAMが心配で辞めました。中で行列を作っているので仕方ないですね。↩︎

  5. 2万行のdfを持っておくのもRAMがもったいない気がします。SimDesignパッケージを使うと、この辺は特に気にせずシミュレーションを回せるんですけどね。↩︎

  6. 参考：https://stackoverflow.com/questions/59686541/ggplot-order-of-multiple-legends-geom-line-and-geom-line-plus-geom-point↩︎