Logistic Regression: The Model

1 이 포스트의 위치

EIV 시리즈(포스트 147–150)가 측정오차 모형을 다루었다면, 이 포스트는 Casella & Berger(2002) Ch.12의 두 번째 주제 §12.3 — 로지스틱 회귀(Logistic Regression) 를 다룬다.

핵심 질문: 반응 변수가 0/1인 베르누이 데이터에서 선형 회귀가 아닌 GLM이 필요한 이유는 무엇이며, 로지스틱 링크는 어떻게 이 문제를 해결하는가?

이 시리즈의 선행 포스트 GLM 프레임워크가 개관을 제공한다. 이 포스트는 수학적 전개와 추정론에 집중한다.

---

2 왜 선형 회귀가 아닌가

\(Y_i \in \{0, 1\}\) 인 이진 반응을 \(\mu_i = E Y_i = P(Y_i=1) = \pi_i\) 로 연결하는 단순한 방법은 선형 확률 모형(Linear Probability Model)이다:

\[\pi_i = \alpha + \beta x_i\]

그러나 이것은 두 가지 근본적인 문제를 갖는다:

1. 범위 위반: \(\pi_i\) 는 확률이므로 반드시 \(0 \leq \pi_i \leq 1\) 이어야 한다. 그러나 \(\alpha + \beta x_i\) 는 임의의 실수 값을 가질 수 있어 예측값이 \([0,1]\) 을 벗어난다.
2. 이분산성 무시: \(Y_i \sim \text{Bernoulli}(\pi_i)\) 이면 \(\text{Var}(Y_i) = \pi_i(1-\pi_i)\) 이다. 이 분산은 \(\pi_i\) 에 따라 달라지므로, OLS의 등분산 가정이 성립하지 않는다.

해결책: \(\pi_i\) 를 직접 선형 예측자로 모형화하는 대신, \(\pi_i\) 의 적절한 함수(링크)를 선형으로 모형화한다.

---

3 GLM의 세 구성요소

로지스틱 회귀를 이해하려면 먼저 GLM(Generalized Linear Model)의 구조를 파악해야 한다. GLM은 세 구성요소로 이루어진다 (Casella & Berger, 2002, §12.3.1).

3.1 랜덤 구성요소 (Random Component)

반응 변수 \(Y_1, \ldots, Y_n\) 이 지수족(exponential family) 분포를 따르는 독립 확률변수라고 가정한다.

• 동일 분포일 필요는 없다: 각 \(Y_i\) 는 같은 분포족에 속하지만 모수 \(\pi_i\) 가 다를 수 있다
• 로지스틱 회귀에서 랜덤 구성요소: \(Y_i \sim \text{Bernoulli}(\pi_i)\)

베르누이 PMF를 지수족 형태로 쓰면:

\[\pi^y (1-\pi)^{1-y} = (1-\pi)\exp\!\left\{y\log\!\left(\frac{\pi}{1-\pi}\right)\right\}\]

자연 모수(natural parameter)는 \(\eta = \log\!\left(\pi/(1-\pi)\right)\), 즉 로그 오즈(log-odds)이다.

3.2 선형 구성요소 (Systematic Component)

예측변수 \(x_i\) 의 선형 결합이 선형 예측자 \(\eta_i\) 를 구성한다:

\[\eta_i = \alpha + \beta x_i\]

이것은 모형의 구조 방정식이다. 회귀계수 \(\alpha, \beta\) 가 이 선형 예측자를 결정한다.

3.3 링크 함수 (Link Function)

링크 함수 \(g(\mu_i)\) 는 랜덤 구성요소의 평균 \(\mu_i = E Y_i\) 와 선형 예측자 \(\eta_i\) 를 연결한다:

\[g(\mu_i) = \eta_i = \alpha + \beta x_i\]

로지스틱 회귀의 링크 함수: 로짓(logit) 함수

\[g(\pi) = \log\!\left(\frac{\pi}{1-\pi}\right)\]

이 링크는 베르누이 분포의 정준 링크(canonical link)이다. 정준 링크란 자연 모수가 직접 선형 예측자가 되는 링크로, 수리적으로 가장 편리한 성질을 갖는다.

GLM 비교표

