医学研究で回帰分析を使うと、係数、p値、決定係数が一度に出てきます。そのため、何を知るための解析だったのかが曖昧なまま、数値の大小だけを説明してしまうことがあります。たとえば出生体重を妊娠前体重で回帰した時、傾きは両者の関連を表せますが、妊娠前体重を変えれば出生体重が変わるという因果効果まで保証しません。さらに、観察された人を説明する目的と、未知の人を予測する目的とでは、同じ回帰式でも評価すべき内容が異なります。
回帰分析を読む中心は、目的変数の何を、どの説明変数に条件づけて表した式なのかを確かめることです。この記事を読み終えると、記述、関連の推定、予測、因果効果の推定を区別し、係数の単位と調整条件を研究目的に沿って説明できるようになります。
回帰式が表す条件付き平均
出生体重を$Y$、妊娠前体重を$X$とします。単回帰では、妊娠前体重が同じ対象者について期待される出生体重を、次の直線で表します。
$$
\mathbb{E}(Y\mid X=x)=\beta_0+\beta_1x. \tag{1}
$$
式（1）は、妊娠前体重$X$が$x$ポンドの対象者における出生体重$Y$の条件付き平均を、切片$\beta_0$と傾き$\beta_1$で表しています。$\beta_1$は妊娠前体重が1ポンド大きい対象者同士で比べた平均出生体重の差です。
ここで大切なのは、回帰式が個々の出生体重を直線上に固定しているわけではないことです。同じ妊娠前体重でも出生体重にはばらつきがあり、式が記述するのは条件付き平均です。また、切片は$X=0$での平均を表しますが、今回観察された妊娠前体重は80〜250ポンドなので、切片を実質的な医学的結果として読む理由はありません。
群間比較との関係も、この式から分かります。説明変数を喫煙なしで0、喫煙ありで1とした二値変数$S$へ置き換えると、傾きは喫煙あり群となし群の平均差になります。一方、相関係数は変数同士の直線的な結びつきを単位のない1つの値に要約するので、出生体重が何グラム異なるかを直接は表しません。回帰係数は目的変数と説明変数の役割を分け、元の単位で差を記述できる点が異なります。
係数の意味を決める調整条件
説明変数が妊娠前体重$X$と喫煙状況$S$の2つなら、重回帰式は次の形になります。
$$
\mathbb{E}(Y\mid X=x,S=s)=\beta_0+\beta_1x+\beta_2s. \tag{2}
$$
式（2）は、妊娠前体重$x$ポンド、喫煙状況$s$の対象者における出生体重の条件付き平均です。$\beta_1$は喫煙状況が同じ対象者同士で妊娠前体重が1ポンド異なる時の平均差、$\beta_2$は妊娠前体重が同じ対象者同士で喫煙ありとなしを比べた平均差を表します。
「妊娠前体重の係数」とだけ報告すると、この条件が抜け落ちます。MASSパッケージのbirthwtデータセットでは、妊娠前体重だけを使った傾きは1ポンド当たり4.43gでしたが、喫煙状況も含めると4.24gでした。10ポンド当たりに直すと、それぞれ44.29gと42.37gです。Leffondréら（2014）が述べるように、連続変数の係数には変化量の単位と、ほかに何を調整した値なのかを添える必要があります。
係数表の解釈には、目的変数、中心となる説明変数、推定したい量、対象集団の定義が必要です。これらが定まっていなければ、変数を追加して係数が変わった理由を評価できません。
4つの目的で異なる推定対象と評価
記述の対象は、そのデータに含まれる観察標本です。birthwtデータの189例について、出生体重の平均や喫煙状況別の平均を示すのがこれに当たります。標本自体の要約なら評価すべきなのは集計定義と欠測の扱いです。その平均からスプリングフィールドのより広い対象集団へ一般化するなら、標本がどう選ばれたかと、標本抽出に伴う不確実性も別に検討します。
関連の推定では、モデルを置いた上で、説明変数の値が異なる対象者間の条件付き平均の差を推定します。式（2）の妊娠前体重係数なら、喫煙状況を同じとした1ポンド当たりの平均差が推定対象です。記述統計量と異なり、係数には回帰式の仮定が入り、95%信頼区間で推定の不確実性を示します。どの変数に条件づけた関連かを外すと、同じ「平均差」でも推定対象が変わります。
Shmueli（2010）は、データがどのような仕組みで生じたかを説明するモデルと、新しい観測値を当てる予測モデルを区別しました。関連の推定が目的なら、係数が研究上の概念にどう対応し、どの条件の下で解釈されるかが中心です。予測では未知の対象者の出生体重が推定対象となり、係数のp値よりも、観察標本と独立したデータでの予測誤差を評価します。
今回の重回帰モデルでは決定係数が0.0678でした。この値だけを見て「回帰分析が失敗した」と判断するのは適切ではありません。birthwtデータで観察された出生体重のばらつきのうち、当てはめた直線が標本内で説明した割合を表すだけで、妊娠前体重や喫煙状況の係数を何のために推定したかまでは答えないためです。
逆方向の読み違いにも注意が要ります。標本内の決定係数が高くても、未知の対象者で予測がよいとは限りません。この記事で計算した当てはめ値は、観察標本における条件付き平均の例であり、標本外予測性能の評価ではありません。説明と予測を同じ物差しで評価すると、係数の意味とモデルの用途が混ざります。
関連から因果効果へは進めない
因果効果の推定対象は、同じ対象集団に異なる介入を行った場合の反事実アウトカムの差です。そのため、回帰係数を自動的に因果効果と呼ぶのは適切ではありません。回帰式は観察された変数間の条件付き関連を表せますが、介入した場合の結果との差を、式を当てはめただけで識別するわけではありません。不確実性の評価も、回帰式の標準誤差だけでなく、交絡や選択について置いた識別仮定を含めて行う必要があります。
Hernánら（2019）は、記述、予測、反事実に基づく因果推論を別の課題として分類しています。たとえば「喫煙状況が同じ人の中で、妊娠前体重と出生体重がどう関連したか」は式（2）で関連として推定できます。一方、「妊娠前体重を変える介入によって出生体重がどれだけ変化するか」は、介入、対象集団、比較条件を定めた別の推定対象であり、この観察データの回帰係数だけからは答えられません。
多変量モデルに入れた変数は全て交絡因子になり、全係数を因果効果として読めると考えてしまいがちですが、これは適切ではなく、WestreichとGreenland（2013）は、1つのモデルの各係数を同じ種類の効果として提示する誤りをTable 2 fallacyと呼びました。さらにGreenlandら（1999）が示したように、調整変数はp値や係数変化だけではなく、研究で想定する因果構造に基づいて選ぶ必要があります。
このbirthwtデータでは、喫煙状況を加えると妊娠前体重の傾きが1ポンド当たり0.192gだけ小さくなりました。この係数差は標本内の省略変数公式と一致しますが、交絡が除かれた証明ではありません。変化が小さいことも、因果解釈を正当化する根拠にはなりません。
出生体重データで読む係数
MASSパッケージのbirthwtデータセットには、1986年に米国マサチューセッツ州スプリングフィールドのBaystate Medical Centerで得られた189例が収録されています。出生体重の平均は2,944.6gで、妊娠前体重の平均は129.8ポンドでした。喫煙なしは115例、喫煙ありは74例です。
喫煙あり群の平均出生体重は2,771.9g、なし群は3,055.7gで、調整しない平均差は−283.8gでした。妊娠前体重を含む重回帰では、喫煙ありとなしの差は−272.1g、95%信頼区間は−480.4〜−63.8gです。この結果は「妊娠前体重が同じ」というモデル上の条件を置いた、観察データでの平均差として読みます。
95%信頼区間が0を含まないことと、差の医学的な重要性は別の判断です。また、喫煙が出生体重差の原因だと、この解析だけから断定することもできません。統計的な不確実性、臨床的な意味、因果解釈に必要な仮定を分けると、係数が言っている範囲を越えずに済みます。
Rで確かめる調整後の関連
ここではbase Rとstatsパッケージに加え、MASSパッケージのbirthwtデータセットを使います。出生体重を目的変数、妊娠前体重と喫煙状況を説明変数としたモデルから、係数と95%信頼区間を確認します。
library(MASS)
data(birthwt, package = "MASS")
analysis_data <- birthwt[, c("bwt", "lwt", "smoke")]
simple_model <- lm(bwt ~ lwt, data = analysis_data)
adjusted_model <- lm(bwt ~ lwt + smoke, data = analysis_data)
group_mean <- tapply(analysis_data$bwt, analysis_data$smoke, mean)
values <- c(n = nrow(analysis_data), bwt_mean = mean(analysis_data$bwt),
lwt_mean = mean(analysis_data$lwt), no_smoking_n = sum(analysis_data$smoke == 0),
smoking_n = sum(analysis_data$smoke == 1), no_smoking_mean = group_mean["0"],
smoking_mean = group_mean["1"], unadjusted_difference = group_mean["1"] - group_mean["0"],
simple_lwt = coef(simple_model)["lwt"], adjusted_lwt = coef(adjusted_model)["lwt"],
r_squared = summary(adjusted_model)$r.squared,
coefficient_change = coef(simple_model)["lwt"] - coef(adjusted_model)["lwt"])
round(values, 4)
round(cbind(estimate = coef(adjusted_model), confint(adjusted_model)), 2)#> n bwt_mean lwt_mean
#> 189.0000 2944.5873 129.8148
#> no_smoking_n smoking_n no_smoking_mean.0
#> 115.0000 74.0000 3055.6957
#> smoking_mean.1 unadjusted_difference.1 simple_lwt.lwt
#> 2771.9189 -283.7767 4.4291
#> adjusted_lwt.lwt r_squared coefficient_change.lwt
#> 4.2367 0.0678 0.1924
#> estimate 2.5 % 97.5 %
#> (Intercept) 2501.13 2045.73 2956.52
#> lwt 4.24 0.90 7.57
#> smoke -272.08 -480.39 -63.77読み方は、妊娠前体重の係数4.24が、喫煙状況を同じにした時の1ポンド当たりの平均出生体重差を表す、というものです。喫煙係数−272.08は、妊娠前体重を同じ値に固定したモデル上での群間差です。
読み違えやすいのは、2つの係数をそのまま介入効果へ置き換えることです。ここで推定したのは189例の観察データにおける調整済み関連であり、未知の対象者への予測精度も評価していません。
par(mar = c(5.1, 4.1, 4.1, 8.1))
plot(analysis_data$lwt, analysis_data$bwt,
col = ifelse(analysis_data$smoke == 1, "#C44E52", "#4C72B0"),
pch = 16, xlab = "Maternal weight before pregnancy (lb)",
ylab = "Birth weight (g)")
abline(adjusted_model$coefficients[c("(Intercept)", "lwt")],
col = "#4C72B0", lwd = 2)
abline(c(adjusted_model$coefficients["(Intercept)"] +
adjusted_model$coefficients["smoke"],
adjusted_model$coefficients["lwt"]), col = "#C44E52", lwd = 2)
par(xpd = NA)
legend("topright", inset = c(-0.34, 0),
legend = c("No smoking", "Smoking"),
col = c("#4C72B0", "#C44E52"), pch = 16, lty = 1, bty = "n")
点は直線の周囲に広く分布しています。したがって、直線は個人の出生体重を決める規則ではなく、各群の条件付き平均を要約したものとして読む必要があります。図の平行な2本線は交互作用を含めない今回のモデルの形であり、実際の関係が必ず平行だと証明した図ではありません。
論文で残す係数の文脈
Leffondréら（2014）の議論に沿えば、論文では変数名とp値だけでなく、係数の単位、参照群、調整した変数、95%信頼区間を一緒に記載する方が、読者が推定量を復元しやすくなります。今回なら「妊娠前体重と喫煙状況を含む線形回帰モデルで、妊娠前体重10ポンド当たりの平均出生体重差は42.37g、95%信頼区間9.03〜75.71g」と書けます。
ただし、この記載から因果効果へ進むには、対象集団、介入としての説明変数、比較条件、調整変数を選んだ根拠が別途必要です。Greenlandら（1999）の因果ダイアグラムは、その根拠を観察された関連の強さだけに任せないための枠組みを与えます。予測モデルとして報告するなら、係数表や標本内の決定係数に加えて、未知データでの性能を確かめる設計が必要になります。
p値だけで回帰分析の目的達成を判断するのも避けた方がよいでしょう。p値は、置いたモデルと帰無仮説の下でデータとの不整合を測るもので、係数が研究目的に対応しているか、差が医学的に重要か、予測が役立つかを一度に判定する値ではありません。
研究目的に戻る係数解釈
回帰分析を始める前に必要なのは、どの変数を投入するかの長い候補表ではなく、目的変数、中心となる説明変数、推定したい量、対象集団がつながった文章です。その文章が関連の推定を求めているなら、係数には単位と調整条件を添えます。予測が目的なら、未知データでの誤差へ評価軸を移します。因果効果が目的なら、回帰式より先に比較したい介入と因果構造を定める必要があります。
係数は計算結果の終点ではありません。研究目的と結果を結ぶ途中の表現です。
実際の研究では、同じ説明変数でも符号化、非線形性、交互作用、欠測への対応によって解釈が変わります。自分のデータで推定量と調整方針を切り分けにくい場合、Dr.データサイエンスでは、研究目的から解析方針と報告内容をつなぐ相談を受け付けています。
参考文献
- Shmueli G. To Explain or to Predict? Statistical Science. 2010;25(3):289-310.
- Hernán MA, Hsu J, Healy B. A Second Chance to Get Causal Inference Right: A Classification of Data Science Tasks. CHANCE. 2019;32(1):42-49.
- Leffondré K, Jager KJ, Boucquemont J, Stel VS, Heinze G. Representation of exposures in regression analysis and interpretation of regression coefficients: basic concepts and pitfalls. Nephrology Dialysis Transplantation. 2014;29(10):1806-1814.
- Greenland S, Pearl J, Robins JM. Causal diagrams for epidemiologic research. Epidemiology. 1999;10(1):37-48.
- Westreich D, Greenland S. The Table 2 Fallacy: Presenting and Interpreting Confounder and Modifier Coefficients. American Journal of Epidemiology. 2013;177(4):292-298.
この記事は回帰分析カテゴリーの新体系No.1です。前の記事と次の記事は、現行のサイト台帳ではまだ設定されていません。
