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

第10章 ロジスティック回帰分析

この章ではロジスティック回帰分析の原理と結果の解釈方法、判別分析との関係、各種のシグモイド曲線、そして変数選択法について解説します。

10.1 ロジスティック回帰分析の原理

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

ロジスティック回帰分析(logistic regression analysis)は疾患のリスクファクターを分析するためによく用いられる多変量解析です。 この手法は1948年にアメリカのフラミンガムで開始されたフラミンガム研究(Framingham study)のために開発されました。 フラミンガム研究は冠状動脈性疾患に関する大規模なコホート研究であり、複数のリスクファクターつまり多重リスクファクター(multiple risk factor)が疾患に及ぼす影響を分析することを目的のひとつにしています。 そしてそのために開発されたのがこのロジスティック回帰分析であり、現在も主として医学分野で用いられています。

この手法は説明変数が計量尺度のデータで、目的変数が名義尺度を計量尺度化したデータである重回帰分析に相当します。 第5章で説明したように、説明変数が計量尺度で目的変数が名義尺度の時、普通は目的変数を0/1のダミー変数で表し、さらにそれを1の出現率に計量尺度化して回帰直線を求め、その回帰係数の検定としてコクラン・アーミテージの傾向検定を行います。 (→5.3 計数値の相関分析と回帰分析)

しかし第5章の図5.3.3を見てわかるように、この場合の回帰直線は原理的には有り得ない出現率0未満の領域と出現率1以上の領域まで入りこんでしまいます。 そこで説明変数と出現率の関係を直線で回帰せず、第9章の図9.6.1のように、出現率0から1の間でS字状に変化するシグモイド曲線で回帰する方が好ましいと考えられます。 例えば第5章の表5.3.5のデータについて、生後日数分類と尾長8cm以上の出現率を直線で回帰した時とシグモイド曲線で回帰した時を比べると次のようになります。

表5.3.5 マウスの生後日数分類と尾長分類
生後日数分類尾長8cm未満尾長8cm以上8cm以上の出現率
4-5(4.5)2020
6-10(8)3250.4
11-14(12.5)1340.75
全体65110.455
直線で回帰した時:p = -0.355 + 0.0899x
シグモイド曲線で回帰した時:
p:尾長8cm以上の出現率  x:生後日数分類の値
図10.1.1 出現率の回帰直線と回帰曲線

この場合のシグモイド曲線は第9章で紹介したロジスティック曲線を利用しています。 シグモイド曲線には色々なものがありますが、ここでロジスティック曲線を利用したのは、この曲線を利用すると判別分析との整合性が良いからです。 上記のシグモイド曲線において、出現率を対数オッズつまりロジットに変換すると次のようになります。

ロジット(対数オッズ):

これは生後日数分類を説明変数にし、ロジットを目的変数にした直線回帰式に相当します。 そこでこの式を一般化し、さらに説明変数を複数にすると次のような重回帰型のモデルになります。 これを線形ロジスティックモデル(linear logistic model)といい、このモデルに基づいた回帰分析のことをロジスティック回帰分析またはロジット回帰分析(logit regression analysis)といいます。

線形ロジスティックモデル:
η:出現率πのロジット  β0:切片  βj:偏回帰係数(j = 1,…,p)  ε:誤差

(2) 調整オッズ比

ロジスティックモデルは重回帰型ですから、偏回帰係数βj他の説明変数が一定で説明変数xjだけが「1」増加した時にロジットがいくつ変化するかを表す値になります。 そしてロジットは対数オッズのため、ロジットの変化量は対数オッズの変化量になります。 そのため対数オッズを指数変換して元のオッズに戻すと、対数オッズの変化量はオッズが何倍になるかを表す値つまりオッズ比になります。 したがって偏回帰係数βjを指数変換すると、他の説明変数が一定で説明変数xjだけが「1」増加した時のオッズ比になります。 このオッズ比は他の説明変数の影響を取り除いたオッズ比になるため、調整オッズ比(adjusted odds ratio)と呼ばれることがあります。

○xj = 0の時のロジットと出現率
η0 = β0 + β1x1 + … + βj-1xj-1 + βj+1xj+1 + … + βpxp + ε

○xj = 1の時のロジットと出現率
η10 + β1x1 + … + βj-1xj-1 + βj + βj+1xj+1 + … + βpxp + ε

