왈드 검정과 가능도비 검정¶

import numpy as np
import statsmodels.api as sm
from scipy import stats

# 설명변수 셋 중 마지막 하나는 반응변수와 무관하게 만든다
rng = np.random.default_rng(0)
n = 500
X = rng.normal(0, 1, size=(n, 3))
logit = -0.5 + 1.2 * X[:, 0] - 0.8 * X[:, 1] + 0.0 * X[:, 2]
y = (rng.random(n) < 1 / (1 + np.exp(-logit))).astype(int)

# 전체 모형을 적합한다. 세 번째 계수는 참으로 0 이다.
X_full = sm.add_constant(X)
model_full = sm.Logit(y, X_full).fit(disp=0)


def print_summary(res):
    """summary()의 Date/Time 칸은 실행할 때마다 달라지므로 비우고 출력한다."""
    lines = []
    for line in str(res.summary()).split("\n"):
        if line.startswith(("Date:", "Time:")):
            lines.append(line[:19].ljust(38) + line[38:])
        else:
            lines.append(line)
    print("\n".join(lines))


# statsmodels 의 출력표는 기본으로 Wald 검정을 보여 준다. 계수를 그
# 표준오차로 나눈 값이다.
print_summary(model_full)

# 가능도비 검정은 변수를 뺀 모형과 넣은 모형의 로그가능도를 견준다.
# Wald 는 한 모형만 적합하면 되지만, 가능도비는 두 번 적합해야 한다.
# 대신 표본이 작거나 계수가 클 때 가능도비 쪽이 더 믿을 만하다.
model_restricted = sm.Logit(y, X_full[:, :-1]).fit(disp=0)
lr_stat = -2 * (model_restricted.llf - model_full.llf)
p_value = stats.chi2.sf(lr_stat, df=1)

# z^2 과 가능도비 통계량이 비슷하게 나오는지 본다. 표본이 크면 둘이
# 같은 값으로 다가간다(점근적으로 동등하다).
wald_z = model_full.tvalues[-1]
print(f"Wald z = {wald_z:.4f}, z^2 = {wald_z**2:.4f}, p = {model_full.pvalues[-1]:.4f}")
print(f"LRT  = {lr_stat:.4f}, p = {p_value:.4f}")

# 요약표의 LLR 을 두 수로 되살린다 (자유도 3).
ll_full, ll_null = model_full.llf, model_full.llnull
big = -2 * (ll_null - ll_full)
print(f"\n-2*({ll_null:.2f} - ({ll_full:.2f})) = {big:.4f}, df = 3, "
      f"p = {stats.chi2.sf(big, 3):.3e}")
print(f"statsmodels 가 적은 값: LLR = {model_full.llr:.4f}, "
      f"p = {model_full.llr_pvalue:.3e}")

# 네 자리의 Lambda / W^2.  뒤의 셋은 연습문제 2 의 모의실험 결과다.
print("\n  beta_hat      W^2       Lambda   Lambda/W^2")
rows = [(model_full.params[-1], wald_z ** 2, lr_stat),
        (0.966, 3.393 ** 2, 14.16),
        (3.033, 4.780 ** 2, 58.06),
        (6.780, 4.100 ** 2, 90.22)]
for b, w2, lam in rows:
    print(f"{b:10.4f} {w2:10.4f} {lam:10.4f} {lam / w2:10.4f}")

                           Logit Regression Results                           
==============================================================================
Dep. Variable:                      y   No. Observations:                  500
Model:                          Logit   Df Residuals:                      496
Method:                           MLE   Df Model:                            3
Date:                                   Pseudo R-squ.:                  0.2078
Time:                                   Log-Likelihood:                -269.48
converged:                       True   LL-Null:                       -340.15
Covariance Type:            nonrobust   LLR p-value:                 1.943e-30
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
const         -0.2445      0.105     -2.329      0.020      -0.450      -0.039
x1             1.2057      0.137      8.776      0.000       0.936       1.475
x2            -0.7111      0.118     -6.016      0.000      -0.943      -0.479
x3            -0.1287      0.103     -1.249      0.212      -0.331       0.073
==============================================================================
Wald z = -1.2493, z^2 = 1.5608, p = 0.2115
LRT  = 1.5719, p = 0.2099

-2*(-340.15 - (-269.48)) = 141.3402, df = 3, p = 1.943e-30
statsmodels 가 적은 값: LLR = 141.3402, p = 1.943e-30

  beta_hat      W^2       Lambda   Lambda/W^2
   -0.1287     1.5608     1.5719     1.0071
    0.9660    11.5124    14.1600     1.2300
    3.0330    22.8484    58.0600     2.5411
    6.7800    16.8100    90.2200     5.3670

import numpy as np
import statsmodels.api as sm
from scipy import stats

rng = np.random.default_rng(7)
n = 100
x = rng.normal(size=n)

for c in [1.0, 3.0, 6.0]:
    p = 1 / (1 + np.exp(-c * x))
    y = (rng.random(n) < p).astype(float)
    X = sm.add_constant(x)
    m = sm.Logit(y, X).fit(disp=0)
    m0 = sm.Logit(y, np.ones((n, 1))).fit(disp=0)
    lr = -2 * (m0.llf - m.llf)
    print(f"c={c}: beta={m.params[1]:.3f} se={m.bse[1]:.3f} "
          f"z={m.tvalues[1]:.3f} wald_p={m.pvalues[1]:.3g} "
          f"LR={lr:.2f} LR_p={stats.chi2.sf(lr, 1):.3g}")

c=1.0: beta=0.966 se=0.285 z=3.393 wald_p=0.000692 LR=14.16 LR_p=0.000168
c=3.0: beta=3.033 se=0.635 z=4.780 wald_p=1.75e-06 LR=58.06 LR_p=2.54e-14
c=6.0: beta=6.780 se=1.654 z=4.100 wald_p=4.13e-05 LR=90.22 LR_p=2.13e-21