베이지안 로지스틱 회귀(Bayesian logistic regression)

단순 로지스틱 회귀(Simple logistic regression) 모델은 다음과 같다.

 

 

▷ 가능도는 베르누이 분포로 정하고, 가능도 모수의 로짓(Logit)을 독립변수의 선형결합으로 정의한다. 이때, 위의 식에서는 하나의 독립변수에 대한 선형결합으로 표현하였지만, 독립변수가 여러개이면 이 독립변수들의 선형결합을 통해 로짓을 정의한다.

 

▷ 로지스틱 회귀의 예측은 개별 데이터의 가능도의 평균, 즉 가능도의 모수를 통해 이루어진다.

 

▶ 베이지안 로지스틱 회귀(Bayesian logistic regression)는 로짓을 나타내는 독립변수의 선형결합된 모수의 사전분포를 정의한다는 점에서 로지스틱 회귀와 차이가 있다.

 

R의 boot 패키지의 urine 데이터를 이용하여 베이지안 로지스틱 회귀 분석을 수행할 것이다. 수행 과정은 다음과 같다.

 

1. 데이터 확인

2. 모델링

3. 모델 확인

4. 모델 예측

 

1. 데이터 확인

 

In:

library(boot)

data('urine')

pairs(urine)

 

Out:

 

▷ 범주형 변수는 r, 이외의 6개의 변수는 연속형 변수로 데이터가 구성되어 있다.

 

▷ 그림에서 보다시피 몇 개의 연속형 변수는 강한 선형관계를 나타내고 있는 것을 확인할 수 있다.

 

2. 모델링

 

r을 종속변수로, 이외의 변수를 독립변수로 두어 베이지안 로지스틱 회귀 분석을 할 것이다. 위의 데이터를 이용하여 모델을 할 경우, 독립변수 간의 관계로부터 다중공선성(Multicollinearity)를 예상할 수 있다. 다중공선성은 모델링을 통한 예측만을 다루고자 한다면 상관이 없지만, 모델의 결과를 통한 해석을 하고자 한다면 문제가 발생하게 된다. 왜냐하면 다중공선성으로 인한 불안정한(Unstable) 추정치를 얻어, 종속변수와의 관계에 대한 잘못된 결론을 이끌어 낼 수 있기 때문이다.

 

다중공선성 문제를 해결하기 위해 모수의 사전분포로 라플라스 분포(Laplace distribution)를 가정할 것이다. 라플라스 분포는 0인 지점을 선호(Favor)하기 때문에 변수 선택이 가능하다는 특징이 있다.

 

In:

library(rjags)

X = scale(data[, -1], center = T, scale = T)

mod_string = " model {
  for (i in 1:length(y)) {
    y[i] ~ dbern(p[i])
    logit(p[i]) = b0 + b[1]*gravity[i] + b[2]*ph[i] + b[3]*osmo[i] + b[4]*cond[i] + b[5]*urea[i] + b[6]*calc[i]
  }

  b0 ~ dnorm(0.0, 1.0/25.0)

  for (j in 1:6) {
    b[j] ~ ddexp(0.0, sqrt(2.0)) # has variance 1.0
  }
} "

data_jags = list(y = data$r, 
                 gravity = X[,'gravity'], 
                 ph = X[,'ph'], 
                 osmo = X[,'osmo'], 
                 cond = X[,'cond'], 
                 urea = X[,'urea'], 
                 calc = X[,'calc'])

params = c('b0', 'b')

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: 77
   Unobserved stochastic nodes: 7
   Total graph size: 1085

Initializing model

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

 

▷ 개별 모수에 대하여 같은 사전분포를 사용하기 위하여 독립변수를 표준화(Standardization)하였다.

 

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

 

JAGS(Just Another Gibbs Sampler) 사용법

 

3. 모델 확인

 

마르코프 체인의 수렴 과정을 trace plot을 통해 확인하여 보자.

 

In:

plot(mod_sim)

 

Out:

 

▷ 모든 마르코프 체인이 추세(Trend)를 보이지 않으며, 수렴하고 있는 것을 확인할 수 있다.

 

▷ b[1], b[4], b[6]의 분포는 0으로부터 벗어난 형태를 보이고 있고, 이외의 모수는 0으로부터 크게 벗어나지 못한 것을 확인할 수 있다. 모수의 분포가 0으로부터 크게 벗어나지 못했다는 의미는 해당 모수의 변수가 종속변수에 크게 영향을 주지 않는 것으로 해석할 수 있다. 따라서 이러한 변수들을 제외한 새로운 모델에 대해 생각해볼 수 있다.

 

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

 

In:

gelman.diag(mod_sim)

 

Out:

Potential scale reduction factors:

     Point est. Upper C.I.
b[1]          1       1.00
b[2]          1       1.00
b[3]          1       1.01
b[4]          1       1.00
b[5]          1       1.00
b[6]          1       1.00
b0            1       1.00

Multivariate psrf

1

 

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

 

4. 모델 예측

 

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

 

In:

post_coef = colMeans(mod_comb_sim)
post_Xb = post_coef['b0'] + X %*% post_coef[-7]

phat = 1.0 / (1.0 + exp(-post_Xb))

plot(phat, jitter(data$r))

sum((data$r[phat >= 0.5] == 1), (data$r[phat < 0.5] == 0))/nrow(data)

 

Out:

[1] 0.8311688

 

▷ 실제 결과가 0인 경우의 phat은 0 근처에 모여있고, 1인 경우의 p hat은 1근처에 모여있는 것을 확인할 수 있다. 따라서 분류가 어느 정도 되고 있음을 확인할 수 있다.

 

▷ 경계값을 0.5로 두고 분류한 결과,  분류의 정확도가 83%인 것으로 나타났다.

 

▶ 경계값을 어떻게 설정하느냐에 따라 모델의 예측 결과가 달라진다. 따라서 모델의 예측 결과에 따른 성능 최적화는 따로 다루지 않겠다.

 

---

 

Reference:

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