베이지안 통계로 보는 2022 대선 결과 예측 (Bayesian Hierarchical Modeling)

해당 게시글은

(1) Bayesian Hierarchical modeling의 개념

(2) Bayesian Hierarchical modeling 활용 사례 - 2022 대선 결과 예측

에 대한 내용을 담고 있습니다.

 

---

Chapter #1 - Bayesian Hierarchical modeling의 개념

Bayesian Hierarchical modeling은 무엇이고, 왜 쓰는 것인가요?

Photo of Bonnie Kittle

Hierarchical modeling 이란, 알고 싶은 파라미터를 여러 계층으로 구분하여 추정하는 방법론을 의미합니다.

베이즈 정리를 활용하여 사전 (prior) 분포를 가정하고, 사후 (posterior) 분포를 추정하는 베이지안 방법론과 같이 쓰일 때 Bayesian Hierarchical modeling이라 합니다. 일반적으로 실험이나 여론조사에서 서로 다른 집단에 대해 여러 차례 진행된 조사 결과를 종합하기 위한 방식으로 사용합니다. 한 번 실시된 조사보다는 여러 차례 실시된 조사를 종합하면 더 신뢰할 수 있는 데이터를 얻을 수 있겠죠.

 

예를 들어 실험실에서 어떤 신약을 개발하였고, 신약이 암의 발병 확률을 낮출 가능성이 있다고 가정해 봅시다.

암은 신약 처치 여부와 관련 없이 발생할 수 있기 때문에, 아무런 처치도 하지 않은 경우에 암이 발병할 확률을 먼저 알아야 합니다.

따라서 아무런 처치도 하지 않았을 때 암이 어느정도 확률로 발생하는지 확인하기 위해, 쥐를 대상으로 여러 차례 실험을 진행 하였습니다.

* 데이터 출처 : Bayesian Data Analysis - Rat Tumor

 

총 71번의 서로 독립적인 실험을 진행 하였고, 실험 결과는 아래와 같이 표시합니다.

 

• $ y_j $ : j번째 실험에서 암이 발병한 쥐의 수
• $ n_j $ : j번째 실험에서 전체 쥐의 수
• $ j = 1,2, ..., 71 $
• 데이터 형태 : $(y_j, n_j)$
  • (0, 20), (0, 18), (4,20), ...

각 실험에서 쥐의 수는 $n_j$ 마리이고, 암이 발병한 쥐의 수는 $y_j$ 마리입니다.

각 실험에서 암이 발병할 확률을 $\theta_j$ 라고 하면, $y_j$ 는 이항분포를 따르기 때문에 아래와 같이 작성할 수 있습니다.

* 이항분포 : 연속된 n번의 독립적인 시행에서 각 시행이 확률 p를 가질 때의 이산 확률 분포

 

$$ y_j | \theta_j \sim Bin(n_j, \theta_j) $$

 

그럼 서로 독립인 여러 번의 실험 결과를 어떻게 종합할 수 있을까요? 

 

Option #1 : 모든 실험이 공통의 파라미터를 갖는다고 가정

첫 번째 방법은 모든 실험에 대해 암이 발병할 확률은 동일하다고 가정하는 것 입니다. 즉 $\theta_j = \theta$ 라고 가정합니다.

이 경우, 모든 실험 결과를 종합해서 아래와 같이 작성할 수 있습니다.

 

$$ y_1 + y_2 + ... + y_{71}  | \theta \sim Bin(n_1 + n_2 + ... + n_71, \theta) $$

 

물론 이렇게 단순하게 가정할 경우, 한계가 명확합니다.

• 실제 데이터에서는 암이 발병한 비율이 실험마다 다름. 즉 확률 모델과 데이터가 불일치
• 서로 다른 여러 실험 간의 variability를 다룰 수가 없음

 

Option #2 : 각 실험마다 서로 독립인 파라미터를 갖는다고 가정

두 번째 방법은 각 실험마다 서로 독립인 발병 확률을 갖는다고 가정하는 것 입니다. 

이 경우 실험마다 암이 발병한 비율이 다르게 나타난 데이터를 설명할 수는 있지만, 마찬가지로 한계가 있습니다.

• 각 실험을 구성하는 샘플 크기가 작음 -> 사후 확률의 분산이 커짐
• 새로운 experiment 에 대한 답을 할 수가 없음 (모든 실험은 다르기 때문에)

 

Option #3 : 계층 구조를 갖는 파라미터를 가정 (Hierarchical modeling)

마지막은 Option #1과 Option #2를 절충하는 안 입니다.

