10.5 順序ロジスティック回帰分析

(1) 累積ロジスティックモデル

今まで説明したデータは目的変数が「0：反応無、1：反応有」という2分類のものでした。 しかし表10.5.1のように、目的変数が3分類以上のグレードデータつまり順序分類尺度のデータという時もあります。 このようなデータはグレードデータをそのまま計量尺度として扱い、重回帰分析を適用するのが普通です。 しかしグレードデータにロジスティックモデルを当てはめ、ロジスティック回帰分析を適用することもできます。

この場合、まずグレード1を反応なしと考え、グレード2と3を反応ありと考えて、1番目のロジスティックモデルを当てはめます。

ロジスティックモデル1：
π₁：グレード2・3の出現率　　η₁：π₁のロジット 　　β₁₀：モデル1の切片　　β_j：偏回帰係数(j=1,…,p) 　　ε₁：モデル1の回帰誤差

次にグレード1と2を反応なしと考え、グレード3を反応ありと考えて、2番目のロジスティックモデルを当てはめます。 この時、偏回帰係数β_jはモデル1と同じ値であり、切片β₂₀と誤差ε₂だけが異なるという仮定をします。

ロジスティックモデル2：
π₂：グレード3の出現率　　η₂：π₂のロジット 　　β₂₀：モデル2の切片　　β_j：偏回帰係数(j=1,…,p) 　　ε₂：モデル2の回帰誤差

モデル1とモデル2の偏回帰係数が同じということは、2つのモデルのロジット(対数オッズ)の違いは切片β₁₀とβ₂₀の差に影響されるだけで、説明変数には影響されないということです。 その結果、2つのモデルのオッズの間には比例関係があり、その比例定数はβ₁₀とβ₂₀の差を指数変換した値になります。 そのためこのようなモデルのことを比例オッズモデル(POM：proportional odds model)または累積ロジスティックモデル(cumulative logistic model)といいます。 累積ロジスティックモデルのロジスティック曲線は、図10.5.1と図10.5.2のように2つのモデルのロジスティック曲線の立ち上がりの位置が異なるだけで傾きは同じです。 _(注1)

モデル1：
モデル2：
l₂ - l₁=β₂₀ - β₁₀
∴

(2) 順序ロジスティック回帰分析

このように目的変数が3分類以上の順序分類尺度のデータの時に、(分類数-1)個の累積ロジスティックモデルを当てはめ、それらが同時に成り立つようなロジスティック回帰式を求める手法を順序ロジスティック回帰分析(ordinal logistic regression analysis)といいます。 また目的変数が3分類以上の分類尺度のデータの時についてもロジスティック回帰式を求める手法が開発されていて、それを多項ロジスティック回帰分析(multinomial logistic regression analysis)といいます。

多項ロジスティック回帰分析に対して、目的変数が2分類の時の普通のロジスティック回帰分析を二項ロジスティック回帰分析(binomial logistic regression analysis)と呼ぶことがあります。 そして二項ロジスティック回帰分析を後ろ向き研究用にすると(線形)判別分析に対応するのと同様に、多項ロジスティック回帰分析を後ろ向き研究用にすると特殊な重判別分析に対応します。 (→第18章　重判別分析)

目的変数が、例えば「心疾患非発症・心疾患による入院・心疾患による死亡」のように3分類あったとします。 このような場合、入院と死亡をまとめて「反応有」として扱い、二項ロジスティック回帰分析を強引に適用することがあります。 しかし医学的には入院よりも死亡の方が重篤な反応ですから、これらを順序分類尺度扱いして順序ロジスティック回帰分析を適用する方が合理的です。 また入院と死亡に優越を付けず、これらを順序が付けられない分類扱いして多項ロジスティック回帰分析を適用することもできます。

順序ロジスティック回帰分析は、最尤法を利用した繰り返し計算によって解を求めるのが普通です。 表10.5.1のデータに順序ロジスティック回帰分析を適用すると次のようになります。 _(注2)

○モデル1
ロジスティック回帰式：l₁=0.0846 - 0.892x₁ + 0.836x₂ - 0.00344x₃
ロジスティック曲線：
○モデル2
ロジスティック回帰式：l₂=-1.875 - 0.892x₁ + 0.836x₂ - 0.00344x₃
ロジスティック曲線：
○偏回帰係数の検定と推定
χ_β1²=2.236(p=0.1348)＜χ²(1,0.05)=3.841 … 有意水準5％で有意ではない
95％信頼区間：下限 β_1L=-2.062　上限 β_1U=0.277
χ_β2²=4.763(p=0.0291)＞χ²(1,0.05)=3.841 … 有意水準5％で有意
95％信頼区間：下限 β_2L=0.085　上限 β_2U=1.587
χ_β3²=0.016(p=0.8980)＜χ²(1,0.05)=3.841 … 有意水準5％で有意ではない
95％信頼区間：下限 β_3L=-0.0561　上限 β_3U=0.0492

