ロジスティック回帰モデル

def logistic(eta: np.ndarray) -> np.ndarray:
    """ロジスティック関数（数値安定版）"""
    return np.where(eta >= 0,
                    1 / (1 + np.exp(-eta)),
                    np.exp(eta) / (1 + np.exp(eta)))


def logistic_regression_irls(
    X: np.ndarray,
    y: np.ndarray,
    max_iter: int = 25,
    tol: float = 1e-8,
) -> dict:
    """
    IRLSによるロジスティック回帰の最尤推定

    Parameters
    ----------
    X : np.ndarray, shape (n, p)
        計画行列（切片列を含む）
    y : np.ndarray, shape (n,)
        目的変数（0/1）
    max_iter : int
        最大反復回数
    tol : float
        収束判定の閾値

    Returns
    -------
    dict
        beta, se, log_likelihood, n_iter
    """
    n, p = X.shape
    beta = np.zeros(p)

    for iteration in range(max_iter):
        eta = X @ beta
        pi = logistic(eta)

        # 重み行列の対角要素
        w = pi * (1 - pi)
        w = np.clip(w, 1e-10, None)  # 数値安定性

        # 作業従属変数
        z = eta + (y - pi) / w

        # 重み付き最小二乗法
        W = np.diag(w)
        XtWX = X.T @ W @ X
        XtWz = X.T @ W @ z

        beta_new = np.linalg.solve(XtWX, XtWz)

        # 収束判定
        if np.max(np.abs(beta_new - beta)) < tol:
            beta = beta_new
            break
        beta = beta_new

    # 最終結果
    eta = X @ beta
    pi = logistic(eta)
    log_lik = np.sum(y * np.log(pi + 1e-15) + (1 - y) * np.log(1 - pi + 1e-15))

    # 標準誤差
    w = pi * (1 - pi)
    XtWX = X.T @ np.diag(w) @ X
    cov = np.linalg.inv(XtWX)
    se = np.sqrt(np.diag(cov))

    return {
        "beta": beta,
        "se": se,
        "log_likelihood": log_lik,
        "n_iter": iteration + 1,
        "cov": cov,
    }

import pandas as pd
import statsmodels.api as sm

# データ生成
rng = np.random.default_rng(42)
n = 500
x1 = rng.normal(0, 1, n)
x2 = rng.normal(0, 1, n)
beta_true = np.array([-0.5, 1.2, -0.8])  # beta_0, beta_1, beta_2
X = np.column_stack([np.ones(n), x1, x2])
eta = X @ beta_true
pi_true = logistic(eta)
y = rng.binomial(1, pi_true)

print(f"サンプルサイズ: {n}")
print(f"真の係数: β₀={beta_true[0]}, β₁={beta_true[1]}, β₂={beta_true[2]}")
print(f"y=1 の割合: {y.mean():.3f}")

# スクラッチ実装
result_irls = logistic_regression_irls(X, y)

# statsmodels
model_sm = sm.GLM(y, X, family=sm.families.Binomial())
result_sm = model_sm.fit()

# 比較
comparison = pd.DataFrame({
    "真の値": beta_true,
    "IRLS (自前)": result_irls["beta"],
    "statsmodels": result_sm.params,
    "SE (自前)": result_irls["se"],
    "SE (statsmodels)": result_sm.bse,
}, index=["β₀", "β₁", "β₂"])

print(f"IRLS 反復回数: {result_irls['n_iter']}")
print(f"対数尤度 (自前): {result_irls['log_likelihood']:.4f}")
print(f"対数尤度 (statsmodels): {result_sm.llf:.4f}")
print()
comparison

print(result_sm.summary())

# オッズ比と95%信頼区間
from scipy import stats

beta_hat = result_sm.params
se = result_sm.bse
z_crit = stats.norm.ppf(0.975)

odds_ratio_df = pd.DataFrame({
    "β": beta_hat,
    "OR (exp(β))": np.exp(beta_hat),
    "OR 95% CI lower": np.exp(beta_hat - z_crit * se),
    "OR 95% CI upper": np.exp(beta_hat + z_crit * se),
    "p-value": result_sm.pvalues,
}, index=["intercept", "x₁", "x₂"])

odds_ratio_df

# statsmodels 組み込みの Titanic データ
titanic = sm.datasets.get_rdataset("Titanic", "datasets").data

# 集計データを個票データに展開
rows = []
for _, row in titanic.iterrows():
    for _ in range(int(row["Freq"])):
        rows.append({
            "Class": row["Class"],
            "Sex": row["Sex"],
            "Age": row["Age"],
            "Survived": 1 if row["Survived"] == "Yes" else 0,
        })
df = pd.DataFrame(rows)
print(f"サンプルサイズ: {len(df)}")
print(f"生存率: {df['Survived'].mean():.3f}")
print()
df.groupby(["Class", "Sex"])["Survived"].mean().unstack()

from statsmodels.formula.api import logit

# ロジスティック回帰モデルの推定
model_titanic = logit("Survived ~ C(Class) + C(Sex) + C(Age)", data=df).fit(disp=0)
print(model_titanic.summary())

# オッズ比の表示
or_df = pd.DataFrame({
    "β": model_titanic.params,
    "OR": np.exp(model_titanic.params),
    "OR 95%CI lower": np.exp(model_titanic.conf_int()[0]),
    "OR 95%CI upper": np.exp(model_titanic.conf_int()[1]),
    "p-value": model_titanic.pvalues,
})
or_df