즉 $\theta_1$ 부터 $\theta_71$은 서로 다르지만, 공통적인 무언가를 가지고 있을 거라는 아이디어 인데요,

이 경우 $\theta_j$ 가 어떤 분포를 가지고 있고, 해당 분포에서 각 실험을 규정하는 파라미터인 $\theta_j$가 뽑혔다고 가정합니다.

 

계층 구조를 가지는 파라미터를 가정하면 위와 같은 모델링이 가능하며,

이 경우 Option #1과 Option #2에서는 없었던 장점이 있습니다.

• 일반화 가능
• 서로 다른 여러 실험 간의 variablility를 다룰 수 있음
• 다른 실험 결과를 활용하여 각 $y_j$ 에 대한 예측 정확도를 높일 수 있음

그럼 $\theta_j$ 가 따르는 분포는 무엇으로 가정하면 좋을까요? 바로 베타분포 입니다.

$\theta_j$ 와 관련된 제약 조건들이 있는데, 베타 분포를 사용하면 해당 제약 조건들을 만족시킬 수 있기 때문입니다.

• $\theta_j$ 는 이항 분포의 확률을 규정하는 값이기 때문에 0에서 1 사이 값을 가져야 함
  • 베타분포는 0에서 1 사이 값을 가짐
•  이항 분포와 결합하여 사후 분포를 예측하기 쉬운 분포여야 함
  
  • 베타분포는 이항 분포의 conjugate prior

* conjugate prior : prior과 posterior 가 같은 probability distribution family를 이루게 하는 prior 분포

 

따라서 아래와 같은 문제 세팅이 가능하며,

이와 같은 문제 세팅을 Hierarchical (=multi-level) model 이라 합니다.

 

• Lower level : $y_j | \theta_j \sim Bin(n_j, \theta_j) $
• Higher level : $\theta_j \sim Beta(\alpha, \beta)$
• $\alpha, \beta \sim iid \; Exp(\lambda)$

DAG Model

 

여기서 Higher level을 구성하는 파라미터인 $\alpha$ 와 $\beta$를 hyperparameters 라고 하며,

$\alpha$ 와 $\beta$ 도 분포를 갖습니다. 이들의 사전 분포를 hyperprior 이라 합니다.

 

위 예시에서는 $\alpha$ 와 $\beta$가 각각 파라미터가 $\lambda$ 인 지수분포를 따르는 것으로 가정 하였는데요,

$\alpha$ 와 $\beta$ 에 대한 정보가 전혀 없으므로 too informative 한 분포가 아니라면 사용이 가능합니다.

사후 분포가 proper distribution 이 된다는 가정 하에 improper distribution 도 사용이 가능한 경우가 있으며, 위의 예시에서는 지수 분포를 가정 했으므로 $\lambda$가 0에 가까울수록 평탄한 (not informative) 분포가 됩니다.

 

---

Chapter #2  - 2022 대선 결과 예측

 

제20대 대통령 선거 벽보 (출처 : 나무위키)

 

그럼 위에서 소개한 Bayesian Hierarchical modeling 을 활용하여 2022년 3월 9일에 치러진 제20대 대통령 선거 결과를 예측해 보겠습니다. 제20대 대통령 선거는 막판까지 접전 양상을 보였고, 최종적으로 1-2위가 0.73%p 득표율 차이로 역대 최소 득표율 차이를 기록한 선거 입니다.

* 데이터 출처 : 제20대 대통령 선거

 

 

1. 분석 개요

• 2022년 2월 마지막 주, 각 조사 기관의 여론조사 결과
  
  • 2022년 2월 마지막 주 결과가 없는 경우, 3월 결과로 대체
  • 지지율 상위 4자 이상을 대상으로 실시된 여론조사
• 예측하고자 하는 것 (y) : 당시 지지율 top2를 기록하던 이재명 후보와 윤석열 후보의 득표율 차이
• 사용한 프로그래밍 툴 : R의 rjags 패키지
  • (참고) 파이썬의 경우 pyjags 라이브러리에서 유사한 기능을 제공합니다.

 

분석에 사용한 여론조사 데이터는 아래와 같습니다.

