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.

Modelos lineales de efectos mixtos

Ver en TensorFlow.org

Ejecutar en Google Colab

Ver fuente en GitHub

Descargar el bloc de notas

Un modelo lineal de efectos mixtos es un enfoque simple para modelar relaciones lineales estructuradas (Harville, 1997; Laird y Ware, 1982). Cada punto de datos consta de entradas de distintos tipos (clasificadas en grupos) y una salida de valor real. Un modelo lineal de efectos mixtos es un modelo jerárquico: comparte solidez estadística entre grupos para mejorar las inferencias sobre cualquier punto de datos individual.

En este tutorial, demostramos modelos lineales de efectos mixtos en TensorFlow Probability con un ejemplo del mundo real. Usaremos los módulos JointDistributionCoroutine y el método de Monte Carlo basado en cadenas de Markov (tfp.mcmc).

Dependencias y requisitos previos

#@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')

¡Aceleremos el proceso!

Antes de sumergirnos, asegurémonos de que estamos usando una GPU para esta demostración.

Para hacer esto, seleccione "Tiempo de ejecución" -> "Cambiar tipo de tiempo de ejecución" -> "Acelerador de hardware" -> "GPU".

El siguiente fragmento verificará si tenemos acceso a una 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.

Nota: Si, por alguna razón, no puede acceder a una GPU, esta colaboración seguirá funcionando. (La capacitación simplemente llevará más tiempo).

Datos

Usamos el conjunto de datos InstEval del popular paquete en R lme4 (Bates et al., 2015). Es un conjunto de datos de cursos y sus calificaciones de evaluación. Cada curso incluye metadatos como students, instructors y departments, y la variable de respuesta de interés es la calificación de evaluación.

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

Cargamos y preprocesamos el conjunto de datos. Guardamos el 20 % de los datos para poder evaluar nuestro modelo ajustado en puntos de datos invisibles. A continuación, visualizamos las primeras filas.

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)

Configuramos el conjunto de datos en términos de un diccionario de features de entradas y una salida labels correspondiente a las calificaciones. Cada característica se codifica como un número entero y cada etiqueta (calificación de evaluación) se codifica como un número de punto flotante.

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

Modelo

Un modelo lineal típico supone independencia, donde cualquier par de puntos de datos tiene una relación lineal constante. En el conjunto de datos InstEval, las observaciones surgen en grupos, cada uno de los cuales puede tener diferentes pendientes e intersecciones. Los modelos lineales de efectos mixtos, también conocidos como modelos lineales jerárquicos o modelos lineales multinivel, capturan este fenómeno (Gelman & Hill, 2006).

Ejemplos de este fenómeno incluyen lo siguiente:

• Estudiantes. Las observaciones de un estudiante no son independientes: algunos estudiantes pueden calificar sistemáticamente las clases con un valor bajo (o alto).

• Instructores. Las observaciones de un instructor no son independientes: esperamos que los buenos profesores generalmente tengan buenas calificaciones y los malos profesores generalmente tengan malas calificaciones.

• Departamentos. Las observaciones de un departamento no son independientes: ciertos departamentos generalmente pueden tener material poco variado o una clasificación más estricta y, por lo tanto, recibir calificaciones más bajas que otros.

Para capturar esto, recuerde que para un conjunto de datos de $N\times D$ características $\mathbf{X}$ y $N$ etiquetas $\mathbf{y}$, la regresión lineal postula el modelo

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

donde hay un vector de pendiente $\beta\in\mathbb{R}^D$, intersección $\alpha\in\mathbb{R}$ y ruido aleatorio $\epsilon\sim\text{Normal}(\mathbf{0}, \mathbf{I})$. Decimos que $\beta$ y $\alpha$ son "efectos fijos": son efectos que se mantienen constantes en toda la población de puntos de datos $(x, y)$. Una formulación equivalente de la ecuación como probabilidad es $\mathbf{y} \sim \text{Normal}(\mathbf{X}\beta + \alpha, \mathbf{I})$. Esta probabilidad se maximiza durante la inferencia para encontrar estimaciones puntuales de $\beta$ y $\alpha$ que se ajusten a los datos.

Un modelo lineal de efectos mixtos extiende la regresión lineal de este modo:

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

