§ 12.3 — ZIP Model 의 깊이: Mixture, Unified PDF, EM vs Newton-Raphson

1 들어가며 — 12-0, 12-1 의 깊이

Ch.12 Overview (12-0) 가 ZIP 의 큰 그림을, § 12.1 ~ 12.2 (12-1) 가 modified Poisson family 의 분류를 다뤘다. 본 sub-post 는 § 12.3 ZIP 모형의 수학적 도출 + 추정 알고리즘 을 깊이 있게.

한 줄 요약

“§ 12.3 = Lambert (1992) ZIP 의 정확한 정의 (식 12.7-12.9) — mixture of perfect state (\(\pi_i\)) + Poisson state (\(1-\pi_i\)). ZIP(τ) variants (식 12.10-12.11) 가 같은 covariate + functional 관계로 모수 절약. Greene (1994) 의 unified pdf (식 12.13) — indicator \(I(y_i)\) 로 두 component 결합 = mixture 의 elegant 표현. Log-likelihood 의 score (식 12.14-12.16) chain rule 도출 — ‘잔차 가중 + perfect zero 보정’ 의 두 항 분해. 추정 알고리즘 3 가지 — Lambert 의 EM (두 separate log-likelihoods), Newton-Raphson (ZIP(τ) 의 직접 최대화), Greene 의 gradient method (빠른 수렴). BHHH (1974) information matrix — score 의 outer product 합. 한계: random effects 없음 → § 12.4 mixed-effects 확장 필요.”

2 § 12.3 — ZIP 모형의 정확한 정의

2.1 Lambert (1992) Mixture 모형

식 (12.7) — ZIP 의 mixture 표현

응답 \(y_i\) 가 두 종류 generation process:

\[ y_i \sim \begin{cases} 0 & \text{with probability } \pi_i \\ \text{Poisson}(\lambda_i) & \text{with probability } 1 - \pi_i \end{cases} \tag{12.7} \]

\(i = 1, \ldots, N\).

식 (12.8-12.9) — Link 함수

\(\pi_i\) 는 logit link 의 logistic regression:

\[ \text{logit}(\pi_i) = w_i^\top \gamma \tag{12.8} \]

\(\lambda_i\) 는 log link 의 Poisson regression:

\[ \log(\lambda_i) = x_i^\top \beta \tag{12.9} \]

표기:

• \(w_i\): \((s+1) \times 1\) — logistic part 의 covariate.
• \(\gamma\): \((s+1) \times 1\) — logistic 회귀 계수.
• \(x_i\): \((p+1) \times 1\) — Poisson part 의 covariate.
• \(\beta\): \((p+1) \times 1\) — Poisson 회귀 계수.

\(x\) 와 \(w\) 가 같을 수도, 다를 수도 — 모형의 유연성.

직관 — Canonical Link 의 자연 선택

식 (12.8) 의 logit link, 식 (12.9) 의 log link 가 각각 Bernoulli, Poisson 의 canonical link (§ 12.1 의 12-1 참조).

McCullagh & Nelder (1989) 의 GLM framework — canonical link 사용 시:

• Score 의 단순한 형태 (‘잔차 × 공변량’).
• Newton-Raphson = Fisher scoring.
• 추정의 안정성.

ZIP 가 두 GLM 의 결합 — 두 canonical link 모두 사용. 자연스러운 framework.

Lambert 의 명명:

• \(y_i = 0\) as perfect state or zero state (Lambert 의 정의 — 제조에서 결함 없음).
• \(y_i \sim \text{Poisson}\) as imperfect state or Poisson state (결함 발생 가능).

Zorn (1996) 의 명명:

• Transition stage: 어떤 state 인지 결정 (perfect vs Poisson).
• Event stage: Poisson state 라면 몇 개 사건 발생.

→ 두 명명이 같은 mathematical structure, 다른 metaphor. Lambert 는 manufacturing, Zorn 은 dual regime.

2.2 ZIP(τ) Variants — 식 (12.10-12.11)

ZIP(τ) — Functional 관계

ZIP 의 두 부분이 같은 covariate \(x_i\) 사용 + functional 관계:

\[ \text{logit}(\pi_i) = -\tau x_i^\top \beta \tag{12.10} \]

\[ \log(\lambda_i) = x_i^\top \beta \tag{12.11} \]

추가 모수 \(\tau\) — 두 부분의 효과 비례 상수.

직관 — ZIP vs ZIP(τ) 의 모수 절약

ZIP (식 12.7-12.9):

• 자유 모수: \(\beta\) (\(p+1\) 개) + \(\gamma\) (\(s+1\) 개) = \(p + s + 2\) 개.
• 가장 유연 — 각 part 별도 covariate.

ZIP(τ) (식 12.10-12.11):

• 자유 모수: \(\beta\) (\(p+1\) 개) + \(\tau\) (1 개) = \(p + 2\) 개.
• 절약적 — 같은 covariate.

\(-\tau\) 부호의 의미:

• \(\tau > 0\): \(x^\top\beta\) ↑ → \(\lambda\) ↑ (Poisson rate ↑) + \(\text{logit}(\pi)\) ↓ (perfect 확률 ↓).
• 즉 covariate 가 Poisson rate 를 증가시킬 때, 동시에 perfect (구조적 zero) 확률을 감소시킴.

임상적 의미 (mental health 예):

• 보험 가입 (\(x \uparrow\)) → 진료 빈도 ↑ (\(\lambda\)) + 진료 미필요 비율 ↓ (\(\pi\)).
• 두 효과가 자연스럽게 같은 방향 (보험이 의료 접근성과 인지 둘 다 영향).

ZIP(τ) 가 적합한 경우:

• 두 부분의 효과 메커니즘이 같음 (예: 같은 socioeconomic 변수).
• 모수 절약 필요 (작은 표본).
• 해석 단순 (single set of coefficients).

ZIP 가 적합한 경우:

• 두 부분의 효과 메커니즘이 다름.
• 별도 covariate 사용 (예: \(w\) = 보험 종류, \(x\) = 거주지).
• 효과 패턴 분리 필요.

2.3 대안 Link 함수

\(\pi_i\) 의 다양한 link options (Lambert 1992)

ZIP(τ) framework 에서 \(\pi_i\) 의 link 다양화:

Logit (default): \(\text{logit}(\pi_i) = -\tau x_i^\top \beta\).

Log-log: \(\log(-\log(\pi_i)) = \tau x_i^\top \beta\).

→ \(\pi_i = \exp(-\exp(\tau x_i^\top \beta))\).

Complementary log-log: \(\log(-\log(1 - \pi_i)) = -\tau x_i^\top \beta\).

→ \(\pi_i = 1 - \exp(-\exp(-\tau x_i^\top \beta))\).

직관 — 비대칭 link 의 의미

Logit 은 0 과 1 모두 대칭으로 다룸. Log-log 와 complementary log-log 는 비대칭:

• Log-log: \(\pi_i\) 가 1 에 가까운 영역에 sensitive (perfect state 가 흔할 때).
• C-log-log: \(\pi_i\) 가 0 에 가까운 영역에 sensitive (Poisson state 가 우세할 때).

§ 9.5.5 의 c-log-log 와 같은 발상 — 비대칭 잠재 분포 (Gumbel) 가 link 함수의 형태 결정.

언제 비대칭 link 사용:

• 응답에 명확한 비대칭이 있을 때.
• Hurdle 모형의 logit vs c-log-log 와 같은 선택.

실무에서는 logit 이 가장 흔함 — 단순성 + 표준성. 비대칭 link 는 특수 응용 (예: 산업 결함 모니터링).

3 Greene (1994) 의 Unified Probability Density

3.1 식 (12.12) 의 자연 도출

식 (12.12) — Re-expressed mixture

식 (12.7) 의 mixture 를 풀면:

\[ P(y_i = 0) = \pi_i + (1 - \pi_i) e^{-\lambda_i} \]

\[ P(y_i = k) = (1 - \pi_i) \frac{e^{-\lambda_i} \lambda_i^k}{k!}, \quad k = 1, 2, \ldots \tag{12.12} \]

직관 — 두 종류의 zero 가 합산

\(P(y_i = 0)\) 의 두 항:

• \(\pi_i\): Perfect zero — perfect state 의 확률 (구조적 zero).
• \((1 - \pi_i) e^{-\lambda_i}\): Poisson zero — Poisson state 인데 우연히 0 인 확률.

두 zero 종류가 분리 안 됨 — 관측에서 어느 종류인지 구분 불가능.

그러나 모형에서는 mixing parameter \(\pi_i\) 로 비율 추정:

• 표본의 zero 비율이 \(e^{-\hat\lambda}\) 보다 크면 \(\hat\pi > 0\).
• 차이가 클수록 \(\hat\pi\) 가 큼 (perfect state 비율 ↑).

WZ (with zeros) 명명 — Johnson et al. (1992):

• 같은 mathematical form 의 또 다른 이름.
• ZAP (Heilbron 1989) 와 동일.
• ZIP (Lambert 1992) 가 가장 널리 사용 — 통일 권고.

Hurdle 과의 결정적 차이 (12-1 에서 다룸):

• ZIP: Poisson part 에 zero 가능 (mixture).
• Hurdle: Positive part 가 truncated Poisson (zero 안 나옴).
• \(\pi_i = 0\) 일 때 ZIP 는 표준 Poisson 으로 reduce.

3.2 식 (12.13) — Greene Unified PDF

Indicator function 활용

Greene (1994) 의 elegant 표기:

\[ P(y_i) = p(y_i) = (1 - \pi_i) f(y_i) + I(y_i) \pi_i \tag{12.13} \]

표기:

• \(\Psi(\cdot)\): 표준 logistic cdf — \(\pi_i = \Psi(w_i^\top \gamma) = e^{w_i^\top \gamma} / (1 + e^{w_i^\top \gamma})\).
• \(f(y_i)\): Poisson density with mean \(\lambda_i\).
• \(I(y_i)\): indicator function — \(I(y_i) = 1\) if \(y_i = 0\), else 0.