∴η1 - η0 = βj
 (調整オッズ比)

例えば表5.3.5のロジスティック回帰分析の結果から、生後日数分類のオッズ比を求めると次のようになります。

OR = exp(0.445) = 1.560

このオッズ比から、生後日数分類の値つまり生後日数が「1」増加するとオッズが約1.6倍になることがわかります。 しかしこの値は生後日数が1日増加した時のオッズ比であり、生後日数分類が1つ上になった時のオッズ比ではありません。 例えば2番目の生後日数分類の値は8であり、3番目の生後日数分類の値は12.5です。 そのため2番目の生後日数分類と3番目の生後日数分類のオッズ比は、生後日数が4.5増加した時のオッズ比として計算する必要があります。 その場合のオッズ比は、次のように偏回帰係数を4.5倍してから指数変換して計算します。

OR = exp(0.445×4.5) = exp(2.0025) = 7.408

表5.3.5のデータから、2番目の生後日数分類と3番目の生後日数分類の実際のオッズ比を計算すると次のようになります。

2番目の生後日数分類のオッズ=
3番目の生後日数分類のオッズ=
実際のOR=

この値は偏回帰係数から計算した値と少し違います。 この違いは、偏回帰係数から計算した値が生後日数分類全体の平均的なオッズ比に相当することに起因します。 特に表5.3.5の場合、1番目の生後日数分類の出現率が0のため1番目の分類と2番目の分類のオッズ比が非常に大きくなり、平均的なオッズ比を大きくしてしまうのです。

医学分野ではオッズ比をよく用いるので、計量尺度の説明変数を特定の境界値で2分し、「0:境界値未満 1:境界値以上」という名義尺度のデータにしてオッズ比を強引に計算することがあります。 これは第2章で説明した尺度合わせの一種であり、データが持っている情報量を落としてしまい、結果の信頼性が低くなるので原則として避けるべきです。 (→2.6 尺度合わせと外れ値)

例えば表5.3.5の元となった第5章の表5.3.4のデータについて10を境界値にして生後日数を2分類にし、オッズ比を計算すると次のようになります。

表10.1.1 マウスの生後日数分類と尾長分類
生後日数分類尾長8cm未満尾長8cm以上8cm以上の出現率
10日未満(平均値=6.5)5160.17
10日以上(平均値=12)1450.8
全体65110.455

この時のオッズ比は、実は生後日数が10日未満の群の平均値6.5である時と10日以上の群の平均値12である時の値、つまり生後日数が5.5増加した時のオッズ比に相当します。 生後日数が5.5日増加した時のオッズ比をロジスティック回帰分析の偏回帰係数から求めると次のようになります。 この値は実際のオッズ比と少し違いますが、その原因は表5.3.5の2番目の生後日数分類と3番目の生後日数分類の場合と同じです。

生後日数が5.5日増加した時のオッズ比:OR = exp(0.445×5.5) = exp(2.4475) = 11.5594

生後日数を境界値で2分類にした時、実際のデータが境界値の周辺に集中していると、境界値未満の群の平均値と境界値以上の群の平均値の差が小さくなりオッズ比も小さくなります。 それに対して実際のデータが境界値を中心にして広く分布していると、境界値未満の群の平均値と境界値以上の群の平均値の差が大きくなりオッズ比も大きくなります。

同一項目について境界値が同じなら別々の研究でも同じ条件に見え、両者のオッズ比を公平に比較することができるように思えてしまいます。 ところがたとえ境界値が同じでもデータの分布状態によってオッズ比が変わるので、実は公平な比較ではありません。 一方、実測値のままロジスティック回帰分析によってオッズ比を求めれば、それはデータが「1」増加した時のオッズ比になります。 そのためデータの分布状態が多少違ってもオッズ比はあまり変わらず、公平な比較をすることができます。

ただし例えば産婦人科領域における分娩回数のように、初回分娩は色々なリスクが高いものの、2回目はリスクがかなり低くなり、3回目以後はリスクがあまり変わらないという項目が有り得ます。 このような場合は分娩回数を実測値のまま用いるよりも、「0:初産婦 1:経産婦」という2分類データにしてオッズ比を求める方が実態をうまく反映します。