모형	랜덤 구성요소	링크	\(g(\mu)\)
선형 회귀	\(N(\mu, \sigma^2)\)	항등(identity)	\(\mu\)
로지스틱 회귀	\(\text{Bernoulli}(\pi)\)	로짓(logit)	\(\log\!\frac{\pi}{1-\pi}\)
포아송 회귀	\(\text{Poisson}(\lambda)\)	로그(log)	\(\log\lambda\)
프로빗 회귀	\(\text{Bernoulli}(\pi)\)	프로빗(probit)	\(\Phi^{-1}(\pi)\)

---

4 로지스틱 회귀 모형

4.1 모형 정의

정의: 로지스틱 회귀 모형 (Casella & Berger, 2002, 식 12.3.1)

\(Y_1, \ldots, Y_n\) 이 독립이고 \(Y_i \sim \text{Bernoulli}(\pi_i)\) 일 때, \(\pi_i\) 와 예측변수 \(x_i\) 의 관계를

\[\log\!\left(\frac{\pi_i}{1-\pi_i}\right) = \alpha + \beta x_i\]

로 가정하는 모형이다.

좌변 \(\log(\pi/(1-\pi))\) 는 성공의 로그 오즈(log-odds) 또는 로짓이라 부른다.

이 모형은 역함수를 취해 \(\pi_i\) 에 대해 직접 쓸 수 있다 (Casella & Berger, 2002, 식 12.3.2):

\[\pi(x) = \frac{e^{\alpha + \beta x}}{1 + e^{\alpha + \beta x}}\]

또는 동치적으로:

\[\pi(x) = \frac{1}{1 + e^{-(\alpha + \beta x)}}\]

이 함수를 로지스틱 함수(logistic function) 또는 시그모이드(sigmoid)라 한다.

4.2 로지스틱 함수와 로지스틱 분포

\(F(w) = e^w / (1 + e^w)\) 는 로지스틱(0,1) 분포의 CDF이다. 따라서 로지스틱 회귀는

\[\pi(x) = F(\alpha + \beta x)\]

로도 쓸 수 있다. 이 관점은 다른 CDF로 교체하여 대안 모형을 만들 수 있음을 시사한다:

• \(F = \Phi\) (표준정규 CDF) → 프로빗 회귀(probit regression)
• \(F =\) Gumbel CDF → 로그-로그 링크(log-log link)

4.3 로지스틱 함수의 수학적 성질

도함수: 로지스틱 함수의 도함수는 아름다운 형태를 갖는다 (Casella & Berger, 2002, 식 12.3.3):

\[\frac{d\pi(x)}{dx} = \beta\,\pi(x)(1-\pi(x))\]

이 공식이 말하는 것:

• \(\pi(x)(1-\pi(x)) > 0\) 이므로 \(d\pi/dx\) 의 부호는 \(\beta\) 의 부호와 같다
• \(\beta > 0\): \(\pi(x)\) 는 \(x\) 의 단조 증가 함수 — \(x\) 가 커질수록 성공 확률이 올라간다
• \(\beta < 0\): \(\pi(x)\) 는 \(x\) 의 단조 감소 함수
• \(\beta = 0\): \(\pi(x) = e^\alpha/(1+e^\alpha)\) 로 상수 — \(x\) 와 \(\pi\) 사이에 관계가 없다 (SLR에서 \(\beta=0\) 과 동일한 의미)

핵심 직관: \(\pi(x)(1-\pi(x))\) 는 베르누이 분포의 분산 \(\text{Var}(Y|x) = \pi(x)(1-\pi(x))\) 이다. 도함수 공식은 “로지스틱 곡선의 기울기가 해당 점에서의 분산에 비례한다”는 것을 말한다. \(\pi = 0.5\) 일 때 분산이 최대(= 0.25)이므로 곡선이 가장 가파르고, \(\pi \to 0\) 또는 \(\pi \to 1\) 로 갈수록 분산이 0에 가까워져 곡선이 평평해진다. 이것이 S자형의 수학적 이유이다.

변곡점과 대칭성: 로지스틱 함수는 \(x = -\alpha/\beta\) 에서 \(\pi(-\alpha/\beta) = 1/2\) 이다.

또한 이 점을 중심으로 대칭이다:

\[\pi\!\left(\frac{-\alpha}{\beta} + c\right) = 1 - \pi\!\left(\frac{-\alpha}{\beta} - c\right)\]

직관: 로지스틱 곡선은 “성공 확률 1/2”인 점을 축으로 S자 모양으로 대칭이다.

---

5 모수의 해석: α와 β

5.1 α의 해석

(12.3.1)에서 \(x=0\) 을 대입하면:

