MLコミュニティデーは11月9日です! TensorFlow、JAXからの更新のために私たちに参加し、より多くの詳細をご覧ください

線形混合効果モデル

TensorFlow.orgで表示 GoogleColabで実行 GitHubでソースを表示 ノートブックをダウンロード

線形混合効果モデルは、構造化された線形関係をモデル化するための単純なアプローチです(Harville、1997; Laird and Ware、1982)。各データポイントは、グループに分類されたさまざまなタイプの入力と実数値の出力で構成されます。線形混合効果モデルは、階層モデルである:それはすべての個々のデータポイントについての推論を改善するために、グループ間で統計的強度を共有しています。

このチュートリアルでは、TensorFlowProbabilityの実際の例を使用して線形混合効果モデルを示します。私たちは、JointDistributionCoroutineとマルコフ連鎖モンテカルロ法(使用しますtfp.mcmc )モジュールを。

依存関係と前提条件

インポートして設定

物事を速くする!

飛び込む前に、このデモにGPUを使用していることを確認しましょう。

これを行うには、「ランタイム」->「ランタイムタイプの変更」->「ハードウェアアクセラレータ」->「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()))
SUCCESS: Found GPU: /device:GPU:0

データ

私たちは、使用InstEval人気からデータセットをlme4 Rパッケージ(ベイツら。、2015)。これは、コースとその評価評価のデータセットです。各コースは、次のようなメタデータが含まれstudentsinstructors 、および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()

私たちは、の観点でデータセットを設定する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ポイント高くなっています」。

このチュートリアルでは、次の効果を想定しています。

  • 固定効果: serviceserviceコースは講師の主要部門に属しているかどうかに対応するバイナリ共変量です。収集する追加データの量に関係なく、値は$ 0 $と$ 1 $のみになります。
  • ランダム効果: studentsinstructors 、およびdepartments 。コース評価評価の母集団からのより多くの観察を考えると、私たちは新しい学生、教師、または学部を検討している可能性があります。

Rのlme4パッケージ(Bates et al。、2015)の構文では、モデルは次のように要約できます。

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

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

以下では、このモデルをJointDistributionとして実装します。パラメータトラッキング(例えば、我々はすべて追跡するためのより良いサポートさせるにはtf.Variableしてmodel.trainable_variables )、我々は、モデルテンプレートを実装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_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),
      # 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 '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>)

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

統計的なツールとして、私たちはより良い参照するために、グラフに見えるかもしれません、例えば、そのintercepteffect_service条件付きで依存している与えられたratings 。プログラムがクラス、モジュール間の相互参照、および/またはサブルーチンを使用して記述されている場合、これをソースコードから確認するのは難しい場合があります。計算ツールとして、我々はまた、潜在変数がに流入気づくかもしれませんratingsを介して可変tf.gatherオプス。インデックス場合、これは、特定のハードウェア・アクセラレータ上のボトルネックになることがありTensor S高価です。グラフを視覚化すると、これがすぐにわかります。

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 et al。、1977; Wei and Tanner、1990)¹。マルコフ連鎖モンテカルロを実行して、ランダム効果(「Eステップ」)、およびパラメーターに関する期待値を最大化するために勾配降下を実行します(「Mステップ」)。

  • Eステップでは、ハミルトニアンモンテカルロ(HMC)を設定しました。現在の状態(学生、インストラクター、および部門の効果)を取得し、新しい状態を返します。新しい状態をTensorFlow変数に割り当てます。これは、HMCチェーンの状態を示します。

  • 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=-528062.5>
# 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

ウォームアップステージを実行します。このステージでは、1つのMCMCチェーンを複数の反復で実行して、トレーニングを事後確率質量内で初期化できるようにします。次に、トレーニングループを実行します。 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.754
Warm-Up Iteration: 999 Acceptance Rate: 0.707
Iteration:    0 Acceptance Rate: 1.000 Loss: 98220.266
Iteration:  500 Acceptance Rate: 0.703 Loss: 96003.969
Iteration: 1000 Acceptance Rate: 0.678 Loss: 95958.609
Iteration: 1499 Acceptance Rate: 0.685 Loss: 95921.891