表10.5.1のデータはグレードが3つあるため、切片だけが異なる2つのロジスティック回帰式とロジスティック曲線が求められます。 これらのロジスティック回帰式とロジスティック曲線から、ある説明変数の値を持つ個体の目的変数がどのグレードになる確率が高いか調べる方法を考えてみましょう。 例えば表10.5.1の13番目の個体について、各説明変数の値を2つのロジスティック回帰式に代入してロジットと出現率を求めると次のようになります。

13番目の個体のデータ：重症度=2　x₁=0　x₂=1　x₃=23
○モデル1
ロジスティック回帰式：l₁=0.0846 - 0.892×0 + 0.836×1 - 0.00344×23≒0.841
ロジスティック曲線：
○モデル2
ロジスティック回帰式：l₂=-1.875 - 0.892×0 + 0.836×1 - 0.00344×23≒-1.118
ロジスティック曲線：

モデル1はグレード1だけが反応なしで、グレード2と3は反応ありとしたモデルです。 そのためこのモデルの出現率p₁はこの個体がグレード2または3になる確率になります。 したがってこの個体がグレード1になる確率は1からp₁を引いた値になります。

モデル2はグレード1と2が反応なしで、グレード3が反応ありとしたモデルです。 そのためこのモデルの出現率p₂はこの個体がグレード3になる確率になります。 したがってこの個体がグレード2になる確率は、この個体がグレード2または3になる確率であるp₁からp₂を引いた値になります。

これら3つの確率値から、この個体はグレード2になる確率が一番高いと予測できます。 そしてこの個体のグレードは2ですから、確率値による予測が当たっていることがわかります。 2つのモデルと3つの確率値の関係を模式図にすると図10.5.3のようになります。

しかしこのような確率値によるグレードの予測は非常に面倒でわかりにくく、あまり良い方法とは言えません。 そのため順序ロジスティック回帰分析は説明変数の値に基づいて反応率を予測したり、目的変数の値を予測するための手法というよりも、説明変数が目的変数に対してどの程度の影響を与えているかを分析するための手法と言えるでしょう。 実際、通常の二項ロジスティック回帰分析も順序ロジスティック回帰分析も、ほとんどの場合はリスクファクターの重要度を分析するために利用されます。

比例オッズモデルは2つのモデルのオッズの間に比例関係があり、その比例定数は説明変数には影響されないという前提で組み立てられたモデルです。 この前提が成り立つ可能性はかなり低く、現実問題として目的変数のグレードデータをそのまま計量値扱いできるという可能性とあまり変わりません。 このように順序ロジスティック回帰分析は、やたらと複雑な計算をするわりには得られる結果の信頼性が低く、しかも結果の使い勝手も良くありません。 そのためこのようなデータは、グレードデータをそのまま計量値扱いして重回帰分析を適用する方が簡単でしかも便利です。

参考までに、表10.5.1のデータに重回帰分析を適用すると次のようになります。

目的変数の単位が異なるため、偏回帰係数の値そのものは順序ロジスティック回帰分析の結果と異なります。 しかし順序ロジスティック回帰分析の検定統計量χ²値と、この場合の検定統計量F値を比較すると、説明変数の相対的な重要性はよく似ていることがわかります。 そして目的変数の予測値はほぼ2であり、実際の目的変数の値とよく近似しています。 また重回帰分析では重回帰式の適合性を表す重寄与率を求めることができます。

これらのことから、このようなデータには順序ロジスティック回帰分析よりも重回帰分析を適用した方が良いことがわかると思います。

(注1)　分類数がa個の順序分類尺度の場合、k番目の分類までは反応なし、(k+1)番目以後の分類を反応ありとした、次のような累積確率とそのオッズを考えることができます。

k番目の分類に属す確率：π_k=P(y=k|ｘ) (k=1,…,a)
(k+1)番目以後の分類に属す累積確率：
累積オッズ：

