Linear_Mixed_Effects_Models.ipynb

Kernel: Python 3

Copyright 2018 The TensorFlow Probability Authors.

Licensed under the Apache License, Version 2.0 (the "License");

#@title Licensed under the Apache License, Version 2.0 (the "License"); { display-mode: "form" }
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
# https://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

선형 혼합 효과 모델

TensorFlow.org에서 보기

Google Colab에서 실행

GitHub에서 소그 보기

노트북 다운로드

선형 혼합 효과 모델은 구조적인 선형 관계를 모델링하기 위한 간단한 접근 방식입니다(Harville, 1997, Laird 및 Ware, 1982). 각 데이터 포인트는 그룹으로 분류된 다양한 유형의 입력과 실제 값 출력으로 구성됩니다. 선형 혼합 효과 모델은 계층적 모델입니다. 개별 데이터 포인트에 대한 추론을 개선하기 위해 그룹 간에 통계적 강도를 공유합니다.

이 튜토리얼에서는 TensorFlow Probability의 실제 예제를 사용하여 선형 혼합 효과 모델을 보여줍니다. JointDistributionCoroutine 및 마르코프 사슬 몬테카를로(tfp.mcmc) 모듈을 사용합니다.

종속성과 전제 조건

#@title Import and set ups{ display-mode: "form" }


import csv
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import requests

import tensorflow as tf
import tensorflow_probability as tfp

tfd = tfp.distributions
tfb = tfp.bijectors

dtype = tf.float64

%config InlineBackend.figure_format = 'retina'
%matplotlib inline
plt.style.use('ggplot')

빠르게 처리하세요!

시작하기 전에 이 데모에 GPU를 사용하고 있는지 확인하겠습니다.

확인하려면 'Runtime'-> 'Change runtime type'-> 'Hardware accelerator'-> 'GPU'를 선택합니다.

다음 코드 조각은 GPU에 대한 액세스 권한이 있는지 확인합니다.

if tf.test.gpu_device_name() != '/device:GPU:0':
  print('WARNING: GPU device not found.')
else:
  print('SUCCESS: Found GPU: {}'.format(tf.test.gpu_device_name()))

WARNING: GPU device not found.

참고: 어떤 이유로 GPU에 액세스할 수 없는 경우에도 이 colab은 계속 동작합니다(훈련은 더 오래 걸립니다).

데이터

R에서 인기 있는 lme4 패키지의 InstEval 데이터세트를 사용합니다(Bates 등, 2015). 과정 및 평가 등급의 데이터세트입니다. 각 과정에는 students, instructors 및 departments와 같은 메타데이터가 포함되며 관심 있는 응답 변수는 평가 등급입니다.

def load_insteval():
  """Loads the InstEval data set.

  It contains 73,421 university lecture evaluations by students at ETH
  Zurich with a total of 2,972 students, 2,160 professors and
  lecturers, and several student, lecture, and lecturer attributes.
  Implementation is built from the `observations` Python package.

  Returns:
    Tuple of np.ndarray `x_train` with 73,421 rows and 7 columns and
    dictionary `metadata` of column headers (feature names).
  """
  url = ('https://raw.github.com/vincentarelbundock/Rdatasets/master/csv/'
         'lme4/InstEval.csv')
  with requests.Session() as s:
    download = s.get(url)
    f = download.content.decode().splitlines()

  iterator = csv.reader(f)
  columns = next(iterator)[1:]
  x_train = np.array([row[1:] for row in iterator], dtype=np.int)
  metadata = {'columns': columns}
  return x_train, metadata

데이터세트를 로드하고 전처리합니다. 보이지 않는 데이터 포인트에 대해 적합 모델을 평가할 수 있도록 데이터의 20%를 보관합니다. 아래에서는 처음 몇 행을 시각화합니다.

data, metadata = load_insteval()
data = pd.DataFrame(data, columns=metadata['columns'])
data = data.rename(columns={'s': 'students',
                            'd': 'instructors',
                            'dept': 'departments',
                            'y': 'ratings'})