あなたはまたにforループウォームアップ書くことができますtf.while_loopに、トレーニング段階をtf.scantf.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()

png

また、特定の潜在的な次元にわたるマルコフ連鎖モンテカルロアルゴリズムの軌跡を示すトレースプロットも使用します。以下では、特定のインストラクター効果が実際に初期状態から有意義に移行し、状態空間を探索していることがわかります。トレースプロットは、効果がインストラクター間で異なるが、混合動作が類似していることも示しています。

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

png

批判

上に、モデルを取り付けました。ここで、データを使用してその適合性を批判することを検討します。これにより、モデルを調査してよりよく理解することができます。そのような手法の1つは、各データポイントのモデルの予測とグラウンドトゥルースの差をプロットする残差プロットです。モデルが正しければ、それらの差は標準正規分布である必要があります。プロット内のこのパターンからの逸脱は、モデルの不適合を示します。

最初に評価に対する事後予測分布を形成することによって残差プロットを作成します。これにより、変量効果の事前分布が、与えられた事後トレーニングデータに置き換えられます。特に、モデルを前方に実行し、推定された事後平均を使用して、以前の変量効果への依存を傍受します。²

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まで実際には順序値であり、明らかになりました順序分布、または相対的な順序を破棄するのに十分なデータがある場合はカテゴリカルです。これは、上記のモデルに対する1行の変更です。同じ推論コードが適用されます。

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

png

モデルがどのように個々の予測を行うかを調べるために、学生、インストラクター、および部門の効果のヒストグラムを調べます。これにより、データポイントの特徴ベクトルの個々の要素が結果にどのように影響するかを理解できます。

当然のことながら、以下では、各学生が通常、インストラクターの評価評価にほとんど影響を与えないことがわかります。興味深いことに、インストラクターが所属する部門が大きな影響を及ぼしていることがわかります。

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

png

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

png

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

png

脚注

¹線形混合効果モデルは、周辺密度を分析的に計算できる特殊なケースです。このチュートリアルの目的のために、モンテカルロEMを示します。これは、可能性が通常ではなくカテゴリに拡張された場合など、非分析的周辺密度に簡単に適用されます。

²簡単にするために、モデルの1つのフォワードパスのみを使用して予測分布の平均を作成します。これは事後平均の条件付けによって行われ、線形混合効果モデルに有効です。ただし、これは一般的には有効ではありません。事後予測分布の平均は通常、扱いにくく、事後サンプルが与えられたモデルの複数のフォワードパスにわたって経験的平均を取る必要があります。

謝辞

このチュートリアルは、もともとエドワード1.0(で書かれたソース)。そのバージョンの作成と改訂に貢献してくれたすべての人に感謝します。

参考文献

  1. ダグラスベイツとマーティンメクラーとベンボルカーとスティーブウォーカー。 lme4を使用した線形混合効果モデルの近似。統計ソフトウェア、67のジャーナル(1):1-48、2015。

  2. アーサー・P・デンプスター、ナン・M・レアード、ドナルド・B・ルービン。 EMアルゴリズムによる不完全なデータからの最尤法。王立統計学会誌、シリーズB(方法論)、1-38、1977。

  3. アンドリューゲルマンとジェニファーヒル。回帰モデルとマルチレベル/階層モデルを使用したデータ分析。ケンブリッジ大学出版局、2006年。

  4. デビッドA.ハービル。分散成分の推定および関連する問題への最尤法。アメリカ統計学会誌、72(358):320から338、1977。

  5. マイケルI.ジョーダン。グラフィカルモデルの紹介。テクニカルレポート、2003年。

  6. ナン・M・レアードとジェームズ・ウェア。縦断的データの変量効果モデル。バイオメトリクス、963から974、1982。

  7. グレッグウェイとマーティンA.タナー。 EMアルゴリズムと貧乏人のデータ拡張アルゴリズムのモンテカルロ実装。アメリカ統計学会誌、699から704、1990。