ClasCat® AI Research

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

Posted on by Sales Information

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

翻訳 : (株)クラスキャット セールスインフォメーション
更新日時 : 10/31/2018
作成日時 : 10/20/2018

* 本ページは、Edward サイトの Tutorials : Linear Mixed Effects Models を翻訳した上で適宜、補足説明したものです：

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

 

線形混合効果モデル

線形混合効果モデルでは、様々なタイプの入力を持ち、サブグループにカテゴリー分けされ、そして実数値出力に関連付けられたデータポイントのために線形関係をモデル化することを望みます。

Edward で例を示します。Jupyter notebook による対話的バージョンは こちら で利用可能です。

 

データ

ポピュラーな lme4 R パッケージ から InstEval データセットを使用します (Bates, Mächler, Bolker, & Walker, 2015)。それは講師 (= instructor) の評価格付け (= ratings) で、そこでは students (生徒) や departments (学部) のようなカテゴリーを含み、そして興味のある応答変数は講師評価格付けです。

import pandas as pd
from observations import insteval

data, metadata = insteval("~/data")
data = pd.DataFrame(data, columns=metadata['columns'])

# s - students - 1:2972
# d - instructors - codes that need to be remapped
# dept also needs to be remapped
data['s'] = data['s'] - 1
data['dcodes'] = data['d'].astype('category').cat.codes
data['deptcodes'] = data['dept'].astype('category').cat.codes
data['y'] = data['y'].values.astype(float)

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

データセットは以下を意味します :

  • students as s
  • instructors as d
  • departments as dept
  • service as service 
s_train = train['s'].values
d_train = train['dcodes'].values
dept_train = train['deptcodes'].values
y_train = train['y'].values
service_train = train['service'].values
n_obs_train = train.shape[0]

s_test = test['s'].values
d_test = test['dcodes'].values
dept_test = test['deptcodes'].values
y_test = test['y'].values
service_test = test['service'].values
n_obs_test = test.shape[0]

n_s = max(s_train) + 1  # number of students
n_d = max(d_train) + 1  # number of instructors
n_dept = max(dept_train) + 1  # number of departments
n_obs = train.shape[0]  # number of observations

print("Number of students: {}".format(n_s))
print("Number of instructors: {}".format(n_d))
print("Number of departments: {}".format(n_dept))
print("Number of observations: {}".format(n_obs))

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

 

モデル

線形回帰では、独立性仮定をします、そこでは各データポイントはお互いに定数傾き (= constant slope) で回帰します。私達の設定では、観測は様々な傾きと切片を持つかもしれないグループに由来します。こうして私達はこの挙動を捕捉できるモデルを構築したいです (Gelman & Hill, 2006)。

この現象の例として :

  • 単一の生徒からの観測は互いに独立ではありません。むしろ、ある生徒達はシステマチックに低い (or 高い) 講義格付けを与えるかもしれません。
  • 単一の教師からの観測は互いに独立ではありません。良い教師は一般に良い格付けを得て悪い教師は一般に悪い格付けを得ることを期待します。
  • 単一の学部からの観測は互いに独立ではありません。一つの学部は一般に無味乾燥なマテリアルを持ちそのため他よりも低い格付けをされるかもしれません。

典型的な線形回帰は次の形式を取ります :

\[
\mathbf{y} = \mathbf{X}\beta + \epsilon,
\]

ここで $\mathbf{X}$ は係数 $\beta$ を持つ固定効果 (= fixed effects) に対応して $\epsilon$ はランダム・ノイズ, $\epsilon\sim\mathcal{N}(\mathbf{0}, \mathbf{I})$ に対応しています。

線形混合効果モデルでは、追加の項 $\mathbf{Z}\eta$ を追加します、ここで $\mathbf{Z}$ は係数 $\eta$ を持つ変量効果 (= random effects, ランダム効果) に対応します。モデルは次の形式を取ります :

\[
\begin{aligned}
\eta &\sim \mathcal{N}(\mathbf{0}, \sigma^2 \mathbf{I}), \\
\mathbf{y} &= \mathbf{X}\beta + \mathbf{Z}\eta + \epsilon.\end{aligned}
\]

データが与えられたとき、目的は $\beta$, $\eta$ と $\sigma^2$ を推論することです、ここで $\beta$ はモデルパラメータ (“fixed effects”) で、$\eta$ は潜在変数 (“random effects”)、そして $\sigma^2$ は分散成分パラメータです。