이재명 지지율 (%)	윤석열 지지율 (%)	오차 범위 (%)	조사 기관	주차
39.8	39.8	2.2	KBS	22년 2월 2차
39.6	41.9	3.1	MBC	22년 2월 2차
34.1	42.4	3.1	JTBC	22년 2월 2차
34.9	36.5	3.1	TV조선	22년 2월 3차
36.4	43.3	3.1	채널A	22년 2월 2차
31.6	36.1	3.1	넥스트리서치	22년 2월
37	39	3.1	NBS	22년 2월 4주
38	37	3.1	한국갤럽	22년 2월 4주
38.3	39	3.1	머니투데이	22년 2월 4주
40.5	41.9	2.2	리얼미터	22년 2월 4주 주중
43.8	36.1	3.1	KSOI	22년 2월 4주
39.4	40.2	3.1	엠브레인(중앙일보)	22년 2월 2차
40.2	42.4	3.1	엠브레인(news1)	22년 2월 4주
42.5	46.5	1.5	PNR(뉴데일리)	22년 3월
41	43.8	3.1	PNR(프라임경제)	22년 2월 4주
42.3	45.4	1.8	여론조사공정	22년 2월 4주
42	44.2	2.6	미디어토마토	22년 2월 4주
41	46	3.1	리서치뷰	22년 2월 4주
40.9	43.6	3.1	한길리서치	22년 3월
42.1	43.6	3.1	조원씨앤아이	22년 2월 3주
42.2	43.2	3.1	미디어리서치	22년 2월 4주
40	40.4	3.1	서던포스트	22년 2월 4주
39.5	44	3.1	코리아정보리서치	22년 2월 4주

 

2. 시각화

여론조사 기관들의 데이터를 그래프로 표현하면 아래와 같습니다.

 

with(d[, c('organization','y', 'ME')],
	 plot(1:n, y, ylim=c(min(y-ME), max(y+ME)), xaxt="n",
        	main="이재명 후보 지지율 (%) - 윤석열 후보 지지율 (%)",
        	xlab="", ylab="차이"))
with(d[, c('organization','y', 'ME')], segments(1:n, y-ME, 1:n, y+ME))
axis(1, at=1:n, labels=d$organization, las=2)
abline(h=0, lty=2, col="blue")

 

제20대 대통령 선거 - 여론조사 결과

 

• KSOI : 이재명 후보의 우세, 오차 범위 밖
• JTBC, 채널A, 리서치뷰, 정보리서치 : 윤석열 후보의 우세, 오차 범위 밖
• 나머지 여론조사들은 윤석열 후보의 우세가 많았으나 오차 범위 안

y = "이재명 후보 지지율 - 윤석열 후보 지지율" 의 평균 값은 -2.139 입니다.

 

 

3. 문제 세팅

• $y_j$ : 이재명 후보가 윤석열 후보를 지지율에서 몇%p 앞서는지
  • "이재명 후보 지지율 - 윤석열 후보 지지율"
• j = 1,2, ... , 23
• 서로 다른 여론조사 끼리는 서로 독립이다.
• $\sigma_j$ : margin of error (오차 범위) 의 1/2
  • $\sigma_j$ 는 고정된 값이고, 알고 있다고 가정합니다.
  • 물론 $\sigma_j$ 에 대해서도 prior distribution을 적용할 수도 있습니다.
• $y_j$ 값은 정규분포를 따른다고 가정 하였습니다.
• $y_j$ 값의 평균인 $\theta_j$ 도 정규 분포를 따른다고 가정 하였습니다.
• $\theta_j$ 의 분포를 규정하는 $\mu$ 와 $\tau$ 에 대해서는 정보가 전혀 없으므로, 최대한 non-informative prior에 가까운 분포를 사용 하였습니다. (범위가 넓은 균일분포)
• Lower level : $y_j | \theta_j \sim N(\theta_j, \sigma_j^2) $
• Higher level : $\theta_j | \mu, \tau \sim N(\mu, \tau^2)$
• $\mu \sim flat \; on \; (-\infty, \infty) \doteq U(-1000, 1000)$
• $\tau \sim flat \; on \; (0, \infty) \doteq U(0, 1000)$

 

DAG Model

4. 문제 풀이

Gibbs sampler를 활용하여 분포를 만족하는 샘플을 생성해 봅시다.

우선 rjags 패키지에서 활용할 수 있도록 위의 문제 세팅을 코드로 구현해 줍니다. 

model {

  for (j in 1:length(y)) {
    y[j] ~ dnorm(theta[j], 1/sigma[j]^2)
    theta[j] ~ dnorm(mu, 1/tau^2)
  }

  mu ~ dunif(-1000,1000)
  tau ~ dunif(0,1000)

}

위 코드를 polls2022_korean.bug 파일로 저장합니다.

 

library(rjags)

initial.vals <- list(
  list(mu=100, tau=0.01),
  list(mu=-100, tau=0.01),
  list(mu=100, tau=100),
  list(mu=-100, tau=100)
)

m1 <- jags.model("polls2022_korean.bug", d, initial.vals, n.chains=4)