직관 — Indicator 의 역할

식 (12.13) 의 두 항:

• \((1 - \pi_i) f(y_i)\): Poisson contribution (모든 \(y_i\)).
• \(I(y_i) \pi_i\): Perfect contribution (only when \(y_i = 0\)).

\(y_i > 0\) 일 때:

• \(I(y_i) = 0\) → 둘째 항 사라짐.
• \(p(y_i) = (1 - \pi_i) f(y_i)\) — 식 (12.12) 의 둘째 줄.

\(y_i = 0\) 일 때:

• \(I(y_i) = 1\) → 둘째 항 = \(\pi_i\).
• \(p(y_i) = (1 - \pi_i) f(0) + \pi_i = (1 - \pi_i) e^{-\lambda_i} + \pi_i\) — 식 (12.12) 의 첫 줄.

Unified 표기의 장점:

• 미분이 단순화 — 두 case 분기 없이 chain rule 직접 적용.
• Score 와 information matrix 도출이 elegant.
• 코드 구현이 간단.

이것이 식 (12.14-12.18) 의 score 와 information matrix 도출의 토대.

4 Log-Likelihood 와 Score — 식 (12.14-12.16)

4.1 Log-Likelihood — 식 (12.14)

표본 log-likelihood

\(N\) 독립 환자의 log-likelihood:

\[ \log L = \sum_{i=1}^N \log p(y_i) = \sum_i \log [(1 - \pi_i) f(y_i) + I(y_i) \pi_i] \tag{12.14} \]

직관 — Mixture 의 log 의 어려움

식 (12.14) 의 log 안에 두 항의 합 — log-sum 의 형태. 미분이 복잡:

• 표준 GLM (Poisson, logistic): log 안에 한 항 — chain rule 단순.
• Mixture (ZIP): log 안에 합 — chain rule 복잡 (분모에 sum 등장).

이 어려움이 EM 알고리즘 (Dempster et al. 1977) 의 자연 적용 동기:

• Latent indicator \(z_i\) (“perfect state 인지”) 로 분해.
• E-step: \(z_i\) 의 사후 확률 계산.
• M-step: \(z_i\) 가 알려진 것처럼 두 GLM 별도 적합.
• → log-sum 회피, 두 simple GLM 으로 분해.

자세한 알고리즘 비교는 다음 절.

4.2 Score 의 Chain Rule 도출 — 식 (12.15-12.16)

\(\beta\) (Poisson part) 의 score — 식 (12.15)

\[ \frac{\partial \log L}{\partial \beta} = \sum_i \frac{1}{p(y_i)} \frac{\partial p(y_i)}{\partial \beta} \]

\(\partial p / \partial \beta\) 의 도출 (식 12.13 미분):

\[ \frac{\partial p(y_i)}{\partial \beta} = (1 - \pi_i) \frac{\partial f(y_i)}{\partial \beta} + (I(y_i) - f(y_i)) \frac{\partial \pi_i}{\partial \beta} \]

(만약 \(\pi_i\) 가 \(\beta\) 의 함수면 — ZIP(τ) 의 경우.)

ZIP 모형 (별도 covariate) 에서는 \(\partial \pi_i / \partial \beta = 0\):

\[ \frac{\partial \log L}{\partial \beta} = \sum_i \frac{(1 - \pi_i) f(y_i)}{p(y_i)} \frac{\partial \log f(y_i)}{\partial \beta} \]

\(\gamma\) (logistic part) 의 score — 식 (12.16)

\(\pi_i = \Psi(w_i^\top \gamma)\) 의 미분 — chain rule:

\[ \frac{\partial \pi_i}{\partial \gamma} = \pi_i (1 - \pi_i) w_i = \pi_i' w_i \]

여기서 \(\pi_i' = \pi_i (1 - \pi_i)\) — logistic pdf (§ 9.5.5 의 \(\Psi' = \Psi(1-\Psi)\)).

식 (12.16):

\[ \frac{\partial \log L}{\partial \gamma} = \sum_i \frac{1}{p(y_i)} [I(y_i) - f(y_i)] \frac{\partial \pi_i}{\partial \gamma} \]

또는

\[ = \sum_i \frac{[I(y_i) - f(y_i)] \pi_i (1 - \pi_i)}{p(y_i)} w_i \]

직관 — Score 의 두 단계 분해

\(\gamma\) score 의 핵심 항 \([I(y_i) - f(y_i)]\):

• \(y_i = 0\): \(I(0) - f(0) = 1 - e^{-\lambda_i}\). 양수 — perfect zero 의 신호.
• \(y_i > 0\): \(I(y_i) - f(y_i) = 0 - f(y_i) < 0\). 음수 — Poisson state 의 신호.

해석:

• 표본의 zero 비율이 \(f(0) = e^{-\lambda}\) 보다 크면 → score 가 양수 → \(\pi\) 증가 방향.
• Zero 비율이 \(f(0)\) 와 비슷하면 → score 가 0 근처 → \(\pi \approx 0\) 안정.
• → Score 가 자동으로 excess zeros 의 강도 측정.

