베이지안 포아송 회귀(Bayesian poisson regression)

단순 포아송 회귀(Simple poisson regression) 모델은 다음과 같다.

 

 

▷ 가능도는 포아송 분포로 정하고, 가능도 모수의 로그를 취한 것에 대하여 독립변수의 선형결합으로 정의한다. 이때, 위의 식에서는 하나의 독립변수에 대한 선형결합으로 표현하였지만, 독립변수가 여러개이면 이 독립변수들의 선형결합을 통해 가능도 모수의 로그를 취한 값에 대하여 정의한다.

 

▷ 포아송 회귀 모델의 예측은 개별 데이터의 가능도의 평균, 즉, 가능도의 모수를 통해 이루어진다.

 

▷ 베이지안 포아송 회귀(Bayesian poisson regression)는 독립변수와 선형결합된 모수 beta의 사전분포를 정의한다는 점에서 포아송 회귀와 차이가 있다.

 

R의 COUNT 패키지의 badhealth 데이터를 이용하여 베이지안 포아송 회귀 분석을 수행할 것이다. 수행 과정은 다음과 같다.

 

1. 데이터 확인

2. 모델링

3. 모델 확인

4. 모델 예측 및 분석

 

1. 데이터 확인

 

In:

library(COUNT)

data('badhealth')

pairs(badhealth)

 

Out:

 

▷ numvisit과 age는 연속형 변수, badh는 범주형 변수로 데이터가 구성되어 있다.

 

2. 모델링

 

numvisit을 종속변수로, age와 badh를 독립변수로 두어 베이지안 포아송 회귀 분석을 할 것이다.

 

In:

library(rjags)

mod_string = " model {
    for (i in 1:length(numvisit)) {
        numvisit[i] ~ dpois(lam[i])
        log(lam[i]) = b0 + b_badh*badh[i] + b_age*age[i] + b_badh_age*age[i]*badh[i]
    }
    
    b0 ~ dnorm(0.0, 1.0/1e6)
    b_badh ~ dnorm(0.0, 1.0/1e4)
    b_age ~ dnorm(0.0, 1.0/1e4)
    b_badh_age ~ dnorm(0.0, 1.0/1e4)
} "

data_jags = as.list(badhealth)

params = c('b0', 'b_badh', 'b_age', 'b_badh_age')

mod = jags.model(textConnection(mod_string), 
                 data = data_jags, 
                 n.chains = 3)

update(mod, 1e3)

mod_sim = coda.samples(model = mod,
                       variable.names = params,
                       n.iter = 5e3)
mod_comb_sim = as.mcmc(do.call(rbind, mod_sim))

 

Out:

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Graph information:
   Observed stochastic nodes: 1127
   Unobserved stochastic nodes: 4
   Total graph size: 3665

Initializing model

  |++++++++++++++++++++++++++++++++++++++++++++++++++| 100%
  |**************************************************| 100%
  |**************************************************| 100%

 

▷ 독립변수로써 b_badh와 b_age의 교호작용을 고려하여 교차항을 추가하였다.

 

▷ 이외의 코드에 대한 내용은 생략하도록 하겠다. JAGS의 사용법에 대해 알아보고자 한다면, 다음의 포스팅을 참고하길 바란다.

 

JAGS(Just Another Gibbs Sampler) 사용법

 

3. 모델 확인

 

Gelman-Rubin diagostic을 통해 마르코프 체인이 수렴하였는지 확인하여 보자.

 

In:

gelman.diag(mod_sim)

 

Out:

Potential scale reduction factors:

           Point est. Upper C.I.
b0               1.02       1.03
b_age            1.01       1.03
b_badh           1.01       1.02
b_badh_age       1.01       1.02

Multivariate psrf

1.01

 

▷ 모든 모수의 결과가 1로 수렴하고 있다는 것을 확인할 수 있다.

 

4. 모델 예측 및 분석

 

위의 모델을 통해 예측을 하고, 실제 결과와 비교해보자.

 

In:

X = as.matrix(badhealth[, -1])
X = cbind(X, with(badhealth, badh*age))

post_coef = apply(mod_comb_sim, 2, median)

