[통계] 베이지안 회귀분석 - 계층적선형모델(HLM) with PyMC3

 계층적 선형 모델(Hierarchical linear mode)은 관측한 표본 데이터에 계층 구조가 있을 경우 적용해볼 수 있는 모델이다.

 예를 들면, 전국에서 학생들의 특징(x)과 성적(y)에 관해 데이터를 모았다고 해보자. 이 때 나타나는 x와 y의 관계만으로는 올바른 설명이 이루어 질 수 없다. 전국에서 추출한 학생들의 특징은 '학생-반-학교-지역'으로 이루어 지는 계층 구조에 속해 있다. 이렇게 보면 전체 데이터에 나타난 집단 수준(반, 학교, 지역)의 차이에 의해 x와 y의 관계를 올바로 예측할 수 없다.

 그래서 등장한 것이 계층적선형모델이다. 특히 오늘 소개할 베이지안 회귀분석을 이용한 계층적선형모델에서 유명한 예는 Gelman과 Hill에 의한 집으로 침투하는 라돈(Radon) 함유량 연구다. Radon 함유량을 집의 지하실 보유 여부와 지역의 우라늄 레벨이라는 두 가지 요인을 계층화하여 분석한 예로, 여기에 아주 자세한 예제가 있다. 

 이 글에서는 계층 구조로 생성한 데이터를 통해 올바른 회귀분석이 이루어지는지 확인하는 예를 정리해 보려 한다.

---

모델

 이전 글에서 예로 든 선형 모델은 다음과 같다.

$y \sim \theta_{1}x + \theta_{2} + \epsilon $

 이번에 제시하는 계층 모델에서는 $\theta_{1}$, $\theta_{2}$가 각 관측 값이 속한 group에 따라 다른 경우를 상정한다.

$y \sim \theta_{1;g}x + \theta_{2;g} + \epsilon $

 이번에 생성할 예는 group별 $\theta_{1;g}$와 $\theta_{2;g}$ 또한 정규분포를 따르게 하고, 이 분포를 계층적선형모델을 통해 추측해 보려 한다.

$\theta_{1;g} \sim N(\mu_{1}, \sigma_{1})$

$\theta_{2;g} \sim N(\mu_{2}, \sigma_{2})$

 만약 윗 식의 확률분포에 대한 정보를 알 수 있다면, 우리는 x와 y의 관계에 group간 차이는 존재하지만, 이를 고려해도 전반적으로 x와 y는 어떤 관계에 있는지 확인할 수 있게 된다.

---

예제 데이터 생성

$y \sim \theta_{1;g}x + \theta_{2;g} + \epsilon $

$\theta_{1;g} \sim N(\mu_{1}, \sigma_{1})$

$\theta_{2;g} \sim N(\mu_{2}, \sigma_{2})$

$\mu_{1} = 2$

$\sigma_{1} = 1$

$\mu_{2} = 5$
$\sigma_{2} = 2$

$\epsilon \sim N(0, \sigma)$

$\sigma = 0.5$

$\theta_{1;g}$, $\theta_{2;g}$는 정규분포에서 추출된 값이며, 이 값이 그룹 별로 model을 이루는 경우에 해당하는 데이터를 생성하는 것이 목표다. 최종적으로는 베이지안 회귀를 통해 $\mu_{1}$, $\sigma_{1}$, $\mu_{2}$, $\sigma_{2}$의 확률분포를 추정하는 것이 목표다.

코드)

import matplotlib.pyplot as plt
import numpy as np


""" 데이터 사이즈 """
g_n = 20        # 그룹 개수
g_size = 30     # 한 그룹에 속한 데이터 개수

""" 파라미터 """
mu_1 = 2
sigma_1 = 1
theta1_g = mu_1 + np.random.normal(loc=0, scale=sigma_1, size=g_n)
mu_2 = 5
sigma_2 = 2
theta2_g = mu_2 + np.random.normal(loc=0, scale=sigma_2, size=g_n)
sigma = 0.5

""" 데이터 생성 """
for i in range(g_n):
    x_g = np.linspace(start=0, stop=10, num=g_size)
    y_g = theta1_g[i] * x_g\
          + theta2_g[i]\
          + np.random.normal(loc=0, scale=sigma, size=g_size)
    g_idx = np.array([i] * g_size)
    if i == 0:
        x = x_g.copy()
        y = y_g.copy()
        g = g_idx.copy()
    else:
        x = np.append(x, x_g)
        y = np.append(y, y_g)
        g = np.append(g, g_idx)

""" 플롯 """
fig = plt.figure()
plt.scatter(x, y, c=g, cmap='jet')
plt.xlabel('x')
plt.ylabel('x')
plt.show()

 그룹의 개수는 20, 각 그룹 당 데이터 개수는 30, 총 600개의 관 값을 생성한다. 그룹 별로 다른 $\theta_{1;g}$, $\theta_{2;g}$를 적용하기 위해 for문으로 데이터를 생성하였다.

 생성한 데이터의 모습은 다음과 같다.

 각 그룹 별로 색을 다르게 표시했고, 데이터들이 그룹 별로 조금씩 다른 model로 나타나는 걸 육안으로 확인할 수 있다.

---

