Aula 10 - Modelos Multinível
IESP-UERJ
Modelos Multinível
Como vimos na última aula dados em painel ou séries temporais de cortes transversais (time series, cross section) permitem, por meio do que se chama efeitos fixos ou por meio de diferenças em diferenças, estabelecer relações de causalidade. Neste caso a estrutura dos dados nos ajuda a eliminar o viés de seleção de atributos que não variam no nível agregado. Nos modelos multinível a estrutura dos dados nos permite combinar as informações nos níveis individual e de grupo, isto é, entender o quanto e como a variação na relação entre variáveis se deve a características individuais e/ou a efeitos ‘contextuais’. As informações nos dois níveis se combinam permitindo maior eficiência ao estimar coeficientes para grupos específicos ou identificar efeitos no nível agregado.
Além de dados longitudinais, que agregam observações individuais no tempo, dados com estrutura multinível envolvem dados aninhados, sendo o exemplo mais usados o de alunos agregados em salas de aula que, por sua vez, também podem estar reunidas em diferentes escolas. Esse tipo de dado tem um caráter hierárquico. Os dados também podem ser não aninhados, como ocorre em modelos que inferem a posição ideológica de deputados por meio de seu padrão de votação nominal e onde temos votos agregados no nível dos deputados e do projeto de lei.
Podemos pensar nos modelos multinível como uma regressão onde incluímos indicadores para cada grupo. Foi o que fizemos na aula passada para estimar os modelos de efeitos fixos. Ao incluir estes indicadores temos um intercepto para cada grupo e, no caso do que chamamos de efeitos fixos, tratamos esses indicadores como constantes, características permanentes e imutáveisde cada indivíduo/unidade. Podemos também tratar esses indicadores como variáveis aleatórias atribuindo uma distribuição a eles, o que se chama de efeitos aleatórios. Se o centro da distribuição for zero os coeficientes de cada grupo seriam variações com relação ao intercepto. Se tivermos outros preditores além do intercepto podemos, também, permitir que as inclinações variem. Neste último caso teríamos modelos com intercepto e inclinação variáveis.
Uma forma de representar modelos multinível, onde o intercepto e a inclinação variam é
\[ Y_{ij} \sim N(\alpha_{j[i]} + \beta_{j[i]}, \sigma_{y})\]
Onde cada grupo \(j\) ao qual pertence o indivíduo \(i\) tem seu próprio intercepto (\(\alpha_{j[i]}\)) e inclinação (\(\beta_{j[i]}\)). Esses parâmetros, por sua vez têm, também, sua própria distribuição
\[ \left( \begin{array}{c} \alpha_j \\ \beta_j \end{array} \right) = N \left(\left(\begin{array}{c} \mu_{\alpha} \\ \mu_{\beta} \end{array} \right), \left(\begin{array}{c} \sigma_{\alpha}^2 & \rho\sigma_{\alpha}\sigma_{\beta} \\ \rho\sigma_{\alpha}\sigma_{\beta} & \sigma_{\beta} \end{array} \right) \right) \]
Embora esta notação torne explícito que em um modelo multinível tratamos os coeficientes de uma regressão como variáveis aleatórias atribuindo uma distribuição para elas, vou recorrer a uma outra notação para tornar o conceito mais intuitivo1. Usando um modelo de regressão clássica teríamos
\[ y_i = X_i \beta + \epsilon_i^{all}, \text{ } \epsilon_i^{all} \sim N(0, \Sigma)\]
Onde \(X\) é uma matriz com preditores no nível individual e de grupo (mais a constante) e \(\epsilon_i^{all}\) inclui os erros nos diferentes níveis. \(\Sigma\) é a matriz de covariância de \(\epsilon_i^{all}\) onde os efeitos aleatórios aparecem como correlações entre os erros. Um exemplo ajudará a deixar esta ideia mais concreta.
Seguindo Fox (2017), vamos usar como exemplo o survey “High School and Beyond” feito em 1982 com 7185 alunos norte americanos distribuídos em 160 escolas das quais 70 eram católicas e 90 eram públicas. Nosso interesse é na relação entre desempenho escolar (mensurado pelo desempenho em um teste de matemática) e o status socioeconômico da família do aluno (SES). No nosso modelo, o desempenho do aluno é afetado por um “efeito escola” além de caracetrísticas individuais. Esse efeito da escola também impacta a relação entre o desempenho e o SES. Seguindo a notação acima teríamos
\[ \text{matteste}_{ij} = \alpha_{j[i]} + \beta_{j[i]}cses_{ij} + \epsilon_{ij}\]
Onde \(cses - ses_{ij} - \overline{ses_i}\) ou seja o SES centralizado por escola. No nível agregado modelamos os coeficientes conforme características das escolas: se são católicas ou públicas e o nível médio do SES
\[ \begin{align} \alpha_{j[i]} &= \gamma_0^{\alpha} + \gamma_1^{\alpha}*\overline{ses_i} + \gamma_2^{\alpha}*setor_i + \omega_{\alpha} \\ \beta_{j[i]} &= \gamma_0^{\beta} + \gamma_1^{\beta}*\overline{ses_i} + \gamma_2^{\beta}*\overline{ses_i}^2 + \gamma_3^{\beta}*setor_i + \omega_{\beta} \end{align} \]
Onde \(sector\) é uma dummy com 1 para escolas católicas e 0 para públicas. Substituindo a equação do nível agregado na do nível individual temos:
\[ \begin{align} \text{matteste}_{ij} &= [\gamma_0^{\alpha} + \gamma_1^{\alpha}*\overline{ses_i} + \gamma_2^{\alpha}*setor_i + \omega_{\alpha}] + [\gamma_0^{\beta} + \gamma_1^{\beta}*\overline{ses_i} + \gamma_2^{\beta}*\overline{ses_i}^2 + \gamma_3^{\beta}*setor_i + \omega_{\beta}]*cses_{ij} + \epsilon_{ij} \\ &= \gamma_0^{\alpha} + \gamma_1^{\alpha}*\overline{ses_i} + \gamma_2^{\alpha}*setor_i + \gamma_0^{\beta}*cses_{ij} + \gamma_1^{\beta}*\overline{ses_i}*cses_{ij} + \gamma_2^{\beta}*\overline{ses_i}^2*cses_{ij} + \gamma_3^{\beta}*setor_i*cses_{ij} + \omega_{\alpha} + \omega_{\beta}*cses_{ij} + \epsilon_{ij} \end{align} (1)\]
Este modelo, temabém conhecido como modelo linear de efeitos mistos (\(LMM\) em inglês) embora complicado, explicta que os coeficientes das variáveis contextuais (do nível de grupo) aparecem como coeficientes de interações de níveis cruzados na equação final (por exemplo \(\gamma_1^{\beta}*\overline{ses_i}*cses_{ij}\)). Esses são os chamado efeitos fixos de um modelo misto. Os efeitos aletórios são os erros na última equação (\(\omega_{\alpha}, \omega_{\beta}\)), ou seja, eles dão os desvios em torno dos efeitos fixos.
Escrevendo \(\delta_{1j} = \omega_{\alpha}\) e \(\delta_{2j} = \omega_{\beta}\) definimos
\[ \delta_{kj} \sim N(0, \psi_k^2), C(\delta_{kj}, \delta_{k'j}) = \psi_{kk'}\]
\(\delta_{kj}\) e \(\delta_{kj}\) são independentes para \(j \neq j'\)
\[ \epsilon{ij} \sim N(0, \sigma_{\epsilon}^2\lambda_{iij}), C(\epsilon_{ij}, \epsilon_{i'j}) = \sigma_{\epsilon}^2\lambda_{ii'j}\] \(\epsilon_{ij}\) e \(\epsilon_{i'j'}\) são independentes para \(j \neq j'\) e \(\epsilon_{ij}\) é independente de \(\delta_{kj}\) para todo \(i\) e \(j\).
Deste modo teríamos
\[ \text{Var(matteste}_{ij}) = \sigma^2 + \psi_1^2 + cses_{ij}^2 \psi_2^2 + 2cses_{ij} \psi_{12} \\ \text{Cov(matteste}_{ij}, \text{matteste}_{i'j}) = \psi_1^2 + cses_{ij}cses_{i'j} \psi_2^2 + (cses_{ij} + cses_{i'j}) \psi_{12} \]
Seguindo Fox(2019, cap. 23) vou estimar o modelo multinível do desempenho em matemática com relação ao SES, SES médio da escola e setor da escola usando o R. Vamos usar as funções lme() do pacote nlme e lmer() do pacote lme4. Os dados que vamos usar estão disponíveis no pacote nlme.
library(car)
library(nlme)##
## Attaching package: 'nlme'## The following object is masked from 'package:lme4':
##
## lmList## The following object is masked from 'package:dplyr':
##
## collapse# os dados estão organizados em um banco com informações de alunos e outrocom informações das escolas
brief(MathAchieve)## 7185 x 6 nfnGroupedData (7180 rows and 1 columns omitted)
## School Minority Sex . . . MathAch MEANSES
## [c] [f] [f] [n] [n]
## 1 1224 No Female 5.876 -0.428
## 2 1224 No Female 19.708 -0.428
## 3 1224 No Male 20.349 -0.428
## . . .
## 7184 9586 No Female 16.241 0.627
## 7185 9586 No Female 22.733 0.627brief(MathAchSchool)## 160 x 7 data.frame (155 rows omitted)
## School Size Sector PRACAD DISCLIM HIMINTY MEANSES
## [f] [n] [f] [n] [n] [f] [n]
## 1224 1224 842 Public 0.35 1.597 0 -0.428
## 1288 1288 1855 Public 0.27 0.174 0 0.128
## 1296 1296 1719 Public 0.32 -0.137 1 -0.420
## . . .
## 9550 9550 1532 Public 0.45 0.791 0 0.059
## 9586 9586 262 Catholic 1.00 -2.416 0 0.627# vamos juntar os dados dos alunos e das escolas.
# primeiro calculando a média do SES por escola na base de alunos e juntando com a base de escolas
MathAchieve %>% group_by(School) %>% summarize(mean.ses = mean(SES)) -> temp
temp <- merge(MathAchSchool, temp, by = "School")
# Agora juntamos as variáveis do nívelda escola na base dos alunos
HSB <- merge(temp[,c("School", "Sector", "mean.ses")],
MathAchieve[,c("School","SES","MathAch")], by = "School")
names(HSB) <- tolower(names(HSB))
# Por fim centralizamos a variável SES com relação à média de cada escola
HSB$cses <- with(HSB, ses - mean.ses)
brief(HSB)## 7185 x 6 data.frame (7180 rows omitted)
## school sector mean.ses ses mathach cses
## [f] [f] [n] [n] [n] [n]
## 1 1224 Public -0.4343830 -1.528 5.876 -1.09361702
## 2 1224 Public -0.4343830 -0.588 19.708 -0.15361702
## 3 1224 Public -0.4343830 -0.528 20.349 -0.09361702
## . . .
## 7184 9586 Catholic 0.6211525 -0.008 16.241 -0.62915254
## 7185 9586 Catholic 0.6211525 0.792 22.733 0.17084746# Vamos agora examinar nossos dados procurando entender
# a) a relação entre desempenhoescolar e SES
# b) se essa relação varia por setor (católica vs pública)
# c) se essa relação varia aleatoriamente entre escolas de um mesmo setor
# Vamos primeiro olhar para a relação em uma amostra de 20 escolas de cada tipo
set.seed(12345)
cat <- with(HSB,
sample(unique(school[sector == "Catholic"]), 20))
cat.20 <- HSB[is.element(HSB$school, cat),]
pub <- with(HSB,
sample(unique(school[sector == "Public"]), 20))
pub.20 <- HSB[is.element(HSB$school, pub),]
# vamos usar a função xyplot() do pacote latice para visualizar a relação
# entre o desempenho em matemática e o SES centralizado nas nossas amostras de escolas
xyplot(mathach ~ cses | school, data = cat.20, main = "Católicas",
panel = function(x,y){
panel.points(x,y)
panel.lmline(x,y, ly=2, lwd =2, col = "darkgray")
panel.loess(x,y, span = 1, lwd = 2)
})
xyplot(mathach ~ cses | school, data = pub.20, main = "Públicas",
panel = function(x,y){
panel.points(x,y)
panel.lmline(x,y, ly=2, lwd =2, col = "darkgray")
panel.loess(x,y, span = 1, lwd = 2)
})
# O próximo passo é ajustar uma regressão de mathach em cses para cada escola em cada setor.
# A função lmList() faz isto juntando os resultadosde cada regressão em uma lista
cat.list <- lmList(mathach ~ cses | school, subset = sector == "Catholic", data = HSB)
pub.list <- lmList(mathach ~ cses | school, subset = sector == "Public", data = HSB)
brief(pub.list[[1]])## (Intercept) cses
## Estimate 9.72 2.51
## Std. Error 1.10 1.77
##
## Residual SD = 7.51 on 45 df, R-squared = 0.043cat.coef <- coef(cat.list)
head(cat.coef)## (Intercept) cses
## 1308 16.25550 0.1260242
## 1317 13.17769 1.2739128
## 1433 19.71914 1.8542944
## 1436 18.11161 1.6005617
## 1462 10.49556 -0.8288086
## 1477 14.22847 1.2306059pub.coef <- coef(pub.list)
# Podemos examinar a diferença na distribuição dos interceptos e das inclinações usando boxplots.
old.par <- par()
par(mfrow = c(1,2))
boxplot(cat.coef[,1], pub.coef[,1], main = "Intercepto", names = c("Católicas","Públicas"))
boxplot(cat.coef[,2], pub.coef[,2], main = "Inclinação", names = c("Católicas","Públicas"))
par <- old.par
# vamos agora examinar os efeitos no nível das escolas, istoé, vamos examinar a relação entre intercepto, inclinação e
# a média do SES nas escolas. Para isso vamos juntar a variável "mean.ses" na base com os coeficientes.
cat.coef <- merge(cat.coef, temp[,c("School", "mean.ses")],
by.x = "row.names", by.y = "School")
pub.coef <- merge(pub.coef, temp[,c("School", "mean.ses")],
by.x = "row.names", by.y = "School")
colnames(pub.coef)[c(1,2)] <- colnames(cat.coef)[c(1,2)] <- c("school", "intercept")
head(cat.coef)## school intercept cses mean.ses
## 1 1308 16.25550 0.1260242 0.5280000
## 2 1317 13.17769 1.2739128 0.3453333
## 3 1433 19.71914 1.8542944 0.7120000
## 4 1436 18.11161 1.6005617 0.5629091
## 5 1462 10.49556 -0.8288086 -0.6694035
## 6 1477 14.22847 1.2306059 0.1595806scatterplot(intercept ~ mean.ses, data =cat.coef, boxplots = FALSE, main = "Católicas")
scatterplot(cses ~ mean.ses, data =cat.coef, boxplots = FALSE, main = "Católicas")
scatterplot(intercept ~ mean.ses, data =pub.coef, boxplots = FALSE, main = "Públicas")
scatterplot(cses ~ mean.ses, data =pub.coef, boxplots = FALSE, main = "Públicas")
# O exame dos dados mostram que de fato parece haver uma relação linear entre mathach e cses nas escolas
# e que o intercepto e a inclinação desta relaçõ varia com as escolas. Também vemos que há diferenças
# no intercepto e inclinação nos diferentes tipos de escolas (católicas e públicas). Isso justifica a
# a especificação da equação (1)
# Primeiro vamos reordenar o fator"sector" atribuindo o valor zero para
# escolas públicas e um para escolas católicas
HSB$sector <- factor(HSB$sector, levels = c("Public","Catholic"))
#Agora ajustamos o modelo na equação (1)
hsb.lme.1 <- lme(mathach ~ mean.ses*cses + sector*cses, random = ~cses | school, data = HSB)
S(hsb.lme.1)## Linear mixed model fit by REML, Data: HSB
##
## Fixed Effects:
## Formula: mathach ~ mean.ses * cses + sector * cses
##
## Estimate Std.Error df t value Pr(>|t|)
## (Intercept) 12.1279 0.1993 7022 60.854 < 2e-16 ***
## mean.ses 5.3329 0.3692 157 14.445 < 2e-16 ***
## cses 2.9450 0.1556 7022 18.928 < 2e-16 ***
## sectorCatholic 1.2266 0.3063 157 4.005 9.55e-05 ***
## mean.ses:cses 1.0392 0.2989 7022 3.477 0.00051 ***
## cses:sectorCatholic -1.6427 0.2398 7022 -6.851 7.94e-12 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Random effects:
## Formula: ~cses | school
## Structure: General positive-definite, Log-Cholesky parametrization
## StdDev Corr
## (Intercept) 1.5426481 (Intr)
## cses 0.3178419 0.391
## Residual 6.0598009
##
##
## Number of Observations: 7185
## Number of Groups: 160
##
## logLik df AIC BIC
## -23251.83 10 46523.66 46592.45# Podemos também analisar esses resultados graficamente com o pacote effects
library(effects)## Use the command
## lattice::trellis.par.set(effectsTheme())
## to customize lattice options for effects plots.
## See ?effectsTheme for details.plot(predictorEffects(hsb.lme.1, ~cses,
xlevel=list(mean.ses=round(seq(-1.2,0.8, length=6),1))),
lines=list(multiline=TRUE, lwd = 4),
confint=list(style="bands"),
axes=list(x=list(rug=FALSE)),
lattice = list(key.args=list(space="right",columns=1)))
# O exame das variâncias e covariâncias nos informam sobre o quanto os efeitos aleatórios
# fazem sentido e quanta informação eles nos dão.
# Comparando um modelo com e outro sem os componentes da variância nos informa se o
# modelo de efeitos aleatórios faz diferença.
# Primeiro vamos comparar um modeo onde só o intercepto varia por grupo e as inclinações não
# isto é, a correlação entre pares no mesmogrupo é a mesma.
hsb.lme.2 <- update(hsb.lme.1, random = ~1 | school)
anova(hsb.lme.1, hsb.lme.2)## Model df AIC BIC logLik Test L.Ratio p-value
## hsb.lme.1 1 10 46523.66 46592.45 -23251.83
## hsb.lme.2 2 8 46520.79 46575.82 -23252.39 1 vs 2 1.124097 0.57# Podemos também comparar um modelo onde o interceptonão varia (mesmamédia de desempenho em cada setor)
# mas relação entre desempenho e SES varia por escola.
hsb.lme.3 <- update(hsb.lme.1, random = ~ cses - 1 | school)
anova(hsb.lme.1, hsb.lme.3)## Model df AIC BIC logLik Test L.Ratio p-value
## hsb.lme.1 1 10 46523.66 46592.45 -23251.83
## hsb.lme.3 2 8 46740.23 46795.26 -23362.11 1 vs 2 220.5634 <.0001# Além dos efeitos fixos, isto é, da relação entre desempenho em matemática e SES para toda a população
# nosso modelo permite ver como éa relação em cada escola. Para isso vamos usar a função lmer() do pacote lme4
hsb.lmer.1 <- lmer(mathach ~ mean.ses*cses + sector*cses + (cses | school), data = HSB)
display(hsb.lmer.1)## lmer(formula = mathach ~ mean.ses * cses + sector * cses + (cses |
## school), data = HSB)
## coef.est coef.se
## (Intercept) 12.13 0.20
## mean.ses 5.33 0.37
## cses 2.95 0.16
## sectorCatholic 1.23 0.31
## mean.ses:cses 1.04 0.30
## cses:sectorCatholic -1.64 0.24
##
## Error terms:
## Groups Name Std.Dev. Corr
## school (Intercept) 1.54
## cses 0.32 0.39
## Residual 6.06
## ---
## number of obs: 7185, groups: school, 160
## AIC = 46523.7, DIC = 46489.2
## deviance = 46496.4# Para estimar os coeficientes em cada grupo usamos:
coef(hsb.lmer.1)$school[1:5,]## (Intercept) mean.ses cses sectorCatholic mean.ses:cses
## 1224 12.05614 5.332872 2.940096 1.22658 1.039251
## 1288 12.58510 5.332872 2.985510 1.22658 1.039251
## 1296 10.41504 5.332872 2.750412 1.22658 1.039251
## 1308 12.15902 5.332872 2.930121 1.22658 1.039251
## 1317 10.59582 5.332872 2.812846 1.22658 1.039251
## cses:sectorCatholic
## 1224 -1.642682
## 1288 -1.642682
## 1296 -1.642682
## 1308 -1.642682
## 1317 -1.642682# Ou paraver apenas os efeitos aleatórios usamos
ranef(hsb.lme.1)[1:5,]## (Intercept) cses
## 1224 -0.07178617 -0.004945505
## 1288 0.45718314 0.040443092
## 1296 -1.71289856 -0.194549974
## 1308 0.03109671 -0.014924010
## 1317 -1.53213431 -0.132125770# obtendo osdesvios com relação aos efeitos fixos.Exemplo Logit
[Ver exemplo Gelman e Hill cap.14]
ver Gelman e Hill (2007:262) sobre as diferentes maneiras de escrever o mesmo modelo↩