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 lineares de efeitos mistos

Ver em TensorFlow.org

Executar no Google Colab

Ver fonte no GitHub

Baixar notebook

Um modelo linear de efeitos mistos é uma estratégia simples para modelar relações lineares estruturadas (Harville, 1997; Laird e Ware, 1982). Cada ponto de dados consiste de entradas de tipo variado, categorizadas em grupos, e uma saída com valor real. Um modelo linear de efeitos mistos é um modelo hierárquico: ele compartilha a força estatística entre os grupos para melhorar as inferências sobre qualquer ponto de dados específico.

Neste tutorial, vamos demonstrar os modelos lineares de efeitos mistos com um exemplo do mundo real no TensorFlow Probability. Usaremos os módulos JointDistributionCoroutine e Monte Carlo via Cadeias de Markov (tfp.mcmc).

Dependências e pré-requisitos

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

Vamos acelerar as coisas!

Antes de prosseguirmos, vamos usar uma GPU para a demonstração.

Para fazer isso, selecione "Runtime" -> "Change runtime type" (Alterar tipo de runtime) -> "Hardware accelerator" (Acelerador de hardware) -> "GPU".

O seguinte trecho de código verificará se temos acesso a uma 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.

Observação: se, por algum motivo, você não conseguir acessar uma GPU, este Colab ainda funcionará (mas o treinamento levará mais tempo).

Dados

Vamos usar o dataset InstEval do pacote lme4 popular do R (Bates et al., 2015). É um dataset de cursos e suas notas de avaliação. Cada curso inclui metadados, como students (alunos), instructors (professores) e departments (departamentos), e a variável de resposta de interesse é a nota de avaliação.

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

Vamos carregar e pré-processar o dataset. Reservamos 20% dos dados para podermos avaliar o modelo ajustado com pontos de dados ainda não vistos. Abaixo, visualizamos as primeiras linhas.

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)

Vamos configurar o dataset em termos de um dicionário de entradas de features e uma saída de labels correspondente às notas. Cada característica é codificada como inteiro, e cada rótulo (nota de avaliação) é codificado como número de ponto flutuante.

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

Um modelo linear típico pressupõe independência, em que qualquer par de pontos de dados têm uma relação linear constante. No dataset InstEval, as observações surgem em grupos, e cada uma pode ter declives e interceptos variáveis. Os modelos lineares de efeitos mistos, também conhecidos como modelos lineares hierárquicos ou modelos lineares multiníveis, capturam esse fenômeno (Gelman & Hill, 2006).

Veja alguns exemplos desse fenômeno:

• Alunos. As observações de um aluno não são independentes: alguns alunos podem dar notas sistematicamente baixas (ou altas) ao curso.

• Professores. As observações de um professor não são independentes: esperamos que bons professores geralmente tenham notas boas e que professores ruins geralmente tenham notas ruins.

• Departamentos. As observações de um departamento não são independentes: determinados departamentos geralmente podem ter materiais explícitos ou avaliações mais rígidas e, portanto, ter notas menores do que outros.

Para capturar esses fenômenos, lembre-se de que, para um dataset de $N\times D$ características $\mathbf{X}$ e $N$ rótulos $\mathbf{y}$, a regressão linear postula o modelo:

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

em que existe um vetor de declive $\beta\in\mathbb{R}^D$, intercepto $\alpha\in\mathbb{R}$ e ruído aleatório $\epsilon\sim\text{Normal}(\mathbf{0}, \mathbf{I})$. Digamos que $\beta$ e $\alpha$ sejam "efeitos fixos": são efeitos mantidos constantes entre a população de pontos de dados $(x, y)$. Uma formulação equivalente da equação como verossimilhança é $\mathbf{y} \sim \text{Normal}(\mathbf{X}\beta + \alpha, \mathbf{I})$. Essa verossimilhança é maximizada durante a inferência para encontrar estimativas de pontos de $\beta$ e $\alpha$ ajustados para os dados.

Um modelo linear de efeitos mistos estende a regressão linear como:

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

em que há um vetor de declive $\beta\in\mathbb{R}^P$, intercepto $\alpha\in\mathbb{R}$ e ruído aleatório $\epsilon\sim\text{Normal}(\mathbf{0}, \mathbf{I})$. Além disso, há um termo $\mathbf{Z}\eta$, em que $\mathbf{Z}$ é uma matriz de características, e $\eta\in\mathbb{R}^Q$ é um vetor de declives aleatórios; geralmente, $\eta$ tem distribuição normal com parâmetro de componente de variância $\sigma^2$. $\mathbf{Z}$ é formado particionando a matriz de características original $N\times D$ em termos de uma nova matriz $N\times P$ $\mathbf{X}$ e matriz $N\times Q$ $\mathbf{Z}$, em que $P + Q=D$: esta partição permite modelar as características separadamente usando os efeitos fixos $\beta$ e a variável latente $\eta$, respectivamente.

