임의절편 모델(Random intercept model)

R의 car 패키지의 Leinhardt 데이터를 이용하여 임의절편 모델(Random intercept model)을 모델을 만들어보자. 수행과정은 다음과 같다.

 

1. 데이터 확인

2. 모델링

3. 모델확인

 

1. 데이터 확인

 

In:

library(car)

data('Leinhardt')

pairs(Leinhardt)

head(Leinhardt)

 

Out:

          income infant   region oil
Australia   3426   26.7     Asia  no
Austria     3350   23.7   Europe  no
Belgium     3346   17.0   Europe  no
Canada      4751   16.8 Americas  no
Denmark     5029   13.5   Europe  no
Finland     3312   10.1   Europe  no

 

▷ 연속형 변수는 income과 infant, 범주형 변수로는 region과 oil로 구성되어 있다.

 

▷ 그림에서 보다시피 infant와 income은 비선형적인 음의 관계, region별로 infant가 달라지는 경향을 확인할 수 있다. 또한 oil이 yes인 경우, infant가 낮아지는 경향을 확인할 수 있다.

 

2. 모델링

 

임의절편 모델의 선택된 변수의 수준에 따라 절편이 달라지는 모형이다. 빈도주의 관점에서의 임의절편 모델은 다음과 같다.

 

 

▷ beta 0는 절편의 추정량, beta 1은 설명변수의 추정량이다. 기존의 선형회귀 분석과는 다르게 수준에 따라 달라지는 절편의 오차가 추가된다.

 

▶ 베이지안 관점에서의 임의절편 모델은 베이지안 선형회귀 분석에서 절편의 오차인 u의 두 모수에 사전분포에 대한 가정을 추가하여 모델링하게 된다.

 

위의 데이터를 이용한 모델에서는 region의 수준에 따라 절편이 달라지도록 모델링할 것이다.

 

In:

library(rjags)

data = na.omit(Leinhardt)

data$log_income = log(data$income)
data$log_infant = log(data$infant)

mod_string = " model {
  for (i in 1:length(y)) {
    y[i] ~ dnorm(mu[i], prec)
    mu[i] = a[region[i]] + b[1]*log_income[i] + b[2]*is_oil[i]
  }
  
  for (j in 1:max(region)) {
    a[j] ~ dnorm(a0, prec_a)
  }
  
  a0 ~ dnorm(0.0, 1.0/1.0e6)
  prec_a ~ dgamma(1/2.0, 1*10.0/2.0)
  tau = sqrt(1.0/prec_a)
  
  for (j in 1:2) {
    b[j] ~ dnorm(0.0, 1.0/1.0e6)
  }
  
  prec ~ dgamma(5/2.0, 5*10.0/2.0)
  sig = sqrt(1.0/prec)
} "

data_jags = list(y = data$log_infant, 
                 log_income = data$log_income,
                 is_oil = as.numeric(data$oil == 'yes'), 
                 region = as.numeric(data$region))

params = c('a0', 'a', 'b', 'sig', 'tau')

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: 101
   Unobserved stochastic nodes: 9
   Total graph size: 622

Initializing model

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

 

▷ 데이터 확인으로부터 얻은 인사이트를 바탕으로 income과 infant에 로그를 취하여 둘의 관계가 선형이 되도록 하였다.

 

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

 

JAGS(Just Another Gibbs Sampler) 사용법

 

3. 모델 확인

 

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

 

In:

gelman.diag(mod_sim)

# plot(mod_sim)

 

Out:

Potential scale reduction factors:

     Point est. Upper C.I.
a[1]       1.04       1.12
a[2]       1.05       1.13
a[3]       1.05       1.13
a[4]       1.05       1.13
a0         1.01       1.03
b[1]       1.05       1.15
b[2]       1.00       1.01
sig        1.00       1.00
tau        1.01       1.02

Multivariate psrf

1.03

 

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

 

모델의 DIC를 구하여 보자.

 

In:

dic.samples(mod, n.iter = 1e3)

 

Out:

  |**************************************************| 100%
Mean deviance:  213.7 
penalty 7.035 
Penalized deviance: 220.7 

 

▷ 모델의 DIC는 221.7로 위의 데이터를 이용하여 베이지안 선형회귀 모델의 결과(230.1)보다 더 낮은 것으로 나타났다. 따라서 DIC를 고려하였 을때, 임의절편 모델이 더 나은 것을 알 수 있다.

 

▷ 위의 모델의 실제 모수는 9개지만, Effective number of parameters는 7.452개인 것으로 나타났다. 이는 모수들간 정보를 공유하고 있기 때문에 실제 모수의 갯수보다 더 적게 나타난 것으로 볼 수 있다.

 

추정된 모수의 결과를 확인하여 보자.

 

In:

summary(mod_sim)

 

Out:

Iterations = 1001:6000
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
a[1]  6.5393 0.55904 0.0045646      0.0430400
a[2]  5.9943 0.70695 0.0057722      0.0559740
a[3]  5.8331 0.62960 0.0051407      0.0490414
a[4]  5.5144 0.86077 0.0070282      0.0682821
a0    5.9842 1.30529 0.0106576      0.0556278
b[1] -0.3383 0.10649 0.0008694      0.0086461
b[2]  0.6352 0.35149 0.0028699      0.0046075
sig   0.9197 0.06615 0.0005401      0.0006532
tau   2.0316 1.00292 0.0081888      0.0107208

2. Quantiles for each variable:

        2.5%     25%     50%     75%   97.5%
a[1]  5.4841  6.1575  6.5237  6.9064  7.6567
a[2]  4.6616  5.5152  5.9707  6.4650  7.4010
a[3]  4.6407  5.3993  5.8127  6.2531  7.1027
a[4]  3.8707  4.9235  5.4843  6.0860  7.2300
a0    3.4451  5.2072  5.9754  6.7482  8.5718
b[1] -0.5529 -0.4096 -0.3348 -0.2664 -0.1363
b[2] -0.0626  0.4034  0.6363  0.8676  1.3372
sig   0.8013  0.8733  0.9152  0.9623  1.0603
tau   0.9855  1.4085  1.7742  2.3407  4.6289

 

▷ 위의 모델의 결과로부터 수준별 절편에 따른 실제의 해석은 불가능하다. 왜냐하면 절편은 모든 설명변수가 0인 상태에서의 결과를 의미하는데, 위의 데이터의 범위에서는 그런 경우는 존재할 수 없기 때문이다.

 

▷ a0와 tau는 절편에 대한 평균과 표준편차로 해석할 수 있다.

 

▷ 이외의 결과에 대한 해석은 베이지안 선형회귀 분석과 같다.

 

---

Reference:

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