ランダム化比較試験の主要解析で全体の治療効果が示された後、バイオマーカーの高値群と低値群に分けて同じ検定を繰り返すと、片方だけがp値0.05を下回る、という結果はよく現れます。ここから効果は高値群に限られると読むと、人数の差が作った検出力の違いを、効果の違いとして受け取ってしまうことになります。
群ごとの検定を並べる代わりに、治療効果そのものが部分集団の間で異なるかどうかを一つの推定量として評価するのが交互作用の考え方です。この記事では、乱数で生成したデータセットと反復シミュレーションを使い、二つの読み方でどこまで結論が変わるかを扱います。
効果の異質性を表す記号と推定量
記号を先に決めておきます。アウトカム$Y$は12週間でのLDLコレステロールの変化量、$T$は治療群で1、対照群で0を取る指示変数、$S$はバイオマーカー高値群で1、低値群で0を取る指示変数とします。
この記号を使うと、部分集団ごとの治療効果は低値群で$\delta_{\text{low}}$、高値群で$\delta_{\text{high}}$と書けます。異質性があるかどうかを決めるのは個々の値ではなく、二つの差$\theta = \delta_{\text{high}} – \delta_{\text{low}}$のほうです。
$$\mathbb{E}[Y \mid T, S] = \beta_0 + \beta_1 T + \beta_2 S + \beta_3 TS \tag{1}$$
式(1)は、治療の有無と部分集団の別で平均を四つに分け、それを切片、治療の主効果、部分集団の主効果、交互作用の係数で表したものです。$\beta_0$は低値群の対照群での平均、$\beta_1$は低値群での治療効果、$\beta_2$は対照群における部分集団の差、$\beta_3$は$\theta$に対応する交互作用で、高値群での治療効果は$\beta_1 + \beta_3$になります。
ここで注意が必要なのは、$\beta_3$が部分集団ごとのp値の差ではなく、効果の差そのものを推定している点です。異質性を論じる対象は、この一つのパラメータに最初から定まります。
群内のp値が食い違うときに違うもの
乱数で生成したデータセットで確かめます。個人差の標準偏差を30mg/dLとし、バイオマーカー高値群900例と低値群150例の合わせて1,050例について、各部分集団の中で対照群と治療群へ1対1に割り付け、12週間でのLDL変化量を生成しました。
この設定で置いた真の治療効果は、両方の部分集団とも−8mg/dLで同じ値です。真の異質性はゼロですから、群ごとの結果が食い違えば、それは偶然と情報量の差が作った見かけの食い違いになります。
実際に当てはめると、高値群での治療効果の推定値は−8.594mg/dL、95%信頼区間は−12.543から−4.646で、p値は0.0000215でした。低値群では推定値が−11.152mg/dL、95%信頼区間が−20.906から−1.399、p値は0.0253です。
この二つを並べると、p値は3桁異なりますが、点推定値はむしろ低値群のほうが大きい低下を示しています。信頼区間の幅も低値群で約19.5mg/dL、高値群で約7.9mg/dLと2倍以上の開きがあり、p値の差はこの幅の差から生じています。
差の差としての異質性の推定
同じデータへ式(1)を当てはめると、交互作用項の推定値は2.558mg/dL、標準誤差は5.324、95%信頼区間は−7.890から13.006、p値は0.631でした。交互作用を落としたモデルとの比較によるF検定でも、残差自由度1,046のもとでF値0.231、p値0.631と同じ結果になります。
四つのセル平均を全て自由に推定する飽和したモデルでは交互作用項と群内の差の差が厳密に一致するため、この2.558は高値群の推定値−8.594から低値群の推定値−11.152を引いた差そのものです。
なお、セル平均から組み立てても同じ値に届き、低値群の対照群−2.286と治療群−13.438、高値群の対照群−5.698と治療群−14.293から群内の差を取ると−11.152と−8.594になり、その差が2.558になります。
そのため、p値が片方で小さく片方で大きいという観察から異質性を語ると、区間が−7.890から13.006まで広がる量について、符号すら定まらない段階で結論を出すことになります。効果の差を知りたいのであれば、報告する数値は群ごとのp値ではなく、この差の推定値と区間のほうです。
真の異質性がない場合の食い違い
一度のデータで起きたことが偶然かどうかは、同じ設定で試験を繰り返せば見えてきます。真の交互作用をゼロに固定したまま10,000回生成したところ、高値群だけが有意になった試験が61.15%、低値群だけが有意になった試験が0.65%を占めました。
合わせて61.8%の試験で群内のp値が食い違った一方、同じ10,000回で交互作用検定が有意になったのは5.22%でした。モンテカルロ標準誤差が0.22%ですから、この値は名目水準の5%に収まっています。
つまり、真の異質性が全くない設定でも、群ごとの検定を並べる読み方では6割の試験で片方だけ効いたという記述が生まれます。同じデータを交互作用として評価すれば、誤って異質性を主張する頻度は設定した水準のまま保たれます。
同じ−8mg/dLの効果でも900例の高値群は150例の低値群よりはるかに有意になりやすく、ここで効いているのは検出力の差です。食い違いの内訳も、高値群だけが有意になる側へ大きく偏ります。
交互作用検定の検出力が低くなる理由
交互作用の推定量は二つの推定値の差ですから、その標準誤差は両方の不確実性を合わせたものになります。部分集団が重ならない場合は、次の形で書けます。
$$\operatorname{SE}(\hat\beta_3) = \sqrt{\operatorname{SE}(\hat\delta_{\text{high}})^2 + \operatorname{SE}(\hat\delta_{\text{low}})^2} \tag{2}$$
式(2)は、交互作用の標準誤差が、部分集団ごとの標準誤差の二乗和の平方根になることを表します。今回の設定では高値群で2.000、低値群で4.899となり、交互作用では5.292と、大きいほうの部分集団の2.646倍に広がります。
この広がりは検出力に直結します。真の効果を高値群で−12mg/dL、低値群で−4mg/dLとし、真の交互作用を−8mg/dLと置いて5,000回生成すると、交互作用検定が有意になったのは33.1%にとどまり、正規近似から計算した理論値32.73%と一致しました。
一方、同じ設定で高値群だけが有意になった試験は86.84%に達しています。本当に異質性がある場面でも、群ごとの検定を並べる方式では高値群で効いたという記述だけが量産され、差の大きさそのものは報告に残りません。
そのため、真に8mg/dLの差があっても3回に2回は拾えない33.1%という検出力を踏まえると、交互作用検定が有意でなかったことを根拠に効果が均一だと述べるのは適切ではありません。非有意という結果は、異質性の否定ではなく情報の不足として読むほうが実態に合います。
こうした関係について、Brookes STら（2004）は、シミュレーションによって交互作用の検出に必要な症例数が主効果の検出に必要な症例数よりはるかに大きくなることを示しています。今回得られた2.646倍という標準誤差の比は、同じ事情を別の角度から見たものにあたります。
部分集団因子の増加と偽陽性
実際の報告で並ぶ部分集団は一つではなく、年齢、性別、重症度、併用薬、地域のように複数の因子が同時に検討されます。真の交互作用がどれもゼロの因子を5つ調べる設定で5,000回生成したところ、少なくとも一つが0.05で有意になった試験は23.7%を占め、独立を仮定した解析解22.62%とほぼ一致しました。
この2割強という水準は、多重性を調整しないまま探索した場合に見込まれる偽陽性の大きさを示しています。Yusuf Sら（1991）が早い時期に指摘したとおり、多数の部分集団を調べれば、そのうちどれかで有意な結果が出ることは偶然だけでも十分に起こります。
ただし、調整方法を決めれば済む話ではなく、事前に規定した部分集団の数と検定の数を報告に書くことが先になります。どれが計画された解析で、どれが結果を見てから追加されたのかを読者が区別できる状態になって、はじめてp値の水準を議論できます。
検証的解析と探索的解析での役割の違い
同じ交互作用検定のp値でも、試験実施計画書に事前規定された解析から出たものと、結果を見た後で部分集団を選んで出したものとでは、支えられる主張が異なります。ICH E9は5.7節で部分集団、交互作用、共変量の扱いを事前に規定するよう求めており、事前規定のない部分集団解析には探索的なものという位置づけを与えています。
この区別について、Wang Rら（2007）は、部分集団解析の報告で事前規定か事後の探索かを明記し、群ごとのp値ではなく交互作用の検定結果と効果推定値を示すよう求めています。検定の数と多重性への対応を述べることも、同じ報告要件に含まれています。
もっとも、探索的な結果は次の研究の仮説を作る材料として十分に役立ちます。ただし、その場で治療方針を変える根拠として使うと、先ほどの23.7%という偽陽性の水準をそのまま持ち込むことになります。
報告と査読で求められる記載
公表論文での扱いを見ると、Rothwell PM（2005）は、部分集団ごとに別々の検定結果を並べる報告が広く行われている一方で、その解釈には事前規定、生物学的な妥当性、効果の方向の一貫性といった条件が要ることを整理しています。部分集団解析を載せること自体ではなく、どの条件のもとでどこまで読めるかが論点になります。
この整理を踏まえ、交互作用が有意だった場合の信頼性については、Sun Xら（2010）が評価基準をまとめています。事前に立てた仮説かどうか、効果の方向まで予測していたか、検定の数、交互作用検定を実際に行ったか、他の研究で再現されているか、生物学的な裏づけがあるかを順に確認する形です。
こうした基準を実務へ落とすと、部分集団ごとの推定値と信頼区間を図で示したうえで、交互作用の推定値、区間、p値を本文に置き、事前規定の有無と検定数を方法の節へ書く形になります。今回のデータで言えば、本文へ置く数値は2.558という差と−7.890から13.006という区間になります。
Rでの交互作用の確認
ここまでの数値は、base Rとstatsだけで再現できます。低値群150例と高値群900例からなる1,050例について、個人差の標準偏差30mg/dLのもとで12週間のLDL変化量を用意しました。
set.seed(42)
subgroup <- rep(c("biomarker_low", "biomarker_high"), c(150, 900))
treatment <- c(rep(c("control", "treatment"), each = 75),
rep(c("control", "treatment"), each = 450))
mean_change <- -3 - 2 * (subgroup == "biomarker_high") - 8 * (treatment == "treatment")
trial <- data.frame(
subgroup = factor(subgroup, levels = c("biomarker_low", "biomarker_high")),
treatment = factor(treatment, levels = c("control", "treatment")),
ldl_change = rnorm(1050, mean = mean_change, sd = 30))
用意したデータに対して、部分集団ごとの当てはめと、交互作用を含む当てはめを並べます。推定値、95%信頼区間、p値を同じ形式で表示すると、三つの推定量の関係が読み取りやすくなります。
fit_high <- lm(ldl_change ~ treatment, data = trial[trial$subgroup == "biomarker_high", ])
fit_low <- lm(ldl_change ~ treatment, data = trial[trial$subgroup == "biomarker_low", ])
model_interaction <- lm(ldl_change ~ treatment * subgroup, data = trial)
show_effect <- function(label, fit, row) {
cat(sprintf("%s %.3f (%.3f, %.3f) p = %.3g\n", label, coef(fit)[row],
confint(fit)[row, 1], confint(fit)[row, 2],
summary(fit)$coefficients[row, 4]))
}
show_effect("high ", fit_high, 2)
show_effect("low ", fit_low, 2)
show_effect("interaction", model_interaction, 4)
#> high -8.594 (-12.543, -4.646) p = 2.15e-05
#> low -11.152 (-20.906, -1.399) p = 0.0253
#> interaction 2.558 (-7.890, 13.006) p = 0.631
この出力の3行目が異質性の推定値で、区間はゼロをまたいでいます。読み違えやすいのは1行目と2行目の対比で、p値が0.0000215と0.0253で大きく違っても、それは効果の差ではなく人数の差を反映しており、差そのものの評価は3行目でしか行えません。
同じ対比を図にすると、区間の重なりまで一度に見えます。部分集団ごとの区間と交互作用の区間を同じ軸へ描いたものが図1です。
png("fig1-subgroup-forest.png", width = 1100, height = 688, res = 200, type = "cairo")
par(mar = c(4.2, 9.5, 2.2, 1.2))
est <- c(coef(fit_high)[2], coef(fit_low)[2], coef(model_interaction)[4])
ci <- rbind(confint(fit_high)[2, ], confint(fit_low)[2, ], confint(model_interaction)[4, ])
plot(est, 3:1, xlim = c(-24, 16), ylim = c(0.5, 3.6), yaxt = "n", pch = 16, ylab = "",
xlab = "Difference in LDL change, treatment - control (mg/dL)", main = "Same true effect in both subgroups")
axis(2, at = 3:1, las = 1, cex.axis = 0.8, labels = c("Biomarker high (n=900)", "Biomarker low (n=150)", "Interaction (high - low)"))
segments(ci[, 1], 3:1, ci[, 2], 3:1, lwd = 2)
abline(v = 0, col = "grey55", lty = 2)
dev.off()
反復シミュレーションも同じ枠組みで行い、平均構造を固定したままアウトカムだけを10,000回生成し直しました。真の交互作用はゼロのままなので、観察される食い違いは全て見かけのものになります。
p_high <- p_low <- numeric(10000)
for (i in seq_len(10000)) {
sim <- trial
sim$ldl_change <- rnorm(1050, mean = mean_change, sd = 30)
p_high[i] <- summary(lm(ldl_change ~ treatment,
data = sim[sim$subgroup == "biomarker_high", ]))$coefficients[2, 4]
p_low[i] <- summary(lm(ldl_change ~ treatment,
data = sim[sim$subgroup == "biomarker_low", ]))$coefficients[2, 4]
}
mean(p_high < 0.05 & p_low >= 0.05)
#> [1] 0.6115
この0.6115という割合が、同じ反復での交互作用検定の有意割合0.0522と大きく離れます。離れるのは検定の作り方によるもので、群内の検定が各部分集団の情報量に左右されるのに対し、交互作用検定は効果の差そのものを対象にします。
研究で判断するときの順序
同じ真の効果を置いた今回の設定で6割の試験が群内のp値の食い違いを生んだ以上、その食い違いを効果の差の科学的根拠として使うことはできません。部分集団ごとの治療効果を論じる研究で報告すべきなのは、群ごとの検定を並べた結果ではなく、交互作用の推定値と区間です。
そこで、異質性を論じる場面では、交互作用の推定値と信頼区間を主たる報告対象に置き、部分集団ごとの推定値は参考として図に添える形にします。事前規定した部分集団と検定数を方法の節に書いておけば、探索的な結果が検証的な結論として読まれる余地も狭まります。
逆の向きの誤りも同じ構造で起こります。今回の検出力が33.1%であったように、非有意という結果は効果が均一であることの科学的根拠には届かないため、交互作用検定が有意でなかった場合は区間の幅を示して不確実性の大きさを伝えるほうが読者の判断を助けます。
実際の研究では、事前規定の有無、部分集団の数、アウトカムの種類、検出力の見積もりまで含めて、どこまで書けるかが決まります。Dr.データサイエンスでは、こうした部分集団解析の方針についてのご相談を承っています。
参考文献
- Rothwell PM. Subgroup analysis in randomised controlled trials: importance, indications, and interpretation. The Lancet. 2005;365(9454):176-186.
- Wang R, Lagakos SW, Ware JH, Hunter DJ, Drazen JM. Statistics in medicine — reporting of subgroup analyses in clinical trials. The New England Journal of Medicine. 2007;357(21):2189-2194.
- Brookes ST, Whitely E, Egger M, Davey Smith G, Mulheran PA, Peters TJ. Subgroup analyses in randomized trials: risks of subgroup-specific analyses; power and sample size for the interaction test. Journal of Clinical Epidemiology. 2004;57(3):229-236.
- Sun X, Briel M, Walter SD, Guyatt GH. Is a subgroup effect believable? Updating criteria to evaluate the credibility of subgroup analyses. BMJ. 2010;340:c117.
- Yusuf S, Wittes J, Probstfield J, Tyroler HA. Analysis and interpretation of treatment effects in subgroups of patients in randomized clinical trials. JAMA. 1991;266(1):93-98.
- International Council for Harmonisation. ICH Harmonised Tripartite Guideline E9: Statistical Principles for Clinical Trials. Step 4 version, 5 February 1998. Section 5.7 Subgroups, Interactions and Covariates.