data['students'] -= 1  # start index by 0
# Remap categories to start from 0 and end at max(category).
data['instructors'] = data['instructors'].astype('category').cat.codes
data['departments'] = data['departments'].astype('category').cat.codes

train = data.sample(frac=0.8)
test = data.drop(train.index)

train.head()

<ipython-input-4-5d7a9eabeea1>:21: DeprecationWarning: `np.int` is a deprecated alias for the builtin `int`. To silence this warning, use `int` by itself. Doing this will not modify any behavior and is safe. When replacing `np.int`, you may wish to use e.g. `np.int64` or `np.int32` to specify the precision. If you wish to review your current use, check the release note link for additional information.
Deprecated in NumPy 1.20; for more details and guidance: https://numpy.org/devdocs/release/1.20.0-notes.html#deprecations
  x_train = np.array([row[1:] for row in iterator], dtype=np.int)

입력의 features 사전과 등급에 해당하는 labels 출력의 관점에서 데이터세트를 설정했습니다. 각 특성은 정수로 인코딩되고 각 레이블(평가 등급)은 부동 소수점 숫자로 인코딩됩니다.

get_value = lambda dataframe, key, dtype: dataframe[key].values.astype(dtype)
features_train = {
    k: get_value(train, key=k, dtype=np.int32)
    for k in ['students', 'instructors', 'departments', 'service']}
labels_train = get_value(train, key='ratings', dtype=np.float32)

features_test = {k: get_value(test, key=k, dtype=np.int32)
                 for k in ['students', 'instructors', 'departments', 'service']}
labels_test = get_value(test, key='ratings', dtype=np.float32)

num_students = max(features_train['students']) + 1
num_instructors = max(features_train['instructors']) + 1
num_departments = max(features_train['departments']) + 1
num_observations = train.shape[0]

print("Number of students:", num_students)
print("Number of instructors:", num_instructors)
print("Number of departments:", num_departments)
print("Number of observations:", num_observations)

Number of students: 2972
Number of instructors: 1128
Number of departments: 14
Number of observations: 58737

모델

일반적인 선형 모델은 모든 데이터 포인트 쌍이 일정한 선형 관계를 갖는 독립성을 가정합니다. InstEval 데이터세트에서 관측치는 그룹으로 표시되면, 각 그룹은 다양한 기울기와 절편을 가질 수 있습니다. 계층적 선형 모델 또는 다층 선형 모델이라고도 하는 선형 혼합 효과 모델은 이 현상을 포착합니다(Gelman 및 Hill, 2006).

이 현상의 예는 다음과 같습니다.

• 학생: 학생의 관찰은 독립적이지 않습니다. 일부 학생은 체계적으로 낮은(또는 높은) 강의 평가 등급을 줄 수 있습니다.

• 교사: 교사의 관찰은 독립적이지 않습니다. 좋은 교사는 일반적으로 좋은 평가를 받고 나쁜 교사는 일반적으로 나쁜 평가를 받을 것으로 기대합니다.

• 부서: 부서의 관찰은 독립적이지 않습니다. 특정 부서는 일반적으로 자료가 빈약하거나 등급 매기기가 더 엄격하게 적용될 수 있으므로, 다른 부서보다 낮은 평가를 받을 수 있습니다.

이를 캡처하기 위해 $N\times D$ 특성 $\mathbf{X}$ 및 $N$ 레이블 $\mathbf{y}$의 데이터세트에 대해 선형 회귀가 모델을 가정한다는 점을 기억하세요.

$$\begin{equation*} \mathbf{y} = \mathbf{X}\beta + \alpha + \epsilon, \end{equation*}$$

여기에 기울기 벡터 $\beta\in\mathbb{R}^D$, 절편 $\alpha\in\mathbb{R}$와 무작위 노이즈 $\epsilon\sim\text{Normal}(\mathbf{0}, \mathbf{I})$가 있습니다. $\beta$와 $\alpha$는 '고정 효과'라고 합니다. 이들은 데이터 포인트 $(x, y)$ 전체에 걸쳐 일정하게 유지되는 효과입니다. 등가식의 가능성 공식은 $\mathbf{y} \sim \text{Normal}(\mathbf{X}\beta + \alpha, \mathbf{I})$입니다. 이 가능성은 데이터에 맞춤 조정된 $\beta$와 $\alpha$의 포인트 추정치를 찾기 위해 추론하는 동안 최대화됩니다.