\(\beta\) score 의 핵심 항 \((1 - \pi_i) f(y_i) / p(y_i)\):

• \(y_i > 0\): \((1 - \pi_i) f(y_i) = p(y_i)\) → 비율 = 1 → 표준 Poisson score 그대로.
• \(y_i = 0\): \((1 - \pi_i) f(0) / p(0) = (1 - \pi_i) e^{-\lambda} / [(1 - \pi_i) e^{-\lambda} + \pi]\) < 1.
• → Poisson zero 의 weight 가 줄어듦 — perfect zero 가 weight 일부 차지.

해석:

• ZIP score 가 Poisson part 에 대해 “weight-down zero” 효과.
• 표준 Poisson 의 score 에서 zero contribution 을 줄임.
• 결과적으로 추정된 \(\hat\lambda\) 가 표준 Poisson 보다 큼 — zero 가 perfect state 로 일부 흡수.

5 Information Matrix — 식 (12.17-12.18)

5.1 ZIP 모형 — 식 (12.17)

ZIP 의 information matrix (별도 covariate)

\[ \widehat I(\widehat\beta, \widehat\gamma) = -\begin{bmatrix} \frac{\partial^2 \log p}{\partial \beta \partial \beta^\top} & \frac{\partial^2 \log p}{\partial \beta \partial \gamma^\top} \\ \frac{\partial^2 \log p}{\partial \gamma \partial \beta^\top} & \frac{\partial^2 \log p}{\partial \gamma \partial \gamma^\top} \end{bmatrix} \tag{12.17} \]

Block matrix:

• 대각: \(\beta\) 와 \(\gamma\) 의 self-information.
• 비대각: 두 모수 사이 cross-information.

직관 — Block 구조의 의미

별도 covariate (\(x \neq w\)) 일 때:

• \(\beta\) 와 \(\gamma\) 가 다른 변수.
• 비대각 항이 작음 (또는 0 근처) — 두 추정이 거의 독립.
• → Block diagonal 에 가까움.

같은 covariate (\(x = w\)) 일 때:

• \(\beta\) 와 \(\gamma\) 가 같은 변수에 대한 효과.
• 비대각 항이 큼 — 두 추정 사이 강한 상관.
• → 추정이 unstable, identification 어려움.

이것이 ZIP(τ) 가 발전한 이유 — 같은 covariate 면 functional 관계 부여로 모수 식별 가능.

비대각 항의 추정 어려움:

• \(y_i = 0\) 의 contribution 이 \(\beta, \gamma\) 모두에 영향.
• Mixture 의 ambiguity — 각 zero 가 어느 source 인지 모름.
• → 큰 표본 + clear excess zeros 가 추정에 필수.

5.2 ZIP(τ) 모형 — 식 (12.18)

ZIP(τ) 의 information matrix

\[ \widehat I(\widehat\beta, \widehat\tau) = \begin{bmatrix} \frac{\partial^2 \log p}{\partial \beta \partial \beta^\top} & \frac{\partial^2 \log p}{\partial \beta \partial \tau} \\ \frac{\partial^2 \log p}{\partial \tau \partial \beta^\top} & \frac{\partial^2 \log p}{\partial \tau^2} \end{bmatrix} \tag{12.18} \]

ZIP vs ZIP(τ) 의 모수 차원

차원 비교:

• ZIP: \((p + s + 2) \times (p + s + 2)\) — 큰 행렬.
• ZIP(τ): \((p + 2) \times (p + 2)\) — 훨씬 작음.

→ ZIP(τ) 의 추정이 더 간단 + 안정.

\(\tau\) 추정의 어려움:

• \(\tau\) 가 두 부분의 비례 상수 — 작은 변화가 두 효과 모두에 영향.
• Hessian 의 \(\tau\) block 이 작은 값 → 표준오차 큼.
• 큰 표본 필요 (수백 ~ 수천).

ZIP(τ) 의 적합도가 ZIP 보다 약간 떨어질 수 있지만 (자유도 적음), 모형 절약과 안정성의 trade-off.

6 추정 알고리즘의 비교

6.1 EM Algorithm — Lambert (1992) 의 ZIP

EM 의 분해 발상

ZIP 의 mixture 구조를 latent indicator \(z_i\) 로 분해:

• \(z_i = 1\): perfect state.
• \(z_i = 0\): Poisson state.

Joint distribution:

\[ P(y_i, z_i) = \pi_i^{z_i} [(1 - \pi_i) f(y_i)]^{1-z_i} \]

Complete log-likelihood (만약 \(z_i\) 알려져 있다면):

\[ \log L_c = \sum_i \{z_i \log \pi_i + (1 - z_i) [\log(1 - \pi_i) + \log f(y_i)]\} \]

→ 두 separate GLM 의 log-likelihood:

• Logistic (with \(z_i\) as outcome).
• Poisson (only for \(z_i = 0\) subset, with weights).