この累積オッズを用いたモデルのことを累積オッズモデル(cumulative odds model)といい、このモデルにロジスティック回帰モデルを適用し、対数累積オッズと説明変数の間に線形関係が成り立つと仮定したモデルのことを累積ロジスティックモデルといいます。 そして累積ロジスティックモデルにおいて、(a-1)個のモデルの切片だけが異なり、偏回帰係数は共通と仮定したものが比例オッズモデルです。

比例オッズモデル：

比例オッズモデルが成り立つ可能性はかなり低いものの、このモデルは計算が比較的簡単で結果の解釈も容易です。 そのため累積ロジスティックモデルといえば、普通はこの比例オッズモデルを意味します。 そこでここでは、累積ロジスティックモデルと比例オッズモデルという用語を厳密には区別せずに使うことにします。

a個の分類が順序尺度ではなく名義尺度の時は、k番目の分類に属す確率と1番目の分類に属す確率の比(π_k/π₁)をオッズとして扱います。 そしてこのオッズにロジスティック回帰モデルを適用したものを多項ロジスティックモデル(multinomial logistic model)といい、次のような式で表されます。

多項ロジスティックモデル：

このモデルは比例オッズモデルに比べて成り立つ可能性は高いものの、計算が煩雑になる上に結果の解釈が複雑になります。 そのため利用頻度は累積ロジスティックモデルよりもさらに低くなります。

(注2)　表10.5.1を一般化し、累積ロジスティックモデルを当てはめると次のようになります。

表10.5.2 説明変数がp個の一般的データ

目的変数y：順序分類尺度	説明変数
	x1	…	xj	…	xp
1	x11	…	x1j	…	x1p
:	:		:		:
1	xn11	…	xn1j	…	xn1p
:	:	…	:	…	:
k	x(∑n+1)1	…	x(∑n+1)j	…	x(∑n+1)p
:	:		:		:
k	x(∑n+nk)1	…	x(∑n+nk)j	…	x(∑n+nk)p
:	:	…	:	…	:
a	x(∑n+1)1	…	x(∑n+1)j	…	x(∑n+1)p
:	:		:		:
a	xn1	…	xnj	…	xnp

(k=1,…,a)
(i=1,…,n に応じて k=1,…,a) 　　 (i=1,…,n) 　　 　　
+∞=β₀₀＞β₀₁＞…＞β_0k＞β_0(a-1)＞β_0a=-∞
k番目の分類に属す確率：π_k=P(y=k|ｘ)
(k+1)番目以後の分類に属す累積確率：
○累積ロジスティックモデル(比例オッズモデル)



E_ik=exp(l_k)=exp(β_0k + ｘ_i'β)=e^{β_0k+ｘ_i'β} 　　E_i0=exp(β₀₀ + ｘ_i'β)=exp(∞ + ｘ_i'β)=∞ 　　E_ia=exp(β_0a + ｘ_i'β)=exp(-∞ + ｘ_i'β)=0
1番目以後の分類に属す累積確率：
(a+1)番目以後の分類に属す累積確率：
π_k=p_k-1 - p_k 　　π₁=p₀ - p₁=1 - p₁ 　　π_a=p_{a - 1}-p_a=p_{a - 1}

最尤法によって切片β₀と偏回帰係数βの最尤推定値ｂ₀とｂを求めるには、次のような計算をします。

全体の尤度： (i=1,…,n に応じて k=1,…,a)
対数尤度：

(k≠1　Σはk=2,…,aについてだけ加算する)
(k≠a　Σはk=1,…,(a-1)についてだけ加算する)
(j=1,…,p　Σはk=1,…,aについて加算する)

(k≠1　Σはk=2,…,aについてだけ加算する)
(k≠a　Σはk=1,…,(a-1)についてだけ加算する)
(k≠1,a　Σはk=2,…,(a-1)についてだけ加算する)
(k≠1　Σはk=2,…,aについてだけ加算する)
(k≠a　Σはk=1,…,(a-1)についてだけ加算する)
(j=1,…,p　Σはk=1,…,aについて加算する)
○k=1の時
　　
　　

○k=2,…,(a-1)の時
　　

　　
　　
　　
○k=aの時



　　

これらの値を用いて、ニュートン・ラプソン法によって切片と偏回帰係数の最尤推定値を求めます。 そして最尤推定値の漸近的正規性を利用して、切片と偏回帰係数についてワルドの検定と推定を行うことができます。

検定：＞χ²(1,α)の時、有意水準100α％で有意
＞χ²(1,α)の時、有意水準100α％で有意
推定：100(1-α)％信頼区間
→ 下限： 　上限：
100(1-α)％信頼区間
→ 下限： 　上限：
[Ｉ_f]_0kk^-1：Ｉ_f^-1=E(-Ｈ^-1)の第0対角要素 　　[Ｉ_f]_jj^-1：Ｉ_f^-1の第j対角要素