\[\log\!\left(\frac{\pi(0)}{1-\pi(0)}\right) = \alpha\]

따라서 \(\alpha\) 는 \(x=0\) 에서의 로그 오즈이다. 이것은 선형 회귀의 절편(\(x=0\) 에서의 평균)에 해당한다.

5.2 β의 해석: 로그 오즈비

\(x\) 와 \(x+1\) 에서의 로그 오즈 차이를 계산하면 (Casella & Berger, 2002, §12.3.1):

\[\log\!\left(\frac{\pi(x+1)}{1-\pi(x+1)}\right) - \log\!\left(\frac{\pi(x)}{1-\pi(x)}\right) = (\alpha + \beta(x+1)) - (\alpha + \beta x) = \beta\]

따라서 \(\beta\) 는 \(x\) 가 1 단위 증가할 때의 로그 오즈 변화량이다.

이것을 지수화하면 (Casella & Berger, 2002, 식 12.3.4):

\[e^\beta = \frac{\pi(x+1)/(1-\pi(x+1))}{\pi(x)/(1-\pi(x))}\]

우변은 \(x+1\) 에서의 오즈를 \(x\) 에서의 오즈로 나눈 오즈비(odds ratio)이다. 로지스틱 회귀에서 이 오즈비는 \(x\) 에 관계없이 일정(constant)하다.

그러므로 (Casella & Berger, 2002, 식 12.3.5):

\[\frac{\pi(x+1)}{1-\pi(x+1)} = e^\beta \cdot \frac{\pi(x)}{1-\pi(x)}\]

즉, \(e^\beta\) 는 \(x\) 가 1 단위 증가할 때 오즈의 곱셈적 변화이다.

해석 요약:

모수	해석	단위
\(\alpha\)	\(x=0\) 에서의 로그 오즈	로그 오즈
\(\beta\)	\(x\) 1단위 증가 시 로그 오즈 변화	로그 오즈/단위
\(e^\beta\)	\(x\) 1단위 증가 시 오즈의 배수 변화	오즈비 (단위 없음)

비교: 선형 회귀에서 \(\beta\) 는 “\(x\) 1단위 증가 시 \(E[Y]\) 의 덧셈적 변화”이다. 로지스틱 회귀에서 \(\beta\) 는 “\(x\) 1단위 증가 시 오즈의 곱셈적 변화 \(e^\beta\)” 를 의미한다.

---

6 최대우도 추정 (MLE)

6.1 왜 최소제곱법을 쓸 수 없는가

선형 회귀에서는 \(Y_i = \alpha + \beta x_i + \varepsilon_i\) 로 쓸 수 있어 잔차를 직접 정의하고 최소화할 수 있었다. 로지스틱 회귀에서는 \(Y_i\) 와 \(\alpha + \beta x_i\) 사이에 직접적인 덧셈적 연결이 없다 — 링크 함수가 개입하기 때문이다. 따라서 최소제곱법은 적용 불가하다.

6.2 베르누이 우도 함수

단일 베르누이 관측의 경우 (각 \(x_i\) 에서 한 번씩 관측):

\[L(\alpha, \beta | \mathbf{y}) = \prod_{i=1}^n \pi(x_i)^{y_i} (1-\pi(x_i))^{1-y_i}\]

\(F_i = F(\alpha + \beta x_i) = \pi(x_i)\) 로 쓰면:

\[L(\alpha, \beta | \mathbf{y}) = \prod_{i=1}^n F_i^{y_i}(1-F_i)^{1-y_i}\]

로그 우도:

\[\log L(\alpha, \beta | \mathbf{y}) = \sum_{i=1}^n \left\{\log(1-F_i) + y_i \log\!\left(\frac{F_i}{1-F_i}\right)\right\}\]

6.3 스코어 방정식 (MLE 조건)

\(F(w)\) 의 도함수를 \(f(w) = dF(w)/dw\) 로 쓰고, \(f_i = f(\alpha + \beta x_i)\) 로 정의한다.

\(\partial/\partial\alpha\) 방향으로 미분하면 (Casella & Berger, 2002, 식 12.3.7):

\[\frac{\partial \log L}{\partial \alpha} = \sum_{i=1}^n (y_i - F_i)\frac{f_i}{F_i(1-F_i)} = 0\]

\(\partial/\partial\beta\) 방향 (Casella & Berger, 2002, 식 12.3.8):

\[\frac{\partial \log L}{\partial \beta} = \sum_{i=1}^n (y_i - F_i)\frac{f_i}{F_i(1-F_i)} x_i = 0\]

