Analogia entre modelos multinível e inferência bayesiana

Em uma regressão linear sabemos, pelo Teorema de Gauss-Markov, que um estimador não-enviesado de \(\beta\) que minimiza o soma do quadrado dos resíduos é \(\hat\beta = (X'X)^{-1}X'Y\) e sua variância é dada por \((X'X)^{-1} \sigma^2\) onde estimamos \(\sigma^2\) utilizando o quadrado dos resíduos, isto é \(\hat\sigma^2 = \frac{1}{n-k} \sum_{i = 1}^n (y_i - X_i \beta)^2\). Para se chegar neste resultado parte-se da premissa de homoscedasticidade e da independência dos erros.

O estimador de mínimos quadrados conta todos os \(n\) pontos de dados da mesma maneira. Se algum dado é considerado mais importante que outros (pense em uma amostra de survey desbalanceada) podemos capturar isso minimizando uma soma de quadrados ponderada, \(WSS = \sum_{i = 1}^n w_i(y_i - X_i \beta)^2\) e o estimador é

\[ \hat\beta^{WLS} = (X'WX)^{-1}X'WY\]

O que é \(W\)? Uma matriz diagonal cujos elementos são os pesos \(i\). O que essa matriz faz é incorporar a informação de que nem todo \(x_i\) é igual e “corrige” nossa matriz de dados para refletir essa informação. Se podemos incorporar a informação dos pesos de nossas observações também podemos incorporar outras informação usando uma generalização dessa matriz.

Podemos, por exemplo, incorporar algum conhecimento prévio que temos sobre o parâmetro. Suponha, por exemplo, que \(\beta\) seja o efeito da altura no log da renda e que acreditamos que esse efeito varie entre 0 e 5%. Podemos iaginar, então, que \(\beta\) tenha uma distribuição normal, com média 2.5% e desvio padrão 2.5%, ou seja \(\beta \sim N(0.025, 0.025^2)\). Podemos incorporar esse prior aumentando nossa matriz \(W\) incluindo \(\sigma_y^2 / 0.025^2\) que é a variância dos dados dividida pela variância do prior que, por sua vez, vai multiplicar uma matrix aumentada \(X\), que ganha um indicador 1 para o prior. Estimamos \(\beta\) por uma regressão de um vetor \(y\) aumentado com a média do prior (0.025).

O que esse peso nos diz? Se \(\sigma_y > 0.025\) então a distribuição do prior é mais informativa que qualquer ponto de dado então esse “ponto” terá mais influência e vice-versa. Se \(\sigma_y = 0.025\) o prior traz a mesma informação que outro ponto de dado, portanto o mesmo peso.

Podemos pensar os parâmetros de grupo \(\alpha_j\) e \(\beta_j\) em um modelo multinível como informações a priori e a relação entre sua variância e a variância dos dados vai nos dar o “peso” dos grupos. Podemos incorporar essa informação aumentando nossos dados.

Imagine um modelo onde só o intercepto varia, mas onde temos vários grupos. Temos então \(y \sim N(X \beta, \Sigma_y)\), onde \(\Sigma_y = \sigma^2 W_y^-1\) é a multiplicação de \(\sigma^2\) por uma matriz diagonal com elementos proporcionais ao inverso da variância. Neste caso \(\hat\beta^{WLS} = (X'W_yX)^{-1}X'W_yY\) e \(V_{\beta} = (X'W_yX)^{-1} \sigma^2\). Essa formuação permite heteroscedasticidade e se reduz ao OLS quando os pesos são todos iguais a um e \(W_y\) é a matriz identidade. Podemos modelar \(\beta\) como \(\beta \sim N(\mu_{\beta}, \Sigma_{\beta})\). Definindo \(W_{\beta} = \sigma_y^2 \Sigma_{\beta}^-1\) e sendo \(W_{\beta}\) uma matriz diagonal com elementos \(\sigma_y^2/\sigma_{\beta}^2\) podemos expressar o modelo multinível como a regressão de \(y*\) em \(X*\) com uma matriz de pesos \(W*\) onde

\[ y* = \binom{y}{\mu_{\beta}}, \text{ }X*= \binom{X}{I_g}, \text{ } W* = \begin{pmatrix} W_y & 0 \\ 0 & W_{\beta} \end{pmatrix} \]

Os dados aumentados correspondem à informação extra no modelo de \(\beta\). O aumento tem o efeito de puxar a estimativa de \(\beta\) em direção ao vetor de médis de grupo \(\mu_{\beta}\)

Partial Pooling

Quando temos dados com estrutura multinível, isto é, com observações no nível individual e de grupo e utilizamos um modelo de regressão que não leva em conta a informação dos grupos temos que se chama de complete pooling. Implicitamente esse modelo nos diz que \(\alpha_1 = \alpha_2 = \dots = \alpha_j = \mu_{\alpha}\) e o modelo se reduz a \(y_i \sim N(\mu_{\alpha}, \sigma_y^2)\). Estimamos \(\mu_{\alpha}\) pela média dos dados \(\hat y\). O contrário disso é chamado de no-pooling, isto é, quando fazemos uma regressão para cada grupo. Neste caso cada \(\alpha_j\) é estimado pela média dos dados no grupo \(\hat y_j\).

Modelos multinível tratam os coeficientes dos grupos como priors e incorporam a informação de grupo no modelo no nível individual. Essa modelagem oferece um compromisso entre esses dois extremos ao estimar \(\beta_j\) como uma média ponderada da estimativa do no-pooling, \(\bar{y_j} - \beta\bar{x_j}\) e a média \(\mu_{\beta}\):

\[ \alpha_j \approx \frac{\frac{n_j}{\sigma_y^2}}{\frac{n_j}{\sigma_y^2}+\frac{1}{\sigma_{\alpha}^2}} (\bar{y_j} - \beta\bar{x_j}) + \frac{\frac{1}{\sigma_{\alpha}^2}}{\frac{n_j}{\sigma_y^2}+\frac{1}{\sigma_{\alpha}^2}} \mu_{\alpha}\]

Para entender isso devemos lembrar que o erro de um modelo multinível é dado por \(e_{ig} = \nu_g + \eta_ig\), isto é, temos uma combinação do erro do nível do grupo com o erro individual. Como,por definição, o erro do nível de grupo não varia no grupo \((e_{ig}, e_{i'g})\) são correlacionados. Essa correlação pode ser levada em conta "corrigindo pelo que se chama de coeficiente de correlação intraclasse:

\[ \rho \equiv \frac{\sigma_{\alpha}^2}{\sigma_{\alpha}^2 + \sigma_{y}^2} \]

Modelos multinível permitem que usemos informações do nível de grupos para estimar com maior eficiência as médias em cada grupo. Quando temos preditores neste nível, isto é \(\beta_j \sim N(U_j\gamma, \sigma_{\beta}^2)\) o parâmetro \(\beta_j\) é empurrado para a média do grupo \(U_j\gamma\) e o erro da regressão de \(\beta_j = U_j\gamma + \eta_j\) no nível do grupo \(\eta_j\) é empurrado para zero. Quanto menor \(\sigma_{\beta}\) maior a pressão em direção à média. As estimativas dos grupos com poucas observações serão “encolhidas” em direção a esta média, emprestando a informação da estrutura dos grupos no que se chama de shrinkage1.

Podemos ver que o peso utilizado na shrinkage é o inverso de \(\rho\) ponderado pelo tamanho dos grupos. Outra relação interessante, com vimos mais acima, é a razão \(\sigma_y/\sigma_{\beta}\). Se o número de observações no grupo \(j\), \(n_j\) for maior que a razão entre as variâncias há menos pooling e vice-versa. Essa razão também nos diz a quantas “observações” equivale a variÂncia explicada por \(\sigma_{beta}\).

Estimando a incerteza usando a verossimilhança

Como vimos nas primerias aulas, no método de máxima verossimilhança estimamos a distribuição dos parâmetros \(Pr(\theta|y)\) de interesse a partir da distribuição de nossos dados \(Pr(y|\theta)\), ou seja

\[Pr(\theta|y) = Pr(y|\theta)Pr(\theta)\]

No lugar de estimar \(Pr(\theta|y)\) estimaremos se nosso valor hipótetico \(\tilde{\theta}\) é verossímel ou plausível diante das evidências/dados. Temos então:

\[\mathcal{L}(\tilde{\theta}|y) = k(\theta)f(y|\tilde\theta) \\ \propto f(y|\tilde\theta)\]

Para isso “sorteamos” valores de \(\theta\) que maximizem a verossimilhança, isto é, minimizem a diferença entre valores observados de \(y\) e os valores preditos \(\hat y\) que obtemos a partir dos \(\theta\) sorteados. Esses valores “sorteados” vão se aproximando de uma distribuição dos parâmetros. Uma vez obtida essa distribuição podemos calcular sua dispersão por meio de bootstrap, por exemplo.

Observem que na equação da verossimilhança temos um prior \(k(\theta)\) que é ignorado. Esse prior é a distribuição de \(\theta\) sem levar em conta \(y\), ou seja, levando em conta o que se sabe sobre \(\theta\) antes de se observar os dados. Na verossimilhança não assumimos nada sobre essa distribuição o que equivale a utilizar uma distribuição uniforme.

Na inferência Bayesiana multiplicamos a verossimilhança \(f(y|\tilde\theta)\) por uma distribuição a priori de \(\theta\) e as inferências são obtidas por valores escolhidos da distribuição resultante desta multiplicação, a distribuição posterior. Isso equivale a incluir informações numéricas em uma regressão aumentando nossos dados, como vimos acima.

A estimação dos parâmetros a partir da posterior é feita por métodos numéricos que percorrem todo o espaço paramétrico possibilitando descreve-lo completamente (identificando o ponto máximo e a dispersão). Na inferência Bayesiana os métodos numéricos mais utilizados vêm de uma família de métodos conhecida como Markov Chain Monte Carlo.

Markov Chain Monte Carlo - MCMC

A estimativa de modelos multinível pode ser feita por por máxima verossimilhança tradicional (ML) ou por máxima verossimilhança dos resíduos (REML) que tira a integral dos efeitos fixos e estima os componenetes da variância e, dadas as estimativas, das variâncias recupera-se o efeitos fixos.

A estimativa por ML ou REML pode se tornar instável quando temos muitos parâmetros, o que costuma ocorrer em modelos multinível. Para lidar com isso, outra maneira de se estimar os parâmetros de um modelo multinível é por meio de processos estocásticos. No lugar de maximizar algo (como a verossimilhança) trata-se de explorar o espaço paramétrico de forma aleatória a partir de um ponto inicial. Processo deste tipo são conhecidos como Markov Chain Monte Carlo (MCMC).

Um algoritmo que implementa MCMC é o Metropolis. Esse algoritmo parte de um valor \(x\) de um parâmetro e sorteia-se um novo valor \(x1\) a partir de uma distribuição arbitrária. Calcula-se uma taxa de aceitação \(a = f(x)/f(x1)\) onde \(f(.)\) é a verossimilhança. Gera-se um número aleatório de uma distribuição uniforme e se esse número for menor que \(a\) atualizamos o valor do parâmetro para \(x = x1\), do contrário mantemos \(x\).Esse processo faz com que se explore todo o espaço paramétrico às vezes aceitando a atualização e às vezes permanecendo no mesmo lugar.

Para ilustrar esse processo, McElreath (2020) faz um analogia com a história de um Rei que tem que visitar as ilhas de seu reino permanecendo em cada uma por um tempo proporcional à população e permancendo o menos possível no mar. O Rei segue uma regra simples: se mover para uma ilha ao lado se sua população for maior do que a da ilha atual ou, sendo menor, mover para essa ilha se a proporção da população destino sobre a população atual for maior que um valor sorteado de uma distribuição uniforme entre 0 e 1. Deste modo o rei viajaria por todas as ilhas e ficaria um tempo proporcional à população em cada uma. As ilhas seriam o equivalente aos parâmetros, a população equivalente às probabilidades de cada parâmetro e o tempo de estadia equivalente às amostras de cada parâmetros.

McElreath sugere o seguinte algoritmo para simular a jornada do Rei

num_weeks <- 1e5
positions <- rep(0, num_weeks)
current <- 10
for(i in 1:num_weeks){
  positions[i] <- current
  proposal <- current + sample(c(-1,1), size =1)
  if(proposal < 1) proposal <- 10
  if(proposal > 10) proposal <- 1
  prob_move <- proposal/current
  current <- ifelse(runif(1) < prob_move, proposal, current)
}

plot(1:100, positions[1:100]) #impossível dizer com antecedência onde o Rei vai estar

plot(table(positions)) # mas ele terá explorado todas as ilhas ficando em cada uma um tmpo proporcional ao tamanho.

Esse algoritmo é o ancestral direto de outros algoritmos da famíla MCMC. Um dos mais conhecidos e utilizados é o Gibbs Sampling onde a distribuição do valor de parâmetro proposto é ela mesma atualizada. O Gibbs sampling estima cada parâmetro por vez como uma ‘retirada’ de uma distribuição condicionando nos dados e nos valores dos outros parâmetros. Um exemlo está no cap. 18 de Gelman e Hill

[Exemplo Gibbs Sampler]

Softwares como o BUGS (Bayesian inference UsingGibbs Sampling) ou JAGS (Just Another Gibbs Sampler) implementam o Gibbs Sampler e é o utilizado por Gelman e Hill para ajustar os modelos multinível do livro de 2007.

Recentemente Gelman e outros parceiros desenvolveram o Stan (implementado no R pelo pacote rstan) que usa outra estratégia de estimação conhecida como Hamiltonian Monte Carlo. Neste método imagina-se que o vetor de parâmetros nos dá a posição de uma partícula no log da distribuição posterior (que como vimos acima pode parecer uma colina). Damos então um impulso nesta partícula, como um peteleco em uma bolinha de gude, e assim que a partícula perde momento e para, registramos o valor do parâmetro e recomeçamos (ver pg. 275 do McElreath). Esse método permite superar um problema com o Gibbs Sampler ou o Metrópolis. Como esses algoritmos exploram aos poucos a vizinhança eles podem ficar presos em certos platôs nos dados, zonas de alta correlação entre os parâmetros. Ao dar saltos o HMC evita esse problema.

No R o pacote rstanarm permite ajustar modelos multinível utilizando o HMC com uma sintaxe muito parecida com a que adotamos ao usar o lme4 e sem termos que nos preocupar como todo o processo de modelagem que o uso do Stan exige, se bem que perdemos com isso o controle e a transparência sobre todas as premissas envolvidas no modelo. Vamos continuar a estimar o modelo multinível para os votos em Bush por estado, mas agora usando o rstanarm e o o rstan.

[Exemplo logit multinível Gelman e Hill]

  1. podemos pensar esta relação entre no-poolinge complete pooling como o trade-off entre viés e variância, ver McElreath (2020)↩︎