ClasCat® AI Research

TensorFlow Probability : Tutorials : 線形混合効果モデル

Posted on by Sales Information

TensorFlow Probability : Tutorials : 線形混合効果モデル (翻訳/解説)

翻訳 : (株)クラスキャット セールスインフォメーション
作成日時 : 11/03/2017

* 本ページは、TensorFlow の本家サイトの TensorFlow Probability – Tutorials の次のページを翻訳した上で
適宜、補足説明したものです：

* サンプルコードの動作確認はしておりますが、必要な場合には適宜、追加改変しています。
* ご自由にリンクを張って頂いてかまいませんが、sales-info@classcat.com までご一報いただけると嬉しいです。

 

TensorFlow Probability : 線形混合効果モデル

線形混合効果モデルは構造化された線形関係をモデル化するための単純なアプローチです (Harville, 1997; Laird and Ware, 1982)。各データポイントは — グループにカテゴリー分けされた — 様々なタイプと実数値の出力から成ります。線形混合効果モデルは階層モデルです: それは任意の個々のデータポイントについて推論を改良するためにグループに渡り統計的強度を共有します。

このチュートリアルでは、TensorFlow Probability で現実世界の例で線形混合効果モデルを示します。Edward2 (tfp.edward2) とマルコフ連鎖モンテカルロ (tfp.mcmc) モジュールを使用します。

%matplotlib inline
from __future__ import absolute_import
from __future__ import division
from __future__ import print_function

import csv
import IPython
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import requests
import tensorflow as tf
import tensorflow_probability as tfp
import warnings

from tensorflow_probability import edward2 as ed

plt.style.use('ggplot')

 

データ