선형 혼합 효과 모델은 선형 회귀를 다음과 같이 확장합니다.

$$\begin{align*} \eta &\sim \text{Normal}(\mathbf{0}, \sigma^2 \mathbf{I}), \\ \mathbf{y} &= \mathbf{X}\beta + \mathbf{Z}\eta + \alpha + \epsilon. \end{align*}$$

여기에 여전히 기울기 벡터 $\beta\in\mathbb{R}^P$, 절편 $\alpha\in\mathbb{R}$ 및 무작위 노이즈 $\epsilon\sim\text{Normal}(\mathbf{0}, \mathbf{I})$가 있습니다. 또한, $\mathbf{Z}\eta$라는 항이 있습니다. 여기서 $\mathbf{Z}$는 특성 행렬이고, $\eta\in\mathbb{R}^Q$는 무작위 기울기의 벡터입니다. $\eta$는 분산 구성 요소 매개변수 $\sigma^2$를 사용하여 정규 분포를 따릅니다. $\mathbf{Z}$는 원래 $N\times D$ 특성 행렬을 새로운 $N\times P$ 행렬 $\mathbf{X}$ 및 $N\times Q$ 행렬 $\mathbf{Z}$로 분할하여 형성됩니다. 여기서 $P + Q=D$, 이 파티션을 사용하면 각각 고정 효과 $\beta$ 및 잠재 변수 $\eta$를 사용하여 특성을 개별적으로 모델링할 수 있습니다.

잠재 변수 $\eta$는 '무선 효과'라고 합니다. 이는 모집단에 따라 달라지는 효과입니다(하위 모집단 간에 일정할 수도 있음). 특히 무선 효과 $\eta$의 평균이 0이므로 데이터 레이블의 평균은 $\mathbf{X}\beta + \alpha$로 캡처됩니다. 무선 효과 구성 요소 $\mathbf{Z}\eta$는 데이터의 변형을 캡처합니다. 예를 들어 '교사 54번은 평균보다 1.4점 더 높습니다.'

이 튜토리얼에서는 다음과 같은 효과를 가정합니다.

• 고정 효과: service. service는 과정이 교사의 주요 부서에 속하는지 여부에 해당하는 바이너리 공변량입니다. 수집하는 추가 데이터의 양과 관계없이 $0$ 및 $1$의 값만 가질 수 있습니다.

• 무선 효과: students , instructors 및 departments. 과정 평가 등급의 모집단으로부터 더 많은 관측치가 주어지면 새로운 학생, 교사 또는 부서에 집중할 수 있습니다.

R의 lme4 패키지(Bates 등, 2015)의 구문에서 모델은 다음과 같이 요약할 수 있습니다.

ratings ~ service + (1|students) + (1|instructors) + (1|departments) + 1

여기서 x는 고정 효과를 나타내고, (1|x)는 x에 대한 무선 효과를, 1은 절편 항을 나타냅니다.

이 모델을 JointDistribution으로 구현합니다. 매개변수 추적을 더 잘 지원하기 위해(예: model.trainable_variables의 모든 tf.Variable을 추적하려는 경우) 모델 템플릿을 tf.Module로 구현합니다.

