[R] 빈도 데이터 분석을 위한 포아송 회귀모형과 대안 회귀모형에 대한 소개

이번 포스팅에서는 STEFANY COXE, STEPHEN G. WEST, AND LEONA S. AIKEN, 2009, "The Analysis of Count Data: A Gentle Introduction to Poisson Regression and Its Alternatives" 논문을 요약해서 한글로 번역한 내용을 소개하겠습니다. 

 

빈도(count data) 데이터에 대한 분석을 할 때 어떤 분포를 가정하고, 어떤 회귀모형을 적합하는 것이 좋을지 분석 방향을 잡는데 도움이 될 것입니다. 

 

 

 

 

---

왜 포아송 회귀모형이 필요한가(Why Poisson Regression)?

 

특정 시간 동안 발생한 사건의 건수에 대한 도수 자료(count data, 음수가 아닌 정수)를 목표변수로, 특히 평균이 10 미만일 경우, 최소제곱법(Ordinary Least Squares) 회귀모형을 적합하면 표준 오차와 유의수준이 편향되는 문제가 발생한다. OLS 회귀모형은 오차의 조건부 정규성(conditional normality), 등분산성(homoscedasiticity), 독립성(independence)을 가정한다. 하지만 평균이 작은 도수 데이터는 0 이상의 정수만 있고, 작은 계급에 많은 관측치가 몰려있으며, 우측으로 길게 치우친 분포를 띠어 회귀모형의 가정을 위배하고 음수값도 결과값으로 반환하는 문제가 있다. 종속변수가 도수 데이터이면서 오차가 정규분포, 등분산이 아닌 경우 포아송 회귀모형을 사용한다.

 

 

• 포아송 회귀모형은 무엇인가? (Overview of Poisson Regression)

포아송 회귀모형은 일반화선형모형(Generalized Linear Model, GLM)의 한 종류로서,

 

-       Random component: Poisson Distribution, 

 

 

-       Systematic component: Mixed linear relationship

 

-       Link function: Log Link (ln)

 

 

 

표현 예: 

 

해석 : 

 이며, X₁을 한 단위 증가시켰을 때 효과를 보면 

  이므로, 다른 변수들을 통제했을 때 X₁이 한단위 증가하면 도수의 추정치

 만큼 승법적으로 변화(multiplicative changes)한다.

 

 

선형회귀모형이 최소제곱법(OLS)로 모수를 추정한다면, 포아송 회귀모형은 최대가능도추정(MLE, Maximum Likelihood Estimation)을 통해 모수를 추정한다.

 

 

선형회귀모형은 모델에 의해 설명되는 잔차 제곱합 SS (Sum of Squares)를 총 잔차제곱합(TSS, Total Sum of Squares)로 나눈 값인  R^2로 모델 적합도(Goodness of Fit of the model) 평가하지만, 포아송 회귀모형은 설명변수의 추가에 따른 이탈도의 감소 비율(proportional reduction in deviance)로 표현되는 pseudo-R^2를 사용하여 완전모델(perfect model)에 얼마나 가깝게 적합이 되었는지를 평가한다. 

 

즉, 

 

선형회귀모형은 F-test를 통해 모델의 유의성을 검정하지만, 포아송 회귀모형은 base model과 complete model 간의 이탈도 차이에 대한 Chi-square test를 통해 모델 유의성을 검정한다.  

선형회귀모형은 모형 타당성 평가 (Assessing model adequacy)를 위해 잔차도(residual plot)를 그려서 확인해보는 반면, 포아송 회귀모형은 이탈도 잔차도(deviance residuals vs. predicted values)를 그려보거나 실측치와 예측치를 비교해보는 것이다.

 

 

* [참고] 포아송 분포 : https://rfriend.tistory.com/101

 

 

과대산포(Over-dispersion) 문제는 무엇이며, 어떻게 해결할 수 있나?

 

포아송분포는 평균과 분산이 동일(equidispersion)하다고 가정하므로, 분산이 평균보다 큰 과대산포 (Overdispersion) 시에는 포아송 회귀모형을 사용할 수 없다. 과대산포 시에 적용할 수 있는 대안 모델로는 과대산포 포아송 회귀(Overdispersed Poisson Regression), 음이항회귀(Negative Binomial Regression)가 있다.

 

 

-   Overdispersed Poisson Regression: 과대산포를 모델링하기 위해 두번째 모수인 조건부 분산 overdispersion scaling parameter

를 추정하며, 과대산포 포아송 회귀모형의 오차 분포는 

 를 가진다.

 

 

-   Negative Binomial Regression Models: 음이항회귀모형은 평균에는 영향이 없으면서 과대산포를 유발하는, 설명이 되지 않는 추가적인 가변성이 있다고 가정한다. 음이항 회귀모형은 똑같은 설명변수 값을 가지는 관측치가 다른 평균 모수를 가지고 포아송 회귀모형에 적합될 수 있도록 해준다. 오차 함수는 포아송 분포와 감마분포의 혼합 분포이다.

 

 

과대산포 검정은 가능도비 검정(Likelihood ratio test) 또는 score 검정을 사용한다.

 

모델 선택을 위해서는 AIC(Akaike Information Criterion), BIC(Bayesian Information Criterion) 을 사용하며, AIC, BIC의 값이 작은 모델을 선택하는데, 샘플 크기가 작을 때는 간단한 모델을 선택하고, 샘플 크기가 클 때는 좀더 복잡한 모델을 선택한다.

 

 

 

과대영(Excess Zeros) 문제는 무엇이며, 어떻게 해결할 수 있나?

 

