1 GLMM (3): 포아송 및 음이항 혼합 모델
1.1 카운트 데이터의 특성
카운트 데이터는 세 가지 특성을 가진다:
- 비음수: 0 이상의 정수
- 오른쪽 꼬리: 대부분 작은 값, 가끔 큰 값
- 평균-분산 관계: 포아송 → \(\text{Var}=\mu\), 실제로는 \(\text{Var}>\mu\) (과산포)
| 결과 변수 | 예시 | 분포 선택 |
|---|---|---|
| 세션당 턴 수 | 0, 1, 2, 3, … | 포아송 (과산포 시 NB) |
| 클릭 수 | 0, 0, 0, 3, 0, … | Zero-inflated 포아송 |
| 하루 방문 횟수 | 1, 2, 5, 10, … | 포아송 또는 NB |
| 에러 발생 횟수 | 0, 0, 1, 0, … | Zero-inflated 포아송 |
1.2 포아송 혼합 모델
1.2.1 수식
\[\log(E[Y_{ij} \mid u_i]) = \beta_0 + \beta_1 X_{ij} + u_i\]
\[Y_{ij} \mid u_i \sim \text{Poisson}(\lambda_{ij}), \quad \lambda_{ij} = \exp(\beta_0 + \beta_1 X_{ij} + u_i)\]
\[u_i \sim N(0, \sigma^2_u)\]
계수 해석: \(e^{\beta_1}\) = 비율비(Rate Ratio, RR) — X가 1 증가할 때 카운트 기대값이 몇 배가 되는가
1.2.2 예시: 세션당 턴 수 분석
import numpy as np
import pandas as pd
np.random.seed(42)
n_users, n_weeks = 500, 8
user_ids = np.repeat(range(n_users), n_weeks)
weeks = np.tile(range(1, n_weeks + 1), n_users)
personalized = (weeks >= 5).astype(int)
# 사용자별 랜덤 절편 (대화 참여도 개인차)
u_i = np.repeat(np.random.normal(0, 0.45, n_users), n_weeks)
log_lambda = 0.70 + 0.20 * personalized + 0.02 * weeks + u_i
turn_count = np.random.poisson(np.exp(log_lambda))
df = pd.DataFrame({
"user_id": user_ids, "week": weeks,
"personalized": personalized, "turn_count": turn_count,
"segment": np.repeat(
np.random.choice(["SI","MIEP","N"], n_users, p=[0.6,0.3,0.1]),
n_weeks
)
})
print(df["turn_count"].describe())
print(f"\n개인화 전 평균 턴: {df[df.personalized==0]['turn_count'].mean():.2f}")
print(f"개인화 후 평균 턴: {df[df.personalized==1]['turn_count'].mean():.2f}")count 4000.000
mean 2.351
std 1.701
min 0.000
25% 1.000
50% 2.000
75% 3.000
max 14.000
dtype: float64
개인화 전 평균 턴: 2.04
개인화 후 평균 턴: 2.72library(lme4)
m_pois <- glmer(
turn_count ~ personalized + week + segment + (1 | user_id),
data=df, family=poisson(link="log")
)
summary(m_pois)Random effects:
Groups Name Variance Std.Dev.
user_id (Intercept) 0.197 0.444
Fixed effects:
Estimate Std.Error z-value Pr(>|z|)
(Intercept) 0.678 0.044 15.4 <0.001
personalized 0.198 0.029 6.8 <0.001 ← log RR
week 0.021 0.007 3.0 0.003
segmentMIEP 0.482 0.052 9.3 <0.001
segmentN -0.162 0.071 -2.3 0.022# 비율비 (Rate Ratio)
exp(fixef(m_pois))
# personalized: exp(0.198) = 1.22
# → 개인화 시 세션당 평균 턴 수 22% 증가
# → (같은 사용자 기준, Conditional RR)1.3 과산포 (Overdispersion) 문제
포아송의 가정: \(E[Y] = \text{Var}[Y] = \lambda\)
실제 데이터는 대부분 \(\text{Var}[Y] > E[Y]\) (과산포).
1.3.1 과산포 탐지
DHARMa dispersion test:
dispersion = 1.42, p = 0.003
→ 유의한 과산포 존재
→ 음이항 모델로 전환 필요과산포 원인: - 측정되지 않은 이질성 (관측되지 않은 개인차) - Zero가 예상보다 많음 (Zero-Inflation) - 군집 내 상관 (랜덤 효과로 일부 해결)
1.4 음이항 혼합 모델 (Negative Binomial GLMM)
과산포를 직접 모델에 포함:
\[Y_{ij} \mid u_i \sim \text{NB}(\mu_{ij}, \theta)\]
\[\log(\mu_{ij}) = \beta_0 + \beta_1 X_{ij} + u_i\]
\(\theta\): 분산 파라미터. \(\text{Var}[Y] = \mu + \mu^2/\theta\) (\(\theta \to \infty\) 이면 포아송)
library(glmmTMB) # lme4보다 NB 모델에 강력
m_nb <- glmmTMB(
turn_count ~ personalized + week + segment + (1 | user_id),
data=df,
family=nbinom2 # NB Type II: Var = μ + μ²/θ
)
summary(m_nb)Random effects:
Groups Name Variance Std.Dev.
user_id (Intercept) 0.192 0.438
Conditional model:
Estimate Std.Error z-value Pr(>|z|)
(Intercept) 0.681 0.045 15.1 <0.001
personalized 0.199 0.031 6.4 <0.001
week 0.020 0.008 2.5 0.012
segmentMIEP 0.479 0.055 8.7 <0.001
segmentN -0.159 0.074 -2.1 0.034
Dispersion:
theta (분산 파라미터) = 4.82
→ 과산포 보정됨1.5 Zero-Inflated 모델
데이터에 0이 지나치게 많은 경우:
일반 포아송 기대 0 비율: 20%
실제 0 비율: 45% → Zero-Inflation원인: “구조적 0” (절대 전환 안 하는 사용자) + “카운트 0” (이번 주에 우연히 0)
library(glmmTMB)
# Zero-Inflated 포아송
m_zip <- glmmTMB(
turn_count ~ personalized + week + segment + (1 | user_id),
ziformula = ~1 + segment, # 구조적 0의 확률 모델
data=df,
family=poisson
)
# Zero-Inflated 음이항
m_zinb <- glmmTMB(
turn_count ~ personalized + week + segment + (1 | user_id),
ziformula = ~1 + segment,
data=df,
family=nbinom2
)
# 모델 비교
AIC(m_pois, m_nb, m_zip, m_zinb)AIC 비교:
m_pois: 8234 ← 기본 포아송
m_nb: 8108 ← 음이항 (과산포 보정)
m_zip: 8142 ← Zero-inflated 포아송
m_zinb: 8095 ← Zero-inflated 음이항 (최선)
→ m_zinb 선택1.6 Offset: 노출량 차이 보정
측정 기간이나 노출량이 사람마다 다를 때 offset으로 보정.
예시: 세션 길이가 사람마다 다른데 턴 수를 비교
# 세션 길이(분)을 offset으로 추가
df$log_duration <- log(df$session_duration_min)
m_offset <- glmer(
turn_count ~ personalized + week + segment +
offset(log_duration) + # 세션 길이 보정
(1 | user_id),
data=df, family=poisson
)
# 이제 계수 해석: "같은 세션 길이 대비 턴 수 비율"1.7 포아송 vs 음이항 vs ZIP 선택 가이드
Step 1: 포아송 혼합 모델 적합
↓
Step 2: 과산포 검정 (DHARMa::testDispersion)
↓ 과산포 있음
Step 3: 음이항 혼합 모델 적합
↓
Step 4: Zero-Inflation 검정 (DHARMa::testZeroInflation)
↓ Zero-Inflation 있음
Step 5: Zero-Inflated 음이항 (ZINB) 적합
↓
Step 6: AIC/BIC로 최종 모델 선택1.8 Python에서의 카운트 GLMM
# Python에서 포아송 GLMM은 제한적
# statsmodels는 Poisson GLMM 미지원
# 대안 1: rpy2로 R lme4 호출
import rpy2.robjects as ro
from rpy2.robjects import pandas2ri
pandas2ri.activate()
ro.globalenv['df_r'] = df
ro.r('''
library(lme4)
m <- glmer(turn_count ~ personalized + week + segment + (1|user_id),
data=df_r, family=poisson)
print(summary(m))
''')
# 대안 2: bambi (베이지안 GLMM, Python)
import bambi as bmb
model = bmb.Model(
"turn_count ~ personalized + week + segment + (1|user_id)",
df,
family="poisson"
)
results = model.fit(draws=1000, chains=2)1.9 결과 보고 형식
[방법]
세션당 턴 수(카운트)의 개인화 효과를 분석하기 위해
사용자별 랜덤 절편을 포함한 음이항 혼합 모델을 적용하였다.
과산포 검정 결과(dispersion=1.42, p=0.003) 음이항 분포를 선택하였다.
[결과]
개인화 적용 시 세션당 평균 턴 수가 1.22배 증가하였다
(Conditional RR = 1.22, 95% CI [1.15, 1.29], p < 0.001).
사용자 간 기본 참여도 차이도 유의하였다(σ²_u = 0.192).
MIEP 세그먼트는 SI 대비 1.61배 더 많은 턴을 보였다(RR = 1.61, p < 0.001).다음: 09-mixed-model-gam-intro.qmd — GAM/GAMM: 비선형 관계와 종단 데이터