로지스틱 링크의 특수성: \(F(w) = e^w/(1+e^w)\) 이면 \(f(w) = F(w)(1-F(w))\) 이므로

\[\frac{f_i}{F_i(1-F_i)} = 1\]

스코어 방정식이 대폭 단순해진다:

\[\sum_{i=1}^n (y_i - \hat{\pi}_i) = 0, \qquad \sum_{i=1}^n (y_i - \hat{\pi}_i) x_i = 0\]

이것은 잔차 \((y_i - \hat{\pi}_i)\) 가 1과 \(x_i\) 에 대해 직교해야 함을 요구한다 — 선형 회귀의 정규방정식과 유사한 구조이다.

6.4 다중 관측: 이항 케이스

\(J\) 개의 고유한 \(x\) 값이 있고, \(x_j\) 에서 \(n_j\) 번 관측하여 \(Y_j^* \sim \text{Binomial}(n_j, \pi(x_j))\) 인 경우:

\[L(\alpha, \beta | \mathbf{y}^*) = \prod_{j=1}^J F_j^{y_j^*}(1-F_j)^{n_j - y_j^*}\]

스코어 방정식:

\[\sum_{j=1}^J (y_j^* - n_j F_j)\frac{f_j}{F_j(1-F_j)} = 0, \qquad \sum_{j=1}^J (y_j^* - n_j F_j)\frac{f_j}{F_j(1-F_j)} x_j = 0\]

6.5 수치 최적화

스코어 방정식은 \(\alpha, \beta\) 에 대해 비선형 방정식이므로 닫힌 형태의 해가 없다. 수치 방법이 필요하다:

• Newton-Raphson: 헤시안(2차 도함수 행렬)을 이용한 반복 최적화
• IRLS (Iteratively Reweighted Least Squares): 가중 최소제곱 반복으로 동치 표현

로그 우도의 순오목성(strict concavity): 로지스틱·프로빗 회귀에서 로그 우도는 \((\alpha, \beta)\) 에 대해 순오목이다. 따라서 스코어 방정식의 해가 존재하면 유일하고 그것이 전역 최댓값(MLE)이다.

6.6 완전 분리 (Complete Separation)

예외 상황: 일부 극단적 데이터에서 스코어 방정식의 해가 존재하지 않는다. 이는 데이터가 완전히 분리되어(어떤 임계값으로 \(Y=0\) 과 \(Y=1\) 을 완벽하게 구분 가능) \(\pi(x) \to 0\) 또는 1 로 가는 극한에서 우도가 최대화되기 때문이다.

직관: 로지스틱 모형은 \(0 < \pi(x) < 1\) 을 전제한다. 데이터가 완전 분리되면 “완벽한 예측”이 가능해져 \(\hat{\beta} \to \pm\infty\) 로 발산한다. 이 데이터가 로지스틱 모형에서 나왔다면, 이런 완전 분리가 발생할 확률은 0으로 수렴한다.

---

7 챌린저 데이터 예시

우주왕복선 챌린저(Challenger)호의 O-링 실패와 발사 온도의 관계 데이터 (Dalal et al., 1989)를 Casella & Berger(2002)가 §12.3.1 예시로 사용한다.

• 반응 변수: O-링 실패 여부 (\(Y_i = 1\): 실패, \(Y_i = 0\): 정상)
• 예측변수: 발사 시점의 기온 (화씨, \(x_i\))
• 데이터: 23회 비행 기록 (Table 12.3.1)

MLE: \(\hat{\alpha} = 15.043\), \(\hat{\beta} = -0.232\)

\(\hat{\beta} < 0\) 이므로 기온이 낮을수록 O-링 실패 확률이 올라간다. 1986년 1월 28일 챌린저 발사 당시 기온은 31°F였다. 모형에서의 MLE 예측:

\[\hat{\pi}(31) = \frac{e^{15.043 + (-0.232)(31)}}{1 + e^{15.043 + (-0.232)(31)}} \approx 0.9996\]

사고 당일 O-링 실패 확률이 99.96%로 추정된다.

오즈비 해석: \(e^{\hat{\beta}} = e^{-0.232} \approx 0.793\) 이다. 기온이 1°F 상승하면 O-링 실패 오즈가 0.793배 — 약 21% 감소한다. 반대로 기온이 10°F 낮아지면 오즈가 \(0.793^{-10} \approx 10.3\) 배가 된다. 31°F는 데이터 범위의 평균 기온(약 70°F)보다 39°F 낮으므로, 이 예측은 데이터 범위를 크게 벗어난 극단적 외삽(extrapolation)이다. 오즈비의 일정성 가정이 이 범위에서 성립하는지 별도로 검증할 필요가 있다.