class LinearMixedEffectModel(tf.Module):
  def __init__(self):
    # Set up fixed effects and other parameters.
    # These are free parameters to be optimized in E-steps
    self._intercept = tf.Variable(0., name="intercept")            # alpha in eq
    self._effect_service = tf.Variable(0., name="effect_service")  #  beta in eq
    self._stddev_students = tfp.util.TransformedVariable(
        1., bijector=tfb.Exp(), name="stddev_students")            # sigma in eq
    self._stddev_instructors = tfp.util.TransformedVariable(
        1., bijector=tfb.Exp(), name="stddev_instructors")         # sigma in eq
    self._stddev_departments = tfp.util.TransformedVariable(
        1., bijector=tfb.Exp(), name="stddev_departments")         # sigma in eq

  def __call__(self, features):
    model = tfd.JointDistributionSequential([
      # Set up random effects.
      tfd.MultivariateNormalDiag(
          loc=tf.zeros(num_students),
          scale_diag=self._stddev_students * tf.ones(num_students)),
      tfd.MultivariateNormalDiag(
          loc=tf.zeros(num_instructors),
          scale_diag=self._stddev_instructors * tf.ones(num_instructors)),
      tfd.MultivariateNormalDiag(
          loc=tf.zeros(num_departments),
          scale_diag=self._stddev_departments * tf.ones(num_departments)),
      # This is the likelihood for the observed.
      lambda effect_departments, effect_instructors, effect_students: tfd.Independent(
          tfd.Normal(
              loc=(self._effect_service * features["service"] +
                  tf.gather(effect_students, features["students"], axis=-1) +
                  tf.gather(effect_instructors, features["instructors"], axis=-1) +
                  tf.gather(effect_departments, features["departments"], axis=-1) +
                  self._intercept),
              scale=1.),
              reinterpreted_batch_ndims=1)
    ])

    # To enable tracking of the trainable variables via the created distribution,
    # we attach a reference to `self`. Since all TFP objects sub-class
    # `tf.Module`, this means that the following is possible:
    # LinearMixedEffectModel()(features_train).trainable_variables
    # ==> tuple of all tf.Variables created by LinearMixedEffectModel.
    model._to_track = self
    return model

lmm_jointdist = LinearMixedEffectModel()
# Conditioned on feature/predictors from the training data
lmm_train = lmm_jointdist(features_train)

lmm_train.trainable_variables

(<tf.Variable 'effect_service:0' shape=() dtype=float32, numpy=0.0>,
 <tf.Variable 'intercept:0' shape=() dtype=float32, numpy=0.0>,
 <tf.Variable 'stddev_departments:0' shape=() dtype=float32, numpy=0.0>,
 <tf.Variable 'stddev_instructors:0' shape=() dtype=float32, numpy=0.0>,
 <tf.Variable 'stddev_students:0' shape=() dtype=float32, numpy=0.0>)

확률적 그래픽 프로그램으로서, 계산 그래프 측면에서 모델의 구조를 시각화할 수도 있습니다. 이 그래프는 프로그램에서 확률 변수 간의 데이터 흐름을 인코딩하여 이들 간의 관계를 명확히 합니다(Jordan, 2003).

예를 들어, intercept와 effect_service가 주어진 ratings에 따라 조건부로 종속된다는 것을 더 잘 확인하기 위해 그래프를 통계 도구로서 볼 수 있습니다. 프로그램이 클래스, 모듈 간 상호 참조 및/또는 서브 루틴으로 작성되었으면, 소스 코드에서 보기가 더 어려울 수 있습니다. 계산 도구로서, 잠재 변수가 tf.gather 연산을 통해 ratings 변수로 유입되는 것을 확인할 수도 있습니다. Tensor를 인덱싱하는 데 비용이 많이 들면 특정 하드웨어 가속기에서 병목 현상이 발생할 수 있습니다. 그래프를 시각화하면 명확하게 확인할 수 있습니다.

lmm_train.resolve_graph()

(('effect_students', ()),
 ('effect_instructors', ()),
 ('effect_departments', ()),
 ('x', ('effect_departments', 'effect_instructors', 'effect_students')))

매개변수 예측

주어진 데이터에서 추론의 목표는 모델의 고정 효과 기울기 $\beta$, 절편 $\alpha$ 및 분산 구성 요소 매개변수 $\sigma^2$를 맞춤 조정하는 것입니다. 최대 가능성 원리는 이 작업을 다음과 같이 공식화합니다.

$$ \max_{\beta, \alpha, \sigma}~\log p(\mathbf{y}\mid \mathbf{X}, \mathbf{Z}; \beta, \alpha, \sigma) = \max_{\beta, \alpha, \sigma}~\log \int p(\eta; \sigma) ~p(\mathbf{y}\mid \mathbf{X}, \mathbf{Z}, \eta; \beta, \alpha)~d\eta. $$