変量効果は平均 0 を持つので、データの平均 は $\mathbf{X}\beta$ で獲得されます。変量効果成分 $\mathbf{Z}\eta$ はデータの変位 (= variations) を獲得します (e.g. 講師 #54 は平均よりも 1.4 ポイント高く格付けされます)。

もっともな質問は固定と変量効果の違いです。固定効果は与えられた母集団のための定数 (= constant) である効果です。変量効果は与えられた母集団に対して変化する効果です (i.e. それは部分母集団内では定数かもしれませんが母集団全体の中では変化します)。私達の例で下に示します :

  • 固定効果として service を選択します。それは講義が講義の主要な学部に属するか否かに対応する二値共変量です。どれだけ追加データを集めても、それは 00 と 11 の値を取るだけです。
  • 変量効果として students, teachers と departments のカテゴリー値を選択します。講師評価格付けの母集団からより多くの観測が与えられたとき、新しい students, teachers, または departments を見ているかもしれません。

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

y ~ 1 + (1|students) + (1|instructor) + (1|dept) + service

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

# Set up placeholders for the data inputs.
s_ph = tf.placeholder(tf.int32, [None])
d_ph = tf.placeholder(tf.int32, [None])
dept_ph = tf.placeholder(tf.int32, [None])
service_ph = tf.placeholder(tf.float32, [None])

# Set up fixed effects.
mu = tf.get_variable("mu", [])
service = tf.get_variable("service", [])

sigma_s = tf.sqrt(tf.exp(tf.get_variable("sigma_s", [])))
sigma_d = tf.sqrt(tf.exp(tf.get_variable("sigma_d", [])))
sigma_dept = tf.sqrt(tf.exp(tf.get_variable("sigma_dept", [])))

# Set up random effects.
eta_s = Normal(loc=tf.zeros(n_s), scale=sigma_s * tf.ones(n_s))
eta_d = Normal(loc=tf.zeros(n_d), scale=sigma_d * tf.ones(n_d))
eta_dept = Normal(loc=tf.zeros(n_dept), scale=sigma_dept * tf.ones(n_dept))

yhat = (tf.gather(eta_s, s_ph) +
        tf.gather(eta_d, d_ph) +
        tf.gather(eta_dept, dept_ph) +
        mu + service * service_ph)
y = Normal(loc=yhat, scale=tf.ones(n_obs))

 

推論

データが与えられたとき、モデルの固定と変量効果を推論することが目的です。この解析では、$\text{KL}(q\|p)$ ダイバージェンス measure を持つ変分推論を使用します。変量効果のために完全に分解された (= fully factorized) 正規近似を指定して推論のために総ての訓練データを渡します。アルゴリズムでは、固定効果は変分 EM スキームで推定されます。

q_eta_s = Normal(
    loc=tf.get_variable("q_eta_s/loc", [n_s]),
    scale=tf.nn.softplus(tf.get_variable("q_eta_s/scale", [n_s])))
q_eta_d = Normal(
    loc=tf.get_variable("q_eta_d/loc", [n_d]),
    scale=tf.nn.softplus(tf.get_variable("q_eta_d/scale", [n_d])))
q_eta_dept = Normal(
    loc=tf.get_variable("q_eta_dept/loc", [n_dept]),
    scale=tf.nn.softplus(tf.get_variable("q_eta_dept/scale", [n_dept])))

latent_vars = {
    eta_s: q_eta_s,
    eta_d: q_eta_d,
    eta_dept: q_eta_dept}
data = {
    y: y_train,
    s_ph: s_train,
    d_ph: d_train,
    dept_ph: dept_train,
    service_ph: service_train}
inference = ed.KLqp(latent_vars, data)

フィットしたモデルを批評する (= critique) 一つの方法は残差プロットです、i.e. 各データポイントに対する予測された値と観測された値の間の差異のプロットです。下で私達はアルゴリズムを初期化してループ内で個々の更新を遂行して、手動で推論を実行します。アルゴリズム進捗として残差プロットを形成します。これはアルゴリズムがデータから変量と固定効果を推論することをどのように進めているか検査する助けとなります。

残差を形成するために、最初にテストデータ上で予測します。モデルで定義された yhat をコピーして変量効果をそれらの推論された mean で置き換えることによってこれを行ないます。アルゴリズムの間、テスト入力を供給して予測を評価します。

私達はまた lme4 で同じモデル (y service + (1|dept) + (1|s) + (1|d), 特に、InstEval データセット全体上でフィット) をフィットしました。変量効果推定をセーブしてそれらを学習したパラメータと比較します。

yhat_test = ed.copy(yhat, {
    eta_s: q_eta_s.mean(),
    eta_d: q_eta_d.mean(),
    eta_dept: q_eta_dept.mean()})

今メインループを書きます。

inference.initialize(n_print=20, n_iter=100)
tf.global_variables_initializer().run()

for _ in range(inference.n_iter):
  # Update and print progress of algorithm.
  info_dict = inference.update()
  inference.print_progress(info_dict)

  t = info_dict['t']
  if t == 1 or t % inference.n_print == 0:
    # Make predictions on test data.
    yhat_vals = yhat_test.eval(feed_dict={
        s_ph: s_test,
        d_ph: d_test,
        dept_ph: dept_test,
        service_ph: service_test})

    # Form residual plot.
    plt.title("Residuals for Predicted Ratings on Test Set")
    plt.xlim(-4, 4)
    plt.ylim(0, 800)
    plt.hist(yhat_vals - y_test, 75)
    plt.show()

 

批評

上で、残差プロットを通してモデルの fit を診断するための方法を記述しました。モデルの post-convergence (収束後) のための残差を示します (中間的な残差プロットのためには Jupyter notebook をチェックしてください)。

残差は mean 0 の正規分布に従って見えます。これはモデルのための良いサニティチェックです。

学習されたパラメータを R の ‘lme4‘ で推定されたそれと比較することもできます。

student_effects_lme4 = pd.read_csv('data/insteval_student_ranefs_r.csv')
instructor_effects_lme4 = pd.read_csv('data/insteval_instructor_ranefs_r.csv')
dept_effects_lme4 = pd.read_csv('data/insteval_dept_ranefs_r.csv')

student_effects_edward = q_eta_s.mean().eval()
instructor_effects_edward = q_eta_d.mean().eval()
dept_effects_edward = q_eta_dept.mean().eval()

plt.title("Student Effects Comparison")
plt.xlim(-1, 1)
plt.ylim(-1, 1)
plt.xlabel("Student Effects from lme4")
plt.ylabel("Student Effects from edward")
plt.scatter(student_effects_lme4["(Intercept)"],
            student_effects_edward,
            alpha=0.25)
plt.show()


plt.title("Instructor Effects Comparison")
plt.xlim(-1.5, 1.5)
plt.ylim(-1.5, 1.5)
plt.xlabel("Instructor Effects from lme4")
plt.ylabel("Instructor Effects from edward")
plt.scatter(instructor_effects_lme4["(Intercept)"],
            instructor_effects_edward,
            alpha=0.25)
plt.show()

Great! student と instructor 効果の両者のための推定は ‘lme4‘ からのそれらと密接に適合しているようです。ここでは少しばかり異なるモデルをセットアップして(例えば、student, department, そして instructor 効果のための分散のように、全体的な mean は正規化され、これは ‘lme4′ のモデルには当てはまりません)、そして異なる推論メソッドを持ちますので、‘lme4‘ と正確に同じパラメータを見つけることは期待すべきではありません。しかしそれらが密接に一致することは再確認できます！

#  Add in the intercept from R and edward
dept_effects_and_intercept_lme4 = 3.28259 + dept_effects_lme4["(Intercept)"]
dept_effects_and_intercept_edward = mu.eval() + dept_effects_edward

plt.title("Departmental Effects Comparison")
plt.xlim(3.0, 3.5)
plt.ylim(3.0, 3.5)
plt.xlabel("Department Effects from lme4")
plt.ylabel("Department Effects from edward")
plt.scatter(dept_effects_and_intercept_lme4,
            dept_effects_and_intercept_edward,
            s=0.01 * train.dept.value_counts())
plt.show()

department 効果は ‘lme4‘ からのそれらと近くには一致しません。このためにはおそらく幾つかの理由があります :

  • 全体平均を正規化しますが、一方で ‘lme4‘ はそうではありません。これは Edward モデルに切片の一部を department 効果に入れることを引き起こします、分散を学習しますのでこれはより広範囲に変化することを可能にします。
  • We are using 80 estimate uses the whole ‘InstEval‘ data set
  • department 効果はモデルで最も弱く推定するのが困難です。

 

References

  • Bates, D., Mächler, M., Bolker, B., & Walker, S. (2015). Fitting Linear Mixed-Effects Models Using lme4. Journal of Statistical Software, 67(1), 1–48.
  • Gelman, A., & Hill, J. L. (2006). Data analysis using regression and multilevel/hierarchical models. Cambridge University Press.
 

以上