Digamos que as variáveis latentes $\eta$ sejam "efeitos aleatórios": são efeitos que variam entre a população (embora possam ser constantes entre subpopulações). Especificamente, como os efeitos aleatórios $\eta$ têm média igual a 0, a média dos rótulos de dados é capturada por $\mathbf{X}\beta + \alpha$. O componente de efeitos aleatórios $\mathbf{Z}\eta$ captura as variações nos dados. Por exemplo: "Instructor #54 is rated 1.4 points higher than the mean." (A nota do professor 54 está 1,4 ponto acima da média).

Neste tutorial, postulamos os seguintes efeitos:

• Efeitos fixos: service. service é uma covariável binária correspondente ao fato de o curso pertencer ou não ao departamento principal do professor. Não importa quantos dados adicionais coletarmos, ela só pode assumir os valores $0$ e $1$.

• Efeitos aleatórios: students, instructors e departments. Dadas mais observações da população de notas de avaliação dos cursos, podemos ter novos alunos, professores ou departamentos.

Na sintaxe do pacote lme4 do R (Bates et al., 2015), o modelo pode ser resumido como:

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

em que x denota um efeito fixo, (1|x) denota um efeito aleatório para x, e 1 denota um termo de intercepto.

Implementamos esse modelo abaixo como uma JointDistribution. Para ter um suporte melhor ao rastreamento de parâmetros (por exemplo: queremos rastrear todas as tf.Variables em model.trainable_variables), implementamos o template do 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 um programa gráfico probabilístico, também podemos visualizar a estrutura do modelo quanto ao seu grafo computacional. Esse grafo codifica o fluxo de dados entre as variáveis aleatórias do programa, explicitando suas relações em relação a um modelo gráfico (Jordan, 2003).

Como ferramenta estatística, podemos avaliar o grafo para ver melhor, por exemplo, que intercept e effect_service são condicionalmente dependentes dadas as notas ratings. Pode ser mais difícil ver isso a partir do código-fonte se o programa tiver sido escrito com classes, referências cruzadas entre módulos e/ou subrotinas. Como ferramenta computacional, também podemos observar o fluxo de variáveis latentes para a variável ratings pela operação tf.gather. Isso pode ser um gargalo em determinados aceleradores de hardware se a indexação de Tensors for cara do ponto de vista computacional; visualizar o grafo deixa isso aparente.

lmm_train.resolve_graph()

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

Estimativa de parâmetros

Fornecidos os dados, o objetivo da inferência é ajustar os efeitos fixos do modelo, declive $\beta$, intercepto $\alpha$ e parâmetro de componente de variância $\sigma^2$. O princípio da máxima verossimilhança formaliza essa tarefa como:

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

Neste tutorial, usamos o algoritmo de maximização de expectativa (EM) de Monte Carlo para maximizar essa densidade marginal (Dempster et al., 1977; Wei e Tanner, 1990).¹ Usamos o Monte Carlo via Cadeias de Markov para computar a expectativa da verossimilhança condicional com relação aos efeitos aleatórios ("passo E") e usamos o método do gradiente para maximizar a expectativa com relação aos parâmetros ("passo M"):

• Para o passo E, configuramos o Monte Carlo Hamiltoniano (HMC), que recebe um estado atual (os efeitos de aluno, professor e departamento) e retorna um novo estado. Atribuímos o novo estado a variáveis do TensorFlow, que denotarão o estado da cadeia do HMC.

• Para o passo M, usamos a amostra posterior do HMC para calcular uma estimativa da verossimilhança marginal sem bias como uma constante. Em seguida, aplicamos seu gradiente com relação aos parâmetros de interesse, o que gera um passo de método estocástico sem bias na verossimilhança marginal. Implementamos isso com o otimizador Adam do TensorFlow e minimizamos o negativo da distribuição 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

Fazemos uma etapa de preparação, que executa uma cadeia de MCMC com um determinado número de iterações para que o treinamento possa ser inicializado dentro da massa de probabilidade da distribuição posterior. Em seguida, executamos um loop de treinamento. Os passos E e M são executados conjuntamente, e os valores são registrados durante o treinamento.

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

Você também pode escrever o loop for de preparação em um tf.while_loop, e o passo de treinamento em um tf.scan ou tf.while_loop, para obter uma inferência ainda mais rápida. Por exemplo:

@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

Acima, não executamos o algoritmo até um limiar de convergência ser detectado. Para conferir se o treinamento foi sensível, verificamos se a perda de função tende à convergência ao longo das iterações de treinamento.

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