---

8 피셔 정보행렬과 근사 분포

8.1 2모수 정보행렬

단일 모수의 경우 정보 수(Fisher information number)를 썼다면, 이제 두 모수 \((\alpha, \beta)\) 에 대한 피셔 정보행렬(Fisher information matrix)이 필요하다 (Casella & Berger, 2002, 식 12.3.9):

\[I(\theta_1, \theta_2) = \begin{pmatrix} -\frac{\partial^2}{\partial\theta_1^2}\log L & -\frac{\partial^2}{\partial\theta_1\partial\theta_2}\log L \\ -\frac{\partial^2}{\partial\theta_1\partial\theta_2}\log L & -\frac{\partial^2}{\partial\theta_2^2}\log L \end{pmatrix}\]

8.2 로지스틱 회귀의 정보행렬

이항 케이스에서의 정보행렬 (Casella & Berger, 2002, 식 12.3.10):

\[I(\alpha, \beta) = \begin{pmatrix} \sum_{j=1}^J n_j F_j(1-F_j) & \sum_{j=1}^J x_j n_j F_j(1-F_j) \\ \sum_{j=1}^J x_j n_j F_j(1-F_j) & \sum_{j=1}^J x_j^2 n_j F_j(1-F_j) \end{pmatrix}\]

이 행렬의 원소들이 \(Y_j^*\) 에 의존하지 않는다는 것을 주목하라 — 관측 정보행렬 = 기대 정보행렬이 된다.

가중치 \(w_j = n_j F_j(1-F_j)\) 로 쓰면 정보행렬이 아래처럼 보인다:

\[I(\alpha, \beta) = \begin{pmatrix} \sum_j w_j & \sum_j x_j w_j \\ \sum_j x_j w_j & \sum_j x_j^2 w_j \end{pmatrix}\]

이것은 가중 최소제곱(WLS)에서의 \(\mathbf{X}^T\mathbf{W}\mathbf{X}\) 와 같은 구조이다. IRLS와의 연결이 여기서 나타난다.

8.3 MLE의 근사 분포와 신뢰구간

MLE 점근 이론에 의해:

\[\begin{pmatrix} \hat{\alpha} \\ \hat{\beta} \end{pmatrix} \xrightarrow{d} N\!\left(\begin{pmatrix} \alpha \\ \beta \end{pmatrix},\, I(\alpha,\beta)^{-1}\right)\]

분산을 추정하려면 \(I(\hat{\alpha}, \hat{\beta})^{-1}\) 를 사용한다. 2×2 역행렬:

\[\begin{pmatrix} a & b \\ b & d \end{pmatrix}^{-1} = \frac{1}{ad-b^2}\begin{pmatrix} d & -b \\ -b & a \end{pmatrix}\]

\([\text{se}(\hat{\alpha})]^2\) 와 \([\text{se}(\hat{\beta})]^2\) 는 역행렬의 대각 원소이다.

\(\beta\) 에 대한 근사 신뢰구간:

\[\hat{\beta} \pm z_{\alpha/2}\,\text{se}(\hat{\beta})\]

챌린저 예시:

\[I(\hat{\alpha}, \hat{\beta})^{-1} = \begin{pmatrix} 54.44 & -0.80 \\ -0.80 & 0.012 \end{pmatrix}\]

\(\text{se}(\hat{\beta}) = \sqrt{0.012} \approx 0.110\) 이므로 \(\beta\) 의 95% 근사 신뢰구간:

\[-0.232 \pm 1.96 \times 0.110 \quad \Rightarrow \quad -0.447 \leq \beta \leq -0.017\]

신뢰구간이 0을 포함하지 않으므로 \(\beta < 0\) 이 지지된다.

수치 해석: \(\text{se}(\hat{\alpha}) = \sqrt{54.44} \approx 7.38\) 로 매우 크다. \(\alpha\) 는 \(x = 0\) (0°F)에서의 로그 오즈이므로, 데이터 범위(53–81°F)에서 멀리 외삽된 추정값이라 불확실성이 클 수밖에 없다. \(\text{se}(\hat{\beta}) = \sqrt{0.012} \approx 0.110\) 은 상대적으로 작아서 기울기 추정이 안정적이다. 비대각 원소 \(-0.80\) 은 \(\hat{\alpha}\) 와 \(\hat{\beta}\) 의 음의 상관을 반영한다 — 기울기가 더 가파르면(더 음수이면) 절편이 더 큰 양수를 취하는 교환 관계가 데이터에 존재한다.

