退院先を自宅、回復期リハビリテーション病院、介護施設の3つに分けて記録した研究では、アウトカムを自宅かそれ以外かの二値へまとめ、ロジスティック回帰を当てる解析がよく見られます。しかし、ある要因がリハビリ病院への退院とは関係せず、介護施設への退院だけを増やしている場合、二値へまとめたオッズ比は両者を混ぜた中間の値になり、どちらの比較についても関連の大きさを表さなくなります。
そこで、ランダムに生成した1,200例の退院先データで、多項ロジスティック回帰の式と係数の読み方を確かめます。アウトカムが二値、順序のあるカテゴリー、順序のないカテゴリーのどれに当たるかで、使えるモデルも係数の意味も変わります。
アウトカムの型とモデルの対応
二値アウトカムには二値ロジスティック回帰、順序のあるカテゴリーには比例オッズモデル、順序のないカテゴリーには多項ロジスティック回帰、というのが基本の対応です。ここでの順序とは、重症度の段階のように隣り合うカテゴリーの間で悪化の方向が一つに定まる並びを指し、自宅、リハビリ病院、介護施設のように一つの軸が決まらない退院先は、順序のないカテゴリーとして扱うのが自然です。
もっとも、型の判定はデータの見た目ではなく、研究目的とカテゴリーの意味から決まります。退院先を在宅から遠い順に並べることもできますが、その並びを置くと、年齢や独居が全ての境目で同じ向き、同じ大きさで効くという仮定まで一緒に持ち込むため、順序を置く前にその仮定が研究目的に照らして受け入れられるかを考えることになります。
なお、回帰の係数が何と何を比べた量なのかは回帰分析の目的で扱っています。また、アウトカムが複数あることを指す多変量解析と、説明変数が複数ある多変数回帰の区別は多変量解析と多変数回帰で扱っており、多項ロジスティック回帰は、カテゴリーが3つ以上あっても1つのアウトカムを説明する多変数回帰に当たります。
基準カテゴリーとの対数オッズで表す式
多項ロジスティック回帰では、カテゴリーの1つを基準に選び、残りのカテゴリーそれぞれについて、基準に対する対数オッズを説明変数の一次式で表します。カテゴリー数を$K$、基準をカテゴリー1とすると、モデルは次の形になります。
$$\log\frac{\Pr(Y=k\mid\mathbf{x})}{\Pr(Y=1\mid\mathbf{x})}=\mathbf{x}^\top\boldsymbol{\beta}_k,\quad k=2,\dots,K \tag{1}$$
この式(1)は、自宅を基準にした時のリハビリ病院の対数オッズと介護施設の対数オッズを、それぞれ別の係数ベクトルで表すという意味で、$Y$は退院先、$\mathbf{x}$は切片を含む説明変数のベクトル、$\boldsymbol{\beta}_k$はカテゴリー$k$に固有の係数ベクトルです。
そこで、各カテゴリーの確率へ戻すには、基準カテゴリーの係数を0に固定し、各カテゴリーの指数を全カテゴリーの合計で割ります。
$$\Pr(Y=k\mid\mathbf{x})=\frac{\exp(\mathbf{x}^\top\boldsymbol{\beta}_k)}{\sum_{j=1}^{K}\exp(\mathbf{x}^\top\boldsymbol{\beta}_j)},\quad \boldsymbol{\beta}_1=\mathbf{0} \tag{2}$$
この式(2)は、全カテゴリーの確率の合計が1になるように割り戻す式であり、1つの係数が動くと全カテゴリーの確率が同時に動きます。そのため、あるカテゴリーのオッズ比が1を超えていても、そのカテゴリーの確率が上がるかどうかは他のカテゴリーの係数と説明変数の値にも左右され、結果を伝える時には予測確率を併せて示す方が読み違いが少なくなります。
そのため、係数の解釈は常に基準カテゴリーと比べてという条件つきになります。基準をどれにするかで個々の係数の値は変わっても、当てはまりや予測確率は変わらないため、研究目的に照らして比較の相手として自然なカテゴリー、退院先であれば自宅を基準に置くのが読みやすい選び方です。
カテゴリー数に応じて係数が増える
説明変数が2つの時、二値ロジスティック回帰の係数は切片を含めて3個ですが、3カテゴリーの多項モデルでは2本の式それぞれに3個ずつ、合計6個になります。de Jongら（2019）はシミュレーション研究で、多項の予測モデルに必要な症例数はカテゴリーごとの症例数とパラメータ数から考える必要があり、二値モデルの目安を流用できないことを示しています。
実際、各式の推定精度を支えているのは全体の症例数ではなく、基準カテゴリーとそのカテゴリーの症例数です。後で扱うデータでも、介護施設237例と自宅615例の組み合わせが介護施設側の式の情報量を決めており、症例の少ないカテゴリーを持つ研究で説明変数を増やすと、その式の係数だけが不安定になるという形で問題が表れます。
二値へまとめると効果が薄まる
冒頭の例のように、独居がリハビリ病院への退院とは関係せず、介護施設への退院だけを増やす状況を考えます。自宅以外をまとめた二値のオッズ比は、効果のないリハビリ病院と効果のある介護施設を混ぜたものになるため、介護施設の効果を小さく、リハビリ病院の効果を大きく見せる方向へ偏り、その混ざり方も各カテゴリーの割合によって変わります。
逆に、基準カテゴリーとの二値比較を別々に行う方法、つまり自宅とリハビリ病院だけ、自宅と介護施設だけのデータで二値ロジスティック回帰を当てる方法もあります。BeggとGray（1984）は、この個別回帰の係数が多項モデルの係数に近い値になる一方で一致はせず、効率も多項モデルより低いことを示しており、説明変数全体の検定や予測確率まで求める場合には多項モデルを1本当てる方が扱いやすくなります。
順序を置けるかは比例オッズの仮定で決まる
順序のあるアウトカムには、McCullagh（1980）が体系化した比例オッズモデルが使えます。カテゴリーを$1<2<\dots<K$と並べ、境目$j$ごとの累積確率に対するロジットを考えると、モデルは次の形です。
$$\log\frac{\Pr(Y\le j\mid\mathbf{x})}{\Pr(Y>j\mid\mathbf{x})}=\alpha_j-\mathbf{x}^\top\boldsymbol{\beta},\quad j=1,\dots,K-1 \tag{3}$$
この式(3)は、境目ごとに切片$\alpha_j$だけが異なり、説明変数の係数$\boldsymbol{\beta}$は全ての境目で共通だという意味で、ここでの$\mathbf{x}$は切片を含まない説明変数のベクトルです。そのため、3カテゴリーで説明変数が2つならパラメータは切片2個と係数2個の合計4個となり、多項モデルの6個より少なく済み、順序を正しく反映できていれば推定の効率も上がります。
一方、共通の係数という仮定が崩れると、比例オッズモデルのオッズ比は、境目ごとに異なる効果を1つにまとめた値になります。Brant（1990）はこの仮定の検討方法を示し、BenderとGrouven（1998）は累積の境目ごとに二値ロジスティック回帰を当ててオッズ比を比べる考え方を示しているので、順序モデルを選ぶ前に境目ごとのオッズ比が近いかどうかを見ておくことが判断材料になります。
この確認で境目ごとのオッズ比が大きく異なり、しかもその差が研究目的にとって意味を持つなら、係数を共通にする仮定がそのまま推定値を歪めます。その場合は順序を置かず、多項ロジスティック回帰で各カテゴリーへの効果を別々に示す方が、推定している量と報告する量が一致します。
医学研究での位置づけと報告事項
Biesheuvelら（2008）は、複数の疾患カテゴリーを扱う診断研究で、アウトカムを二値へまとめたり二値モデルを繰り返したりする解析が多いことを指摘し、多項ロジスティック回帰はもっと使われてよいと論じています。また、Roozenbeekら（2011）は外傷性脳損傷の臨床試験データの再解析で、順序のある転帰を二値にまとめずに順序解析すると統計的効率が上がることを報告しています。
こうした文献を踏まえると、論文で多項ロジスティック回帰を報告する時は、基準カテゴリー、カテゴリーごとの症例数、各式のオッズ比と95%信頼区間を示し、説明変数全体の効果を尤度比検定で評価したかどうかを書くことが、読者の再現と解釈を助けます。順序モデルを使った場合は、比例オッズ仮定をどの方法で確認し、境目ごとのオッズ比がどの程度異なったかも併せて記載します。
退院先1,200例のデータを用意する
年齢を平均75歳、標準偏差8歳の正規乱数、独居の割合を0.3として、1,200例のランダムに生成したデータセットを用意しました。自宅を基準に、年齢10歳あたりの対数オッズをリハビリ病院0.3、介護施設0.9、独居の対数オッズをリハビリ病院0、介護施設0.9と置き、多項ロジスティック回帰の当てはめにはnnetパッケージを使っています。
library(nnet)
set.seed(42)
n_patients <- 1200
age <- round(rnorm(n_patients, mean = 75, sd = 8))
age10 <- (age - 75) / 10
alone <- rbinom(n_patients, size = 1, prob = 0.3)
# 基準は自宅。独居はリハビリ病院に効果を置かず、介護施設だけに0.9を置く
true_beta <- rbind(rehab = c(-0.6, 0.3, 0.0), facility = c(-1.4, 0.9, 0.9))
eta <- cbind(0, cbind(1, age10, alone) %*% t(true_beta))
true_prob <- exp(eta) / rowSums(exp(eta))
dest_index <- apply(true_prob, 1, function(p) sample(1:3, size = 1, prob = p))
destination <- factor(c("home", "rehab", "facility")[dest_index],
levels = c("home", "rehab", "facility"))
discharge <- data.frame(destination, age10, alone)
table(discharge$destination)#> home rehab facility
#> 615 348 237 この設定で生成された退院先は自宅615例（51.2%）、リハビリ病院348例（29.0%）、介護施設237例（19.8%）で、独居は383例でした。独居の真の効果はリハビリ病院でオッズ比1、介護施設でオッズ比約2.46と置いてあり、この差が二値へまとめた時に見えなくなるかどうかを後で比べます。
オッズ比と尤度比検定の読み方
このデータに年齢と独居を説明変数とする多項ロジスティック回帰を当てはめ、独居を除いたモデルとの尤度比検定も計算しました。
fit_multi <- multinom(destination ~ age10 + alone, data = discharge,
trace = FALSE, reltol = 1e-12, maxit = 1000)
round(exp(coef(fit_multi)), 3)
fit_no_alone <- multinom(destination ~ age10, data = discharge,
trace = FALSE, reltol = 1e-12, maxit = 1000)
lr_stat <- deviance(fit_no_alone) - deviance(fit_multi)
cat("LR test alone:", round(lr_stat, 2), "df", 2, "p",
signif(pchisq(lr_stat, df = 2, lower.tail = FALSE), 3), "\n")#> (Intercept) age10 alone
#> rehab 0.555 1.201 1.129
#> facility 0.253 2.113 2.681
#> LR test alone: 37.93 df 2 p 5.8e-09 この出力の1行目はリハビリ病院対自宅、2行目は介護施設対自宅のオッズ比です。独居のオッズ比は介護施設で2.681（95%信頼区間1.945〜3.697）、リハビリ病院で1.129（0.843〜1.513、Wald検定p=0.415）となり、置いた真の値2.46と1に近い推定がカテゴリーごとに別々に得られています。
ただし、リハビリ病院側の信頼区間が1をまたぐことから言えるのは、このデータでは独居とリハビリ病院への退院との関連が明確には示されなかった、というところまでです。年齢10歳あたりのオッズ比はリハビリ病院1.201（1.014〜1.422）、介護施設2.113（1.721〜2.596）で、同じ説明変数でもカテゴリーによって効果の大きさが異なることが、多項モデルでは別々の係数として表れます。
独居が退院先全体に関係するかを1つの検定で評価する場合は、2本の式の独居の係数を同時に0と置いたモデルと比べる尤度比検定を使い、今回は統計量37.93、自由度2、$p=5.8\times10^{-9}$でした。このp値は2つの比較をまとめた関連の有無を示すもので、どのカテゴリーで効果があり、その大きさが臨床的に意味を持つかどうかは、個々のオッズ比と信頼区間から読むことになります。
二値へのまとめと累積の境目を比べる
同じデータで、自宅以外を1とする二値と、介護施設だけを1とする二値のロジスティック回帰も当てはめました。後者は自宅とリハビリ病院を合わせた群と介護施設を比べる境目であり、前者と合わせると、退院先に順序を置いた時の累積の境目2つに当たります。
not_home <- as.integer(discharge$destination != "home")
facility <- as.integer(discharge$destination == "facility")
round(exp(coef(glm(not_home ~ age10 + alone, family = binomial, data = discharge))), 3)
round(exp(coef(glm(facility ~ age10 + alone, family = binomial, data = discharge))), 3)#> (Intercept) age10 alone
#> 0.824 1.486 1.606
#> (Intercept) age10 alone
#> 0.162 1.972 2.561 この二値へのまとめでは、独居のオッズ比は1.606（1.253〜2.058）でした。この値は効果のないリハビリ病院と効果のある介護施設を混ぜた結果であり、多項モデルの介護施設2.681より小さく、リハビリ病院1.129より大きいため、どちらの退院先についても独居との関連の大きさを表していません。
これに対して、累積の境目として見直すと、独居のオッズ比は1つ目の境目で1.606、2つ目の境目で2.561（1.896〜3.459）となり、信頼区間の重なりもわずかです。年齢でも1.486と1.972で差があり、このデータに比例オッズモデルを当てると、共通のオッズ比は1つ目の境目では過大、2つ目の境目では過小な要約になる方向です。
なお、自宅とそれぞれのカテゴリーだけで別々に当てた個別の二値回帰では、係数が多項モデルとほぼ同じで、標準誤差の比は1.000〜1.019でした。BeggとGray（1984）が示した効率の差はこのデータでは2%未満と小さかったものの、その大きさは説明変数やカテゴリーの割合で変わるため、多項モデルを当てる手間を省く根拠にはしにくいと考えています。
予測確率で結果を伝える
オッズ比は基準カテゴリーとの比較の量なので、臨床研究の読者には確率で示す方が伝わりやすくなります。年齢60〜90歳について独居の有無ごとに各退院先の予測確率を計算し、図1に示しました。
curve_grid <- expand.grid(age = 60:90, alone = c(0, 1))
curve_grid$age10 <- (curve_grid$age - 75) / 10
pred_curve <- predict(fit_multi, newdata = curve_grid, type = "probs")
par(mfrow = c(1, 2))
for (status in 0:1) {
rows <- curve_grid$alone == status
matplot(60:90, pred_curve[rows, ], type = "l", lty = 1, lwd = 2,
col = c("grey40", "steelblue", "firebrick"), ylim = c(0, 0.7),
xlab = "Age (years)", ylab = "Predicted probability",
main = if (status == 0) "Not living alone" else "Living alone")
legend("topleft", legend = c("Home", "Rehabilitation", "Care facility"),
col = c("grey40", "steelblue", "firebrick"), lwd = 2, bty = "n")
}
この図1の右側では、独居の人の介護施設への予測確率が年齢とともに上がり、85歳で0.449と自宅の0.314を上回ります。一方、リハビリ病院の予測確率は表示した年齢と独居の組み合わせで0.236〜0.307の範囲にとどまり、介護施設の確率が増えた分は主に自宅の確率が下がることで釣り合っています。
この形は置いた係数から生じたもので、年齢のオッズ比がリハビリ病院で1.201と1を超えていても、介護施設の確率が大きく増える高齢の領域ではリハビリ病院の確率が横ばいから低下へ向かいます。式(2)の割り戻しがそのまま表れた結果であり、オッズ比だけを報告すると、この動きは読者に伝わりません。
アウトカムの型からモデルを決める時の判断
アウトカムが3つ以上のカテゴリーを持つ時は、カテゴリー間に研究目的に沿った順序があるかを考え、順序を置けるなら境目ごとのオッズ比を比べて比例オッズの仮定を確かめます。仮定が崩れる時や順序がない時は、多項ロジスティック回帰で基準カテゴリーを明示し、オッズ比、尤度比検定、予測確率を合わせて報告します。
一方、多項モデルはカテゴリーごとに係数を持つため、症例の少ないカテゴリーがある時は説明変数の数を事前に絞る必要があります。p値を見ながらカテゴリーをまとめ直したり、モデルを切り替えたりする進め方は、推定対象そのものを結果に合わせて変えてしまうため推奨しません。
こうした選択を自分の研究で行う時は、カテゴリーの意味、症例数、境目ごとのオッズ比を合わせて判断することになります。Dr.データサイエンスでは、アウトカムの型に応じたモデル選択を含む解析方針のご相談を承っています。
参考文献
- Begg CB, Gray R. Calculation of polychotomous logistic regression parameters using individualized regressions. Biometrika. 1984;71(1):11-18.
- McCullagh P. Regression models for ordinal data. Journal of the Royal Statistical Society: Series B (Methodological). 1980;42(2):109-127.
- Brant R. Assessing proportionality in the proportional odds model for ordinal logistic regression. Biometrics. 1990;46(4):1171-1178.
- Bender R, Grouven U. Using binary logistic regression models for ordinal data with non-proportional odds. Journal of Clinical Epidemiology. 1998;51(10):809-816.
- Biesheuvel CJ, Vergouwe Y, Steyerberg EW, Grobbee DE, Moons KG. Polytomous logistic regression analysis could be applied more often in diagnostic research. Journal of Clinical Epidemiology. 2008;61(2):125-134.
- de Jong VMT, Eijkemans MJC, van Calster B, Timmerman D, Moons KGM, Steyerberg EW, van Smeden M. Sample size considerations and predictive performance of multinomial logistic prediction models. Statistics in Medicine. 2019;38(9):1601-1619.
- Roozenbeek B, Lingsma HF, Perel P, et al. The added value of ordinal analysis in clinical trials: an example in traumatic brain injury. Critical Care. 2011;15(3):R127.
