玄関雑学の部屋雑学コーナー統計学入門

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

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

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

表10.5.1 目的変数が順序分類尺度のデータ
No.重症度
1:軽症 2:中等症 3:重症
説明変数
x1x2x3
110121
210130
310137
410246
511124
611156
711158
811224
911238
1011258
1111326
1211341
1320123
1420143
1520147
1620235
1720241
1820245
1920253
2020340
2121122
2221139
2321152
2421223
2521228
2621232
2721243
2821324
2921327
3021342
3130120
3230144
3330235
3430237
3530341
3630355
3730351
3830336
3931134
4031242
4131251
4231321
4331335
4431336

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

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

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

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

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

モデル1:
モデル2:
l2 - l120 - β10
図10.5.1 累積ロジスティックモデル1 図10.5.2 累積ロジスティックモデル2

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

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

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

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

○モデル1
ロジスティック回帰式:l1=0.0846 - 0.892x1 + 0.836x2 - 0.00344x3
ロジスティック曲線:
○モデル2
ロジスティック回帰式:l2=-1.875 - 0.892x1 + 0.836x2 - 0.00344x3
ロジスティック曲線:
○偏回帰係数の検定と推定
χβ12=2.236(p=0.1348)<χ2(1,0.05)=3.841 … 有意水準5%で有意ではない
95%信頼区間:下限 β1L=-2.062 上限 β1U=0.277
χβ22=4.763(p=0.0291)>χ2(1,0.05)=3.841 … 有意水準5%で有意
95%信頼区間:下限 β2L=0.085 上限 β2U=1.587
χβ32=0.016(p=0.8980)<χ2(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 x1=0 x2=1 x3=23
○モデル1
ロジスティック回帰式:l1=0.0846 - 0.892×0 + 0.836×1 - 0.00344×23≒0.841
ロジスティック曲線:
○モデル2
ロジスティック回帰式:l2=-1.875 - 0.892×0 + 0.836×1 - 0.00344×23≒-1.118
ロジスティック曲線:

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

グレード1になる確率=1 - p1=1 - 0.699=0.301

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

グレード2になる確率=p1 - p2=0.699 - 0.246=0.453
グレード3になる確率=p2=0.246
グレード1になる確率+グレード2になる確率+グレード3になる確率=(1 - p1) + (p1 - p2) + p2=1

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

図10.5.3 各グレードになる確率

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

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

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

○重回帰式:y=1.659 - 0.347x1 + 0.316x2 - 0.00111x3
重寄与率:R2=0.140(14%)  重相関係数:R=0.374
○偏回帰係数の検定と推定
Fβ1=2.288(p=0.1382)<F(1,41,0.05)=4.079 … 有意水準5%で有意ではない
95%信頼区間:下限 β1L=-0.810 上限 β1L=0.116
Fβ2=4.926(p=0.0322)>F(1,41,0.05)=4.079 … 有意水準5%で有意
95%信頼区間:下限 β2L=0.028 上限 β2L=0.603
Fβ3=0.0115(p=0.9153)<F(1,41,0.05)=4.079 … 有意水準5%で有意ではない
95%信頼区間:下限 β3L=-0.022 上限 β3L=0.020
○13番目の個体の目的変数予測値:y=1.659 - 0.347×0 + 0.316×1 - 0.00111×23≒1.949

目的変数の単位が異なるため、偏回帰係数の値そのものは順序ロジスティック回帰分析の結果と異なります。 しかし順序ロジスティック回帰分析の検定統計量χ2値と、この場合の検定統計量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番目の分類に属す確率の比(πk1)をオッズとして扱います。 そしてこのオッズにロジスティック回帰モデルを適用したものを多項ロジスティックモデル(multinomial logistic model)といい、次のような式で表されます。

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

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

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

表10.5.2 説明変数がp個の一般的データ
目的変数y:順序分類尺度説明変数
x1xjxp
1x11x1jx1p
::::
1xn11xn1jxn1p
::::
kx(∑n+1)1x(∑n+1)jx(∑n+1)p
::::
kx(∑n+nk)1x(∑n+nk)jx(∑n+nk)p
::::
ax(∑n+1)1x(∑n+1)jx(∑n+1)p
::::
axn1xnjxnp
(k=1,…,a)
(i=1,…,n に応じて k=1,…,a)    (i=1,…,n)      
+∞=β00>β01>…>β0k>β0(a-1)>β0a=-∞
k番目の分類に属す確率:πk=P(y=k|)
(k+1)番目以後の分類に属す累積確率:
○累積ロジスティックモデル(比例オッズモデル)



Eik=exp(lk)=exp(β0k + i'β)=eβ0k+i'β   Ei0=exp(β00 + i'β)=exp(∞ + i'β)=∞   Eia=exp(β0a + i'β)=exp(-∞ + i'β)=0
1番目以後の分類に属す累積確率:
(a+1)番目以後の分類に属す累積確率:
πk=pk-1 - pk   π1=p0 - p1=1 - p1   πa=pa - 1-pa=pa - 1

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

全体の尤度: (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の時



  

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

検定:>χ2(1,α)の時、有意水準100α%で有意
>χ2(1,α)の時、有意水準100α%で有意
推定:100(1-α)%信頼区間
→ 下限:  上限:
100(1-α)%信頼区間
→ 下限:  上限:
[f]0kk-1f-1=E(--1)の第0対角要素   [f]jj-1f-1の第j対角要素

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

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

定数初期値:   b0k=b0(k-1) - 1 (k=2,…,a-1)
偏回帰係数初期値:bj=aj

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

分類1(軽症)を群2、分類2(中等症)と3(重症)を群1とした時の判別分析の結果
a0=-0.37475  a1=-0.86694   a2=0.71491  a3=-0.011963
定数と偏回帰係数の初期値

b02=0.60608-1=-0.39392  b1=-0.86694   b2=0.71491  b3=-0.011963
傾斜ベクトルとヘスの行列の要素
Ei1=exp(0.60608-0.86694xi1+0.71491xi2-0.011963xi3)   Ei2=exp(-0.39392-0.86694xi1+0.71491xi2-0.011963xi3)


















これらの値を利用してニュートン・ラプソン法を行うと、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
オッズ比:OR1=exp(-0.892)=0.410
オッズ比の95%信頼区間 下限:OR1L=exp(-2.062)=0.127  上限:OR1U=exp(0.277)=1.319
>χ2(1,0.05)=3.841 … 有意水準5%で有意
95%信頼区間=0.836±0.751 → 下限:β2L=0.085 上限:β2U=1.587
オッズ比:OR2=exp(0.836)=2.307
オッズ比の95%信頼区間 下限:OR2L=exp(0.085)=1.089  上限:OR2U=exp(1.587)=4.887
<χ2(1,0.05)=3.841 … 有意水準5%で有意ではない
95%信頼区間=-0.00344±0.0526 → 下限:β3L=-0.0561 上限:β3U=0.0492
オッズ比:OR3=exp(-0.00344)=0.997
オッズ比の95%信頼区間 下限:OR3L=exp(-0.0561)=0.945  上限:OR3U=exp(0.0492)=1.050

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

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