はじめに
大学の演習で R による統計処理とベイズ統計モデリングを学んだので、備忘録も兼ねてまとめてみました。「p < 0.05 だから有意差あり!」から「95%の確率で差がある」と言えるようになるまでの道のりです。
正直、最初は「Excel でいいじゃん...」「なんで検定じゃダメなの?」と思っていましたが、実際に使ってみると確率で結果を語れるのがこんなに直感的で便利だとは思いませんでした。
💭 心の声:統計って数学というより「道具」なんだなって実感しました
なぜ R を選ぶのか?
R 言語の特徴
メリット:
- 無料で使える:GNU GPL ライセンスで完全フリー
- 統計解析に特化:統計処理が標準で豊富に用意されている
- 強力な可視化機能:グラフ作成がとにかく強い
- 豊富なライブラリ:CRAN(The Comprehensive R Archive Network)に膨大なパッケージ
- 再現性:スクリプトが残るので分析手順が明確
デメリット:
- 学習コストが高い:プログラミング未経験だと最初はしんどい
- メモリ使用量が多い:大きなデータセットでは注意が必要
Excel との比較
機能 | Excel | R |
---|---|---|
簡単な集計 | ⭕ 直感的 | △ コード必要 |
大量データ処理 | ❌ 重い | ⭕ 高速 |
高度な統計解析 | ❌ 限定的 | ⭕ 豊富 |
グラフ品質 | △ 普通 | ⭕ 論文レベル |
再現性 | ❌ 手順不明 | ⭕ スクリプト残る |
💭 心の声:最初は「どっちか一つ覚えればいい」と思ってましたが、適材適所で使い分けるのが正解ですね
R の基本操作をマスターする
電卓として使う
# 基本的な四則演算
1 + 1 # 足し算
3 * 3 # かけ算
5 %% 5 # 剰余(注意:C言語と記号が違う!)
6 ^ 6 # べき乗
# 数学関数
sqrt(4) # 平方根
log(5) # 自然対数
sin(pi/2) # 三角関数(πも定義済み!)
重要:添え字の違いに注意!
# Rは1スタート(これ重要!)
a <- c(1, 2, 3, 4, 5)
a[1] # 最初の要素(Pythonなら a[0])
a[1:3] # 1番目から3番目まで(3番目も含む!)
💭 心の声:Python 経験者は絶対にここで躓きます。終端の番号も範囲に含まれるのも慣れが必要でした
ベクトル操作の威力
# ベクトル作成と操作
x <- 1:10 # 1から10まで
y <- seq(0, 1, 0.1) # 0から1まで0.1刻み
# ベクトル同士の計算(これが強力!)
x * 2 # 全要素に2をかける
x + y # 要素ごとの足し算
x^2 # 全要素の二乗
これができるようになったこと:
- 一気に大量の計算ができる
- Excel のように一つ一つセルに数式を入れる必要がない
データの可視化 〜まずは見て理解する〜
基本的なプロット
# irisデータセットを使用
data(iris)
# ヒストグラム
hist(iris$Sepal.Length, breaks=15, main="萼片の長さの分布")
# 箱ひげ図
boxplot(Sepal.Length ~ Species, data=iris,
main="種類別萼片長さの比較")
# 散布図
plot(iris$Sepal.Length, iris$Petal.Length,
xlab="萼片の長さ", ylab="花弁の長さ")
これができるようになったこと:
- データの分布を一目で把握できる
- 外れ値や分布の偏りをすぐに発見できる
- 複数グループの比較が視覚的に可能
💭 心の声:グラフの美しさは本当に Excel とは次元が違います。学会発表で使うグラフは絶対 R で作りたくなります!
統計的検定の基礎
t 検定による群間比較
# 2群間の平均値比較
result <- t.test(iris[iris$Species=="setosa", "Sepal.Length"],
iris[iris$Species=="versicolor", "Sepal.Length"])
print(result)
ここで重要な発見: R のデフォルトは Welch's t 検定!これは等分散を仮定しないバージョンです。
メリット・デメリット:
- メリット: 実装が簡単、理解しやすい
- デメリット: 多重比較問題が発生する
よくある勘違い:
- 「p < 0.05 だから有意差あり!」で終わること → 効果量も重要
- 事前に正規性検定が必要と思うこと → 実は不要(むしろ多重検定問題を増やす)
多重検定問題への対処
複数の検定を同時に行うと、偽陽性(本当は差がないのに「差がある」と判定)の確率が上がってしまいます。
# Bonferroni補正
pairwise.t.test(iris$Sepal.Length, iris$Species,
p.adjust.method="bonferroni")
# Holm法(より検出力が高い)
pairwise.t.test(iris$Sepal.Length, iris$Species,
p.adjust.method="holm")
💭 心の声:最初は「補正しすぎて何も有意にならない...」と思ったけど、それが正しい姿勢なんですね
Games-Howell 検定 〜最強の多重比較〜
個人的に一番気に入った方法です。Welch's t 検定ベースで、サンプルサイズが違っても、分散が違っても OK!
# PMCMRplusパッケージが必要
install.packages("PMCMRplus")
library("PMCMRplus")
# Games-Howell検定
gamesHowellTest(g=iris$Species, x=iris$Sepal.Length)
これができるようになったこと:
- 分散が異なるグループでも安心して多重比較できる
- Bonferroni 法より検出力が高い場合がある
- 事前の ANOVA が不要
💭 心の声:「なんでみんなこれ使わないの?」と思うくらい便利です
従来の検定の限界と問題点
p 値の恣意性
# 極小な差でもサンプルサイズを大きくすれば「有意」になる
set.seed(1)
large_n <- t.test(rnorm(10000, mean=0.01, sd=0.1),
rnorm(10000, mean=0.00, sd=0.1))
small_n <- t.test(rnorm(100, mean=0.01, sd=0.1),
rnorm(100, mean=0.00, sd=0.1))
print(large_n$p.value) # 有意になる
print(small_n$p.value) # 有意にならない
💭 心の声:N=10000 の場合、平均値の差が 0.01 でも有意差が出るけど、この 0.01 に実際的な意味はあるのか...?
従来検定の問題点:
- p 値の恣意性: サンプルサイズで結果が変わる
- 多重比較問題: 検定を増やすと偽陽性率が上昇
- 解釈の困難さ: 「帰無仮説が棄却される確率」は分からない
ベイズ統計モデリングという新しいアプローチ
ベイズ統計の基本的な考え方
ベイズ統計モデリングでは、データが特定の分布から生成されていると仮定し、そのパラメータの事後確率分布を求めます。
事後確率 ∝ 尤度 × 事前確率
例: 母平均 μ、母標準偏差 σ の正規分布からデータが生成されていると仮定
- 観測データ x = (x₁, x₂, x₃, ...) を取得
- 事後確率 p(μ, σ | x) の確率密度関数を推定
- 「μ > 0 の確率」などを計算して活用
💭 心の声:最初は「事後確率って何?」状態でしたが、要は「データを見た後で、パラメータがこの値である確率」のことですね
MCMC(Markov Chain Monte Carlo)による近似
真面目に事後確率を計算するのは困難なので、MCMC 法で近似します。
MCMC とは: どんな分布にでも高速で擬態できる乱数生成器(乱暴な説明)
メリット:
- 複雑な分布でも近似可能
- 計算結果の精度が確保されている
- 様々な統計ソフトで実装済み
デメリット:
- 収束判定が必要
- 計算時間がかかる
- 理論的理解が必要
Stan と RStan の環境構築
Stan は確率的プログラミング言語で、RStan はその R ラッパーです。
# RStanのインストール
install.packages("rstan")
library(rstan)
# 診断用パッケージ
install.packages(c("ggmcmc", "ggplot2"))
library(ggmcmc)
実践例 1: 単一グループの平均推定
データ準備
# 推定対象のデータ作成
set.seed(1)
N_for_norm <- 15
mu <- 5
sigma <- 4
norm_dist <- rnorm(N_for_norm, mean=mu, sd=sigma)
Stan コード(normal.stan)
data{
int N;
vector[N] data_1;
}
parameters{
real mu_1;
real<lower=0> sigma_1;
}
model{
data_1 ~ normal(mu_1, sigma_1);
}
R 側の実行コード
# Stanの引数構築
stan.in <- list(N=length(norm_dist), data_1=norm_dist)
# Stan実行
stan_res <- stan(file="normal.stan", data=stan.in, seed=1)
# 結果確認
print(stan_res)
結果の解釈
# MCMCサンプルの取得
mu_samples <- rstan::extract(stan_res)$mu_1
# 確率計算
f <- function(N){
sum(mu_samples > N)/length(mu_samples)
}
# 3以上の確率
f(3) # 0.9775
# 95%確信区間
quantile(mu_samples, probs=c(0.025, 0.975))
これができるようになったこと:
- 「μ が 3 以上である確率は 97.75%」のような直感的な解釈
- パラメータの不確実性を確率分布として表現
- 信頼区間ではなく確信区間(Credible Interval)が得られる
💭 心の声:「95%以上の確率で差がある」って言えるのが、p 値よりもずっと分かりやすい!
実践例 2: 2 群間の平均値の差
Stan コード(normal.mean.diff.stan)
data{
int N;
int M;
vector[N] data_1;
vector[M] data_2;
}
parameters{
real mu_1;
real<lower=0> sigma_1;
real mu_2;
real<lower=0> sigma_2;
}
model{
data_1 ~ normal(mu_1, sigma_1);
data_2 ~ normal(mu_2, sigma_2);
}
generated quantities{
real diff;
diff = mu_2 - mu_1;
}
実行例
# データ準備
set.seed(1)
mu <- c(0, 5)
sigma <- c(3, 4)
normA <- rnorm(15, mean=mu[1], sd=sigma[1])
normB <- rnorm(15, mean=mu[2], sd=sigma[2])
# Stan実行
stan.in <- list(N=length(normA), M=length(normB),
data_1=normA, data_2=normB)
stan_res <- stan(file="normal.mean.diff.stan", data=stan.in, seed=1)
結果の解釈
# 差の分布を取得
diff_samples <- rstan::extract(stan_res)$diff
# 差が正である確率
sum(diff_samples > 0) / length(diff_samples)
# 95%確信区間
quantile(diff_samples, probs=c(0.025, 0.975))
収束診断の重要性
MCMC は数値計算なので、必ず収束判定が必要です。
# 診断プロット作成
ggmcmc(ggs(stan_res, inc_warmup=TRUE), file="convergence_check.pdf")
チェックポイント:
- 全ての分布が十分に重なっているか?
- warmup 期間中に収束しているか?
- 変なジャンプが起きていないか?
- Rhat < 1.1 になっているか?
よくある勘違い:
- 「とりあえず動いたから大丈夫」→ 必ず収束診断を行う
- 「Rhat が 1.1 を少し超えてるけど大丈夫でしょ」→ 再実行するか設定を見直す
従来手法との比較
手法 | メリット | デメリット | 使用場面 |
---|---|---|---|
t 検定 | 実装簡単、計算高速 | p 値の解釈が困難、多重比較問題 | 簡単な 2 群比較 |
Games-Howell 検定 | 多重比較対応、分散の仮定不要 | 依然として p 値ベース | 多群比較 |
ベイズモデリング | 確率的解釈、柔軟なモデル | 計算コスト、理論的理解必要 | 複雑な推論、確率的判断 |
ベイズ統計のメリット・デメリット
メリット
これができるようになったこと:
- 「差がある確率は 95%」のように直感的な解釈
- 複数比較でも多重比較問題が発生しない
- パラメータの確率分布が得られるため、より豊富な情報が取得可能
- 事前情報を組み込める柔軟性
デメリット
- 計算コストが高い:MCMC サンプリングが必要
- モデル設定が複雑:Stan コードの記述が必要
- 事前分布の設定:主観が入る可能性
- 理論的理解が必要:確率論の基礎知識が必要
よくある勘違い
-
「ベイズは主観的だから科学的でない」
→ 事前分布を無情報的にすれば客観的な分析も可能 -
「計算が複雑すぎて実用的でない」
→ Stan などのツールで実装は簡単になった -
「従来の統計学を完全に置き換える」
→ 目的に応じて使い分けることが重要
まとめ
個人的ベストプラクティス
- まず可視化: データの分布を必ず確認
-
適切な手法選択:
- 簡単な比較 → Games-Howell 検定
- 確率的解釈が必要 → ベイズモデリング
- 収束診断は必須: Rhat とトレースプロットのチェック
- 効果量も重視: 統計的差異と実用的差異は別
学習して良かったこと
- 直感的な解釈: 確率で結果を語れる
- 再現性のあるデータ分析: スクリプトが残る
- 柔軟性: 複雑なモデルに対応可能
- 多重比較問題からの解放: 複数の比較でも問題なし
今後の課題
- より複雑な階層モデルの理解
- 事前分布の適切な設定
- 計算効率の改善
最後に: 統計は完璧を求めるより、まず手を動かすことが大切だと感じました。ベイズ統計は最初のハードルは高いですが、一度慣れると検定よりもずっと直感的で便利です。特に「確率で語れる」というのが最大の魅力だと思います。
💭 心の声:「Excel 卒業」というより「適材適所」が正解。統計は道具なので、目的に応じて使い分けることが重要ですね!
参考文献・情報源
- 実践 Data Science シリーズ R と Stan ではじめる ベイズ統計モデリングによるデータ分析入門 (ISBN-10: 4065165369)
- Stan と R でベイズ統計モデリング (ISBN-10: 4320112423)
- The R Tips: データ解析環境 R の基本技・グラフィックス活用集 第 3 版 (ISBN-10: 4274219585)
- Stan 公式サイト
- R 公式サイト
- 奥村先生のウェブページ