はじめてのシミュレーション

2026年4月8日 R with Pharma Lab はじめてのシミュレーション 中川 雄貴

Introduction R初心者 • 簡単なプログラムでRを学ぶ • SASとの違いを知る シミュレーション • 視覚化による統計手法の理解促進 • 実務での活用 ✓ 症例数設計の妥当性確認 ✓ 推定値の精度および不確実性の評価 Copyright©EPS All rights reserved. 2

2標本t検定 Student or Welch ? 2標本t検定 等分散の仮定 Yes Studentのt検定 Copyright©EPS All rights reserved. No Welchのt検定 3

2標本t検定 Student or Welch ? サンプルサイズ，母平均が同じで，母標準偏差が異なる場合… 本当は２群間の母平均には差がないはず！ ⇒ 誤って有意と判定されるαエラーがどうなるかを確認． サンプルサイズ 群A 25 群B 25 母平均 0 0 母標準偏差 1 0.1, 0.5, 1, 2, 5, 10 有意水準 0.05 シミュレーション回数 10000 Copyright©EPS All rights reserved. 4

2標本t検定 Student or Welch ?
n <- c(25, 25)
mu <- c(0, 0)
sigma_a <- 1
sigma_b <- c(0.1, 0.5, 1, 2, 5, 10)
事前準備

l <- length(sigma_b)
alpha <- 0.05
iter <- 10000

pre-allocation
pvalue <- array(NA, dim = c(l, iter))
set.seed(1234)
for (i in 1 : l) {
for (j in 1 : iter) {
y_a <- rnorm(n[1], mu[1], sigma_a)
y_b <- rnorm(n[2], mu[2], sigma_b[i])
result <- t.test(y_a, y_b, var.equal = TRUE)
pvalue[i, j] <- result$p.value
}
}

ループ等で結果を溜める入れ物（ベクトル・行列・リスト等）を，
最初に必要サイズで作っておく」 というパフォーマンス最適化
2群の平均差のt検定を、繰り返しシミュレーションを実施．
p値をpvalueに蓄積していく．

var.equal = TRUE： Studentのt検定
var.equal = FALSE： Welchのt検定（デフォルト）

type1error <- rep(0, l)
for (i in 1 : l) {
type1error[i] <- ifelse(pvalue[i, ] < 0.05, 1, 0) |> mean()
}
Copyright©EPS All rights reserved.

p値が0.05未満かどうかを判定

5

0.10 2標本t検定 Student or Welch ? 0.06 0.04 0.00 0.02 Alpha Error 0.08 Student Welch 0.1 0.5 1 2 5 10 Population SD in Group B 思ったよりもαエラーは増加しない…? Copyright©EPS All rights reserved. 6

2標本t検定 Student or Welch ? 母平均が同じで，母標準偏差が異なる場合… 本当は２群間の母平均には差がないはず！ ⇒サンプルサイズが異なる場合，誤って有意と判定されるαエラーが どうなるか確認 サンプルサイズ 群A 25 群B 10, 25, 50 母平均 0 0 母標準偏差 1 0.1, 0.5, 1, 2, 5, 10 有意水準 0.05 シミュレーション回数 10000 Copyright©EPS All rights reserved. 7

2標本t検定 Student or Welch ? Group B = 10 0.30 Group B = 25 0.15 Alpha Error 0.20 0.25 Student Welch 0.1 0.5 1 2 5 10 0.00 0.00 0.05 0.10 0.20 0.15 0.10 0.05 Alpha Error 0.25 0.30 Student Welch 0.1 Population SD in Group B 0.5 1 2 5 10 Population SD in Group B 0.30 Group B = 50 Student Welch 0.20 0.15 0.10 0.00 0.05 Alpha Error 0.25 Studentのt検定 𝑥1 − 𝑥2 𝑡= 1 1 𝑠 + 𝑛1 𝑛2 0.1 0.5 1 2 Population SD in Group B Copyright©EPS All rights reserved. 5 10 Welchのt検定 𝑥1 − 𝑥2 𝑡= 𝑠1 2 𝑠2 2 + 𝑛1 𝑛2 𝑥𝑘 : 群𝑘の標本平均 𝑛𝑘 : 群𝑘の標本数 𝑠 :プールされた標本標準偏差 𝑠𝑘 : 群𝑘の標本標準偏差 8

ノンパラメトリック・ブートストラップ法 -パーセンタイル法母集団分布の明確な仮定が置けないけど，何とか母数の推論をしたい… 平均値の95%信頼区間を求める． ① サンプルサイズnの標本実現値から，n個の実現値を復元抽出し，平均値を計算． ② ①をB回（ブートブートストラップ標本数）繰り返す． ③ 得られた平均値の分布から95%信頼区間を算出する． 昇順に並べ，2.5パーセンタイルと97.5パーセンタイルの値を95%信頼区間の下限と上限とす る． ブートストラップ標本数：B = 4 母集団分布 ？ サンプルサイズ5の 標本実現値 1 2 4 3 2 3 5 3 5 5 3 4 2 Copyright©EPS All rights reserved. 2 1 1 3 3 1 3 4 1 3 5 2 9

ノンパラメトリック・ブートストラップ法 -パーセンタイル法data <- c(55, 24, 97, 77, 63, 86, 5, 42, 30, 14)
B <- 10000
boot_data <- rep(0, each = B)
set.seed(123)
for (i in 1:B) {
boot_data[i] <- mean(sample(data, size=length(data), replace=TRUE))
}

mean(boot_data)
quantile(boot_data, prob=c(0.025, 0.975))

replace=TRUE：復元抽出

信頼区間用のパーセンタイルを取得

hist(boot_data, breaks = 10, xlim = c(0, 100), ylim = c(0, 2500), col = "lightblue")
abline(v=quantile(boot_data, prob=0.025), col = "red", lty = 2, lwd = 3)
abline(v=quantile(boot_data, prob=0.975), col = "red", lty = 2, lwd = 3)

Copyright©EPS All rights reserved.

10

ノンパラメトリック・ブートストラップ法 -パーセンタイル法mean(boot_data) [1] 49.15858 quantile(boot_data, prob=c(0.025, 0.975)) 2.5% 97.5% 30.9975 67.7000 1500 1000 0 500 Frequency 2000 2500 Histogram of boot_data 0 20 40 60 80 100 boot_data Copyright©EPS All rights reserved. 11

参考文献 1. 小杉考司，紀ノ定保礼，清水裕士（2023）．数値シミュレーションで読み解く統計のしくみ -Rでためし てわかる心理統計-．技術評論社． 2. 日本統計学会（編）（2020）．統計学実践ワークブック．学術図書出版社． 3. 東京大学教養学部統計学教室（編）（1991）．統計学入門（基礎統計学 I）．東京：東京大学出版会． 4. ブートストラップ法で平均値の95%信頼区間を求めよう！統計ER．https://best-biostatistics.com/toukeier/entry/bootstrap-method-calculate-95-percent-confidence-interval-mean/（最終参照日：2026年3月30 日） Copyright©EPS All rights reserved. 12

• https://best-biostatistics.com/toukei-

経営理念 Copyright©EPS All rights reserved. 13