donde todavía hay un vector de pendiente $\beta\in\mathbb{R}^P$, intersección $\alpha\in\mathbb{R}$ y ruido aleatorio $\epsilon\sim\text{Normal}(\mathbf{0}, \mathbf{I})$. Además, existe un término $\mathbf{Z}\eta$, donde $\mathbf{Z}$ es una matriz de características y $\eta\in\mathbb{R}^Q$ es un vector de pendientes aleatorias; $\eta$ se distribuye normalmente con el parámetro del componente de varianza $\sigma^2$. $\mathbf{Z}$ se forma dividiendo la matriz de características original $N\times D$ en términos de una nueva matriz $N\times P$ $\mathbf{X}$ y una matriz $N\times Q$ $\mathbf{Z}$, donde $P + Q=D$: esta partición nos permite modelar las características por separado usando los efectos fijos $\beta$ y la variable latente $\eta$ respectivamente.

Decimos que las variables latentes $\eta$ son "efectos aleatorios": son efectos que varían a lo largo de la población (aunque pueden ser constantes a lo largo de las subpoblaciones). En particular, debido a que los efectos aleatorios $\eta$ tienen media 0, la media de la etiqueta de datos es capturada por $\mathbf{X}\beta + \alpha$. El componente de efectos aleatorios $\mathbf{Z}\eta$ captura variaciones en los datos: por ejemplo, "El instructor n.º 54 tiene una calificación de 1.4 puntos por encima de la media".

En este tutorial, planteamos los siguientes efectos:

• Efectos fijos: service. service es una covariable binaria que corresponde a si el curso pertenece al departamento principal del instructor. No importa cuántos datos adicionales recopilemos, solo pueden tomar valores $0$ y $1$.

• Efectos aleatorios: students, instructors y departments. Con más observaciones de la población de calificaciones de evaluación de cursos, podríamos estar observando nuevos estudiantes, profesores o departamentos.

En la sintaxis del paquete de R lme4 (Bates et al., 2015), el modelo se puede resumir de la siguiente manera:

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

donde x denota un efecto fijo, (1|x) denota un efecto aleatorio para x y 1 denota un término de intersección.

A continuación, implementamos este modelo como JointDistribution. Para tener un mejor soporte para el seguimiento de parámetros (por ejemplo, queremos hacer un seguimiento de todos los tf.Variable en model.trainable_variables), implementamos la plantilla del modelo como 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>)

Como programa gráfico probabilístico, también podemos visualizar la estructura del modelo en términos de su grafo computacional. Este grafo codifica el flujo de datos a través de las variables aleatorias del programa, haciendo explícitas sus relaciones en términos de un modelo gráfico (Jordan, 2003).

Como herramienta estadística, podríamos mirar el grafo para ver con mas claridad, por ejemplo, que intercept y effect_service dependen condicionalmente de ratings que se den; esto puede ser más difícil de ver en el código fuente si el programa está escrito con clases, referencias cruzadas entre módulos o subrutinas. Como herramienta computacional, también podríamos notar que las variables latentes fluyen hacia la variable ratings a través de las operaciones tf.gather. Esto puede ser un cuello de botella en ciertos aceleradores de hardware si la indexación de Tensor es costosa; visualizar el grafo hace que esto sea evidente.

lmm_train.resolve_graph()

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

Estimación de parámetros

Dados los datos, el objetivo de la inferencia es ajustar la pendiente de efectos fijos $\beta$, la intersección $\alpha$ y el parámetro del componente de varianza $\sigma^2$ modelo. El principio de máxima verosimilitud formaliza esta tarea de este modo:

$$\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.$$

En este tutorial, usamos el algoritmo de EM de Monte Carlo para maximizar esta densidad marginal (Dempster et al., 1977; Wei y Tanner, 1990)¹. Ejecutamos el método de Monte Carlo basado en cadenas de Markov para calcular la expectativa de la probabilidad condicional con respecto a los efectos aleatorios ("paso E"), y llevamos a cabo un descenso de gradiente para maximizar la expectativa con respecto a los parámetros ("paso M"):

• Para el paso E, configuramos el Hamiltoniano de Monte Carlo (HMC). Toma un estado actual (los efectos de estudiante, instructor y departamento) y devuelve un nuevo estado. Asignamos el nuevo estado a las variables de TensorFlow, que denotarán el estado de la cadena HMC.

• Para el paso M, usamos la muestra posterior del HMC para calcular una estimación no sesgada de la probabilidad marginal hasta una constante. Luego, aplicamos su gradiente con respecto a los parámetros de interés. Esto produce un paso de descenso estocástico imparcial sobre la probabilidad marginal. Lo implementamos con el optimizador Adam de TensorFlow y minimizamos el negativo del marginal.

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