이 튜토리얼에서는 이 한계 밀도를 최대화하기 위해 몬테카를로 EM 알고리즘을 사용합니다(Dempster 등, 1977, Wei 및 Tanner, 1990).¹ 무선 효과('E 스텝')에 관한 조건부 가능성의 기대치를 계산하기 위해 마르코프 사슬 몬테카를로를 수행하고, 매개변수('M 스텝')에 대한 기대치를 최대화하기 위해 경사 하강법을 수행합니다.

• E 스텝의 경우 해밀턴 몬테카를로(HMC)를 설정했습니다. 현재 상태(학생, 교사 및 부서 효과)를 사용하고 새 상태를 반환합니다. HMC 사슬의 상태를 나타내는 TensorFlow 변수에 새 상태를 할당합니다.

• M 스텝의 경우 HMC의 사후 확률 샘플을 사용하여 한계 가능성의 편향되지 않은 추정치를 상수까지 계산합니다. 그런 다음 관심 있는 매개변수와 관련하여 그래디언트를 적용합니다. 이것은 한계 가능성에 대해 편향되지 않은 확률적 하강 스텝을 생성합니다. Adam TensorFlow 옵티마이저로 이를 구현하고 한계의 음수를 최소화합니다.

target_log_prob_fn = lambda *x: lmm_train.log_prob(x + (labels_train,))
trainable_variables = lmm_train.trainable_variables
current_state = lmm_train.sample()[:-1]

# For debugging
target_log_prob_fn(*current_state)

<tf.Tensor: shape=(), dtype=float32, numpy=-485996.53>

# Set up E-step (MCMC).
hmc = tfp.mcmc.HamiltonianMonteCarlo(
    target_log_prob_fn=target_log_prob_fn,
    step_size=0.015,
    num_leapfrog_steps=3)
kernel_results = hmc.bootstrap_results(current_state)

@tf.function(autograph=False, jit_compile=True)
def one_e_step(current_state, kernel_results):
  next_state, next_kernel_results = hmc.one_step(
      current_state=current_state,
      previous_kernel_results=kernel_results)
  return next_state, next_kernel_results

optimizer = tf.optimizers.Adam(learning_rate=.01)

# Set up M-step (gradient descent).
@tf.function(autograph=False, jit_compile=True)
def one_m_step(current_state):
  with tf.GradientTape() as tape:
    loss = -target_log_prob_fn(*current_state)
  grads = tape.gradient(loss, trainable_variables)
  optimizer.apply_gradients(zip(grads, trainable_variables))
  return loss

훈련이 사후 확률 질량 내에서 초기화될 수 있도록 여러 번의 반복에서 하나의 MCMC chain을 실행하는 워밍업 단계를 수행합니다. 그런 다음 훈련 루프를 실행합니다. E 및 M 스텝을 공동으로 실행하고 훈련 중에 값을 기록합니다.

num_warmup_iters = 1000
num_iters = 1500
num_accepted = 0
effect_students_samples = np.zeros([num_iters, num_students])
effect_instructors_samples = np.zeros([num_iters, num_instructors])
effect_departments_samples = np.zeros([num_iters, num_departments])
loss_history = np.zeros([num_iters])

# Run warm-up stage.
for t in range(num_warmup_iters):
  current_state, kernel_results = one_e_step(current_state, kernel_results)
  num_accepted += kernel_results.is_accepted.numpy()
  if t % 500 == 0 or t == num_warmup_iters - 1:
    print("Warm-Up Iteration: {:>3} Acceptance Rate: {:.3f}".format(
        t, num_accepted / (t + 1)))

num_accepted = 0  # reset acceptance rate counter

# Run training.
for t in range(num_iters):
  # run 5 MCMC iterations before every joint EM update
  for _ in range(5):
    current_state, kernel_results = one_e_step(current_state, kernel_results)
  loss = one_m_step(current_state)
  effect_students_samples[t, :] = current_state[0].numpy()
  effect_instructors_samples[t, :] = current_state[1].numpy()
  effect_departments_samples[t, :] = current_state[2].numpy()
  num_accepted += kernel_results.is_accepted.numpy()
  loss_history[t] = loss.numpy()
  if t % 500 == 0 or t == num_iters - 1:
    print("Iteration: {:>4} Acceptance Rate: {:.3f} Loss: {:.3f}".format(
        t, num_accepted / (t + 1), loss_history[t]))

