계층적 모델(Hierarchical model)

다음의 예제를 통해 계층적 모델(Hierarchical model)의 특징을 알아보고, 모델링 결과에 대해 분석해 보자.

 

문제)

 

칙촉을 생산하는 5개의 공장이 있다. 각 공장에서 생산된 칙촉 과자 1개에 박혀 있는 초코칩 개수가 포아송 분포를 따르고, 포아송 분포의 모수는 감마분포를 따른다. cookies 데이터를 이용하여 칙촉 과자가 생산될 때, 박혀 있는 초코칩 개수에 대한 모델링을 수행한 후, 분석하시오.

 

풀이)

 

 

▷ 위의 문제에 대해 크게 3가지 모델링 방법으로 접근이 가능하다.

 

(1) Fully independent model: 모든 데이터가 독립이라 가정하고, 하나의 포아송 모델을 만드는 것이다. 이는 각 공장별 차이와 같은 공장에서 생산된 칙촉의 비슷한 특성을 고려하지 못한다는 한계가 있다.

 

(2) Seperate model: 각 공장에서 생산된 칙촉에 대한 개별 모델을 만드는 것이다. 이는 각 모델의 모수를 추정하는데, 다른 공장에서 생산된 칙촉의 데이터를 이용하지 못한다는 한계가 있다.

 

(3) Location dependent model: 각 공장별 포아송 분포의 모수를 따로 추정하고, 각 모수는 감마분포를 따른다고 하였을 때, 감마분포의 두 초모수(Hyperparameter)가 공통적으로 같은 분포를 따른다고 정하여 하나의 모델로 나타내는 것이다. 따라서 하나의 모델에서 모든 데이터를 활용하며, 공장별 생산된 칙촉의 특성을 반영할 수 있다.

 

▶ 위의 Location dependent model과 One-way ANOVA는 비슷한 형태를 가진 것을 알 수 있는데, One-way ANOVA는 모수의 분포가 고정된(Fixed) 초모수를 가지고 있다는 점에서 Location dependent model과 다르다고 할 수 있다.

 

Location dependent model의 그래프 표현(Graphical representation)은 다음과 같다.

 

 

세 모델링 방법 중 계층적 포아송-감마 모델을 통해 분석을 할 것이다. 수행 과정은 다음과 같다.

 

1. 데이터 확인

2. 모델링

3. 모델 확인

4. 모델 예측 및 분석

 

1. 데이터 확인

 

In:

data = read.table(file = '../input/cookies.txt', header = T)

boxplot(chips ~ location, data = data)

 

Out:

 

▷ 데이터는 칙촉 1개에 포함된 초코칩의 갯수인 chips와 생산공장을 의미하는 location으로 구성되어 있다.

 

▷ 공장별 생상된 칙촉의 초코칩 갯수가 서로 다른 분포를 가지고 있음을 확인할 수 있다.

 

2. 모델링

 

In:

library(rjags)

mod_string = " model {
  for (i in 1:length(chips)) {
    chips[i] ~ dpois(lam[location[i]])
  }
  
  for (j in 1:max(location)) {
    lam[j] ~ dgamma(alpha, beta)
  }
  
  mu ~ dgamma(2.0, 1.0/5.0)
  sig ~ dexp(1.0)
  
  alpha = mu^2 / sig^2
  beta = mu / sig^2
} "

data_jags = as.list(data)

params = c('lam', 'mu', 'sig')

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: 150
   Unobserved stochastic nodes: 7
   Total graph size: 315

Initializing model

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

 

▷ alpha와 beta에 대하여 Re-parameterization을 하여 분포를 가정하였다.

 

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

 

JAGS(Just Another Gibbs Sampler) 사용법

 

3. 모델 확인

 

In:

gelman.diag(mod_sim)

 

Out:

Potential scale reduction factors:

       Point est. Upper C.I.
