『医学のための因果推論』のための統計学入門：11.2値データの回帰モデル
『医学のための因果推論Ⅰ 一般化線形モデル』を勉強した時のまとめです。
結構難しい本でしたので自分なりの解釈を入れながらまとめました。
2値データの回帰モデルについて詳しく解説します。
具体的にはプロビット回帰、ロジスティック回帰、積2項分布モデルとその応用（用量反応関数の推定、判別分析、ランダム化臨床試験の解析）について説明します。
モデルの構造
2値データは個々の観測値が「0」または「1」のいずれかをとるデータを指します。
例えばある薬を投与した患者が「治った（1）」か「治らなかった（0）」かを表すデータを2値データといいます。
対象となる個人 $${i}$$ のアウトカム（結果）を $${Y_i}$$ とし、その値が0または1をとるとします。
このとき、$${Y_i}$$ は確率 $${\pi_i = \operatorname{Pr}(Y_i = 1)}$$ で「1」をとる独立なベルヌーイ型の確率変数と考えられます。
N 人 の対象者の同時分布を考えると、これは各 $${Y_i}$$ の確率質量関数を掛け合わせたものです。
$$
P(\boldsymbol{y} \mid \boldsymbol{\pi}) = \prod_{i=1}^N \pi_i^{y_i}(1 - \pi_i)^{1 - y_i}
$$
指数型分布族の形式に合わせるため、対数を取って右辺を展開します。
$$
\log P(\boldsymbol{y} \mid \boldsymbol{\pi}) = \sum_{i=1}^N \left[ y_i \log(\pi_i) + (1 - y_i) \log(1 - \pi_i) \right]
$$
この式を整理します。
$$
\begin{aligned} \log P(\boldsymbol{y} \mid \boldsymbol{\pi}) &= \sum_{i=1}^N \left[ y_i \log(\pi_i) + (1 - y_i) \log(1 - \pi_i) \right] \\ &= \sum_{i=1}^N \left[ y_i \log(\pi_i) + \log(1 - \pi_i) - y_i \log(1 - \pi_i) \right] \\ &= \sum_{i=1}^N \left[ y_i \left( \log(\pi_i) - \log(1 - \pi_i) \right) + \log(1 - \pi_i) \right] \\ &= \sum_{i=1}^N \left[ y_i \log\left( \frac{\pi_i}{1 - \pi_i} \right) + \log(1 - \pi_i) \right] \end{aligned}
$$
以上をまとめると、$${N}$$ 人 の対象者の同時分布は以下のように表されます。
$$
\prod_{i=1}^N \pi_i^{y_i}(1 - \pi_i)^{1 - y_i} = \exp\left[\sum_{i=1}^N y_i \log\left(\frac{\pi_i}{1 - \pi_i}\right) + \sum_{i=1}^N \log(1 - \pi_i)\right]
$$
この式は指数型分布族の一つです。
指数型分布族は多くの確率分布を統一的に扱うための重要な概念であり、一般的に以下の形式で表されます。
$$
p(y; \theta) = \exp\left[ a(y) b(\theta) + c(\theta) + d(y) \right]
$$
ここで
$${a(y)}$$：十分統計量（sufficient statistic）
$$
a(y)=y
$$
$${b(\theta)}$$：自然パラメータ（natural parameter）
$$
b(\theta) = \log \left( \frac{\pi}{1 - \pi} \right)
$$
$${c(\theta)}$$：対数分配関数（log-partition function）の負の値
$$
A(\theta) = - c(\theta) = - \log (1 - \pi)
$$
対数分配関数は正規化を行い、分布の特性を決定します。
$${d(y)}$$：規定測度（base measure）の対数 $${h(y) = \exp[ d(y) ] = 1}$$
一般的な指数型分布族の表現にはいくつかの形式がありますが、この形式では対数分配関数 $${A(\theta)}$$ と $${c(\theta)}$$ の関係が $${A(\theta) = - c(\theta)}$$ となります。
指数型分布族は統計学における重要な分布のクラスで、解析が容易な性質を持っています。
一般化線型モデル（Generalized Linear Model, GLM）では、リンク関数と呼ばれる関数 $${g(\cdot)}$$ を用いて確率パラメータ $${\pi_i}$$ と共変量（説明変数） $${X_1, X_2, \ldots, X_p}$$ の関係を次のようにモデル化します。
$$
g(\pi_i) = \boldsymbol{X_i}^\top \boldsymbol{\beta}
$$
ここで、$${\boldsymbol{X_i}}$$ は共変量のベクトル、$${\boldsymbol{\beta}}$$ は回帰係数のベクトルです。
推定
このモデルはパラメトリックモデルの一種であり、最尤法を用いてパラメータ $${\boldsymbol{\beta}}$$ を推定できます。
しかしこの場合、データから必ずしも最尤推定量が求まるとは限りません。
特に、計算上の問題として完全分離（complete separation）という現象が知られています。
これは、共変量とアウトカムが完全に分離されており、モデルがデータを100%の精度で分類できてしまう状況を指します。
例えば、ある共変量 $${X}$$ に対して、$${X}$$ がある値以上のときは必ず $${Y = 1}$$、それ以外のときは $${Y = 0}$$ となる場合です。
最尤法は、観測データが得られる確率（尤度）が最大になるようにパラメータを選びます。
例としてロジスティック回帰(後述)においてロジットリンク関数を使い、次のように尤度を表現します。
$$
\log P(\boldsymbol{y} \mid \boldsymbol{\pi}) = \sum_{i=1}^N \left[ y_i \log(\pi_i) + (1 - y_i) \log(1 - \pi_i) \right]
$$
ここで、$${{\pi_i = \frac{1}{1 + \exp(-(\boldsymbol{X_i}^\top \boldsymbol{\beta}))} }}$$ です。
目的は、この尤度を最大化する $${\boldsymbol{\beta}}$$ を求めることです。
完全分離の状況では、観測データが「共変量の値がある範囲に入ると必ず $${Y=1}$$ になる」といった構造を持つため、$${Y=1}$$ を予測する確率 $${\pi_i}$$ が完全に1に近づくようにパラメータ $${\boldsymbol{\beta}}$$ を選ぶことができるようになります。
このような場合、以下の問題が発生します。
パラメータが無限大に発散する
尤度を最大化するために、モデルは$${\pi_i}$$をできる限り「1または0に近づける」方向に調整されます。
この結果、回帰係数 $${\boldsymbol{\beta}}$$ のいくつかの要素が無限大に発散してしまいます。
言い換えると、ある閾値を越えると確実に $${Y=1}$$ となる状況では、尤度を最大化するために$${\beta}$$の一部が発散しても尤度が増加し続けます。最尤推定量が存在しない
パラメータが無限大に発散するため、尤度を最大にする有限のパラメータ値が存在しません。
このため最尤推定量が存在せず、通常の方法でパラメータを適切に推定することができなくなります。
完全分離の状況下では尤度が最大となる特定の有限の値には収束せず、推定が不可能な状態に陥ります。
ロジスティック回帰モデルの係数 β を増加させるごとに予測確率 $${\pi_i}$$ の曲線が急峻になり、βが大きくなるにつれて $${\pi_i}$$ が1または0に収束する様子を視覚化。
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LogisticRegression
# 1. 完全分離のデータを生成
np.random.seed(0)
n_samples = 50
X0 = np.random.normal(2, 0.5, n_samples) # Y=0のデータ
X1 = np.random.normal(5, 0.5, n_samples) # Y=1のデータ
X = np.concatenate([X0, X1]).reshape(-1, 1)
Y = np.concatenate([np.zeros(n_samples), np.ones(n_samples)])
# 2. 異なる値の係数βでロジスティック回帰モデルを学習し、発散の様子を確認
# βの値を手動で設定
beta_values = [0.5, 1, 5, 10, 20, 50]
colors = plt.cm.viridis(np.linspace(0, 1, len(beta_values))) # カラーマップで色分け
# プロットの準備
plt.figure(figsize=(10, 6))
plt.scatter(X0, np.zeros(n_samples), color='blue', label='Y=0')
plt.scatter(X1, np.ones(n_samples), color='red', label='Y=1')
plt.axvline(x=3.5, color='grey', linestyle='--', label='Separation threshold')
# 異なるβのロジスティック回帰モデルの予測曲線
x_range = np.linspace(0, 8, 300).reshape(-1, 1)
for beta, color in zip(beta_values, colors):
# ロジスティック回帰モデルをβに比例するようにシミュレーション
model = LogisticRegression(C=1e6, solver='lbfgs') # 高いC値でペナルティを無効化
model.coef_ = np.array([[beta]]) # βを直接指定
model.intercept_ = np.array([-beta * 3.5]) # 分離閾値にあわせた切片
model.classes_ = np.array([0, 1]) # 必要なクラスを指定
# 予測確率の計算
prob = model.predict_proba(x_range)[:, 1]
plt.plot(x_range, prob, color=color, label=f'β = {beta}')
# グラフの設定
plt.xlabel('Covariate X')
plt.ylabel('Predicted Probability P(Y=1)')
plt.title('Coefficient Divergence in Complete Separation Scenario')
plt.legend()
plt.grid(True)
plt.show()完全分離が発生した場合は以下のような対処法が取られることが一般的です。
ペナルティ付き尤度（正則化）
完全分離の影響を軽減するために、リッジ回帰やラッソ回帰のような正則化を導入します。
例えばリッジ回帰は回帰係数の大きさにペナルティを課し、無限大に発散することを防ぎます。
これにより、適切な範囲でパラメータの推定が可能になります。
データ再調整
完全分離の原因となる共変量の数を減らす、あるいはデータを再サンプリングすることで完全分離を避けることができます。
モデルの複雑さを低減させることで、分離を回避する場合もあります。
ベイズ推定
ベイズ推定を用いることで、事前分布に基づく推定が可能です。
正則化を行うと同様に、分離状態に対して適切な推定が可能になる場合があります。
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LogisticRegression
# 1. リッジ正則化の効果を示す図
np.random.seed(0)
n_samples = 50
X0 = np.random.normal(2, 0.5, n_samples) # Y=0のデータ
X1 = np.random.normal(5, 0.5, n_samples) # Y=1のデータ
X = np.concatenate([X0, X1]).reshape(-1, 1)
Y = np.concatenate([np.zeros(n_samples), np.ones(n_samples)])
# 正則化パラメータ C の異なる値でロジスティック回帰モデルを学習
C_values = [1e6, 1, 0.1, 0.01]
colors = plt.cm.viridis(np.linspace(0, 1, len(C_values)))
plt.figure(figsize=(12, 6))
# 元のデータの完全分離状況
plt.subplot(1, 2, 1)
plt.scatter(X0, np.zeros(n_samples), color='blue', label='Y=0')
plt.scatter(X1, np.ones(n_samples), color='red', label='Y=1')
plt.axvline(x=3.5, color='grey', linestyle='--', label='Separation threshold')
# 異なるCのロジスティック回帰モデルの予測曲線
x_range = np.linspace(0, 8, 300).reshape(-1, 1)
for C, color in zip(C_values, colors):
model = LogisticRegression(C=C, solver='lbfgs')
model.fit(X, Y)
prob = model.predict_proba(x_range)[:, 1]
plt.plot(x_range, prob, color=color, label=f'C = {C}')
# グラフの設定
plt.xlabel('Covariate X')
plt.ylabel('Predicted Probability P(Y=1)')
plt.title('Figure 1: Effect of Ridge Regularization on Complete Separation')
plt.legend()
plt.grid(True)用量反応関係
用量反応関係とは、ある薬剤や毒物の用量 $${X}$$ に対して、特定の反応（例えば死亡）の確率 $${\pi}$$ がどのように変化するかを示す関係です。
歴史的に、バイオアッセイの分野で2値データの回帰モデルが用いられてきました。
バイオアッセイでは、異なる用量の物質を用いて生物の反応（例えば死亡率）を調べます。
用量反応関数のモデル化
用量反応関数は、用量 $${X}$$ を変化させたときに、反応確率 $${\pi}$$ が0から1の範囲をとる必要があります。
これを満たすために、以下のようなモデルを用います。
$$
\pi = \int_{-\infty}^X p(y) \, dy
$$
ここで、$${p(y)}$$ は何らかの確率密度関数です。
プロビット関数
初期のバイオアッセイでは、$${p(y)}$$ に正規分布を用いたプロビット関数が使われました。
$$
\pi = \Phi\left(\frac{X - \mu}{\sigma}\right)
$$
ここで、$${\Phi}$$ は正規分布の累積分布関数、$${\mu}$$ は平均、$${\sigma}$$ は標準偏差です。
このモデルは以下のように一般化線型モデルとして表現できます。
$$
\Phi^{-1}(\pi) = \beta_0 + \beta_1 X
$$
ここで$${\beta_0 = -\mu/\sigma, \beta_1 = 1/\sigma}$$ です。
ロジット関数とロジスティック回帰
プロビット関数に代わって現在広く用いられているのがロジット関数です。
$$
g(x) = \log\left(\frac{x}{1 - x}\right)
$$
このときのモデルはロジスティック回帰と呼ばれます。
$$
\log\left(\frac{\pi}{1 - \pi}\right) = \beta_0 + \beta_1 X
$$
ロジット関数の特徴
2項分布の正準リンク関数である。
S字型の曲線を描き、0から1の範囲をとる。
横軸は用量 $${X}$$、縦軸は反応確率 $${\pi}$$。 曲線は $${\exp(\beta_1)}$$ の値を変化させたときのロジット関数を示している。
$${\beta_1}$$ の値が大きいほど、曲線の傾きが急になる。
import numpy as np
import matplotlib.pyplot as plt
def logistic_function(x, beta0, beta1):
return 1 / (1 + np.exp(-(beta0 + beta1 * x)))
x = np.linspace(-10, 10, 400)
beta0 = 0
beta1_values = [0.5, 1, 2, 5]
plt.figure(figsize=(8, 6))
for beta1 in beta1_values:
pi = logistic_function(x, beta0, beta1)
plt.plot(x, pi, label=f'β₁ = {beta1}')
plt.xlabel('Dose X')
plt.ylabel('Response Probability π')
plt.title('Figure 12-1: Example of Logit Function')
plt.legend()
plt.grid(True)
plt.show()ロジット関数とプロビット関数の比較
2値データの回帰モデルとしてプロビット関数とロジット関数はどちらもS字型の曲線を描き、0から1の範囲で反応確率を表現できるため、似た性質を持っています。
しかし現代ではプロビット関数よりもロジット関数が一般的に使用される傾向があります。
その理由には以下のような要因が挙げられます。
計算の簡便さ
ロジット関数は、確率パラメータ $${{\pi}}$$ と共変量の関係を直線的にモデル化するために、リンク関数として単純で使いやすい性質を持っています。
ロジット関数の導関数や計算上の複雑さがプロビット関数よりも低く、計算が効率的です。
また、ロジット関数は正準リンク関数であるため、計算上も統計的にも解析が容易になります。
データの分布に依存しない適用性
プロビット関数は正規分布を前提としており、反応の閾値が正規分布に従う場合には最適ですが、実際のデータが必ずしも正規分布に従うとは限りません。
一方でロジット関数はデータの分布に依存せず、正規分布以外のデータにも広く適用できます。
このためロジット関数はより汎用的で、多様なデータに対応する柔軟性があります。
解釈のしやすさ
ロジット関数を用いた場合、ロジスティック回帰の回帰係数は「オッズ比」という形で解釈されます。
オッズ比は特に臨床研究や疫学で広く利用されており、「特定の因子が1単位増加するごとに、結果が発生する確率が何倍になるか」を明確にします。
この直感的な解釈のしやすさがロジット関数が現代で広く採用される理由の一つです。
判別分析
ロジット関数の導出
判別分析は、データから個体がどのグループに属するかを判別する方法です。
2つのグループ（例えば、陰性と陽性）からのデータ $${X_i}$$ を用いて、その個体がどちらのグループに属するかを推定します。
判別分析において、ある個体がグループ $${Y_i = 1}$$ に属する確率 $${\operatorname{Pr}(Y_i = 1 \mid X_i = x)}$$ と、グループ $${Y_i = 0}$$ に属する確率 $${\operatorname{Pr}(Y_i = 0 \mid X_i = x)}$$ を考えます。
各グループの事前確率を次のように定義します。
$$
\operatorname{Pr}(Y_i = 0) = 1 - \pi, \quad \operatorname{Pr}(Y_i = 1) = \pi
$$
条件付き確率密度関数を $${p(x \mid Y_i)}$$ とします。
条件付き確率を $${\operatorname{Pr}(Y_i = k \mid X_i = x)}$$ のように表します。
尤度比の定義
各グループに属する確率の比を対数で表現するため、まずは尤度比 $${l(x)}$$ を次のように定義します。
$$
l(x) = \log\left[\frac{p(x \mid Y_i = 1)}{p(x \mid Y_i = 0)}\right]
$$
ベイズの定理を適用
ベイズの定理を用いて、$${\operatorname{Pr}(Y_i = 1 \mid X_i = x)}$$ と $${\operatorname{Pr}(Y_i = 0 \mid X_i = x)}$$を次のように表します。
$$
\operatorname{Pr}(Y_i = 1 \mid X_i = x) = \frac{\operatorname{Pr}(Y_i = 1) \cdot p(x \mid Y_i = 1)}{p(x)}
$$
$$
\operatorname{Pr}(Y_i = 0 \mid X_i = x) = \frac{\operatorname{Pr}(Y_i = 0) \cdot p(x \mid Y_i = 0)}{p(x)}
$$
ここで、$${p(x)}$$ は $${X_i = x}$$ が観測される周辺確率で、以下のように定義されます。
$$
p(x) = \operatorname{Pr}(Y_i = 1) \cdot p(x \mid Y_i = 1) + \operatorname{Pr}(Y_i = 0) \cdot p(x \mid Y_i = 0)
$$
対数オッズ比の計算
ベイズの定理の結果を用いて、対数オッズ比を計算します。
$$
\log\left(\frac{\operatorname{Pr}(Y_i = 1 \mid X_i = x)}{\operatorname{Pr}(Y_i = 0 \mid X_i = x)}\right) = \log\left(\frac{\operatorname{Pr}(Y_i = 1) \cdot p(x \mid Y_i = 1)}{\operatorname{Pr}(Y_i = 0) \cdot p(x \mid Y_i = 0)}\right)
$$
事前確率と尤度比を代入
ここで、事前確率 $${\operatorname{Pr}(Y_i = 1) = \pi}$$ と $${\operatorname{Pr}(Y_i = 0) = 1 - \pi}$$ を代入し、$${p(x \mid Y_i = 1) / p(x \mid Y_i = 0)}$$ を $${\exp(l(x))}$$ に置き換えます。
$$
\log\left(\frac{\operatorname{Pr}(Y_i = 1 \mid X_i = x)}{\operatorname{Pr}(Y_i = 0 \mid X_i = x)}\right) = \log\left(\frac{\pi}{1 - \pi}\right) + \log\left(\frac{p(x \mid Y_i = 1)}{p(x \mid Y_i = 0)}\right)
$$
よって、
$$
\log\left(\frac{\operatorname{Pr}(Y_i = 1 \mid X_i = x)}{\operatorname{Pr}(Y_i = 0 \mid X_i = x)}\right) = \log\left(\frac{\pi}{1 - \pi}\right) + l(x)
$$
もし $${X_i}$$ が正規分布に従うと仮定すると、$${p(x \mid Y_i)}$$ は正規分布の確率密度関数で表されるため、尤度比 $${l(x)}$$ は $${X}$$ の1次または2次の関数として表現されます。
この場合、モデルは線形または二次の判別分析として表現することができます。
ROC曲線とC統計量
感度と特異度を以下のように定めます。
感度：実際に陽性である個体を正しく陽性と判別する確率。
特異度：実際に陰性である個体を正しく陰性と判別する確率。
カットオフ値 $${c}$$ を設定して、$${X_i \geq c}$$ なら陽性と判別します。
ROC曲線（Receiver Operating Characteristic曲線）は、カットオフ値を変化させたときの感度と 1− 特異度をプロットしたもので、曲線が左上に近づくほど判別性能が高いとされています。
ROC曲線の下の面積を表し、判別性能の総合的な指標としてC統計量があります。
$${C=1}$$ なら完璧な判別、$${C=0.5}$$ ならランダムな判別です。
C統計量はデータの順位情報にのみ依存すること、Wilcoxon-Mann-WhitneyのU統計量と関連があることに注意が必要です。
ロジスティック回帰モデルで、オッズ比を変化させたときのROC曲線。
オッズ比が大きいほど、ROC曲線が左上に近づき、C統計量が大きくなる。
重要なポイント： オッズ比が20のとき、C統計量は約0.93となり、高い判別性能を示す。
疫学研究では、オッズ比が1.5～5程度のことが多く、その場合C統計量はそれほど高くならない。
import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import roc_curve, auc
np.random.seed(0)
n_samples = 1000
X = np.random.normal(0, 1, n_samples).reshape(-1, 1)
beta0 = 0
beta1_values = [0.2, 0.5, 1, 2]
plt.figure(figsize=(8, 6))
for beta1 in beta1_values:
logits = beta0 + beta1 * X
probabilities = 1 / (1 + np.exp(-logits))
Y = np.random.binomial(1, probabilities).ravel()
model = LogisticRegression(solver='lbfgs')
model.fit(X, Y)
Y_score = model.decision_function(X)
fpr, tpr, thresholds = roc_curve(Y, Y_score)
roc_auc = auc(fpr, tpr)
plt.plot(fpr, tpr, label=f'Odds Ratio = {np.exp(beta1):.2f}, AUC = {roc_auc:.2f}')
plt.plot([0, 1], [0, 1], 'k--')
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.05])
plt.xlabel('1 - Specificity')
plt.ylabel('Sensitivity')
plt.title('ROC Curve of Logistic Regression Model with Different Odds Ratios')
plt.legend(loc='lower right')
plt.grid(True)
plt.show()$${2 \times 2}$$表と積2項分布
リスク差、リスク比、オッズ比
臨床試験や疫学研究では、2つのグループ間での2値アウトカムの比較が重要です。
効果の指標としては以下を用いることが一般的です。
リスク差：$${\pi^1 - \pi^0}$$
リスク比：$${\frac{\pi^1}{\pi^0}}$$
オッズ比：$${\frac{\pi^1 / (1 - \pi^1)}{\pi^0 / (1 - \pi^0)}}$$
ここで、$${\pi^0}$$ はコントロール群のイベント発生確率、$${\pi^1}$$ は処置群のイベント発生確率です。
積2項分布モデル
コントロール群と処置群のイベント数は、それぞれ独立な2項分布に従うと仮定。
全体の対数尤度は、各群の2項分布の対数尤度の和として表される。
リンク関数と効果指標の関係
恒等リンク（$${g(x)=x}$$）：リスク差を直接モデル化。
対数リンク（$${g(x) = \log(x)}$$）：リスク比をモデル化。
ロジットリンク（$${g(x) = \log\left(\frac{x}{1 - x}\right)}$$）：オッズ比をモデル化。
対数尤度の再表現
パラメータの表現を変えることで、対数尤度関数を回帰係数 $${\boldsymbol{\beta}}$$ の関数として表現できる。
これにより、一般化線型モデルとして最尤推定を行うことが可能。
横軸と縦軸は回帰係数 $${\beta_0β}$$ と $${\beta_1}$$。
等高線は対数尤度の値が同じ点を結んだもの。
楕円形の等高線は、$${\beta_0}$$ と $${\beta_1}$$ の間に負の相関があることを示す。
import numpy as np
import matplotlib.pyplot as plt
x = np.linspace(-6, 6, 400)
def threshold_function(x, threshold):
return np.where(x >= threshold, 1, 0)
def gompertz_function(x, a, b):
return np.exp(-a * np.exp(-b * x))
threshold = 0
a = 5
b = 1
y_logistic = 1 / (1 + np.exp(-x))
y_probit = norm.cdf(x)
y_threshold = threshold_function(x, threshold)
y_gompertz = gompertz_function(x, a, b)
plt.figure(figsize=(8, 6))
plt.plot(x, y_logistic, label='Logistic Function')
plt.plot(x, y_probit, label='Probit Function', linestyle='--')
plt.plot(x, y_threshold, label='Threshold Model', linestyle='-.')
plt.plot(x, y_gompertz, label='Gompertz Function', linestyle=':')
plt.xlabel('Dose X')
plt.ylabel('Response Probability π')
plt.title('Comparison of S-Shaped Dose-Response Curves')
plt.legend()
plt.grid(True)
plt.show()参考文献
『医学のための因果推論Ⅰ 一般化線形モデル』