---

9 가설 검정

9.1 H₀: β = 0 (예측변수의 유효성)

선형 회귀에서와 마찬가지로, \(\beta = 0\) 은 예측변수와 반응변수 사이에 관계가 없음을 의미한다.

Wald 검정: 검정 통계량

\[Z = \frac{\hat{\beta}}{\text{se}(\hat{\beta})}\]

은 \(H_0\) 하에 근사적으로 \(N(0,1)\) 을 따른다. 유의수준 \(\alpha\) 에서 \(|Z| \geq z_{\alpha/2}\) 이면 \(H_0\) 를 기각한다.

우도비 검정(LRT): 검정 통계량

\[-2\log\lambda(\mathbf{y}^*) = 2\left[\log L(\hat{\alpha}, \hat{\beta} | \mathbf{y}^*) - \log L(\hat{\alpha}_0, 0 | \mathbf{y}^*)\right]\]

여기서 \(\hat{\alpha}_0\) 는 \(\beta = 0\) 으로 고정했을 때의 \(\alpha\) MLE이다. 이항 인수만 남으면

\[\hat{\alpha}_0 = \log\!\left(\frac{\bar{y}}{1-\bar{y}}\right), \quad \bar{y} = \frac{\sum_j y_j^*}{\sum_j n_j}\]

(전체 성공률의 로짓, Exercise 12.20). \(H_0\) 하에 \(-2\log\lambda \xrightarrow{d} \chi^2_1\) 이므로, \(-2\log\lambda \geq \chi^2_{1,\alpha}\) 이면 기각한다.

Wald vs LRT: Wald 검정이 계산은 더 쉽지만, 표본이 작거나 분리 문제가 있을 때 LRT가 더 신뢰할 수 있다. 이유: Wald 검정은 \(\hat{\beta}\) 의 점근 정규 분포 근사에 의존하는데, \(|\hat{\beta}|\) 가 크거나 표본이 작으면 이 근사가 나빠진다. Hauck & Donner(1977)는 \(|\hat{\beta}|\) 가 매우 클 때 Wald 검정의 p-값이 오히려 커지는 — 즉, 강한 효과가 있음에도 귀무가설을 기각하지 못하는 — 역설적 현상을 보였다. LRT는 두 모형의 로그 우도 차이를 직접 비교하므로 \(\hat{\beta}\) 의 점근 정밀도에 덜 의존하고 이 역설이 발생하지 않는다.

---

10 코드 예시

10.1 Step 1: 순수 Python 구현 (원리 이해)

import numpy as np
from scipy.optimize import minimize

# ──────────────────────────────────────
# 로지스틱 함수와 로그 우도
# ──────────────────────────────────────
def logistic(eta):
    """시그모이드 함수 π(x) = e^η / (1 + e^η)"""
    return np.where(eta >= 0,
                    1 / (1 + np.exp(-eta)),
                    np.exp(eta) / (1 + np.exp(eta)))

def log_likelihood(params, x, y):
    """베르누이 로그 우도"""
    alpha, beta = params
    eta = alpha + beta * x
    pi = logistic(eta)
    # 수치 안정성: clip 처리
    pi = np.clip(pi, 1e-12, 1 - 1e-12)
    return np.sum(y * np.log(pi) + (1 - y) * np.log(1 - pi))

def score(params, x, y):
    """스코어 벡터 (∂ℓ/∂α, ∂ℓ/∂β)"""
    alpha, beta = params
    pi = logistic(alpha + beta * x)
    residual = y - pi
    return np.array([residual.sum(), (residual * x).sum()])

def fisher_info(params, x):
    """피셔 정보행렬 I(α,β) — 로지스틱 링크 (n_j = 1)"""
    alpha, beta = params
    pi = logistic(alpha + beta * x)
    w = pi * (1 - pi)  # 가중치
    I = np.array([
        [w.sum(),       (w * x).sum()],
        [(w * x).sum(), (w * x**2).sum()]
    ])
    return I

