베이지안 선형 회귀(Bayesian linear regression)

R의 car 패키지의 Leihardt 데이터를 이용하여 베이지안 선형 회귀(Bayesian linear regression) 분석을 수행할 것이다. 수행 과정은 다음과 같다.

 

1. 데이터 확인

2. 모델링

3. 모델 확인

4. 잔차 분석

 

1. 데이터 확인

 

In:

library(car)

data('Leinhardt')

pairs(Leinhardt)

 

Out:

 

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

 

▷ income과 infant가 비선형적인 관계를 나타내고 있는 것을 확인할 수 있다. 선형모델은 기본적으로 변수간 선형관계를 가정하기 때문에 두 변수에 로그를 취하여 변형할 것이다.

 

In:

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

plot(log_infant ~ log_income, data = Leinhardt)

 

Out:

 

▷ 위의 그림은 income과 infant에 로그를 취한 후, 시각화한 모습이다. 선형관계가 나타나는 것을 확인할 수 있다.

 

2. 모델링

 

infant를 종속변수로, income을 독립변수로 두어 베이지안 선형회귀 분석을 할 것이다. 모델에 대한 그래프  표현(Graphical representation)은 다음과 같다.

 

 

▷ beta0, ..., betak와 sigma^2는 특정 분포를 따르고 있고, 개별 관측치 yi는 독립이다.

 

▷ beta의 분포를 라플라스 분포(Laplace distribution)로 두면 LASSO(Least Absolute Shrinkage and Selection Operator) 회귀를 모델링할 수 있다. LASSO는 모수의 사전분포로 라플라스 분포를 사용하기 때문에 모수가 0인 지점을 선호(Favor)한다. 따라서 모델링 과정에서 변수 선택을 할 수 있다는 특징이 있다.

 

beta의 사전분포를 정규분포로, sigma^2를 역감마분포로 두고, R의 JAGS를 이용하여 모델링을 할 것이다.

 

In:

library(rjags)

data = na.omit(Leinhardt)

mod_string = " model {
  for (i in 1:n) {
    y[i] ~ dnorm(mu[i], prec)
    mu[i] = b[1] + b[2]*log_income[i]
  }
  
  for (j in 1:2) {
    b[j] ~ dnorm(0.0, 1.0/1.0e6)
  }
  
  prec ~ dgamma(5.0/2.0, 5.0*10.0/2.0)
  sig2 = 1.0/prec
  sig = sqrt(sig2)
} "

data_jags = list(y = data$log_infant, 
                 n = nrow(data), 
                 log_income = data$log_income)

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

inits = function() {
  inits = list('b' = rnorm(2, 0.0, 100.0), 
               'prec' = rgamma(1, 1.0, 1.0))
}

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

update(mod, 1e3)

mod_sim = coda.samples(model = mod, 
                       variable.names = params, 
                       n.iter = 1e4)

mod_comb_sim = do.call(rbind, mod_sim)

 

Out:

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

 

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

 

JAGS(Just Another Gibbs Sampler) 사용법

 

3. 모델 확인

 

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

 

In:

plot(mod_sim)

 

Out:

 

▷ b[1], b[2], sig의 마르코프 체인이 수렴하고 있는 것을 확인할 수 있다.

 

▶ trace plot의 색깔이 다른 선은 서로 다른 초기값을 통한 수렴 과정을 나타내기 위한 것이다.

 

▶ 오른쪽의 그래프는 각 모수의 생성된 데이터로부터 사후분포를 추정한 결과이다.

 

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

 

In:

gelman.diag(mod_sim)

 

Out:

Potential scale reduction factors:

     Point est. Upper C.I.
b[1]          1       1.01
b[2]          1       1.01
sig           1       1.00

Multivariate psrf

1

 

▷ 모든 결과가 1로 마르코프 체인이 수렴한다는 것을 알 수 있다.

 

생성된 데이터의 자기상관성을 확인하여 보자.

 