Ejecutamos una etapa de calentamiento, que ejecuta una cadena del MCMC durante varias iteraciones para que el entrenamiento pueda inicializarse dentro de la masa de probabilidad posterior. Luego, ejecutamos un ciclo de entrenamiento. Ejecuta conjuntamente los pasos E y M, y registra los valores durante el entrenamiento.

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

También puede escribir el bucle for de preparación en un tf.while_loop y el paso de entrenamiento en un tf.scan o tf.while_loop para lograr una inferencia aún más rápida. Por ejemplo:

@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

Arriba, no ejecutamos el algoritmo hasta que se detectó un umbral de convergencia. Para comprobar si el entrenamiento fue válido, verificamos si la función de pérdida tiende a converger durante las iteraciones de entrenamiento.

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

También utilizamos un gráfico de seguimiento, que muestra la trayectoria del algoritmo de Monte Carlo basado en cadenas de Markov a través de dimensiones latentes específicas. A continuación, vemos que los efectos específicos del instructor efectivamente se alejan significativamente de su estado inicial y exploran el espacio de estados. El gráfico de seguimiento también indica que los efectos difieren entre los instructores pero con un comportamiento de mezcla similar.

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

Crítica

Arriba, instalamos el modelo. Ahora nos disponemos a analizar su ajuste a partir de los datos, lo que nos permitirá explorar y comprender mejor el modelo. Una de esas técnicas es un gráfico residual, que traza la diferencia entre las predicciones del modelo y la verdad fundamental para cada punto de datos. Si el modelo fuera correcto, entonces su diferencia debería tener una distribución normal estándar; cualquier desviación de este patrón en el gráfico indica un desajuste del modelo.

Para construir el gráfico residual, primero formamos la distribución predictiva posterior sobre las puntuaciones, que sustituye la distribución previa sobre los efectos aleatorios por su posterior dados los datos de entrenamiento. En particular, ejecutamos el modelo e interseccionamos su dependencia de efectos aleatorios previos con sus medias posteriores inferidas².

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

Tras una inspección visual, los residuales parecen tener una distribución normal. Sin embargo, el ajuste no es perfecto: hay una masa de probabilidad mayor en las colas que en una distribución normal, lo que indica que el modelo podría mejorar su ajuste si se flexibilizaran sus supuestos de normalidad.

En particular, aunque lo más común es usar una distribución normal para modelar calificaciones en el conjunto de datos InstEval, una mirada más cercana a los datos revela que las calificaciones de evaluación de cursos son en realidad valores ordinales del 1 al 5. Esto sugiere que deberíamos usar una distribución ordinal, o incluso categórica, si tenemos suficientes datos para descartar el orden relativo. Este es un cambio de una línea al modelo anterior; se aplica el mismo código de inferencia.

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

Para explorar cómo hace predicciones individuales el modelo, observamos el histograma de efectos para estudiantes, profesores y departamentos. Esto nos permite comprender cómo los elementos individuales en el vector de características de un punto de datos tienden a influir en el resultado.

No resulta sorpresivo comprobar a continuación que, cada estudiante normalmente tiene poco efecto en la calificación de evaluación de un instructor. Curiosamente, vemos que el departamento al que pertenece un instructor tiene un gran efecto.

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

Notas a pie de página

¹ Los modelos lineales de efectos mixtos son un caso especial en el que podemos calcular analíticamente su densidad marginal. Para los propósitos de este tutorial, demostramos el algoritmo de EM de Monte Carlo, que se aplica más fácilmente a densidades marginales no analíticas, como si la probabilidad se extendiera para ser categórica en lugar de normal.

² Para simplificar, formamos la media de la distribución predictiva con solo un paso hacia delante del modelo. Esto se hace mediante el condicionamiento de la media posterior y es válido para modelos lineales de efectos mixtos. Sin embargo, esto no es válido en términos generales: la media de la distribución predictiva posterior suele ser intratable y requiere tomar la media empírica en múltiples pasos hacia adelante del modelo dadas muestras posteriores.

Agradecimientos

Este tutorial fue escrito originalmente en Edward 1.0 (fuente). Agradecemos a todos los contribuyentes por escribir y revisar esa versión.

Referencias

1. Douglas Bates, Martin Machler, Ben Bolker y 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 y 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 y 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 y James Ware. Random-effects models for longitudinal data. Biometrics, 963-974, 1982.

7. Greg Wei y 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.