# ──────────────────────────────────────
# 챌린저 데이터 (Casella & Berger Table 12.3.1)
# ──────────────────────────────────────
temp = np.array([53, 57, 58, 63, 66, 67, 67, 67, 68, 69, 70, 70,
                 70, 70, 72, 73, 75, 75, 76, 76, 78, 79, 81])
failure = np.array([1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0,
                    1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0])

# MLE: 음의 로그 우도 최소화
n = len(temp)
result = minimize(
    fun=lambda p: -log_likelihood(p, temp, failure),
    x0=[0.0, 0.0],
    method='BFGS'
)
alpha_hat, beta_hat = result.x
print(f"α̂ = {alpha_hat:.3f}  (Casella 예시: 15.043)")
print(f"β̂ = {beta_hat:.3f}  (Casella 예시: -0.232)")

# 스코어가 0에 수렴 확인
s = score([alpha_hat, beta_hat], temp, failure)
print(f"\n스코어 벡터 (≈ 0 이어야): {s}")

# 피셔 정보행렬 및 표준오차
I_hat = fisher_info([alpha_hat, beta_hat], temp)
I_inv = np.linalg.inv(I_hat)
se_alpha = np.sqrt(I_inv[0, 0])
se_beta  = np.sqrt(I_inv[1, 1])
print(f"\nI(α̂,β̂)⁻¹ = \n{I_inv.round(3)}")
print(f"se(α̂) = {se_alpha:.3f}")
print(f"se(β̂) = {se_beta:.3f}")

# 95% 신뢰구간
from scipy.stats import norm
z = norm.ppf(0.975)
ci_beta = (beta_hat - z * se_beta, beta_hat + z * se_beta)
print(f"\nβ의 95% CI: [{ci_beta[0]:.3f}, {ci_beta[1]:.3f}]")
print("(Casella 예시: [-0.447, -0.017])")

# 31°F에서의 예측 확률
pi_31 = logistic(alpha_hat + beta_hat * 31)
print(f"\n31°F에서 O-링 실패 확률: {pi_31:.4f}")

# Wald 검정
Z_wald = beta_hat / se_beta
p_wald = 2 * norm.sf(abs(Z_wald))
print(f"\nWald Z = {Z_wald:.3f}, p-value = {p_wald:.4f}")

# LRT
alpha0 = np.log(failure.mean() / (1 - failure.mean()))
ll_full = log_likelihood([alpha_hat, beta_hat], temp, failure)
ll_null = log_likelihood([alpha0, 0.0], temp, failure)
lrt_stat = 2 * (ll_full - ll_null)
from scipy.stats import chi2
p_lrt = chi2.sf(lrt_stat, df=1)
print(f"LRT 통계량 = {lrt_stat:.3f}, p-value = {p_lrt:.4f}")

결과 해석: \(\hat{\alpha} \approx 15.043\), \(\hat{\beta} \approx -0.232\) 는 Casella 예시와 일치한다. 스코어 벡터가 \(\approx (0, 0)\) 임을 확인하면 수렴이 정상이다. \(I^{-1}\) 대각 원소에서 \(\text{se}(\hat{\beta}) \approx 0.110\) 이므로 Wald 통계량 \(Z = -0.232/0.110 \approx -2.11\), p-값 \(\approx 0.035\) 로 유의수준 5%에서 유의하다. LRT 통계량 역시 이와 유사한 결론을 줄 것으로 예상되며, 두 검정 결과가 크게 다르다면 분리 문제나 소표본 문제를 의심해야 한다.

10.2 Step 2: statsmodels 구현 (실무 활용)

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import statsmodels.api as sm
from scipy.stats import norm

# 챌린저 데이터
temp = np.array([53, 57, 58, 63, 66, 67, 67, 67, 68, 69, 70, 70,
                 70, 70, 72, 73, 75, 75, 76, 76, 78, 79, 81])
failure = np.array([1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0,
                    1, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0])

# statsmodels 로지스틱 회귀
X = sm.add_constant(temp)
model = sm.Logit(failure, X)
result = model.fit(disp=False)

print(result.summary())
print(f"\n오즈비 e^β: {np.exp(result.params[1]):.4f}")
print(f"  → 기온 1°F 상승 시 O-링 실패 오즈가 {np.exp(result.params[1]):.4f}배")

# 예측 곡선 시각화
x_grid = np.linspace(30, 90, 300)
X_grid = sm.add_constant(x_grid)
pi_pred = result.predict(X_grid)

