Aula 11 - Modelos Multinível
IESP-UERJ
Modelos Multinível
Todos os modelos paramétricos que vimos até agora partiam da premissa de que uma observação não trazia nenhuma informação sobre outras observações de nossa amostra. As observações seriam independentes uma das outras. Isso faz todo sentido quando as observações se referiam a medidas feitas sobre um mesmo objeto de forma independente por diferentes investigadores. Quando selecionamos uma amostra da população e tratamos cada observação como uma medida feita de modo independente estamos forçando estas observações a se comportarem como uma mesma observação medida por diversos observadores. Quando falamos em erro queremos dizer exatamente isso: diferenças nas obervações se devem a diferenças na forma de observá-las e não em diferenças inerentes a elas.
Na realidade nossas observações não são de fato independentes. Elas sempre tem alguma relação, direta ou indireta, uma com as outras. Isto quer dizer que quando tratamos as observações como independentes estamos ignorando as razões das relações entre as observaçãoes, estamos perdendo informação. Um tipo comum de relação entre as observações é o pertencimento comum a algum grupo (no espaço, tempo ou apenas conceitual). 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 usado 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 modelar a estrutura de clusters como uma regressão onde incluímos indicadores para cada grupo, captando as similaridades das observações em cada grupo e sua diferença com as que pertencem aos demais grupos. É o que se chama de modelos de efeitos fixos. Ao incluir estes indicadores temos um intercepto para cada grupo e tratamos esses indicadores como constantes, características próprias de cada grupo. No caso de um painel, por exemplo, o grupo é o indivíduo e essas características são permanentes e imutáveis o que possibilita identificar efeitos causais. Mas não estamos levando em conta a interdependência no erro das unidades no mesmo cluster.
Podemos também tratar essas características dos grupos, representadas por indicadores, como variáveis aleatórias atribuindo uma distribuição a eles. A isso chamamos de efeitos aleatórios ou efeitos mistos. Isso sigifica que além da distribuição de \(Y\) assumimos uma distribuição para o coeficiente do grupo, isto é adicionamos mais um erro aleatório ao modelo. Ao estipular uma distribuição do erro por cluster nosso modelo incorpora a informação de que nossas observações são oriundas de grupos diferentes que podem ter uma relação diferente com nossa variável resposta segundo características desses cluster.
Vamos examinar um modelo de efeitos mistos linear para clarear a ideia.
\[\begin{align*} & y_ij = \beta_1 + \beta_{2}X_{2ij} + \dots + \beta_{p}X_{pij} + \delta{1i}Z_{1ij} + \delta{qi}Z_{qij} + \epsilon_{ij} \\ &\delta_{ki} \sim N(0, \psi_{k}^2 \\ & C(\delta_{ki}, \delta_{k'i}) = \psi_{kk'} \\ & \delta_{ki}, \delta_{k'i} \text{ são independentes para }i \neq i' \\ & \epsilon_{ij} \sim N(0, \sigma_{\epsilon}^2 \lambda_{ijj}) \\ & C(\epsilon_{ij}. \epsilon_{i'j'} \text{ são independentes para }i \neq i') \\ & \delta_{ki}, \epsilon_{i'j} \text{ são independentes para todo } i,i',k,j (\text{incluindo }i = i') \end{align*}\]
Onde:
\(Y_{ij}\) é o valor da variável resposta para a \(j^a\) de \(n_i\) observações no \(i^o\) de \(m\) clusters.
\(\beta\)s são coeficientes de efeitos fixos e \(X\)s são os regressores cujos efeitos são fixos.
\(\delta\)s são coeficiente de efeitos aleatórios e \(Z\)s são os regresores cujos efeitos são aleatórios (têm uma distribuição). Podemos pensar em \(\delta\) não como parâmetros, mas como variáveis aleatórias não observadas.
\(\psi_{k}^2\) são as variâncias e \(\psi_{kk'}\) as covariâncias entre os efeitos aleatórios, que assumimos ser constantes dentros dos clusters.
\(\epsilon_{ij}\) são erros para a observação \(j\) dentro do cluster \(i\).
\(\sigma_{\epsilon}^2 \lambda_{ijj}\) são as covariâncias entre os erros no grupo \(i\). Se as observações em cada grupo são escolhidas de maneira independente então \(\lambda_{ijj} = 1\) e \(\lambda_{ijj} = 0\). Se não forem, como acontece em estudos longitudinais, quando o grupo é o indivíduo, cada observação e esse indivíduo em um dado momento do tempo e cada observação é serialmente correlacionada com a anterior, a estrutura de \(\lambda\) deve ser especificada.
Existem duas propriedades essenciais que distinguem o modelo linear padrão do modelo linear de efeitos mistos: (1) Existem efeitos aleatórios estruturados no nível do cluster no modelo misto (os \(\delta_{ki}\)), em adição ao erro no nível individual (os \(\epsilon_{ij}\)), e (2) o modelo misto pode acomodar certas formas de variância de erro que não são constantes além de dependências entre or erros.
Exemplo
Seguindo Fox (2019, cap. 23), 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). Suspeitamos que o desempenho do aluno é afetado por um “efeito escola” além de caracetrísticas individuais e que esse efeito se deve ao tipo de escola (católica vs pública).
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# Analisando os dados e identificando sua natureza multinível.
# os dados estão organizados em um banco com informações de alunos e outro com 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)
})
# Os gráficos indicam um impacto do tipo de escola no desempenho. Para refinar mais a análise o próximo passo é ajustar uma regressão linear simples 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"),
cex.main = 0.7,
cex.axis = 0.7)
boxplot(cat.coef[,2], pub.coef[,2],
main = "Inclinação",
names = c("Católicas","Públicas"),
cex.main = 0.7,
cex.axis = 0.7)
par <- old.par
# 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ção 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 de um modelo multinível.No nível individual temos o seguinte modelo:
\[\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.
Agora vamos modelar as variações entre os grupos. 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 i} \\ \beta_{j[i]} &= \gamma_0^{\beta} + \gamma_1^{\beta} \overline{ses_i} + \gamma_3^{\beta} setor_i + \omega_{\beta i} \end{align}\]
Onde \(sector\) é um indicador com 1 para escolas católicas e 2 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} \\ \\ \text{matteste}_{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_3^{\beta}setor_i cses_{ij} \\ &+ \omega_{\alpha} + \omega_{\beta} cses_{ij} + \epsilon_{ij} \end{align*}\]
Este modelo, també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})\)
\(Cov(\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} \]
Continuando com Fox 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.
# Primeiro vamos reordenar o fator"sector" atribuindo o valor um para
# escolas públicas e dois para escolas católicas
HSB$sector <- factor(HSB$sector, levels = c("Public","Catholic"))
# Agora ajustamos o modelo na equação composta
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,cex = 0.7),y=list(cex = 0.7)),
lattice = list(key.args=list(space="right",columns=1,cex.title=0.8),strip=list(factor.names=TRUE,values=TRUE,cex=0.6)))
# 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 modelo onde só o intercepto varia por grupo e as inclinações não
# isto é, a correlação entre pares no mesmo grupo é a mesma.
hsb.lme.2 <- update(hsb.lme.1, random = ~cses - 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 46740.23 46795.26 -23362.11 1 vs 2 220.5634 <.0001# Podemos também comparar um modelo onde o intercepto não varia (mesma média de desempenho em cada setor)
# mas relação entre desempenho e SES varia por escola.
hsb.lme.3 <- update(hsb.lme.1, random = ~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 46520.79 46575.82 -23252.39 1 vs 2 1.124097 0.57# 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.