EM 알고리즘의 단계

E-step: \(z_i\) 의 사후 확률 계산:

\[ \hat z_i = E(z_i \mid y_i) = \begin{cases} \frac{\pi_i}{p(y_i)} = \frac{\pi_i}{\pi_i + (1 - \pi_i) e^{-\lambda_i}} & \text{if } y_i = 0 \\ 0 & \text{if } y_i > 0 \end{cases} \]

M-step: \(\hat z_i\) 가 알려진 것처럼 두 GLM 별도 적합:

1. Logistic regression: \(z_i\) vs \(w_i\) → \(\hat\gamma\) 갱신.
2. Weighted Poisson regression: \(y_i\) vs \(x_i\) with weight \((1 - \hat z_i)\) → \(\hat\beta\) 갱신.

수렴까지 반복.

EM 의 장점:

• 각 step 이 표준 GLM — 단순.
• 결정론적 수렴 (likelihood 매 iteration 증가).
• 코드 구현 쉬움.

EM 의 단점:

• 수렴이 느림 (특히 ZIP 의 모호한 zero 분리).
• 표준오차 직접 산출 안 됨 (별도 계산 필요).
• 매 step 의 GLM 적합도 자체 iteration → 전체 비용 큼.

6.2 Newton-Raphson — ZIP(τ) 직접 최대화

Newton-Raphson update

Lambert 의 권고 (ZIP(τ) 의 경우):

\[ (\beta, \tau)^{new} = (\beta, \tau)^{old} + I^{-1}(\beta, \tau) \nabla \log L \]

식 (12.15-12.16) 의 score + 식 (12.18) 의 information matrix 사용.

매 iteration 에서:

1. 현재 \((\hat\beta, \hat\tau)\) 에서 \(\hat\lambda_i, \hat\pi_i\) 계산.
2. Score 계산.
3. Hessian 또는 Fisher information 계산.
4. Update.

ZIP(τ) 가 Newton-Raphson 에 적합한 이유

ZIP (별도 covariate) 의 EM 이 자연스러운 반면, ZIP(τ) 는 Newton-Raphson 이 더 자연:

• ZIP(τ) 의 모수가 적음 (p + 2) — Newton-Raphson 의 \(O(p^2)\) 부담 작음.
• \(\tau\) 의 EM 분해가 어려움 — 두 부분이 functional 관계로 묶여 있음.
• Newton-Raphson 의 직접 최대화가 더 단순.

수렴 속도:

• Newton-Raphson: 2 차 수렴 (수렴 시 매우 빠름).
• EM: 1 차 수렴 (느림).

→ ZIP(τ) 에서 Newton-Raphson 권고. ZIP 에서는 EM (또는 hybrid).

6.3 Greene (1994) Gradient Method

더 빠른 수렴을 위한 alternative

Greene (1994) 의 권고 — ZIP 와 ZIP(τ) 모두에:

• BHHH (Berndt, Hall, Hall, Hausman 1974) algorithm: Gradient 만 사용 (Hessian 없음).
• Information matrix 추정: score 의 outer product 합.

\[ \widehat I_{\text{BHHH}} = \sum_i \frac{\partial \log p_i}{\partial \theta} \frac{\partial \log p_i}{\partial \theta^\top} \]

여기서 \(\theta = (\beta, \gamma)\) 또는 \((\beta, \tau)\).

직관 — BHHH 의 특별한 성질

BHHH 의 두 장점:

1. Hessian 계산 불필요: 2 차 도함수 도출이 어려운 모형 (ZIP) 에 유리.
2. 항상 양정치: outer product 의 합 → eigendecomposition 안정 — Newton-Raphson 의 negative Hessian 문제 회피.

근거 (Information matrix equality, 정칙성 조건 하):

\[ E\left[\frac{\partial \log L}{\partial \theta} \frac{\partial \log L}{\partial \theta^\top}\right] = -E\left[\frac{\partial^2 \log L}{\partial \theta \partial \theta^\top}\right] = I(\theta) \]

→ Score 의 outer product 의 기댓값 = 정보 행렬. BHHH 는 표본 평균으로 추정.

ZIP 에서의 장점:

• ZIP 의 mixture 형태로 Hessian 계산이 복잡.
• BHHH 가 score 만으로 정보 추정 → 알고리즘 단순.
• 빠른 수렴 (Greene 의 시뮬레이션 검증).

단점:

• 수렴이 진정한 Newton-Raphson 보다 느릴 수 있음 (BHHH 는 BFGS-like quasi-Newton).
• 표준오차가 Hessian-based 보다 약간 부정확.

현재 (2026) 의 표준:

• R pscl::zeroinfl: BFGS 기본 (BHHH 와 비슷).
• SAS PROC NLMIXED: Newton-Raphson 또는 quasi-Newton.
• Python statsmodels.discrete.count_model.ZeroInflatedPoisson: BFGS.

7 한계 — Random Effects 부재

ZIP 의 결정적 한계

저자 본문 인용:

“A major limitation of these models is that they do not accommodate random effects.”