R のポピュラーな lme4 パッケージ からの InstEval データセットを使用します (Bates et al., 2015)。それは課程とそれらの評価格付けのデータセットです。各課程は生徒、講師そして学部のようなメタデータを含みそして関心のある応答変数は評価格付けです。

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.darray `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()

生徒
(students)		 講師
(instructors)		 studage lectage service 学部
(departments)		 格付け
(ratings)
40736 1651 17 4 1 1 12 5
23953 967 226 6 2 0 11 4
4252 159 897 6 5 1 1 5
64541 2616 561 8 6 1 7 1
59310 2404 838 4 2 1 1 3

入力の特徴辞書と格付けに相当するラベル出力の観点からデータセットをセットアップします。各特徴は整数としてエンコードされて各ラベル (評価格付け) は浮動小数点数としてエンコードされます。


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)。

現象の例は次を含みます :

  • 生徒 (= Students): 生徒からの観測は独立ではありません: 何人かの生徒はシステマチックに低い (or 高い) 講義格付けを与えるかもしれません。
  • 講師 (= Instructors): 講師からの観測は独立ではありません: 良い講師は一般に良い格付けを持ちそして悪い講師は一般に悪い格付けを持つことを期待します。
  • 学部 (= Departments): 学部からの観測は独立ではありません: ある学部は一般に無味乾燥なマテリアルやより厳しい等級付け (= grading) を持ち、そしてそのため他よりも低く格付けされるかもしれません。

これを捕捉するために、$N\times D$ 特徴 $\mathbf{X}$ と $N$ ラベル $\mathbf{y}$ のデータセットについて、線形回帰は次のモデルを仮定することを思い出してください :

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

ここで傾きベクトル $\beta\in\mathbb{R}^D$、切片 intercept $\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 P$ 行列 $\mathbf{X}$ と $N\times Q$ 行列 $\mathbf{Z}$ の観点で元の $N\times D$ 特徴行列を分割することにより形成されます、ここで $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 (学部)。課程評価格付けの母集団からより多くの観測が与えられたとき、新しい生徒、先生 (= teachers) または学部を見ているかもしれません。

R の lme4 パッケージ (Bates et al., 2015) のシンタックスでは、モデルは次のように要約できます :

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

ここで x は固定効果を表し、(1|x) は x のためのランダム効果を表し、そして 1 は切片項を表します。

このモデルを下で Edward プログラムとして実装します。

def linear_mixed_effects_model(features):
  # Set up fixed effects and other parameters.
  intercept = tf.get_variable("intercept", [])  # alpha in eq
  effect_service = tf.get_variable("effect_service", [])  # beta in eq
  stddev_students = tf.exp(
      tf.get_variable("stddev_unconstrained_students", []))  # sigma in eq
  stddev_instructors = tf.exp(
      tf.get_variable("stddev_unconstrained_instructors", [])) # sigma in eq
  stddev_departments = tf.exp(
      tf.get_variable("stddev_unconstrained_departments", [])) # sigma in eq

  # Set up random effects.
  effect_students = ed.MultivariateNormalDiag(
      loc=tf.zeros(num_students),
      scale_identity_multiplier=stddev_students,
      name="effect_students")
  effect_instructors = ed.MultivariateNormalDiag(
      loc=tf.zeros(num_instructors),
      scale_identity_multiplier=stddev_instructors,
      name="effect_instructors")
  effect_departments = ed.MultivariateNormalDiag(
      loc=tf.zeros(num_departments),
      scale_identity_multiplier=stddev_departments,
      name="effect_departments")

  # Set up likelihood given fixed and random effects.
  # Note we use `tf.gather` instead of matrix-multiplying a design matrix of
  # one-hot vectors. The latter is memory-intensive if there are many groups.
  ratings = ed.Normal(
      loc=(effect_service * features["service"] +
           tf.gather(effect_students, features["students"]) +
           tf.gather(effect_instructors, features["instructors"]) +
           tf.gather(effect_departments, features["departments"]) +
           intercept),
      scale=1.,
      name="ratings")
  return ratings

# Wrap model in a template. All calls to the model template will use the same
# TensorFlow variables.
model_template = tf.make_template("model", linear_mixed_effects_model)

Edward プログラムとして、その計算グラフの視点からモデルの構造も可視化できます。このグラフはプログラムの確率変数に渡るデータフローを、グラフィカルモデルの視点からそれらの関係を明確にしてエンコードします (Jordan, 2003)。

統計ツールとして、例えば、intercept と effect_service が与えられた格付けに条件付きで依存していることをより良く見るためにグラフを見るかもしれません; これはプログラムがクラス、モジュールに渡る交差参照、and/or サブルーチンで書かれている場合ソースコードから見ることはより困難かもしれません。計算ツールとして、tf.gather ops を通して ratings 変数になる潜在変数に気付くかもしれません。これは indexing Tensor がコスト高であればあるハードウェア・アクセラレータではボトルネックかもしれません; グラフの可視化はこれを容易に明白にします。

def strip_consts(graph_def, max_const_size=32):
  """Strip large constant values from graph_def."""
  strip_def = tf.GraphDef()
  for n0 in graph_def.node:
    n = strip_def.node.add()
    n.MergeFrom(n0)
    if n.op == 'Const':
      tensor = n.attr['value'].tensor
      size = len(tensor.tensor_content)
      if size > max_const_size:
        tensor.tensor_content = bytes("<stripped %d bytes>"%size, 'utf-8')
  return strip_def

def draw_graph(model, *args, **kwargs):
  """Visualize TensorFlow graph."""
  graph = tf.Graph()
  with graph.as_default():
    model(*args, **kwargs)
  graph_def = graph.as_graph_def()
  strip_def = strip_consts(graph_def, max_const_size=32)
  code = """
      <script>
        function load() {{
          document.getElementById("{id}").pbtxt = {data};
        }}
      </script>
      <link rel="import" href="https://tensorboard.appspot.com/tf-graph-basic.build.html" onload=load()>
      <div style="height:600px">
        <tf-graph-basic id="{id}"></tf-graph-basic>
      </div>
  """.format(data=repr(str(strip_def)), id='graph'+str(np.random.rand()))

  iframe = """
      <iframe seamless style="width:1200px;height:620px;border:0" srcdoc="{}"></iframe>
  """.format(code.replace('"', '"'))
  IPython.display.display(IPython.display.HTML(iframe))


draw_graph(linear_mixed_effects_model, features_train)

 

パラメータ推定

データが与えられたとき、推論のゴールはモデルの固定効果傾き $\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 et al., 1977; Wei and Tanner, 1990)。ランダム効果 (“E-step”) に関する条件付き尤度の期待値を計算するためにマルコフ連鎖モンテカルロを遂行し、そしてパラメータ (“M-step”) に関する期待値を最大化するために勾配降下を遂行します :

  • E-step については、ハミルトニアン・モンテカルロ (HMC) をセットアップします。それは現在の状態— student, instructor と department 効果 — を取りそして新しい状態を返します。その新しい状態を TensorFlow 変数に割り当て、これは HMC chain の状態を表します。
  • M-step については、定数までの周辺尤度の公平な (= unbiased) 推定を計算するために HMC からの事後分布サンプルを使用します。それから関心のあるパラメータに関する勾配を適用します。これは周辺尤度上の公平な確率的勾配ステップを生成します。それを Adam TensorFlow optimizer で実装して周辺の負を最小化します。 
log_joint = ed.make_log_joint_fn(model_template)

def target_log_prob_fn(effect_students, effect_instructors, effect_departments):
  """Unnormalized target density as a function of states."""
  return log_joint(  # fix `features` and `ratings` to the training data
    features=features_train,
    effect_students=effect_students,
    effect_instructors=effect_instructors,
    effect_departments=effect_departments,
    ratings=labels_train)

tf.reset_default_graph()

# Set up E-step (MCMC).
effect_students = tf.get_variable(  # `trainable=False` so unaffected by M-step
    "effect_students", [num_students], trainable=False)
effect_instructors = tf.get_variable(
    "effect_instructors", [num_instructors], trainable=False)
effect_departments = tf.get_variable(
    "effect_departments", [num_departments], trainable=False)

hmc = tfp.mcmc.HamiltonianMonteCarlo(
    target_log_prob_fn=target_log_prob_fn,
    step_size=0.015,
    num_leapfrog_steps=3)

current_state = [effect_students, effect_instructors, effect_departments]
next_state, kernel_results = hmc.one_step(
      current_state=current_state,
      previous_kernel_results=hmc.bootstrap_results(current_state))

expectation_update = tf.group(
    effect_students.assign(next_state[0]),
    effect_instructors.assign(next_state[1]),
    effect_departments.assign(next_state[2]))

# Set up M-step (gradient descent).
# The following should work. However, TensorFlow raises an error about taking
# gradients through IndexedSlices tensors. This may be a TF bug. For now,
# we recompute the target's log probability at the current state.
# loss = -kernel_results.accepted_results.target_log_prob
with tf.control_dependencies([expectation_update]):
  loss = -target_log_prob_fn(effect_students,
                             effect_instructors,
                             effect_departments)
  optimizer = tf.train.AdamOptimizer(learning_rate=0.01)
  minimization_update = optimizer.minimize(loss)

warm-up ステージを遂行します、これは訓練が事後分布の確率質量内で初期化されるように多くの反復のために一つの MCMC chain を実行します。それから訓練ループを実行します。それは結合的に E と M-steps を実行して訓練の間に値を記録します。

init = tf.global_variables_initializer()

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])

sess = tf.Session()
sess.run(init)

# Run warm-up stage.
for t in range(num_warmup_iters):
  _, is_accepted_val = sess.run(
      [expectation_update, kernel_results.is_accepted])
  num_accepted += is_accepted_val
  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):
  for _ in range(5):  # run 5 MCMC iterations before every joint EM update
    _ = sess.run(expectation_update)
  [
      _,
      _,
      effect_students_val,
      effect_instructors_val,
      effect_departments_val,
      is_accepted_val,
      loss_val,
  ] = sess.run([
      expectation_update,
      minimization_update,
      effect_students,
      effect_instructors,
      effect_departments,
      kernel_results.is_accepted,
      loss,
  ])
  effect_students_samples[t, :] = effect_students_val
  effect_instructors_samples[t, :] = effect_instructors_val
  effect_departments_samples[t, :] = effect_departments_val
  num_accepted += is_accepted_val
  loss_history[t] = loss_val
  if t % 500 == 0 or t == num_iters - 1:
    print("Iteration: {:>4} Acceptance Rate: {:.3f} Loss: {:.3f}".format(
        t, num_accepted / (t + 1), loss_val))

Warm-Up Iteration:   0 Acceptance Rate: 1.000
Warm-Up Iteration: 500 Acceptance Rate: 0.806
Warm-Up Iteration: 999 Acceptance Rate: 0.736
Iteration:    0 Acceptance Rate: 1.000 Loss: 103901.820
Iteration:  500 Acceptance Rate: 0.647 Loss: 96190.906
Iteration: 1000 Acceptance Rate: 0.671 Loss: 96115.094
Iteration: 1499 Acceptance Rate: 0.674 Loss: 96050.664

上で、収束閾値が検出されるまでアルゴリズムを実行しませんでした。訓練に道理があったかを確認するために、損失関数が訓練反復に渡り実際に収束する傾向があるかを検証します。

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

trace プロットもまた使用します、これは特定の潜在次元に渡るマルコフ連鎖モンテカルロ・アルゴリズムの軌道を示します。下では特定の講師効果が実際にそれらの初期状態から離れて意味深く推移して状態空間を探究することを見ます。trace プロットはまた (類似の混合した挙動を持ちながら)効果が講師に渡り異なることを示しています。

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()

 

批評

上で、モデルを適合させました。今はデータを使用してその適合を批評する調査をします、これはモデルを探究してより良く理解させてくれます。そのようなテクニックの一つは残差プロットです、これは各データポイントに対するモデルの予測と正解の差異をプロットします。モデルが正しければ、それらの差異は標準正規分布に従います; プロットのこのパターンからの任意の偏差はモデルの不適合を示します。

最初に訓練データが与えられたときの事後分布を持つランダム効果上の格付けに渡り事後予想分布を形成することにより残差プロットを構築します。特にモデルを forward に実行して推論された事後平均² を持つ事前ランダム効果上の依存性をインターセプトします。

def interceptor(rv_constructor, *rv_args, **rv_kwargs):
  """Replaces prior on effects with empirical posterior mean from MCMC."""
  name = rv_kwargs.pop("name")
  if name == "effect_students":
    rv_kwargs["value"] = np.mean(effect_students_samples, 0)
  elif name == "effect_instructors":
    rv_kwargs["value"] = np.mean(effect_instructors_samples, 0)
  elif name == "effect_departments":
    rv_kwargs["value"] = np.mean(effect_departments_samples, 0)
  return rv_constructor(*rv_args, **rv_kwargs)

with ed.interception(interceptor):
  ratings_posterior = model_template(features=features_test)

ratings_prediction = ratings_posterior.distribution.mean()

視覚的な調査において、残差はある程度標準正規分布に従うように見えます。けれども適合は完全ではありません; 終わりの方で正規分布よりもより大きな確率質量があり、これはモデルがその normality assumption (訳注: 正規分布に従うという仮定) をその適合を改善するかもしれないことを示します。

特に、InstEval データセットで格付けをモデル化するために正規分布を使用することは殆ど一般的ですが、データを良く見ると課程評価格付けは実際に 1 から 5 の元の値であることを明らかにします。これは、十分なデータがあれば元の分布か Categorical でさえも使用するべきであることを提案します。これが上の Edward プログラムへの 1-行の変更です; 同じ推論コードが適用可能です。

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

モデルが個々の予測をどのように行なうかを調査するために、students, instructors と departments のための効果のヒストグラムを見ます。これはデータポイントの特徴ベクトルの個々の成分がどのように出力結果に影響する傾向があるかを理解させてくれます。

驚くほどのことではありませんが、各生徒は講師の評価格付け上に典型的には殆ど効果を持たないことを下で見ます。興味深いことに、講師が所属する学部は大きな効果を持つことを見ます。

[
    effect_students_mean,
    effect_instructors_mean,
    effect_departments_mean,
] = sess.run([
    tf.reduce_mean(effect_students_samples, 0),
    tf.reduce_mean(effect_instructors_samples, 0),
    tf.reduce_mean(effect_departments_samples, 0),
])

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()

 

References

  • Douglas Bates and Martin Machler and Ben Bolker and Steve Walker. Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software, 67(1):1-48, 2015.
  • Arthur P. Dempster, Nan M. Laird, and Donald B. Rubin. Maximum likelihood from incomplete data via the EM algorithm. Journal of the Royal Statistical Society, Series B (Methodological), 1-38, 1977.
  • Andrew Gelman and Jennifer Hill. Data analysis using regression and multilevel/hierarchical models. Cambridge University Press, 2006.
  • 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.
  • Michael I. Jordan. An Introduction to Graphical Models. Technical Report, 2003.
  • Nan M. Laird and James Ware. Random-effects models for longitudinal data. Biometrics, 963-974, 1982.
  • Greg Wei and 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.
 

以上