---
lang: ja-jp
tags: TensorFlow Probability, TensorFlow, Python
---
==WIP==
# 【翻訳】 Linear Mixed Effects Models
- [GitHub: tensorflow/probability Linear_Mixed_Effects_Models.ipynb](https://github.com/tensorflow/probability/blob/master/tensorflow_probability/examples/jupyter_notebooks/Linear_Mixed_Effects_Models.ipynb)
- [TensorFlow: Linear Mixed Effects Models](https://www.tensorflow.org/probability/examples/Linear_Mixed_Effects_Models) ちょっと綺麗に表示される
## Data
R の著名なパッケージである `lme4` の `InstEval` データセットを題材にします。このデータセットは講義とその評価（レーティング）を取り扱うものです。各講義には学生や講師、専攻（部門）といったメタデータが含まれており、評価の値を応答変数として取り上げます。
```python=
def load_insteval():
"""InstEval データセットを読み込む
73421 件もの学生の大学講義に対する評価が含まれる
総学生数は 2972 講師の数は 2160 となっている
Returns:
73421 行 x 7 列の `x_train` という名称の np.ndarray と
カラムのヘッダーとして `metadata` という名称の辞書を返す
"""
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% をホールドアウト検証用のデータとして扱い、学習させたモデルに対する当てはまりの良さを評価するのに利用します。最初の数行を具体的に見ていくことにしましょう。
```python=
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 # インデックスは 0 から始める
# カテゴリ変数を 0 ~ max( n_categories ) に番号付ける
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()
```
入力特徴量を辞書型オブジェクト、対応するラベルとしてのレーティング情報を NumPy 配列として取得します。各特徴量は整数型としてエンコードされ、ラベルは浮動小数点数として表現されます。
```python=
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)
```
```python=
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
```
## Model
典型的には、線形モデルはどの点の組を選んでも線形性の保たれた関係を持つといった独立性を仮定しています。 InstEval データセットでは、観測結果としてグループごとに異なる傾きと切片を持ちうるかもしれません。線形混合モデル（階層線形モデルあるいは多変量線形モデル）はこのような事象に当てはまりの良いモデルでしょう（ Gelman & Hill, 2006 ）。
- **Students** 学生からの評価は互いに独立とは限らない。一部の学生は組織的に悪い（または良い）評価を講義に与えているかもしれない
- **Instructors** 講師からの評価もまた互いに独立とは限らない。良い講師は良い評価を得るだろうし、不得手な講師は悪い評価を得ると考えるのが自然である
- **Departments** 専攻によって評価はなんらか構造を持つはずだ。専攻によってはシビアな課題あるいは厳しい指導を課す傾向にあるだろうから、そのような専攻において評価は悪化するであろう
このことを理解するために、 $N \times D$ の特徴量 $\mathbf{X}$ と $N$ 個のラベル $\mathbf{y}$ で線形回帰モデルを考えてみましょう：
$$
\mathbf{y} = \mathbf{X} \beta + \alpha + \epsilon
$$
$\beta \in \mathbb{R}^{P}$ は傾きベクトルであり、$\alpha \in \mathbb{R}$ は切片、$\epsilon \sim \text{Normal} \left(\mathbf{0}, \mathbf{I} \right)$ は誤差項です。 $\beta$ と $\alpha$ は `fixed effects` と呼ばれています：データ点 $\left( x, y \right)$ の母集団全体に渡って固定的な影響を及ぼします。尤度関数とした場合の等価な表現は $\mathbf{y} \sim \text{Normal} \left( \mathbf{X} \beta + \alpha, \mathbf{I} \right)$ となります。この尤度関数はデータに適合する $\beta$ と $\alpha$ を推定する際に最大化される量となります。
線形回帰の拡張版としての線形混合モデルは次のように表現されます：
$$
\eta \sim \text{Normal} \left( \mathbf{0}, \sigma^{2} \mathbf{I} \right) \\
\mathbf{y} = \mathbf{X} \beta + \mathbf{Z} \eta + \alpha + \epsilon
$$
傾きベクトル $\beta \in \mathbb{R}^{P}$ と切片 $\alpha \in \mathbb{R}$ 、誤差項 $\epsilon \sim \text{Normal} \left( \mathbf{0}, \mathbf{I} \right)$ はこちらでも共通です。追加されたのは $\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$ となっています。このように分割することにより、 `fixed effect` を表現する $\beta$ と 潜在変数を表現する $\eta$ をそれぞれ別々に取り扱うことが可能になります。
:::success
【訳者注】日本語では `fixed effect` は「固定効果」、 `random effect` は「変量効果」と呼んでいるようです。
:::
ここでの閃亜愛変数のベクトル $\eta$ を `random effect` と呼んでいます：母集団全体で異なる（とはいえ部分的には共通しているかもしれない）効果です。特筆すべきは $\eta$ の平均値は $0$ であるので、データセットのラベルの平均値は `fixed effect` $\mathbf{X}\beta + \alpha$ により捉えられるということです。 `random effect` としての $\mathbf{Z}\eta$ はデータの「ばらつき」を捉える役割をになっています：『講師 #54 は平均より $1.4$ 点高く評価されている』といった事象を表現できます。
R 言語のパッケージ `lme4` の文法では、このモデルは以下のように表現されます：
```r=
ratings ~ service + (1|students) + (1|instructors) + (1|departments) + 1
```
`x` は `fixed effect` を、 `(1|x)` は `x` に対する `random effect` を、そして `1` は切片を表現しています。
:::info
【訳者注】 `service` という特徴量と、学生・講師・専攻ごとの個体差項を定式化していると思われます。
:::
これからこのモデルを JointDistribution として実装します。パラメータを効率的に記録する（例えば、 `model.trainable_variables` に含まれる全ての `tf.Variable` を記録したい場合など）ために、 `tf.Module` としてモデルのテンプレートを実装します。
```python=
class LinearMixedEffectModel(tf.Module):
def __init__(self):
# fixed effect とその他のパラメータを設定する
# これらは E ステップを通じて最適化される free parameter として扱う
self._intercept = tf.Variable(0., name='intercept') # 式の \alpha
self._effect_service = tf.Variable(0., name='effect_service') # 式の \beta
self._stddev_students = tfp.util.TransformedVariable( # 式の \sigma
1., bijector=tfb.Exp(), name='stddev_students',
)
self._stddev_instructors = tfp.util.TransformedVariabel( # 式の \sigma
1., bijector=tfb.Exp(), name='stddev_instructor',
)
self._stddev_departments = tfp.util.TransformedVariable( # 式の \sigma
1., bijector=tfb.Exp(), name='stddev_departments',
)
def __call__(self, features):
model = tfd.JointDistributionSequential([
# random effect を設定する
tfd.MultivariateNormalDiag(
loc=tf.zeros(num_students),
scale_identity_multiplier=self._stddev_students,
),
tfd.MultivariateNormalDiag(
loc=tf.zeros(num_instructors),
scale_identity_multiplier=self._stddev_instructors,
),
tfd.MultivariateNormalDiag(
loc=tf.zeros(num_departments),
scale_identity_multiplier=self._stddev_departments,
),
# 観測されたデータに対する尤度関数
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,
)
])
# 分布の生成を通じて学習可能な変数の記録を可能にするために
# `self` に対する参照を渡す。全ての TFP オブジェクトは `tf.Module` のサブクラス
# であるので、 LinearMixedEffectModel()(features_train).trainable_variables
# のようにして全ての tf.Variables を含むタプルを取得できます
model._to_track = self
return model
lmm_jointdist = LinearMixedEffectModel()
# 特徴量 / 予測変数により条件付ける
lmm_train = lmm_jointdist(features_train)
```
```python=
lmm_train.trainable_variables
# (<tf.Variable 'stddev_students:0' shape=() dtype=float32, numpy=0.0>,
# <tf.Variable 'stddev_instructors:0' shape=() dtype=float32, numpy=0.0>,
# <tf.Variable 'stddev_departments:0' shape=() dtype=float32, numpy=0.0>,
# <tf.Variable 'effect_service:0' shape=() dtype=float32, numpy=0.0>,
# <tf.Variable 'intercept:0' shape=() dtype=float32, numpy=0.0>)
```
確率的グラフプログラムでは、モデルの構造を計算グラフの観点から可視化することも可能です。プログラム中の確率変数のデータの流れを、グラフィカルモデルに基づいて（確率変数どうしの）関係を明示したものとなります。
より良い洞察を得るために、モデルのグラフによる表現を見てみましょう。例えば、 `intercept` と `effect_service` は `ratings` が与えられた上では条件付き独立です。このような関係性は、複数のクラスを利用していたり、モジュール全体に渡って定義が分散されていたり、条件分岐により呼び出されるサブモジュールが異なったりする場合には、見抜くことは困難です。また、 `tf.gather` を利用して `ratings` への潜在的な影響をモデリングしており、 `Tensor` に対するインデクシングのコストが大きい `accelerators` （ CPU, GPU, TPU etc ）ではボトルネックになり得ます。
このような統計ツール・計算ツールとしての困難を取り除く一つの方法として、グラフの可視化が取り上げられます。可視化により、様々なボトルネックの可能性を徐々に明らかにすることが可能です。
```python=
lmm_train.resolve_graph()
# (('effect_students', ()),
# ('effect_instructors', ()),
# ('effect_departments', ()),
# ('x', ('effect_departments', 'effect_instructors', 'effect_students')))
```
## Parameter Estimation
データの生成過程を記述できたので、次の目標は `fixed effect` の傾きベクトル $\beta$ と切片 $\alpha$ 、分散を表現するパラメータ $\sigma^{2}$ の推定を行うことです。これを尤度最大化問題として捉えた場合の最適化タスクは次の通りです（[セミコロン記法の説明](https://math.stackexchange.com/questions/2559085/what-does-a-semicolon-denote-in-the-context-of-probability-and-statistics/2559092)）：
$$
\max_{\beta, \alpha, \sigma} \log p \left( \mathbf{y} | \mathbf{X}, \mathbf{Z}; \beta, \alpha, \sigma \right) = \max_{\beta, \alpha, \sigma} \log \int p \left( \eta; \sigma \right) p \left( \mathbf{y} | \mathbf{X}, \mathbf{Z}, \eta; \beta, \alpha \right) \text{d} \eta
$$
このチュートリアルでは、 Monte Carlo EM アルゴリズムを使用して周辺尤度を最大化します（ Depster et al 1997, Wei and Tanner 1990 ）。 E-step として `random effects` についての条件付き尤度の期待値を MCMC により計算し、 M-step としてパラメータについての期待値の最大化を勾配降下法により実現します。
- E-step では Hamiltonian Monte Carlo （ HMC ）を設定する。学生、講師、専攻による影響（ random effect ）といった **current state** を受け取り、新たな状態を返す。 HMC によるサンプリング結果のチェーンを保存しておくための新たな TensorFlow 変数を用意することになる
## Appendix
### EM アルゴリズム
- [SlideShare: 数式を使わずイメージで理解するEMアルゴリズム](https://www.slideshare.net/yag_ays/em-algorithm-animation)
- [PDF: EM アルゴリズムの幾何学](https://staff.aist.go.jp/s.akaho/papers/josho-main.pdf)
- Expectation Step
- Maximization Step
- パラメータが決まった時の「下界」はどのように決めるのか
- 対数尤度関数の偏微分係数を何かと比較していることは分かる
Sign in with Wallet
Connect another wallet