ZIP 의 한계:

• Cross-sectional 데이터에만 적용.
• 같은 환자의 반복 측정 → 독립 가정 위반.
• 클러스터 효과 (지역, 가족, 학교) → 처리 불가.
• → §12.4 mixed-effects ZIP 으로 확장 필요.

직관 — Random Effects 의 두 부분 모두 필요

ZIP 의 두 component 에 모두 random effects 추가가 자연:

• Logistic part 의 random effect (\(\sigma_2 \theta_{2i}\)): 환자별 perfect state 확률의 차이.
  • “어떤 환자는 일관되게 진료 필요 없는 그룹, 어떤 환자는 항상 필요한 그룹.”
• Poisson part 의 random effect (\(\sigma_1 \theta_{1i}\)): 환자별 Poisson rate 의 차이.
  • “필요 그룹 안에서도 어떤 환자는 자주, 어떤 환자는 가끔 방문.”

두 종류의 환자 이질성 — § 12.4.3 의 mixed-effects ZIP 가 처리 (식 12.30-12.33, 12-0 에서 다룸).

Hall (2000) 의 부분적 해결:

• Poisson part 에만 random effects.
• Logistic part 는 fixed.
• 한계: logistic 부분의 환자 차이 무시.

Hur et al. (2002) 의 단일 random effect:

• 두 부분에 같은 random effect 사용.
• 가정: 두 부분의 환자 효과가 완전 상관.
• 단순하지만 가정 강함.

Hedeker & Gibbons (2006, § 12.4.3) 의 두 random effects:

• 가장 일반적 — 두 부분에 별도 random effects.
• 두 random effects 간 상관 가능 (이변량 정규).
• → 다음 sub-post 에서 자세히.

8 응용 분야

분야	ZIP vs ZIP(τ)	적용 이유
의료 이용	ZIP (별도 covariate)	“필요성” 과 “빈도” 의 결정 요인 다름
흡연 횟수	ZIP(τ) (같은 covariate)	같은 socioeconomic 변수가 두 부분에 영향
보험 청구	ZIP	청구 의사 vs 청구 횟수 분리
사고 발생	ZIP(τ)	위험 노출이 두 부분에 비례적 효과
산업 결함	ZIP (Lambert 1992 원전)	정상 공정 vs 결함 발생 공정
자살 예방	ZIP	자살 발생 가능성 vs 발생 횟수
마케팅	ZIP	비구매자 vs 구매 횟수

9 코드 예시

9.1 Step 1: ZIP Density 의 직접 구현

import numpy as np
from scipy.special import expit, factorial


def zip_density(y: np.ndarray, lam: np.ndarray, pi: np.ndarray) -> np.ndarray:
    """식 (12.13) — Greene 의 unified ZIP density.

    P(y) = (1 - pi) * Poisson(y; lam) + I(y=0) * pi
    """
    poisson_pmf = np.exp(-lam) * lam ** y / factorial(y)
    indicator = (y == 0).astype(float)
    return (1 - pi) * poisson_pmf + indicator * pi


# 단순 사례
y_test = np.array([0, 1, 2, 5, 0, 0])
lam_test = np.array([2.0, 2.0, 2.0, 2.0, 2.0, 2.0])
pi_test = np.array([0.3, 0.3, 0.3, 0.3, 0.3, 0.3])

p = zip_density(y_test, lam_test, pi_test)
print(f"ZIP probabilities (lambda=2, pi=0.3):")
for y, prob in zip(y_test, p):
    print(f"  P(Y = {y}) = {prob:.4f}")

# 표준 Poisson 과 비교
poisson_p = np.exp(-2.0) * 2.0 ** y_test / factorial(y_test)
print(f"\nPoisson probabilities (lambda=2):")
for y, prob in zip(y_test, poisson_p):
    print(f"  P(Y = {y}) = {prob:.4f}")

# Excess zeros 검증
print(f"\nP(Y = 0):")
print(f"  Poisson: {poisson_p[0]:.4f}")
print(f"  ZIP:     {p[0]:.4f}  (excess by {p[0] - poisson_p[0]:.4f})")

검증 포인트

ZIP \(P(y = 0) = \pi + (1 - \pi) e^{-\lambda}\) = 0.3 + 0.7 * 0.135 = 0.395. Poisson \(P(y = 0) = e^{-\lambda}\) = 0.135.

→ ZIP 가 표준 Poisson 보다 zero 확률 0.260 만큼 많음. Excess zeros 의 mathematical expression.

다른 \(y > 0\) 에서:

• ZIP 의 \(P(y = k)\) 가 Poisson 보다 약간 작음 (perfect state 로 weight 일부 빠짐).
• 합산 = 1 보존.

9.2 Step 2: ZIP MLE 직접 구현 (BFGS)

import numpy as np
from scipy.special import expit, factorial
from scipy.optimize import minimize