切片と偏回帰係数の初期値は、普通のロジスティック回帰分析と同様に判別分析で求めた切片と判別係数を利用します。 その場合、最初の分類を群2とし、2番目以後の分類を群1として判別分析を行います。 これは累積ロジスティックモデルの1番目のモデルに判別分析を適用したものになります。 (→10.1 ロジスティック回帰分析の原理)

また順序ロジスティック回帰分析の切片は複数あり、2番目以後の切片は次第に小さくなります。 そこで判別分析の切片を利用して最初の切片初期値を求め、以後の切片初期値は前の値から定数——例えば1——を引いた値にします。 判別分析の切片と判別係数をa₀とa_jとすると、順序ロジスティック回帰分析の切片と偏回帰係数の初期値は次のようになります。

切片初期値： 　　b_0k=b_0(k-1) - 1 (k=2,…,a-1)
偏回帰係数初期値：b_j=a_j

表10.5.1の例題について実際に計算してみましょう。

分類1(軽症)を群2、分類2(中等症)と3(重症)を群1とした時の判別分析の結果
a₀=-0.37475　　a₁=-0.86694 　　a₂=0.71491　　a₃=-0.011963
切片と偏回帰係数の初期値

b₀₂=0.60608-1=-0.39392　　b₁=-0.86694 　　b₂=0.71491　　b₃=-0.011963
傾斜ベクトルとヘスの行列の要素
E_i1=exp(0.60608-0.86694x_i1+0.71491xi2-0.011963x_i3) 　　E_i2=exp(-0.39392-0.86694x_i1+0.71491xi2-0.011963x_i3)

これらの値を利用してニュートン・ラプソン法を行うと、5回目で繰り返し計算が収束します。

　　 　　


切片と偏回帰係数に関するワルドの検定と推定
＜χ2(1,0.05)=3.841 … 有意水準5％で有意ではない
95％信頼区間=0.0846±2.623 → 下限：β_01L=-2.538　上限：β_01U=2.708
＜χ2(1,0.05)=3.841 … 有意水準5％で有意ではない
95％信頼区間=-1.875±2.690 → 下限：β_02L=-4.565　上限：β_02U=0.815
＜χ2(1,0.05)=3.841 … 有意水準5％で有意ではない
95％信頼区間=-0.892±1.170 → 下限：β_1L=-2.062　上限：β_1U=0.277
オッズ比：OR₁=exp(-0.892)=0.410
オッズ比の95％信頼区間　下限：OR_1L=exp(-2.062)=0.127 　上限：OR_1U=exp(0.277)=1.319
＞χ2(1,0.05)=3.841 … 有意水準5％で有意
95％信頼区間=0.836±0.751 → 下限：β_2L=0.085　上限：β_2U=1.587
オッズ比：OR₂=exp(0.836)=2.307
オッズ比の95％信頼区間　下限：OR_2L=exp(0.085)=1.089 　上限：OR_2U=exp(1.587)=4.887
＜χ2(1,0.05)=3.841 … 有意水準5％で有意ではない
95％信頼区間=-0.00344±0.0526 → 下限：β_3L=-0.0561　上限：β_3U=0.0492
オッズ比：OR₃=exp(-0.00344)=0.997
オッズ比の95％信頼区間　下限：OR_3L=exp(-0.0561)=0.945 　上限：OR_3U=exp(0.0492)=1.050

この累積ロジスティックモデルはk番目までの分類を反応なし、(k+1)番目以後の分類を反応ありとしたものです。 このモデルとは反対にk番目までの分類を反応あり、(k+1)番目以後の分類を反応なしとしたモデルもあります。 その場合の累積ロジスティックモデルは次のようになります。 このモデルに基づいた順序ロジスティック回帰式は、切片の符号が反対になるだけで偏回帰係数は変わりません。

また分類の順序を逆転し、数字が小さいほどグレードが高い、つまり大きな順位になるとしたモデルもあります。 有名な統計ソフトSASはそのようなモデルを採用しています。 そのモデルに基づいた順序ロジスティック回帰式は切片と偏回帰係数の符号が反対になります。

最終更新日：2016年2月25日 　第4節へ 第11章へ

　webmaster@snap-tck.com 　Copyleft (C) 2000 SNAP(Sugimoto Norio Art Production)