Warm-Up Iteration:   0 Acceptance Rate: 1.000
Warm-Up Iteration: 500 Acceptance Rate: 0.758
Warm-Up Iteration: 999 Acceptance Rate: 0.729
Iteration:    0 Acceptance Rate: 1.000 Loss: 98200.422
Iteration:  500 Acceptance Rate: 0.649 Loss: 98190.469
Iteration: 1000 Acceptance Rate: 0.656 Loss: 98068.664
Iteration: 1499 Acceptance Rate: 0.663 Loss: 98155.070

더 빠른 추론을 위해 tf.while_loop에 루프를 위한 워밍업을 작성하고 tf.scan 또는 tf.while_loop에 훈련 스텝을 작성할 수도 있습니다. 예제는 다음과 같습니다.

@tf.function(autograph=False, jit_compile=True)
def run_k_e_steps(k, current_state, kernel_results):
  _, next_state, next_kernel_results = tf.while_loop(
      cond=lambda i, state, pkr: i < k,
      body=lambda i, state, pkr: (i+1, *one_e_step(state, pkr)),
      loop_vars=(tf.constant(0), current_state, kernel_results)
  )
  return next_state, next_kernel_results

위에서 수렴 임계값이 감지될 때까지 알고리즘을 실행하지 않았습니다. 훈련이 합리적인지 확인하기 위해 손실 함수가 실제로 훈련 반복 동안 수렴하는 경향이 있는지 확인합니다.

plt.plot(loss_history)
plt.ylabel(r'Loss $-\log$ $p(y\mid\mathbf{x})$')
plt.xlabel('Iteration')
plt.show()

또한 특정 잠재 차원에서 마르코프 사슬 몬테카를로 알고리즘의 궤적을 보여주는 추적 플롯을 사용합니다. 아래에서 특정 교사 효과가 실제로 초기 상태에서 의미 있게 전환되고 상태 공간을 탐색하는 것을 볼 수 있습니다. 추적 플롯은 또한 효과가 교사마다 다르지만, 혼합 동작이 유사함을 나타냅니다.

for i in range(7):
  plt.plot(effect_instructors_samples[:, i])

plt.legend([i for i in range(7)], loc='lower right')
plt.ylabel('Instructor Effects')
plt.xlabel('Iteration')
plt.show()

비평

위에서 모델을 맞춤 조정했습니다. 이제 모델을 탐색하고 더 잘 이해할 수 있도록 데이터를 사용하여 적합성을 비평하는 방법을 살펴봅니다. 이러한 기술 중 하나는 각 데이터 포인트에 대한 모델의 예측과 실측 자료의 차이를 나타내는 잔차 플롯입니다. 모델이 정확하다면 그 차이는 표준 정규 분포여야 합니다. 플롯에서 이 패턴의 편차는 모델이 적합하지 않음을 나타냅니다.

먼저 등급에 대한 사후 확률 예측 분포를 형성하여 잔차 플롯을 작성합니다.이 분포는 무선 효과에 대한 사전 확률 분포를 주어진 훈련 데이터의 사후확률로 대체합니다. 특히, 모델을 전달 실행하고 추론된 사후 확률 평균을 사용하여 사전 확률 무선 효과에 대한 의존성을 절편으로 표현합니다.²

lmm_test = lmm_jointdist(features_test)

[
    effect_students_mean,
    effect_instructors_mean,
    effect_departments_mean,
] = [
     np.mean(x, axis=0).astype(np.float32) for x in [
       effect_students_samples,
       effect_instructors_samples,
       effect_departments_samples
       ]
]

# Get the posterior predictive distribution
(*posterior_conditionals, ratings_posterior), _ = lmm_test.sample_distributions(
    value=(
        effect_students_mean,
        effect_instructors_mean,
        effect_departments_mean,
))