したがって原則として計量尺度のデータは実測値のまま説明変数にし、ロジスティック回帰分析によって求められたオッズ比そのものではなく、説明変数が医学的に意義のある値だけ増加した時のオッズ比を偏回帰係数から求め、それを医学的に検討するのが合理的です。 ただし分娩回数のように特定の境界値の前後で非連続的な変化をする時だけは、特定の境界値で2分してオッズ比を求めるのが合理的です。

それからオッズ比は2種類の値の比なので、出現率の差が同じでも例数が異なるとオッズ比は異なる値になります。 普通の重回帰分析では2種類の説明変数x1とx2が無相関なら、x2を含めた時のx1の偏回帰係数とx2を除外した時のx1の偏回帰係数の値は変わりません。 ところがロジスティック回帰分析ではx1とx2が無相関でも、x2を含めた時のx1の調整オッズ比(指数変換した偏回帰係数)とx2を除外した時のx1の調整オッズ比は値が変わることがあります。

この現象は第4章のロジット検定の解説で詳しく説明したので、詳細はそちらを参照してください。 ここでは第4章の表4.4.6のデータにロジスティック回帰分析を適用した結果と、調整リスク比を求めることができるポアソン回帰分析を適用した結果を示しておきます。 (→4.4 繰り返しのある多標本・多時期の計数値第15章 ポアソン回帰分析)

表4.4.6 薬剤と改善度の性別2×2分割表
(男女別に層別無作為割付けしたもの)
性別薬剤非改善(%)改善(%)計(%)
標準薬25(50)25(50)50(100)
新薬5(10)45(90)50(100)
30(30)70(70)100(100)
標準薬45(90)5(10)50(100)
新薬25(50)25(50)50(100)
70(70)30(30)100(100)
全体標準薬70(70)30(30)100(100)
新薬30(30)70(70)200(100)
100(50)100(50)200(100)
○層別解析
・男
 リスク差(改善率の差):RD = 0.4
 リスク比(改善率の比):RR = 1.8  ln(RR) = 0.587787
 オッズ比:OR = 9  ln(OR) = 2.19722
・女
 リスク差(改善率の差):RD = 0.4
 リスク比(改善率の比):RR = 5  ln(RR) = 1.60944
 オッズ比:OR = 9  ln(OR) = 2.19722
・全体
 リスク差(改善率の差):RD = 0.4
 リスク比(改善率の比):RR = 2.33333  ln(RR) = 0.847298
 オッズ比:OR = 5.44444  ln(OR) = 1.6946
○ポアソン回帰分析(p:改善する確率)
・性別を説明変数に含めた時
 ポアソン回帰式:p = exp(-0.867501 + 0.847298・x1 - 0.847298・x2)
 x1:薬剤(0:標準薬 1:新薬)  調整リスク比:RR1 = exp(0.847298) = 2.33333
 x2:性別(0:男 1:女)  調整リスク比:RR2 = exp(-0.847298) = 0.428571
 ※x1とx2の相関係数:r12 = 0
・性別を説明変数に含めない時
 ポアソン回帰式:p = exp(-1.20397 + 0.847298・x1)
 x1:薬剤(0:標準薬 1:新薬)  調整リスク比:RR1 = exp(0.847298) = 2.33333
○ロジスティック回帰分析(p:改善する確率)
・性別を説明変数に含めた時
 ロジスティック回帰式:l = ln{p/(1-p)} = 2.19722・x1 - 2.19722・x2
 x1:薬剤(0:標準薬 1:新薬)  調整オッズ比:OR1 = exp(2.19722) = 9
 x2:性別(0:男 1:女)  調整オッズ比:OR2 = exp(-2.19722) = 0.111111
 ※x1とx2の相関係数:r12 = 0
・性別を説明変数に含めない時
 ロジスティック回帰式:l = ln{p/(1-p)} = -0.847298 + 1.6946・x1
 x1:薬剤(0:標準薬 1:新薬)  調整オッズ比:OR1 = exp(1.6946) = 5.44444
図4.4.1 男女別と全体のロジスティック曲線