log_lam_hat = post_coef['b0'] + X %*% post_coef[c('b_badh', 'b_age', 'b_badh_age')]
lam_hat = exp(log_lam_hat)

plot(lam_hat, badhealth$numvisit)

 

Out:

 

▷ 모델의 예측결과가 데이터의 변동을 잘 설명하지 못하는 것을 확인할 수 있다. 하지만 lamda hat이 증가할 수록 실제 결과도 어느 정도 증가하는 것을 확인할 수 있다.

 

▷ 포아송 분포의 특성으로 평균과 분산이 모수로 표현된다. 즉, 예측값인 평균이 크면 클수록 분산도 증가하게 되고, 위의 결과에서도 이러한 특성을 확인할 수 있다.

 

▶ 위의 경우와 같이 포아송 회귀가 데이터의 변동을 잘 설명하지 못할 때, 대안으로 음이항 회귀를 사용할 수 있다.

 

모델의 결과에 대하여 분석하여보자.

 

In:

summary(mod_sim)

 

Out:

Iterations = 2001:7000
Thinning interval = 1 
Number of chains = 3 
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
b0          0.340722 0.082551 6.740e-04      0.0054368
b_age       0.008652 0.002115 1.727e-05      0.0001386
b_badh      1.566322 0.177987 1.453e-03      0.0127482
b_badh_age -0.010900 0.004126 3.369e-05      0.0003004

2. Quantiles for each variable:

                2.5%       25%       50%      75%     97.5%
b0          0.178886  0.283849  0.340609  0.39707  0.498929
b_age       0.004571  0.007236  0.008664  0.01009  0.012775
b_badh      1.209814  1.445276  1.567056  1.68999  1.904278
b_badh_age -0.018683 -0.013763 -0.010927 -0.00807 -0.002601

 

▷ b_age와 b_badh의 추정치의 평균과 백분위수 모두 양의 값인 것을 확인할 수 있다. 이는 numvisit의 평균과 이들 변수간의 관계가 양의 관계를 나타낸다고 볼 수 있다. 즉, b_age이나 b_badh가 증가한다면 평균적으로 numvisit도 증가한다고 해석할 수 있다.

 

▷ b_age와 b_badh의 교차항은 음수인 것으로 나타났다. 이는 b_age와 b_badh에 대한 보정값으로 해석할 수 있다.

 

X1(age: 35, badh: 0)과 X2(age: 35, badh: 1)의 사후 예측 분포를 구하고, 둘을 비교하여 보자.

 

In:

x_1 = c(0, 35, 0)
x_2 = c(1, 35, 35)

log_lam_1 = mod_comb_sim[, 'b0'] + mod_comb_sim[, c('b_badh', 'b_age', 'b_badh_age')] %*% x_1
log_lam_2 = mod_comb_sim[, 'b0'] + mod_comb_sim[, c('b_badh', 'b_age', 'b_badh_age')] %*% x_2

lam_1 = exp(log_lam_1)
lam_2 = exp(log_lam_2)

n_sim = length(log_lam_1)

y_1 = rpois(n = n_sim, lambda = lam_1)
y_2 = rpois(n = n_sim, lambda = lam_2)

plot(table(factor(y_1, levels = 0:18))/n_sim, pch = 2, ylab = 'Posterior Prob.', xlab = 'Num. of Visits')
points(table(y_2 + 0.1)/n_sim, col = 'red')

 

Out:

 

▷ 검정색의 분포는 X1, 빨강색의 분포는 X2의 사후예측분포를 의미한다. 두 사후예측분포에서 보다시피 X2의 사후예측평균이 X1의 사후예측평균보다 큰 것을 알 수 있다.

 

X2의 numvisit이 X1의 numvisit보다 클 확률을 구하여 보자.

 

In:

mean(y_2 > y_1)

 

Out:

[1] 0.9157333

 

▷ 사후예측분포로부터 X2의 numvisit이 X1의 numvisit보다 클 확률은 약 91%인 것으로 나타났다.

 

---

Reference:

"Bayesian Statistics: From Concept to Data AnalysisTechniques and Models," Coursera, https://www.coursera.org/learn/bayesian-statistics/.