ratings_prediction = ratings_posterior.mean()

육안 검사에서 잔차는 어느 정도 표준 정규 분포로 보입니다. 그러나 적합도가 완벽하지 않습니다. 정규 분포보다 꼬리에 더 큰 확률 질량이 있으며, 이는 모델이 정규성 가정을 완화하여 적합도를 향상할 수 있음을 나타냅니다.

특히 InstEval 데이터세트의 등급을 모델링하기 위해 정규 분포를 사용하는 것이 가장 일반적이지만, 데이터를 자세히 살펴보면 과정 평가 등급이 실제로 1에서 5까지의 서수 값임을 알 수 있습니다. 이는 순위 분포를 사용하거나 상대 순위를 버리기에 충분한 데이터가 있으면 Categorical을 사용해야 함을 의미합니다. 위 모델에서 한 줄만 변경하면 되며, 같은 추론 코드가 적용됩니다.

plt.title("Residuals for Predicted Ratings on Test Set")
plt.xlim(-4, 4)
plt.ylim(0, 800)
plt.hist(ratings_prediction - labels_test, 75)
plt.show()

모델이 개별 예측을 수행하는 방법을 탐색하기 위해 학생, 교사 및 부서에 대한 효과 히스토그램을 살펴봅니다. 이를 통해 데이터 포인트의 특성 벡터에 있는 개별 요소가 결과에 어떤 영향을 미치는지 이해할 수 있습니다.

당연히 아래에서 각 학생은 일반적으로 교사의 평가 등급에 거의 영향을 미치지 않습니다. 흥미롭게도 교사가 속한 부서의 효과가 크다는 것을 알 수 있습니다.

plt.title("Histogram of Student Effects")
plt.hist(effect_students_mean, 75)
plt.show()

plt.title("Histogram of Instructor Effects")
plt.hist(effect_instructors_mean, 75)
plt.show()

plt.title("Histogram of Department Effects")
plt.hist(effect_departments_mean, 75)
plt.show()

각주

¹ 선형 혼합 효과 모델은 한계 밀도를 분석적으로 계산할 수 있는 특별한 경우입니다. 이 튜토리얼의 목적을 위해, 가능성이 Normal이 아닌 Categorical로 확장된 경우와 같은 비 분석적 한계 밀도에 더 쉽게 적용되는 몬테카를로 EM을 보여줍니다.

² 단순화를 위해 모델의 순방향 전달을 하나만 사용하여 예측 분포의 평균을 형성합니다. 이는 사후 확률 평균에 대한 조건화에 의해 수행되며 선형 혼합 효과 모델에 유효합니다. 그렇지만 이는 일반적으로 유효하지 않습니다. 사후 확률 예측 분포의 평균은 일반적으로 다루기 어렵고, 사후 확률 샘플이 주어진 모델의 여러 순방향 전달 간의 경험적 평균을 사용해야 합니다.

감사의 말

이 튜토리얼은 원래 Edward 1.0(소스)로 작성되었습니다. 해당 버전을 작성하고 수정하는데 기여해 주신 모든 분께 감사드립니다.

참고 문헌

1. Douglas Bates, Martin Machler, Ben Bolker 및 Steve Walker. Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software, 67(1):1-48, 2015.

2. Arthur P. Dempster, Nan M. Laird 및 Donald B. Rubin. Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society, Series B (Methodological), 1-38, 1977.

3. Andrew Gelman 및 Jennifer Hill. Data analysis using regression and multilevel/hierarchical models. Cambridge University Press, 2006.

4. David A. Harville. Maximum likelihood approaches to variance component estimation and to related problems. Journal of the American Statistical Association, 72(358):320-338, 1977.

5. Michael I. Jordan. An Introduction to Graphical Models. Technical Report, 2003.

6. Nan M. Laird 및 James Ware. Random-effects models for longitudinal data. Biometrics, 963-974, 1982.

7. Greg Wei 및 Martin A. Tanner. A Monte Carlo implementation of the EM algorithm and the poor man's data augmentation algorithms. Journal of the American Statistical Association, 699-704, 1990.