Também usamos um gráfico de traçado, que mostra a trajetória do algoritmo de Monte Carlo via Cadeias de Markov entre as dimensões latentes específicas. Podemos ver abaixo que os efeitos específicos de professor são transicionados do estado inicial de forma significativa e exploram o espaço de estados. O gráfico de traçado também indica que os efeitos diferem entre os professores, mas com comportamento de mistura 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

Acima, fizemos o ajuste do modelo. Agora, vamos criticar o ajuste usando os dados, que nos permite explorar e entender melhor o modelo. Essa técnica é um gráfico residual, que plota a diferença entre as previsões do modelo e a verdade fundamental para cada ponto de dados. Se o modelo estivesse correto, então a diferença deveria ter uma distribuição normal padrão; qualquer desvio desse padrão no gráfico indica um modelo mal ajustado.

Criamos o gráfico residual primeiro formando a distribuição preditiva posterior para as notas, que substitui a distribuição anterior para os efeitos aleatórios pela distribuição posterior após o fornecimento dos dados de treinamento. Especificamente, executamos o modelo para a frente e interceptamos sua dependência dos efeitos aleatórios anteriores com as médias 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()

Após uma inspeção visual, os residuais parecem ter uma distribuição normal padrão. Entretanto, o ajuste não está perfeito: há uma massa de probabilidade maior nas caudas do que em uma distribuição normal, indicando que o modelo pode melhorar seu ajuste ao abrandar as suposições de normalidade.

Especificamente, embora seja mais comum usar uma distribuição normal para modelar as notas no dataset InstEval, uma análise mais detalhada dos dados revela que as notas de avaliação dos cursos são, na verdade, valores ordinais de 1 a 5, o que indica que devemos usar uma distribuição ordinal ou até mesmo de categoria se tivermos dados suficientes para descartar a ordenação relativa. É necessário alterar somente uma linha do modelo acima; o mesmo código de inferência é aplicável.

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 ver como o modelo faz previsões específicas, avaliamos o histograma dos efeitos para alunos, professores e departamentos, o que nos permite compreender como os elementos individuais em um vetor de características de um ponto de dados tendem a influenciar o resultado.

Sem surpresa nenhuma, podemos ver abaixo que, tipicamente, cada aluno tem pouco efeito na nota de avaliação de um professor. É interessante notar que o departamento do professor tem um grande efeito.

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 de rodapé

¹ Os modelos lineares de efeitos mistos são um caso especial em que podemos computar analiticamente sua densidade marginal. Para a finalidade deste tutorial, demonstramos o EM de Monte Carlo, que se aplica melhor a densidades marginais não analíticas se a verossimilhança for estendida para ser de categoria em vez de normal, por exemplo.

² Por questões de simplicidade, formamos a média da distribuição preditiva usando somente um passo para frente do modelo, o que é feito por meio do condicionamento da média posterior e é válido para modelos lineares de efeitos mistos. Porém, isso não é válido de forma geral: a média da distribuição preditiva posterior é tipicamente intratável e requer o cálculo da média empírica entre diversos passos para frente do modelo dadas amostras posteriores.

Agradecimentos

Este tutorial foi escrito originalmente em Edward 1.0 (fonte). Agradecemos a todos os contribuidores por escreverem e revisarem essa versão.

Referências

1. Douglas Bates e Martin Machler, Ben Bolker e Steve Walker. Fitting Linear Mixed-Effects Models Using lme4 (Ajuste de modelos lineares de efeitos mistos usando o lme4). Journal of Statistical Software, 67(1):1-48, 2015.

2. Arthur P. Dempster, Nan M. Laird e Donald B. Rubin. Maximum likelihood from incomplete data via the EM algorithm (Máxima verossimilhança a partir de dados incompletos pelo algoritmo de EM). Journal of the Royal Statistical Society, Series B (Methodological), 1-38, 1977.

3. Andrew Gelman e Jennifer Hill. Data analysis using regression and multilevel/hierarchical models (Análise de dados usando modelos de regressão e hierárquicos/multiníveis). Cambridge University Press, 2006.

4. David A. Harville. Maximum likelihood approaches to variance component estimation and to related problems (Estratégias de máxima verossimilhança para estimativa de componente de variância e problemas relacionados). Journal of the American Statistical Association, 72(358):320-338, 1977.

5. Michael I. Jordan. An Introduction to Graphical Models (Introdução a modelos gráficos). Technical Report, 2003.

6. Nan M. Laird e James Ware. Random-effects models for longitudinal data (Modelos de efeitos aleatórios para dados longitudinais). Biometrics, 963-974, 1982.

7. Greg Wei e Martin A. Tanner. A Monte Carlo implementation of the EM algorithm and the poor man's data augmentation algorithms (Implementação de Monte Carlo do algoritmo de EM e algoritmos de dados com restrição de recursos). Journal of the American Statistical Association, 699-704, 1990.