lam[1]          1          1
lam[2]          1          1
lam[3]          1          1
lam[4]          1          1
lam[5]          1          1
mu              1          1
sig             1          1

Multivariate psrf

1

 

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

 

4. 모델 예측 및 분석

 

In:

post_params = colMeans(mod_comb_sim)

y_hat = rep(post_params[1:5], each = 30)
resid = data$chips - y_hat

plot(jitter(y_hat), resid)

var(resid[y_hat < 7])
var(resid[y_hat > 11])

 

Out:

[1] 6.447126
[1] 13.72414

 

▷ 예측값인 y 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
lam[1]  9.283 0.5359 0.004376       0.004595
lam[2]  6.227 0.4580 0.003740       0.004430
lam[3]  9.523 0.5440 0.004442       0.004543
lam[4]  8.939 0.5327 0.004349       0.004337
lam[5] 11.762 0.6196 0.005059       0.005693
mu      9.106 0.9866 0.008055       0.012361
sig     2.078 0.7131 0.005823       0.012066

2. Quantiles for each variable:

         2.5%    25%    50%    75%  97.5%
lam[1]  8.271  8.921  9.269  9.633 10.375
lam[2]  5.374  5.907  6.216  6.530  7.166
lam[3]  8.495  9.143  9.511  9.887 10.601
lam[4]  7.934  8.571  8.923  9.292 10.000
lam[5] 10.593 11.337 11.752 12.177 13.018
mu      7.239  8.492  9.066  9.694 11.165
sig     1.093  1.574  1.949  2.423  3.798

 

▷ 5번째 공장의 lambda가 가장 큰 것으로 나타났다.

 

사후예측분포를 확인하여 보자.

 

In:

n_sim = nrow(mod_comb_sim)

lam_pred = rgamma(n = n_sim, 
                  shape = mod_comb_sim[, 'mu']^2/mod_comb_sim[, 'sig']^2, 
                  rate = mod_comb_sim[, 'mu']/mod_comb_sim[, 'sig']^2)
hist(lam_pred)

 

Out:

 

각 공장별 사후예측분포를 확인하여 보자.

 

In:

y_pred_1 = rpois(n_sim, lambda = mod_comb_sim[, 'lam[1]'])

plot(density(y_pred_1, bw = 1), col = 'red', ylim = c(0, 0.15), main = NA)

y_pred_2 = rpois(n_sim, lambda = mod_comb_sim[, 'lam[2]'])

lines(density(y_pred_2, bw = 1), col = 'blue')

y_pred_3 = rpois(n_sim, lambda = mod_comb_sim[, 'lam[3]'])

lines(density(y_pred_3, bw = 1), col = 'green')

y_pred_4 = rpois(n_sim, lambda = mod_comb_sim[, 'lam[4]'])

lines(density(y_pred_4, bw = 1), col = 'yellow')

y_pred_5 = rpois(n_sim, lambda = mod_comb_sim[, 'lam[5]'])

lines(density(y_pred_5, bw = 1), col = 'black')

 

Out:

 

▷ 각 공장별 사후예측분포가 서로 다른 것을 확인할 수 있고, 5번째 공장의 평균이 가장 큰 것을 확인할 수 있다.

 

▷ 위의 그림을 그리기 위해 생성한 사후예측분포에 샘플을 이용하여 다양한 가설 검증이 가능하다. 예를 들어, 1번째 공장에서 생산한 칙촉의 초코칩의 갯수가 2번째 공장에서 생산한 초코칩보다 많을 가능성, 1번째 공장에서 생산한 칙촉의 초코칩의 갯수가 7개보다 작을 확률 등이 있다.

 

▶ 계층적 모델의 장점은 하나의 모델로부터 다양한 가설에 대한 분석이 가능하다는 것이다. 따라서 다른 소스 또는 시간에서 얻은 데이터를 이용하여 분석하는 메타 분석(Meta analysis)에 활용된다.

 

---

Reference:

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