fig, ax = plt.subplots(figsize=(10, 5))
ax.scatter(temp, failure, color='k', zorder=5, label='관측 데이터')
ax.plot(x_grid, pi_pred, 'b-', lw=2, label='로지스틱 적합 곡선')
ax.axvline(31, color='r', ls='--', lw=1.5, label='챌린저 발사 온도 (31°F)')
ax.axhline(0.5, color='gray', ls=':', lw=1)
ax.set_xlabel('발사 온도 (°F)')
ax.set_ylabel('O-링 실패 확률 π(x)')
ax.set_title('챌린저 O-링 데이터 — 로지스틱 회귀 적합')
ax.legend()
plt.tight_layout()
plt.show()

# 로지스틱 vs 프로빗 비교
model_probit = sm.Probit(failure, X)
result_probit = model_probit.fit(disp=False)
pi_probit = result_probit.predict(X_grid)

fig, ax = plt.subplots(figsize=(10, 5))
ax.plot(x_grid, pi_pred,   'b-',  lw=2, label='로지스틱 (logit)')
ax.plot(x_grid, pi_probit, 'r--', lw=2, label='프로빗 (probit)')
ax.scatter(temp, failure, color='k', zorder=5)
ax.set_xlabel('발사 온도 (°F)')
ax.set_ylabel('π(x)')
ax.set_title('로지스틱 vs 프로빗 회귀 비교')
ax.legend()
plt.tight_layout()
plt.show()

# 모수 해석: 오즈비 테이블
params = result.params
conf = result.conf_int()
odds_table = pd.DataFrame({
    '계수 (log-odds)': params,
    '오즈비 (e^β)': np.exp(params),
    'CI 하한 (OR)': np.exp(conf[0]),
    'CI 상한 (OR)': np.exp(conf[1])
})
print("\n모수 해석 (오즈비 기준):")
print(odds_table.round(4))

---

11 요약

항목	내용
GLM 3요소	랜덤(베르누이), 선형(\(\alpha+\beta x\)), 링크(로짓)
정준 링크	\(g(\pi) = \log[\pi/(1-\pi)]\) — 베르누이 자연 모수가 선형 예측자
로지스틱 함수	\(\pi(x) = e^{\alpha+\beta x}/(1+e^{\alpha+\beta x})\), 범위 \((0,1)\)
도함수	\(d\pi/dx = \beta\pi(1-\pi)\) — \(\beta\) 의 부호가 단조성 결정
α 해석	\(x=0\) 에서의 로그 오즈
β 해석	\(x\) 1단위 증가 시 로그 오즈 변화; \(e^\beta\) = 오즈비
스코어 방정식	\(\sum(y_i - \hat{\pi}_i) = 0\), \(\sum(y_i-\hat{\pi}_i)x_i=0\) (비선형, 수치 해법)
완전 분리	\(\hat{\beta} \to \pm\infty\), 스코어 해 없음 — 극단 데이터에서 발생
정보행렬	\(I(\alpha,\beta)_{jk}\): 가중 설계 행렬 \(\mathbf{X}^T\mathbf{W}\mathbf{X}\) 구조
근사 CI	\(\hat{\beta} \pm z_{\alpha/2}\text{se}(\hat{\beta})\), \(\text{se}\) 는 정보행렬 역행렬 대각 원소
Wald 검정	\(Z = \hat{\beta}/\text{se}(\hat{\beta}) \sim N(0,1)\) under \(H_0:\beta=0\)
LRT	\(-2\log\lambda \sim \chi^2_1\) under \(H_0\)
챌린저 결과	\(\hat{\alpha}=15.04\), \(\hat{\beta}=-0.232\), 31°F에서 실패 확률 \(\approx 0.9996\)

---

12 관련 주제

선행 지식

• GLM 프레임워크 — GLM 개관, 설계 행렬, 실무 예시
• 최대우도추정량 — MLE 일반론
• 점근 평가 — MLE 점근 분포, 정보 행렬

후속 주제

• Logistic Regression: Estimation Algorithm — IRLS 상세, Newton-Raphson 유도
• Robust Regression — Casella §12.4, M-추정, 영향함수
• Multiple Linear Regression — 행렬 OLS, Gauss-Markov
• GLM 확장 — 포아송 회귀, 감마 회귀, 과산포(overdispersion)

관련 개념

• 지수족 분포 — 정준 링크의 이론적 기반
• EIV 모형 — 같은 Ch.12의 §12.2
• 회귀 모형 개관 — EIV·로지스틱·로버스트의 계보