【検定力分析】R で仮説検定を行う

今回は, 仮説検定の中でも基本的な t検定とカイ二乗検定 , 加えて検定力分析について書きました。

仮説検定

一般的に仮説検定 (hypothesis testing) では, 主張したい仮説である対立仮説 H1 の逆の仮説を帰無仮説 H0 として, 標本から H0 を棄却できるかを検証するアプローチを取ります。仮説検定には様々な種類があり, 検定の目的と母集団の仮定により使い分ける必要があります。ここで挙げている以外にも多くの検定手法があります。

2群に関する検定手法	正規性	等分散性	特徴
Student’s t	要	要	平均の差
Welch’s t	要	不要	平均の差
Mann-Whitney	不要	要	平均の差
Wilcoxon rank sum	不要	不要	中央値の差
Brunner-Munzel	不要	不要	中央値の差 [1]

一般的に特定の確率分布を仮定する parametric な検定の方が non-parametric な検定より検出力は高いとされています。
正規性の検定には, Shapiro-Wilkの正規性検定また, 等分散の検定は F検定が代表的です。

仮説検定におけるエラー

仮説検定の議論では, 第一種の過誤（Type I error）と第二種の過誤 (Type II error）の話が付いて回ります。

第一種の過誤（Type I error, α-error） : 帰無仮説が正しいにも関わらず, 帰無仮説を棄却してしまう
第二種の過誤 (Type II error, β-error） : 帰無仮説が正しくない (対立仮説が正しい) にも関わらず, 帰無仮説を棄却できない

仮説検定における有意水準とは, α-error を起こす確率 P(α) のことです。
この有意水準 P(α) の設定で大事なのが, 第一種の過誤を抑えたいがために有意水準を 1%, 0.1%… と厳しくすると, 第二種の過誤が増えてしまうことです。つまり, 基本的には第一種の過誤と第二種の過誤はトレードオフの関係にあります。

t検定

t検定は基本的には母集団を正規分布と仮定したパラメトリック検定で, 2郡の平均値に差があるかどうかの検定を行います。

以下の x, y の2群に平均差があるかの検定を行ってみます。

x <- c(1,2,1,1,1,1,1,1,1,1,2,4,1,1)
y <- c(3,3,4,3,1,2,3,1,1,5,4,6)

Rの t.test()は Welch’s t test で等分散かわからない場合にも使えます。

t.test(x, y, paired = F)

結果は下記のようになりました。
t統計量が -3.2, p値が 0.0054 となったので2群の平均値に差がないという帰無仮説 H0 は有意水準 1% で棄却され対立仮説 H1 を採用し, 2群の平均値に差がありそうだと言えます。

	Welch Two Sample t-test

data:  x and y
t = -3.205, df = 16.101, p-value = 0.005484
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -2.7289451 -0.5567692
sample estimates:
mean of x mean of y 
 1.357143  3.000000

Mann–Whitney U testは, 正規性を仮定しないノンパラメトリック検定ですが, 等分散でない場合は問題があるようです。ノンパラメトリック検定で等分散性も仮定しない Brunner-Munzel test は {lawstat} パッケージに含まれており, 中央値の検定に使うと良さそうです。

対応のあるt検定

対応のあるt検定 (paired t-test) はデータがペアとして対応している2群における各観測値の差に関する検定です。t検定の枠組みなので正規性を仮定します。
正規性の仮定が怪しい場合には, 事前にシャピロ・ウィルク検定 (Shapiro-Wilk test) で確認します。
帰無仮説 H0 は2群のデータに差がないこととします。従って, ペアのデータの差に着目するので, データ間の分散は気にする必要がありません。

例で扱うデータセット Student’s Sleep Data は, 2つの睡眠薬の効果の比較データです。

シャピロ・ウィルク検定で正規性の検定を行います。帰無仮説は標本分布が正規分布に従うことなので, 帰無仮説が棄却されないことで正規分布に従うと考えることができます。

shapiro.test(grp1$extra)
# Shapiro-Wilk normality test
#
# data:  grp1$extra
# W = 0.92581, p-value = 0.4079

shapiro.test(grp2$extra)
# Shapiro-Wilk normality test
#
# data:  grp2$extra
# W = 0.9193, p-value = 0.3511

t.test で paired を TRUE にします。

t.test(grp1$extra, grp2$extra, paired = T)

p値が 0.0028 となり, 2つの睡眠薬の効果には差がある可能性が高いと言えそうです。

Paired t-test

data:  grp1$extra and grp2$extra
t = -4.0621, df = 9, p-value = 0.002833
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-2.4598858 -0.7001142
sample estimates:
mean of the differences
              -1.58

補足になりますが, t.test() で両側検定, 片側検定は alternative = c(“two.sided”, “less”, “greater”) で設定できます。

カイ二乗検定

カイ二乗検定 (Chi-squared test) は適合度検定, 独立性検定で独立な標本間の比率の差を検定します。
帰無仮説 H0 が正しい場合に検定統計量がカイ二乗分布に従います。頻度分布自体は特定のものに限定されないので, ノンパラメトリック検定に分類されます。

Rで Pearson’s Chi-squared testを行ってみます。

例として以下の架空のデータを使います。
交通手段として自転車と徒歩を男女別に集計したデータで, 交通手段に性別差があるか否かの検定を行います。

         transportation
gender   bicycle walk
male         2   15
female      10    3

chisq.test()がPearson’s Chi-squared testになります。
X-squared統計量が 10.4, p値が 0.001 になりました。今回(架空)のデータだと交通手段に性別差があると言えそうです。

d <- matrix(
  c(2, 10, 15, 3),
  nrow = 2,
  dimnames = list("gender" = c("male", "female"),
                  "transportation" = c("bicycle", "walk"))
)

# Chi-squared test
chisq.test(d)
# X-squared = 10.4581, df = 1, p-value = 0.001221

サンプルサイズが小さく, 分割表の期待度数に 10 以下の値がある場合や数値の偏りが大きい場合には正確確率検定 (Fisher’s Exact Test) の方が良い場合もありそうです。サンプルサイズがよほど大きくない限り, 常に Fisher’s Exact Test を使った方が良いという声もあります。

# Fisher exact test
fisher.test(d)
# p-value = 0.0005367
# alternative hypothesis: true odds ratio is not equal to 1
# 95 percent confidence interval:
#  0.003325764 0.363182271
# sample estimates:
# odds ratio
# 0.04693661

検定力分析

検定力分析を考える上でもう一度, 仮説検定におけるエラーに触れておきます。

第一種の過誤 (Type I error, α-error): 帰無仮説が正しい場合に, 帰無仮説を誤って棄却してしまう事
第二種の過誤 (Type II error, β-error): 帰無仮説が正しくない場合に, 帰無仮説を誤って採択してしまう事

β は P(α) を 0.01 や 0.001 と厳しくする程, 起きやすくなる場合があります。
検定力 (Power) または検出力は, 対立仮説 H1 が正しい場合に帰無仮説 H0 が棄却される確率です。検定力は 1 – P(β) となります。
例えば, 検定力が 0.8 の時は, 80% の確率で対立仮説が正しい場合に帰無仮説が棄却できることを意味します。

サンプルサイズが小さすぎる場合, 検定力が下がってしまう場合があり, 逆にサンプルサイズが大きすぎる場合には α が発生する可能性が高くなります。

臨床研究では, 多くの場合 P(α)=0.05, Power=0.8 または 0.9 に設定しているようです。

stats::power.t.test() で検定力分析を行ってみます。下記の中で, 求めたいものに NULL を指定します。

検定力
有意水準
サンプルサイズ

例えば, サンプルサイズを決めたい時は, n = NULL を指定します。

# calc power
power.t.test(
  n = 10, # sample size
  delta = 0.5, # true difference in means
  sd = 1, # standard deviation
  sig.level = 0.05, # significance level
  power = NULL, # power of test
  strict = TRUE # two-side case
)
# Two-sample t test power calculation 
#
#          n = 10
#      delta = 0.5
#         sd = 1
#  sig.level = 0.05
#      power = 0.1850957
# alternative = two.sided

# calc sample size
power.t.test(
  n = NULL,
  delta = 0.2,
  sd = 1,
  sig.level = 0.05,
  type = "paired",
  power = 0.8,
  alternative = "one.sided"
)
# Paired t test power calculation
#
#          n = 155.9257
#      delta = 0.2
#         sd = 1
#  sig.level = 0.05
#      power = 0.8
# alternative = one.sided

検定力分析は {stats} 以外にも {pwr} を使う方法もあり, {pwr} の場合, 有意な差があった場合のどれだけ実験的操作の実質的に意味があるかの指標, 変数間の関係の強さを表す指標である効果量 (Effect size) を指定できます。標準化しているため, 観測データの単位の影響を受けない特徴があります。

require(pwr)

pwr.t.test(
    n = NULL,
    d = 0.5, # Effect size
    sig.level = 0.05,
    power = 0.8
)

A Power Primer [7] によると, t検定の効果量 d の目安として小さな差が0.2, 中程度の差が0.5, 大きな差が0.8 となっています。対応ありの時の効果量 r は, 小さな差が0.2, 中程度の差が0.3, 大きな差が0.5 を目安としている場合もあります。 (水本篤, 竹内理 2008)

おわりに

Code は GitHub にあります。

仮説検定の概念の基本的な理解は, 『実証分析入門データから「因果関係」を読み解く作法』を参考にしました。

また, 検定力の部分は『新版入門医療統計学 -Evidenceを見いだすために』を参考にしました。

[1] マイナーだけど最強の統計的検定 Brunner-Munzel 検定
[2] カイ二乗分布
[3] Fisherの正確検定かカイ2乗検定か
[4] １３章対応のある２群間の量的データの検定
[5] 対応のあるt検定
[6] シャピロ・ウィルク検定
[7] A Power Primer [Cohen 1992]
[8] Rによるやさしい統計学第20章「検定力分析によるサンプルサイズの決定」
[9] 第５５回　検定力分析
[10] 統計的検定における検定力と効果量