def negative_log_likelihood_zip(theta: np.ndarray, y: np.ndarray,
                                  X: np.ndarray, W: np.ndarray) -> float:
    """식 (12.14) 의 negative log-likelihood (BFGS 최소화용).

    theta: (beta, gamma) concatenated.
    X: Poisson part 공변량 행렬.
    W: Logistic part 공변량 행렬.
    """
    p = X.shape[1]
    beta = theta[:p]
    gamma = theta[p:]

    lam = np.exp(X @ beta)
    pi = expit(W @ gamma)

    # 식 (12.13) 의 unified density
    poisson_pmf = np.exp(-lam) * lam ** y / factorial(y)
    indicator = (y == 0).astype(float)
    p_y = (1 - pi) * poisson_pmf + indicator * pi

    # log L
    return -np.sum(np.log(p_y + 1e-300))


# 시뮬레이션
np.random.seed(2026)
n = 1000
x = np.random.normal(size=n)
X = np.column_stack([np.ones(n), x])
W = X.copy()  # 같은 covariate

beta_true = np.array([0.5, 0.8])  # log lambda
gamma_true = np.array([-0.5, 1.0])  # logit pi

lam_true = np.exp(X @ beta_true)
pi_true = expit(W @ gamma_true)
is_perfect = np.random.binomial(1, pi_true)
y = np.where(is_perfect == 1, 0, np.random.poisson(lam_true))

# BFGS 적합
theta0 = np.concatenate([np.zeros(2), np.zeros(2)])  # 초기값
result = minimize(negative_log_likelihood_zip, theta0,
                  args=(y, X, W), method='BFGS')

print(f"Converged: {result.success}")
print(f"Estimated beta:  {result.x[:2].round(3)}  (true: {beta_true})")
print(f"Estimated gamma: {result.x[2:].round(3)}  (true: {gamma_true})")

# Hessian-based standard errors
hessian_inv = result.hess_inv
se = np.sqrt(np.diag(hessian_inv))
print(f"SE: {se.round(3)}")

검증

scipy 의 BFGS 가 Greene 의 quasi-Newton 알고리즘과 같은 family. 추정값이 진짜 모수에 가까움 + 표준오차 안정.

R pscl::zeroinfl(y ~ x | x, dist = "poisson") 와 비교 — 같은 결과 나오면 검증 성공.

일반 권고:

• 작은 n: 직접 BFGS (간단한 코드).
• 중간/큰 n: R pscl::zeroinfl 또는 Python statsmodels.discrete.count_model.ZeroInflatedPoisson.
• 매우 큰 n: SAS PROC NLMIXED (parallelization).

9.3 Step 3: EM Algorithm 직접 구현

import numpy as np
from scipy.special import expit, factorial
from scipy.optimize import minimize
from sklearn.linear_model import LogisticRegression
import statsmodels.api as sm


def em_zip(y: np.ndarray, X: np.ndarray, W: np.ndarray,
           max_iter: int = 100, tol: float = 1e-6) -> dict:
    """ZIP 의 EM 알고리즘 (Lambert 1992).

    E-step: z_i (perfect state indicator) 의 사후 확률.
    M-step: weighted GLMs.
    """
    n = len(y)
    p = X.shape[1]
    s = W.shape[1]

    # 초기값
    beta = np.zeros(p)
    beta[0] = np.log(y.mean() + 0.01)
    gamma = np.zeros(s)
    gamma[0] = -1.0  # 초기 pi 작게

    log_lik_old = -np.inf
    for it in range(max_iter):
        # ===== E-step =====
        lam = np.exp(X @ beta)
        pi = expit(W @ gamma)
        poisson_pmf = np.exp(-lam) * lam ** y / factorial(y)
        # z_i 의 사후 확률
        z_post = np.where(y == 0,
                          pi / (pi + (1 - pi) * np.exp(-lam)),
                          0.0)

        # ===== M-step =====
        # 1. Logistic regression: z_post vs W (treating z_post as continuous)
        # → 정확히는 weighted logistic, 근사로 sklearn 사용
        # 더 정확한 구현은 statsmodels GLM with family Binomial
        glm_logistic = sm.GLM(z_post, W, family=sm.families.Binomial())
        gamma_new = glm_logistic.fit(start_params=gamma).params

        # 2. Weighted Poisson regression: y vs X with weight (1 - z_post)
        weight = 1 - z_post
        glm_poisson = sm.GLM(y, X, family=sm.families.Poisson(),
                             freq_weights=weight)
        beta_new = glm_poisson.fit(start_params=beta).params

        # Convergence check
        log_lik = np.sum(np.log((1 - pi) * poisson_pmf
                                + (y == 0) * pi + 1e-300))
        if abs(log_lik - log_lik_old) < tol:
            return {"beta": beta_new, "gamma": gamma_new,
                    "converged": True, "iterations": it + 1,
                    "log_lik": log_lik}
        log_lik_old = log_lik
        beta = beta_new
        gamma = gamma_new

    return {"beta": beta, "gamma": gamma, "converged": False,
            "iterations": max_iter, "log_lik": log_lik}