# burn-in for 2,500
update(m1, 2500)
x1 <- coda.samples(m1, c("mu", "tau"), n.iter=5000)

rjags 패키지를 불러 오고, $\mu$와 $\tau$의 초깃값을 잡아 줍니다.

샘플 추출을 위한 chain은 4개를 사용 하겠습니다.

mixing이 잘 되도록 하기 위해 최초 샘플 2500개는 burn-in 해 줍니다.  

여론조사들의 분포를 정의하는 평균 $\mu$ 와 표준편차인 $\tau$ 가 가장 중요한 값이기 때문에, 

위 분포 정의에 맞게 뽑은 $\mu$, $\tau$ 샘플을 각각 5천개*4 = 2만개씩 저장합니다.

 

traceplot 과 $\mu$, $\tau$ 의 분포를 그려 보면 아래와 같습니다.

 

 

 

summary(x1)

Iterations = 3501:8500
Thinning interval = 1 
Number of chains = 4 
Sample size per chain = 5000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

      Mean     SD Naive SE Time-series SE
mu  -2.153 0.6695 0.004734       0.007539
tau  2.780 0.5847 0.004134       0.008289

2. Quantiles for each variable:

      2.5%    25%    50%    75%   97.5%
mu  -3.497 -2.582 -2.147 -1.717 -0.8376
tau  1.793  2.366  2.726  3.126  4.1035

참고로 Iterations 가 3501 부터 시작하는 이유는,

아까 2500개 burn-in 에 이어서 auto-correlation을 줄이기 위한 adaptation 과정이 기본 옵션으로 들어갔기 때문입니다. 

adaptation 과정에서 1000개 sample이 추가적으로 버려집니다.

 

$\mu$ 값의 평균은 -2.153, 중간값은 -2.147 입니다.

이는 최초 여론조사 기관들의 "이재명 지지율 - 윤석열 지지율" 의 평균이 -2.139 였던 것과 일맥상통하는 결과입니다.

평균을 보면 윤석열 후보가 우세할 확률이 가능성이 높지만, $\tau$ 값이 꽤 크기 때문에 실제 결과는 변동성이 클 것으로 예측이 됩니다.

 

$\mu$의 95% 사후 확률 구간은 -3.497 ~ -0.8376 입니다. 

이제 뽑은 2만개 $\mu$, $\tau$ 값들을 가지고 대선 결과를 예측해 보겠습니다.

 

post <- as.matrix(x1)
Nsim <- dim(post)[1]
ys <- vector(length=Nsim)

for (s in 1:Nsim){
  now_mu <- post[s, 1]
  now_tau <- post[s, 2]
  now_theta <- rnorm(1, now_mu, now_tau)
  now_y <- rnorm(1, now_theta, 1/2*(3.1)) # 각 여론조사 기관의 오차범위는 3.1%로 가정
  ys[s] <- now_y
}

# 이재명 후보 당선 확률
> print(sum(ys > 0) / length(ys))
[1] 0.25025
# 윤석열 후보 당선 확률
> print(sum(ys < 0) / length(ys))
[1] 0.74975

이재명 후보의 당선 확률은 약 0.25 (25%), 윤석열 후보의 당선 확률은 약 0.75 (75%) 입니다.

 

---

Summary

이번 포스팅에서는 Bayesian Hierarchical modeling 에 대해 알아 보았고, 예시로 제20대 대통령선거 결과를 예측해 보았습니다.

 

Bayesian Hierarchical modeling은 실험이나 여론조사와 같은, 서로 다른 집단에 대해 여러 차례 진행된 조사 결과를 종합할 때 유용한 도구입니다. 여러 조사들의 평균만 구매 보는 것 보다 훨씬 풍부한 분석 결과를 얻으실 수 있으니, 관심 분야의 데이터를 모아 직접 분석을 수행해 보셔도 좋겠습니다.

 

잘못된 내용에 대한 지적이나 문의사항은 편하게 댓글로 남겨주세요! 읽어 주셔서 감사합니다.

 

---

References

[1] Gelman, A., Carlin, J.B., Stern, H.S., Dunson, D.B., Vehtari, A., & Rubin, D.B. (2013). Bayesian Data Analysis (3rd ed.). Chapman and Hall/CRC. https://doi.org/10.1201/b16018

Bayesian Data Analysis | Andrew Gelman, John B. Carlin, Hal S. Stern,

[2] Park, Trevor H. “Hierarchical modeling fundementals”, Advanced Bayesian Modeling, 2022 spring, University of Illinois Urbana Champaign, Lecture.

 

[3] 제20대 대통령 선거 (나무위키)