上記の結果から、薬剤と性別が無相関つまり標準薬群と新薬群の男女の割合が同じならリクス差(改善率の差)は性別の影響を受けないことがわかります。 それに対してリスク比(改善率の比)は、層別解析では性別の影響を受けるものの多変量解析(ポアソン回帰分析)では性別の影響を受けず、性別の影響を補正してもしなくても調整リスク比は変わらないことがわかります。 ところがオッズ比は、層別解析でも多変量解析(ロジスティック回帰分析)でも、性別の影響を補正した時の調整オッズ比と性別の影響を補正しない時の調整オッズ比の値が変わります。

これらのことから、たとえ特定の背景因子で層別無作為割り付けしても、また薬剤群の背景因子が揃っていて薬剤と背景因子が無相関でも、オッズ比を指標にした薬剤比較はロジスティック回帰分析を用いて背景因子の影響を補正しなければならないことがわかると思います。

(3) リスク比

出現率が非常に低い時はオッズと出現率が近似するため、オッズ比と出現率比も近似します。 表10.1.1では「尾長8cm以上」というのが「反応有」に相当しますが、反応有が疾患の発症の場合、出現率比はリスク比RR(相対危険度)になります。 そのためオッズ比のことを近似的にリスク比と解釈することがよくあります。 しかしオッズ比のことをリスク比と解釈できるのは出現率がだいたい10%未満の時です。 しかも出現率が低い時はリスク比の値の信頼性が低くなり、相対的な危険性という意味が怪しくなります。

そのため出現率が低い時はオッズ比がリスク比に近似すると解釈するよりも、むしろリスク比がオッズ比に近似し、どちらも単なる関連性の指標にすぎなくなると解釈した方が妥当です。 したがってオッズ比のことを近似的にリスク比と解釈するのは本当は好ましくありません。 (→1.9 科学的研究のデザイン)

例えば表10.1.1の生後日数分類のリスク比を計算すると次のようになります。

生後日数10日未満の出現率=
生後日数10日以上の出現率=
リスク比:

出現率が比較的高いため、実際のデータから計算したオッズ比20も、偏回帰係数から計算したオッズ比11.5594もこの値とはかなり違います。 したがってこのような時はオッズ比をリスク比と解釈するのは難しいことがよくわかると思います。 このような時はロジスティック曲線を利用して正確なリスク比を求めることができます。 例えば表5.3.5から求めたロジスティック曲線を利用して、生後日数が6.5日(生後日数10日未満の平均値)の時と12日(生後日数10日以上の平均値)の時のリスク比を求めると次のようになります。 オッズ比と違い、このリスク比は実際のリスク比に比較的近いことがわかると思います。

生後日数6.5日の時の出現率:
生後日数12日の時の出現率:
リスク比:

表5.3.5の尾長分類は、目的変数が名義尺度の時の回帰分析の説明をするために、あえて計量尺度である尾長のデータを8cmを境界値にして2分類にしています。 しかし目的変数が計量尺度の時は、それをわざわざ名義尺度にせず、計量尺度のまま回帰分析を適用する方が合理的です。 表5.3.5の元になった表5.1.3は生後日数も尾長も計量尺度であり、これに直線回帰分析を適用した結果は次のようなものでした。 (→5.1 相関係数と回帰直線)

表5.1.3 マウスの生後日数と尾長
個体No.生後日数(x)尾長(y)
144.26
255.68
367.24
474.82
586.95
698.81
7108.04
8118.33
91210.84
10137.58
11149.96
平均97.50
回帰直線:y=3 - 0.5x  寄与率:r2=0.667(66.7%)

この回帰直線を利用すると生後日数が特定の値の時の尾長を推測できますし、尾長が8cmになる時の生後日数を逆算することもできます。 さらに回帰直線の予測限界を利用して、尾長が8cm以上になる確率を求めることができます。 そこでその確率を利用してリスク比を求めることができます。 例えば生後日数が6.5日の時と12日の時のリスク比を求めると次のようになります。 (注1)

○尾長が8cmの時の生後日数:8=3 + 0.5x → x=10
○生後日数が6.5日の時
 y=3 + 0.5x=3 + 0.5×6.5=6.25<8
 yの片側予測限界上限が8になる時の信頼係数=0.890442 → yが8以上になる確率=1-0.890442=0.109558(0.109558%)
○生後日数が12日の時
 y=3 + 0.5x=3 + 0.5×12=9>8
 yの片側予測限界下限が8になる時の信頼係数=0.7628625(76.28625%)=yが8以上になる確率
リスク比:

