手術を受けた患者について、術後30日以内に合併症が起きたかどうかを解析する場面がありますが、アウトカムが0か1の場合、ロジスティック回帰から得られる代表的な指標はオッズ比です。たとえば2.33という値が出ても「合併症が2.33倍起きやすい」とそのまま書くことはできず、2.33倍になっているのは確率ではなく、イベントが起きる確率を起きない確率で割ったオッズです。
この違いはイベントがまれな時には目立ちませんが、合併症が20%、30%と起きるデータでは、オッズ比とリスク比の差を無視しにくくなります。そこで、予定手術後の合併症を扱う2,000例のランダムに生成したデータセットを使い、確率からオッズ、ロジスティック回帰の係数、調整オッズ比、予測確率までをつなげて確認します。
ここではアウトカムを$Y$とし、合併症が起きた場合を$Y=1$、起きなかった場合を$Y=0$、説明変数のベクトルを$\mathbf{x}$、係数を$\boldsymbol{\beta}$、合併症が起きる条件付き確率$\Pr(Y=1\mid\mathbf{x})$を$p$と書きます。曝露として糖尿病の既往を置き、年齢を共変量として扱います。
二値アウトカムをロジットで扱う理由
連続アウトカムなら、説明変数とアウトカムの関係を直線で表し、係数を単位あたりの平均値の変化として読むことができ、この考え方は回帰分析の目的で扱いました。一方、合併症の有無は0か1なので、同じ考え方をそのまま確率へ持ち込むと扱いが少し変わります。
ここで「二値アウトカムには線形回帰を使えない」と覚えると単純化しすぎており、0と1を目的変数にした線形確率モデルも計算できます。ただし、予測値が0未満や1を超えることがあるうえ、誤差分散も一定にならないため、二値データの確率モデルとしては扱いにくい性質があります。
そこでロジスティック回帰では、確率をそのまま線形式へ結びつけず、まずオッズへ変換します。
$$\mathrm{odds}=\frac{p}{1-p} \tag{1}$$
式(1)ではイベントが起きる確率$p$を起きない確率$1-p$で割っており、たとえば$p=0.05$ならオッズは0.053、$p=0.50$なら1、$p=0.80$なら4になります。つまり、確率とオッズは同じ尺度ではありません。
さらに、確率は0から1の範囲にありますが、オッズは0から正の無限大まで動き、オッズの対数をとると負の無限大から正の無限大まで広がります。この変換によって、二値アウトカムの確率を説明変数の線形結合と結びつけられるようになります。
probability <- seq(0.01, 0.95, by = 0.01)
odds <- probability / (1 - probability)
plot(probability, odds, type = "l",
xlab = "Probability", ylab = "Odds", ylim = c(0, 5))
abline(a = 0, b = 1, lty = 2)
この区別について、BlandとAltman（2000）もオッズと確率を別の量として扱う必要性を説明しており、「イベントがまれならオッズ比をリスク比に近い値として読める」という近似があります。ただし、どこから突然近似が十分になるという境界があるわけではなく、両者の差は基準リスクと効果の大きさの両方によって変わります。
オッズ比は2つのオッズを比べている
糖尿病の既往ありを$p_1$、既往なしを$p_0$とすると、オッズ比は次のように書けます。
$$\mathrm{OR}=\frac{p_1/(1-p_1)}{p_0/(1-p_0)} \tag{2}$$
式(2)は糖尿病既往あり群のオッズを既往なし群のオッズで割った値であり、リスク$p_1$を$p_0$で割るリスク比とは分子も分母も異なります。
実際に2,000例をランダムに生成すると、糖尿病既往ありでは257例に合併症が起き、316例には起きなかったのに対し、既往なしでは合併症ありが323例、なしが1,104例でした。その結果、オッズは既往ありで0.813、既往なしで0.293となり、粗オッズ比は2.780でした。
set.seed(42)
age <- round(rnorm(2000, mean = 65, sd = 10))
diabetes <- rbinom(2000, size = 1,
prob = plogis(-1.0 + 0.05 * (age - 65)))
complication <- rbinom(2000, size = 1,
prob = plogis(-1.2 + 0.85 * diabetes + 0.035 * (age - 65)))
cohort <- data.frame(age, diabetes, complication)
tab <- table(cohort$diabetes, cohort$complication)
odds_by_group <- tab[, "1"] / tab[, "0"]
or_crude <- unname(odds_by_group["1"] / odds_by_group["0"])
round(c(odds_by_group, or_crude = or_crude), 3)#> 0 1 or_crude
#> 0.293 0.813 2.780ここで2.780を「合併症リスクが2.78倍」と読むと尺度を取り違えており、同じ集計表から求めたリスクは糖尿病既往ありで0.449、既往なしで0.226なので、リスク比は1.98です。したがって、オッズ比2.78が表しているのは合併症が起きる確率そのものの倍率ではありません。
係数を指数化すると調整オッズ比になる
ロジスティック回帰では、オッズの対数である対数オッズを説明変数の線形結合として表します。
$$\log\left\{\frac{\Pr(Y=1\mid\mathbf{x})}{1-\Pr(Y=1\mid\mathbf{x})}\right\}=\beta_0+\boldsymbol{\beta}^{\top}\mathbf{x} \tag{3}$$
式(3)の左辺が対数オッズ、$\beta_0$が切片、$\boldsymbol{\beta}$が説明変数に対応する係数であり、回帰係数そのものは確率の変化ではなく対数オッズの変化を表しています。
また、交互作用項や非線形項がないモデルで、説明変数$x_j$だけを1単位増やして他の説明変数を同じ値に固定すると、対数オッズの差は$\beta_j$となります。これを指数化すると、次の関係になります。
$$\exp(\beta_j)=\frac{\mathrm{odds}(x_j+1)}{\mathrm{odds}(x_j)} \tag{4}$$
式(4)から、$\exp(\beta_j)$は他の説明変数を同じ値にした条件で$x_j$が1単位異なる2者のオッズ比となり、糖尿病のような0/1変数なら既往あり対なし、年齢なら1歳増加あたりのオッズ比です。
実際に年齢を調整して糖尿病の係数を求めると、調整オッズ比は2.331、95%信頼区間は1.884から2.884でした。データ生成時に置いた糖尿病の係数は0.85なので、真の条件付きオッズ比は$\exp(0.85)=2.34$です。
fit_adjusted <- glm(complication ~ diabetes + age,
family = binomial(link = "logit"), data = cohort)
or_adj <- exp(cbind(OR = coef(fit_adjusted),
confint.default(fit_adjusted)))
round(or_adj["diabetes", ], 3)#> OR 2.5 % 97.5 %
#> 2.331 1.884 2.884このデータでは年齢が高いほど糖尿病既往を持つ確率が高く、同時に合併症の確率も高くなるよう生成しているため、年齢による交絡が生じます。その結果、粗オッズ比2.780と年齢調整後の2.331にはかなりの差が出ました。
ただし、2.780−2.331=0.449をそのまま「年齢による交絡の量」と呼ぶことはできず、オッズ比には交絡がなくても条件付きオッズ比と周辺オッズ比が一致しないnoncollapsibilityという性質があります。この点は後で改めて確かめます。
さらに、年齢のような連続変数では報告単位も考える必要があり、1歳あたりのオッズ比が1に近く読みにくい場合は、10歳あたりのオッズ比として$\exp(10\beta_{\mathrm{age}})$を報告できます。これは単位を変えているだけなので、モデルを当てはめ直す必要はありません。
オッズ比2.33はリスク比2.33ではない
「糖尿病既往ありでは合併症が2.3倍起きやすかった」と書くとオッズの倍率をリスクの倍率へ置き換えることになり、今回のようにイベントが珍しくないデータでは両者の差が実際に大きくなります。この点について、ZhangとYu（1998）はアウトカム頻度が高くなるにつれてオッズ比とリスク比の差が広がることを示しています。
そこで、同じ条件の患者について非曝露時のリスク$p_0$とオッズ比$\mathrm{OR}$が分かれば、対応するリスク比は次の式から求められます。
$$\mathrm{RR}=\frac{\mathrm{OR}}{(1-p_0)+p_0\mathrm{OR}} \tag{5}$$
式(5)は同じ条件における$p_0$とオッズ比からリスク比へ変換する関係で、ZhangとYu（1998）はこの関係を調整オッズ比からリスク比を近似する方法として提案しました。ただし、調整ORへ集団全体の非曝露リスクを機械的に代入すれば常に調整RRが得られる、という使い方はできません。
今回のモデルでは65歳で糖尿病既往なしの予測リスクは約0.23、既往ありでは約0.41なので、同じ65歳で比べたリスク比は約1.8となり、調整オッズ比2.331よりかなり1に近くなります。この差があるため、私はイベント頻度が高い研究でオッズ比だけを日本語の「○倍起きやすい」へ置き換える書き方は適切ではありません。
また、一定なのはオッズ比であってリスク差ではなく、ロジスティック回帰で糖尿病の係数を一定にしても、年齢によって基準リスクが変われば糖尿病ありとなしの予測確率の差も動きます。
age_grid <- 40:90
risk_no <- predict(fit_adjusted,
newdata = data.frame(diabetes = 0, age = age_grid), type = "response")
risk_yes <- predict(fit_adjusted,
newdata = data.frame(diabetes = 1, age = age_grid), type = "response")
plot(age_grid, risk_yes, type = "l", lwd = 2.2, ylim = c(0, 0.9),
xlab = "Age (years)", ylab = "Predicted probability of complication")
lines(age_grid, risk_no, lwd = 2.2, lty = 2)
legend("topleft", c("Diabetes: yes", "Diabetes: no"),
lty = c(1, 2), lwd = 2.2, bty = "n")
この予測確率との関係について、Nortonら（2019）はロジスティック回帰の結果を予測確率の変化として表すmarginal effectsを解説しています。コホート研究などリスク自体を推定できる研究では、オッズ比だけでは大きさを直感的に読み取りにくい場合に予測確率やリスク差を併記すると、どの程度の確率差があるのかを伝えやすくなります。
粗ORと調整ORの差は交絡だけではない
ここはロジスティック回帰で特に読み違えやすいところで、粗オッズ比と調整オッズ比の差をそのまま共変量による交絡の大きさと解釈すると、尺度の性質を取りこぼします。オッズ比では、交絡がない状況でも周辺オッズ比と条件付きオッズ比が一致するとは限りません。
そこで、糖尿病を年齢とは独立に割り付けた20,000例を別に生成し、この条件では少なくとも年齢による糖尿病と合併症の交絡が生じないようにしました。それでも、この1回のデータでは年齢を含むモデルのオッズ比が2.292、年齢を含まないモデルでは2.265となりました。
age_nc <- round(rnorm(20000, mean = 65, sd = 10))
diabetes_nc <- rbinom(20000, size = 1, prob = 0.35)
complication_nc <- rbinom(20000, size = 1,
prob = plogis(-1.2 + 0.85 * diabetes_nc + 0.035 * (age_nc - 65)))
nc <- data.frame(age = age_nc, diabetes = diabetes_nc,
complication = complication_nc)
or_nc_cond <- exp(coef(glm(complication ~ diabetes + age,
family = binomial, data = nc))["diabetes"])
or_nc_marg <- exp(coef(glm(complication ~ diabetes,
family = binomial, data = nc))["diabetes"])
round(c(conditional = or_nc_cond, marginal = or_nc_marg), 3)#> conditional marginal
#> 2.292 2.265ただし、この0.027を「noncollapsibilityの大きさ」とそのまま読むこともできず、1回のランダム標本なので標本変動が含まれています。ここで確認したいのは、交絡を作っていないデータでも、共変量を含む条件付きORと集団をまとめた周辺ORが必ず同じ値になるわけではないことです。
この区別についてGreenland、Robins、Pearl（1999）はnoncollapsibilityと交絡を分けて論じており、最初の2,000例で粗ORが2.780から2.331へ動いたという事実だけから、0.449全体を交絡の大きさへ変換することはできません。このデータ生成には年齢による交絡が実際に入っていますが、粗ORと調整ORの差には尺度の性質と標本変動も関係します。
そのため、私なら粗ORと調整ORのどちらを使うかを差の大きさから決めず、年齢を同じ値にした患者同士の関連を知りたいなら条件付きORを使い、対象集団全体で糖尿病ありとなしの予測リスクを比較したいなら調整モデルから標準化した予測確率や周辺効果を求めます。ただし、交絡のある観察データで未調整の粗ORをそのまま「周辺効果」と呼ぶこととは区別が必要です。
論文ではORと予測確率を分けて書く
症例対照研究では研究対象者を症例と対照から抽出するため、その標本だけから発症リスクを直接求めることは通常できませんが、標準的な症例対照研究の抽出条件のもとではオッズ比を推定できます。BlandとAltman（2000）が説明しているように、オッズ比が中心的な関連指標となるのは、リスク比の代わりに何となく使われているからではなく、研究デザインと指標の性質によるものです。
一方、コホート研究ではリスクそのものを推定できるため、ロジスティック回帰を使った場合でも結果の説明をオッズ比だけに限定する必要はありません。今回なら「年齢を調整した糖尿病既往のオッズ比は2.331、95%信頼区間1.884〜2.884」とモデルの結果を書いたうえで、代表的な年齢での予測確率や、研究目的に合うなら標準化したリスク差・リスク比を併記できます。
また、モデルを安定して推定できるだけの情報量があるかも別に確認する必要があり、Peduzziら（1996）のシミュレーションではEPVが10未満の条件で係数の偏りや信頼区間の被覆などに問題が増えました。ただし、「10イベントあれば1変数を安全に入れられる」という普遍的な境界ではなく、van Smedenら（2016）は推定の問題がEPVだけでなく総標本サイズや分離などにも依存することを示しています。
今回の2,000例では合併症が580例あり、説明変数は糖尿病と年齢の2つなので、少なくともイベント数が極端に少ない例ではありません。ただし、「EPV=290だからモデルは妥当」と結論するのではなく、連続変数の関数形、分離、強い相関、欠測、クラスタリングなどは別に確認します。
さらに、施設ごとに患者がまとまっているデータでは観測が独立という前提が崩れるため、通常のロジスティック回帰の標準誤差をそのまま使う根拠が弱くなります。研究目的に応じてクラスターロバスト標準誤差、GEE、混合効果ロジスティック回帰などを考えることになり、標準誤差だけを扱う場面はロバスト標準誤差ともつながります。
そのため論文の一文へ移す時は「2.33倍」という数字だけを取り出さず、今回のモデルなら「年齢を調整すると、糖尿病既往ありの患者では既往なしに比べて術後30日以内の合併症オッズが2.33倍だった（95%信頼区間1.88〜2.88）」までがオッズ比から直接言える範囲です。そのうえで予測確率を示すなら、「65歳では予測確率が約0.23と0.41だった」と別の尺度として書きます。
実際の研究では、条件付きの関連を報告したいのか、対象集団全体でのリスク差やリスク比を知りたいのかによって、ロジスティック回帰の結果から次に計算する量が変わります。Dr.データサイエンスでは、こうした推定対象と解析モデル、結果の報告方法を含めた解析方針のご相談を承っています。
参考文献
- Peduzzi P, Concato J, Kemper E, Holford TR, Feinstein AR. A simulation study of the number of events per variable in logistic regression analysis. Journal of Clinical Epidemiology. 1996;49(12):1373-1379.
- Zhang J, Yu KF. What’s the relative risk? A method of correcting the odds ratio in cohort studies of common outcomes. JAMA. 1998;280(19):1690-1691.
- Greenland S, Robins JM, Pearl J. Confounding and collapsibility in causal inference. Statistical Science. 1999;14(1):29-46.
- Bland JM, Altman DG. Statistics notes. The odds ratio. BMJ. 2000;320(7247):1468.
- van Smeden M, de Groot JAH, Moons KGM, Collins GS, Altman DG, Eijkemans MJC, Reitsma JB. No rationale for 1 variable per 10 events criterion for binary logistic regression analysis. BMC Medical Research Methodology. 2016;16:163.
- Norton EC, Dowd BE, Maciejewski ML. Marginal Effects—Quantifying the Effect of Changes in Risk Factors in Logistic Regression Models. JAMA. 2019;321(13):1304-1305.