# 시뮬레이션 데이터에 EM 적용 (Step 2 의 데이터 활용)
em_result = em_zip(y, X, W)
print(f"EM Converged: {em_result['converged']} in {em_result['iterations']} iterations")
print(f"Estimated beta:  {em_result['beta'].round(3)}")
print(f"Estimated gamma: {em_result['gamma'].round(3)}")

EM vs BFGS 의 비교

같은 데이터, 두 알고리즘:

• BFGS (Step 2): 빠름 (수십 iteration), 표준오차 자동.
• EM (Step 3): 느림 (수백 iteration), 표준오차 별도 계산 필요.

EM 의 장점:

• 각 step 이 표준 GLM — 디버깅 쉬움.
• 결정론적 수렴 (likelihood 매 iteration 증가 — 검증 용이).
• 작은 표본 / 어려운 case 에 안정.

BFGS 의 장점:

• 빠른 수렴 (quasi-Newton).
• 표준오차 자동.
• 구현 단순 (scipy minimize).

현대 표준: BFGS (또는 quasi-Newton 변형) — pscl, statsmodels, gamlss 모두. EM 은 mixed-effects ZIP (적분 추가) 에서 더 자주 사용.

9.4 Step 4: ZIP vs ZIP(τ) 비교 (R)

library(pscl)


# 시뮬레이션 데이터
set.seed(2026)
n <- 1000
x <- rnorm(n)

# True ZIP(τ): logit(pi) = -tau * (intercept + beta * x), log(lambda) = intercept + beta * x
intercept_true <- 0.5
beta_true <- 0.8
tau_true <- 0.7

lam <- exp(intercept_true + beta_true * x)
pi <- plogis(-tau_true * (intercept_true + beta_true * x))
is_perfect <- rbinom(n, 1, pi)
y <- ifelse(is_perfect == 1, 0, rpois(n, lam))

df <- data.frame(y = y, x = x)

# ZIP 적합 (별도 covariate)
fit_zip <- zeroinfl(y ~ x | x, data = df, dist = "poisson")
summary(fit_zip)

# ZIP(τ) 는 zeroinfl 에 직접 옵션 없음 — 직접 구현 또는 brms 사용
# 대안: gamlss::gamlss(y ~ x, sigma.formula = ~x, family = ZIP)

# AIC 비교
cat("ZIP AIC:", AIC(fit_zip), "\n")
# ZIP(τ) AIC 가 자유도 적어 약간 더 작을 가능성 (모수 절약)

# Newton-Raphson vs BFGS (zeroinfl 의 method 옵션)
fit_zip_nr <- zeroinfl(y ~ x | x, data = df, dist = "poisson",
                       control = zeroinfl.control(method = "L-BFGS-B"))
fit_zip_bfgs <- zeroinfl(y ~ x | x, data = df, dist = "poisson",
                          control = zeroinfl.control(method = "BFGS"))

# 두 method 의 추정값 비교 (거의 같음)
print(coef(fit_zip_nr))
print(coef(fit_zip_bfgs))

알고리즘 선택의 실무

R pscl::zeroinfl 의 method 옵션:

• "BFGS" (기본): 표준 BFGS.
• "L-BFGS-B": 메모리 효율 BFGS (큰 모형).
• "CG": Conjugate gradient (작은 모형, 빠름).

ZIP(τ) 적합:

• pscl 직접 지원 안 함.
• gamlss 패키지 (Stasinopoulos & Rigby) — family = ZIP, sigma.formula = ~x 로 functional 관계.
• brms::brm(y ~ x, family = zero_inflated_poisson(link = "log", link_zi = "logit"), ...) — Bayesian.

실무 권고:

• 별도 covariate ZIP: pscl::zeroinfl (frequentist 표준).
• ZIP(τ): gamlss 또는 직접 구현.
• 큰 표본 + 여러 random effects: Bayesian (brms, Stan).

10 관련 주제

선행 지식

• Ch.12 Overview — Counts GLMM 의 큰 그림 (식 12.7-12.12)
• § 12.1 ~ 12.2 — Poisson + Modified Poisson family 분류
• § 9.5.5 c-log-log — 비대칭 link 함수
• § 9.6 Marginal MLE — Mixed-effects 추정 framework

후속 주제 (Ch.12 sub-posts)

• § 12.4 — Mixed-effects Poisson + ZIP (식 12.19-12.33)
• § 12.5 — Suicide rate illustration (Gibbons et al. 2005)
• § 12.6 — Ch.12 Summary

관련 개념

• Lambert (1992) — ZIP / ZIP(τ) 정립
• Greene (1994) — Unified pdf + ZINB + gradient method
• Dempster, Laird & Rubin (1977) — EM 알고리즘 원전
• Berndt, Hall, Hall & Hausman (1974) — BHHH estimator
• Zorn (1996) — Dual regime data generating process 통합
• McCullagh & Nelder (1989) — GLM canonical link
• Hall (2000), Hur et al. (2002) — Mixed-effects ZIP 부분적 확장
• Hedeker & Gibbons (2006, § 12.4.3) — 완전한 mixed-effects ZIP
• Statistics 음이항 분포 — ZINB 의 NB 부분