このリスク比はロジスティック曲線を利用して求めたリスク比よりも大きな値です。 ロジスティック曲線と比べて回帰直線は信頼性が高く、回帰直線によって求めた目的変数の値の信頼性も高くなります。 つまりx=6.5の時のy=6.25と、x=12の時のy=9の信頼性が高いので、「yが8以上になるリスク比」も大きな値になるわけです。 ただしリスク比は目的変数が名義尺度の時のラフな指標です。 目的変数が計量尺度の時は、回帰直線を利用して目的変数の値そのものを求めることができます。 したがって通常はラフな指標であるリスク比をわざわざ求める必要はないでしょう。

(4) ロジスティック回帰分析と判別分析の関係

線形ロジスティックモデルと第9章のロジスティック曲線を比べると、切片が異なるだけで同じ形式の式であることがわかります。 試しに表5.3.5のデータで尾長8cm以上を群1に、8cm未満を群2にして判別分析を適用し、群1に属す確率をロジスティック曲線で表すと次のようになります。 (→9.6 ロジスティック曲線)

判別関数:z=-3.770 + 0.412x
ロジスティック曲線:
ロジット変換:
z:判別スコア  x:生後日数分類の値  p:群1に属す確率  π1:群1の事前確率  l:ロジット

ロジット変換の式において、群1の事前確率π1に実際の群1の出現率5/11を代入すると次のようになります。


l=-0.182+z=-3.952 + 0.412x

このように判別分析の結果をロジスティック曲線で表し、群1の事前確率として実際の群1の出現率を代入したものはロジスティック回帰分析の結果とよく似たものになります。 これがシグモイド曲線としてロジスティック曲線を利用する理由であり、これによってロジスティック回帰分析は判別分析の親類筋に当たる手法になります。

ただしロジスティック回帰分析は疾患の発症に影響するリスクファクターを分析し、疾患が発症する前に、ある被験者が疾患を発症するかどうかを予測するための手法です。 そのためこの手法は原則として前向き研究で得られたデータ(目的変数にだけ誤差がある)に適用し、説明変数としてはリスクファクターだけを用います。 そして上式からわかるようにロジスティック回帰式の切片には事前確率の情報は含まれておらず、実際のデータの出現率つまり発症率が切片に反映されます。 (→1.9 科学的研究のデザイン)

それに対して判別分析は診断率の分析を多変量に拡張した手法に相当し、疾患の診断に有用な診断指標を分析し、疾患が発症した後で、ある被験者が疾患であるかどうかを診断するための手法です。 そのためこの手法は原則として後ろ向き研究で得られたデータ(説明変数にだけ誤差がある)に適用し、説明変数としては診断指標が本来ですが、リスクファクターを混合してもかまいません。 そして判別分析におけるロジスティック曲線の切片には事前確率の情報が含まれていて、それがなければ正しい出現率を求められません。 これは診断率の分析において、疾患の一般的な有病率つまり事前確率を使わなければ正しい陽性予測値を求められないのと全く同じ理屈です。

判別分析におけるロジスティック曲線に上記の例のように事前確率として群1の実際の出現率を代入すれば、見かけ上の出現率は計算できます。 しかし後ろ向き研究では群1の例数を任意に決定することができるため、群1の出現率も任意に決定できます。 そのためロジスティック曲線によって計算する出現率も任意の値にすることができてしまいます。 したがってロジスティック曲線を利用して正確な出現率を計算するためには正確な事前確率が必要になります。 (→9.2 群の判別と診断率)

このようにロジスティック回帰分析と判別分析は見かけ上はよく似た手法ですが、目的と適用すべきデータが明確に異なる手法です。 したがって目的とデータに応じて両者をうまく使い分ける必要があります。

判別分析とロジスティック回帰分析の比較として、判別分析は多変量正規分布を前提にしているのに対して、ロジスティック回帰分析はそのような前提を必要としないのでロジスティック回帰分析の方が実際のデータに適している、という説明をしている解説書があります。 しかしロジスティック回帰分析はロジットつまり対数オッズ比と説明変数の間に線形関係があるという前提を必要としています。 これは説明変数が多変量正規分布をするという前提と同じくらい、現実にはかなり無理がある前提です。