특정 평균을 모수로 가지는 포아송 분포에서 나타나는 ‘0’보다 많은 ‘0’을 가지는 샘플의 경우 과대영(excess zeros)이 발생했다고 말한다. 과대영은 (a) 특정 행동을 나타내지 않은 그룹의 포함 (예: 비흡연자 그룹), (b) 특정 행동을 나타내지만 샘플링 계획에서 제외 (예: 병원 입원한 사람만 연구) 하는 경우에 발생할 수 있다. 과대영일 경우 ZIP 회귀모형을 대안으로 사용할 수 있다.

 

-   Zero-Inflated Poisson(ZIP) Regression: 두 부분으로 나뉘는데, (1) 로지스틱 회귀모형을 사용하여 항상 ‘0’인 집단(structural zeros) 또는 항상 ‘0’은 아니지만 조사 시점에 ‘0’이라고 응답한 집단에 속할 확률을 구하고, (2) 포아송 회귀 또는 음이항 회귀모형으로 structural zeros를 포함하지 않은 나머지 관측치를 추정하는 모델을 적합한다.

 

 

---

음이항분포(Negative Binomial Distribution)

 

* Reference : https://en.wikipedia.org/wiki/Negative_binomial_distribution

 확률 이론과 통계에서, 음 이항 분포는 특정 (무작위) 실패 횟수 (r로 표시)가 발생하기 전에 독립적이고 동일하게 분산 된 베르누이 시행 순서에서 성공 횟수의 이산 확률 분포이다. 예를 들어 1을 실패로, 1이 아닌 모든 것을 성공으로 정의하고 1이 세 번째로 나타날 때까지 반복적으로 던지면 (r = 3 번의 실패), 출현한 1이 아닌 수(non-1)의 확률 분포는 음 이항 분포가 된다.

 

 파스칼 분포 (Blaise Pascal 이후)와 Polya 분포 (George Pólya의 경우)는 음 이항 분포의 특수한 경우이다. 엔지니어, 기후 학자 및 다른 사람들 사이의 협약은 정수 값 정지 시간 매개 변수 r의 경우 "음 이항" 또는 "파스칼"을 사용하고, 실수 값의 경우에는 "Polya"를 사용한다. 토네이도 발생과 같은 "전염성 있는" 이산 이벤트가 발생하는 경우 Polya 분포를 사용하여 Poisson과 달리 평균 및 분산을 다르게 하여 Poisson 분포보다 더 정확한 모델을 제공 할 수 있다. "전염성 있는" 사건은 양의 공분산으로 인해 사건이 독립적이었던 경우보다 양의 상관관계가 있어 Poisson 분포에 의한 분산보다 더 큰 분산을 일으킨다. Poisson 분포는 평균과 분산이 동일하다고 가정하기 때문에 평균보다 분산이 큰 분포(Overdispersed)의 경우 적절하지 않으며, 2개의 모수를 가져서 과대산포분포를 적합할 수 있는 음 이항 분포를 사용할 수 있다.

 

 일련의 독립적 인 베르누이 재판이 있다고 가정하자. 따라서 각 임상 시험은 "성공"과 "실패"라는 두 가지 결과를 나타낸다. 각 시험에서 성공 확률은 p이고 실패 확률은 (1 - p)이다. 미리 정의 된 실패 횟수 r이 발생할 때까지의 서열을 관찰한다. 그러면 우리가 본 성공의 무작위 수 X는 음 이항 (또는 파스칼) 분포를 가질 것이다. 

 

 

---

R을 활용한 Poisson model, Negative binomial model, Zero-Inflated model

 

Poisson model, Negative binomial model, Zero-Inflated model 의 유형, 분포, 모수 추정 방법, R패키지와 함수를 표로 정리하면 아래 [표1]을 참고하세요. 

 

* Reference: http://cran.at.r-project.org/web/packages/pscl/vignettes/countreg.pdf

 

 

[표1] Overview of count regression models

 

유형	분포	추정 방법	Description	R package & function
GLM	Poisson	ML	Poisson Regression: classical GLM, estimated by maximum likelihood(ML)	{stats} 패키지의 glm()
	NB	ML	NB Regression: exteded GLM, estimated by ML including additional shape parameter	{stats} 패키지의 glm(), {MASS} 패키지의 glm.nb() 함수
Zero-augmented	Poisson	ML	Zero-Inflated Poission(ZIP)	{pscl} 패키지의 zeroinfl() 함수

 

각 회귀모형별로 R 패키지와 함수의 활용 방법은 아래와 같다.

-       (1) R을 활용한 Poisson model : glm() in the {stats} package

 

glm(formula, data, subset, na.action, weights, offset,     family = poisson, start = NULL, control = glm.control(...),     model = TRUE, y = TRUE, x = FALSE, ...)

 

 

-       (2) R을 활용한 Negative binomial model
   : glm() in the {stats} package

 

glm(formula, data, subset, na.action, weights, offset,     family = negative.binomial, start = NULL, control = glm.control(...),     model = TRUE, y = TRUE, x = FALSE, ...)

 

 

   : glm.nb() in the {MASS} package

 

glm.nb(formula, data, weights, subset, na.action,        start = NULL, etastart, mustart,        control = glm.control(...), method = "glm.fit",        model = TRUE, x = FALSE, y = TRUE, contrasts = NULL, ...,        init.theta, link = log)

 

-       (3) R을 활용한 Zero-inflated regression model: zeroinfl() in the {pscl} package

 

zeroinfl(formula, data, subset, na.action, weights, offset, dist = "poisson", link = "logit", control = zeroinfl.control(...), model = TRUE, y = TRUE, x = FALSE, ...)

 

 

많은 도움이 되었기를 바랍니다.