선형 회귀 모델 (pooled)

 계층적 선형모델을 적용하기 전에 선형 회귀모델 결과를 먼저 살펴보자. 선형 회귀모델은 데이터를 한꺼번에 몰아서 때문에 pooling이라고 불린다. 선형 회귀모델은 이전 글과 같은 방법으로 분석한다.

코드)

import pymc3 as pm
import arviz as az


""" model 생성 """
with pm.Model() as pooled:

    # Prior probability (proposal distribution)
    a = pm.Normal('theta1', mu=0, sd=10)
    b = pm.Normal('theta2', mu=0, sd=10)
    eps = pm.HalfCauchy('sigma', 1)

    # Model
    y_est = a * x + b

    # Likelihood distribution
    likelihood = pm.Normal('y', mu=y_est, sd=eps, observed=y)

""" MCMC 샘플링 """
with pooled:
    pooled_trace = pm.sample(3000, progressbar=True, chains=2, cores=1)

""" 사후분포 plot """
pm.traceplot(pooled_trace)

사후분포)

 파라미터 탐색을 통해 얻은 사후분포에서 $\theta_{1}$ 은 평균 2 부근, $\theta_{2}$는 평균 5 부근으로  데이터 생성 시 설정한 값 부근으로 탐색이 이루어 진 것을 확인할 수 있다. 하지만, 오차 eps의 표준편차가 상당히 높게 측정 되므로, 현 모델의 예측력이 높다고 할 수 없다는 것을 알 수 있다.

---

선형 회귀 모델 (unpooled)

 Unpooled model은 데이터 전체를 한꺼번에 처리하지 않고, 그룹 별로 다른 $\boldsymbol{\theta}$를 갖는 모델이다. 계층적 선형모델과 다른점은 그룹별로 다른 $\boldsymbol{\theta}$의 확률분포에 대해 고려하지 않는다.

코드)

""" model 생성 """
with pm.Model() as unpooled:

    # Prior probability (proposal distribution)
    a = pm.Normal('theta1_g', mu=0, sd=100, shape=g_n)
    b = pm.Normal('theta2_g', mu=0, sd=100, shape=g_n)
    eps = pm.HalfCauchy('eps', 5)

    # Model
    y_est = a[g] * x + b[g]

    # Likelihood distribution
    likelihood = pm.Normal('y', mu=y_est, sd=eps, observed=y)

""" MCMC 샘플링 """
with unpooled:
    unpooled_trace = pm.sample(3000, progressbar=True, chains=2, cores=1)

""" 사후분포 plot """
pm.traceplot(unpooled_trace)

사후분포)

 사후분포 결과로부터 20개의 조직에 대한 파라미터들의 확률분포의 예측이 이루어진 것을 확인할 수 있다. 그리고 이 때, eps의 표준편차는 0.5로 우리가 처음 데이터 생성 시 설정한 값에 가까운 오차의 분산을 확인할 수 있다.

---

선형적 계층 모델 (Hierarchical linear model)

 선형적 계층 모델은 위의 unpooled model에서 추정한 theta1_g와 theta2_g가 하나의 확률분포에서 추출되었다는 것을 상정한다. 즉, 파라미터 $\boldsymbol{\theta}$는 그룹에 따라 오차가 존재할 수 있지만, 크게 보면 하나의 분포를 따른다는 관점의 모델이다.

코드)

""" model 생성 """
with pm.Model() as HLM:

    # Prior probability (proposal distribution) - level2
    mu_a = pm.Normal('mu_1', mu=0, sd=100)
    sigma_a = pm.HalfCauchy('sigma_1', 5)
    mu_b = pm.Normal('mu_2', mu=0, sd=100)
    sigma_b = pm.HalfCauchy('sigma_2', 5)

    # Prior probability (proposal distribution) - level1
    a = pm.Normal('theta1_g', mu=mu_a, sd=sigma_a, shape=g_n)
    b = pm.Normal('theta2_g', mu=mu_b, sd=sigma_b, shape=g_n)
    eps = pm.HalfCauchy('eps', 5)

    # Model
    y_est = a[g] * x + b[g]

    # Likelihood distribution
    likelihood = pm.Normal('y', mu=y_est, sd=eps, observed=y)

""" MCMC 샘플링 """
with HLM:
    HLM_trace = pm.sample(3000, progressbar=True, chains=2, cores=1)

""" 사후분포 plot """
pm.traceplot(HLM_trace)

사후분포)

 사후분포확률에서 mu_1, mu_2, sigma_1, sigma2라는 파라미터가 탐색 된 것을 확인할 수 있다. 예제 데이터 생성 시, theta_1은 평균 mu_1에 표준편차 sigma_1의 분포에서, theta_2는 평균 mu_2에 표준편차 sigma_2의 분포로 부터 추출했고, 각각의 값은 다음과 같았다.

$\mu_{1} = 2$

$\sigma_{1} = 1$

$\mu_{2} = 5$
$\sigma_{2} = 2$

 베이지안 모델로부터 탐색한 파라미터의 확률분포 값이 우리가 처음 설정한 파라미터 값에 가깝게 탐색 되었다는 것을 확인할 수 있다.

참고 사이트)

https://docs.pymc.io/notebooks/multilevel_modeling.html?highlight=sampler

https://docs.pymc.io/notebooks/GLM-hierarchical.html