そして第9章で説明したように、説明変数が近似的に多変量正規分布をする時、対数オッズ比と説明変数の間に近似的な線形関係が成り立ちます。 つまり2つの手法が正当性を持つためには、全く同じ前提を必要としているのです。 したがって2つの手法の使い分けはデータの正規性といったことではなく、あくまでも分析の目的を主眼にすべきです。 (→9.6 ロジスティック曲線 (注1))

基本的に疾患が発症する前に、疾患が発症するかどうかをリスクファクターから予測するための手法がロジスティック回帰分析であり、疾患が発症した後で、疾患であるかどうかを診断指標から診断するための手法が判別分析です。 したがってこの基本を踏まえた上で、2つの手法をうまく使い分けることが大切です。

(5) 一般化線形モデル

ロジスティック回帰分析において、出現率をロジット変換するのは説明変数と目的変数の関係を直線状つまり線形にするためです。 このように目的変数を線形にするための変換関数のことをリンク関数といいます。 そして適当なリンク関数を利用すれば色々な回帰曲線を線形にすることができ、重回帰型モデルつまり線形モデルにすることができます。 そのようにして一般化した線形モデルのことを一般化線形モデル(GLM:generalized linear model)といいます。

普通の線形モデルの場合、目的変数の回帰誤差が近似的に正規分布することを前提にして回帰分析を行います。 しかし一般化線形モデルの場合、目的変数の回帰誤差が正規分布ではなく特殊な分布になることがよくあります。 そのような場合は第9章で説明した最尤法を利用して回帰分析を行うことがあります。 (→9.3 1変量の場合)

ロジスティック回帰分析の場合もロジットの回帰誤差が特殊な分布になります。 そのため回帰誤差が近似的に正規分布するとみなして重回帰分析と同じように最小2乗法を利用して回帰分析を行う方法と、最尤法を利用した繰り返し近似計算によって回帰分析を行う方法の2種類の計算法があります。 そしてコンピュータが発達した現在では最尤法を利用する方法が主流になっています。

図10.1.1の回帰曲線は最尤法を利用した繰り返し近似計算で求めたロジスティック曲線です。 この回帰曲線の切片と回帰係数は、判別分析から求めたロジスティック曲線の切片と回帰係数とは少し異なっています。 そのため2本のロジスティック曲線を描くと図10.1.2のように微妙にずれます。 この図の「ロジスティック回帰分析のロジスティック曲線」は図10.1.1の回帰曲線と同じものです。

○ロジスティック回帰分析のロジスティック曲線:l=-4.258 + 0.445x

○判別分析のロジスティック曲線:l=-3.952 + 0.412x
図10.1.2 2本のロジスティック曲線

判別分析は説明変数が近似的に多変量正規分布すると仮定して、最尤法を利用して計算します。 これは目的変数を0/1のダミー変数にし、その回帰誤差が近似的に正規分布すると仮定して最小2乗法を利用して計算した結果と原理的に一致します。 そのためロジスティック回帰分析の計算結果とよく似たものになりますが、正確に一致するとは限りません。 (→9.5 変数の選択 (注2))


(注1) 通常、回帰直線の信頼限界や予測限界は両側です。 しかし平均値の信頼区間と同様に、片側信頼限界や片側予測限界を求めることも可能です。 (→5.5 各種手法の相互関係 (注3))

両側予測限界:
片側予測限界: , yU=∞
 または yL=-∞ , 
mx:xの標本平均  Sxx:xの平方和   VR:残差分散
t(n-2,α):自由度(n-2)のt分布における両側100α%点   t(n-2,2α):自由度(n-2)のt分布における片側100α%点
表5.1.3の場合
両側予測限界:
片側予測限界:
 または

これらの式を利用してx=6.5の時の片側予測限界上限が8になる時の信頼係数と、x=12の時の片側予測限界下限が8になる時の信頼係数を求めてみましょう。

○x=6.5
 片側予測限界上限 yU=8の時:
 t(9,2α)=3.6819の時の片側確率 α=0.219116÷2=0.109558
 ∴信頼係数(1-α)=1-0.109558=0.890442(89.0442%)
○x=12
 片側予測限界下限 yL=8の時:
 t(9,2α)=0.8467401の時の片側確率 α=0.474275÷2=0.2371375
 ∴信頼係数(1-α)=1-0.2371375=0.7628625(76.28625%)