In:

autocorr.diag(mod_sim)

 

Out:

             b[1]       b[2]          sig
Lag 0  1.00000000 1.00000000  1.000000000
Lag 1  0.95097734 0.95103631  0.030863049
Lag 5  0.78000619 0.77872439 -0.003679779
Lag 10 0.61439905 0.61215120  0.007809461
Lag 50 0.07377375 0.07599907  0.003517834

 

▷ 세 모수 모두 시차가 50인 지점에서는 자기상관성이 거의 나타나지 않는 것을 확인할 수 있다.

 

생성된 데이터의 ESS(Effective Sample Size)를 알아보자.

 

In:

effectiveSize(mod_sim)

 

Out:

      b[1]       b[2]        sig 
  723.2895   722.9372 27469.2788 

 

▷ b[1]과 b[2]는 ESS가 약 700여개, sig는 27,000여개인 것으로 나타났다. 총 생성된 데이터의 수가 각 모수별로 30,000개인 것을 감안하면 b[1]과 b[2]의 ESS는 상대적으로 sig에 비해 작다는 것을 알 수 있다.

 

▷ b[1]과 b[2]의 생성된 데이터를 이용하여 각 모수의 평균을 구한다면 괜찮지만, 95% 신용구간을 구하기에는 ESS가 부족한 것을 알 수 있다. 이에 대한 기준은 raftery.diag()를 이용하여 확인할 수 있다.

 

생성된 데이터의 요약된 결과를 확인하여 보자.

 

In:

summary(mod_sim)

 

Out:

Iterations = 1001:11000
Thinning interval = 1 
Number of chains = 3 
Sample size per chain = 10000 

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

        Mean      SD  Naive SE Time-series SE
b[1]  7.1775 0.44692 0.0025803      0.0166242
b[2] -0.5170 0.07236 0.0004178      0.0026928
sig   0.9709 0.06848 0.0003953      0.0004135

2. Quantiles for each variable:

        2.5%     25%     50%     75%   97.5%
b[1]  6.3027  6.8773  7.1805  7.4734  8.0663
b[2] -0.6600 -0.5649 -0.5177 -0.4683 -0.3749
sig   0.8485  0.9227  0.9669  1.0142  1.1145

 

4. 잔차(Residual) 분석

 

잔차 분석을 통해 모델의 가정을 만족하는지 확인하여 보자.

 

In:

X = cbind(rep(1.0, data_jags$n), data_jags$log_income)
post_params = colMeans(mod_comb_sim)

yhat = drop(X %*% post_params[1:2])
resid = data_jags$y - yhat

plot(resid)
plot(yhat, resid)
qqnorm(resid)

 

Out:

 

▷ 첫 번째 그림은 잔차의 분포를 나타낸 것이다. 추세(Trend)나 패턴을 확인할 수 없다. 다만, 잔차의 표준편차의 2배 이상 큰 두 개의 이상치(Outlier)가 존재하는 것을 확인할 수 있다.

 

▷ 두 번째 그림은 예측한 결과와 잔차의 관계를 나타낸 것이다. 추세나 패턴을 확인할 수 없고, 잔차의 대략적인 평균이 0임을 확인할 수 있다. 하지만 예측 결과가 증가함에 따라 분산이 증가하는 것과 앞의 결과와 마찬가지로 두 개의 이상치가 존재하는 것을 확인할 수 있으므로 완벽히 선형회귀의 가정을 만족한다 보기 어렵다.

 

▷ 세 번째 그림은 잔차가 정규분포를 따르는지 확인하기 위한 정규확률 그림(Q-Q plot)이다. 대체적으로 직선의 형태를 띄지만, 마지막 두 점은 직선에서 벗어나는 것을 확인할 수 있다. 이는 앞의 그림에서 보았다시피 두 개의 이상점에 따른 결과라는 것을 추론할 수 있다.

 

---

Reference:

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