1 Modelos Mistos

1.1 O que são Modelos Mistos

Um pressuposto importante para regressões lineares é o pressuposto da independência. Isso implica que os nossos preditores são independentes entre si; ou seja, um valor não impacta o outro. No entanto, existem situações em que esse pressuposto é quebrado, resultando em correlação forte entre preditores (ou seja, X1 está correlacionado/depende de uma outra variável o que causa essa correlação, que pode ser causada por fatores espaciais, temporais ou pela presença de clusters.

1.2 O que são Clusters

Clusters referem-se a agrupamentos de dados que apresentam um nível de organização mais elevado em um estudo. Esses agrupamentos podem ser categóricos ou contínuos e resultam em quebras na variabilidade dos dados.

1.3 Efeitos Aleatórios e Mistos

Quando lidamos com dados que apresentam correlações fortes ou dependência entre preditores, empregamos modelos multivariados de análise ou modelos mistos. Esses modelos permitem separar os efeitos fixos e aleatórios decorrentes dos clusters.

  • Os efeitos fixos são aquelas variáveis que acreditamos não serem significativas para explicar a variância entre diferentes clusters. Em outras palavras, não há evidências de que as diferenças entre os clusters tenham um impacto significativo nas variáveis dependentes.

  • Por outro lado, os efeitos aleatórios (random effects) são variáveis que, de fato, explicam significativamente a variância dentro dos clusters. Isso significa que as variâncias para cada grupo são diferentes, e essas diferenças são importantes para a modelagem.

Os modelos mistos, portanto, produz uma regressão para cada cluster com interceptos ou slopes aleatórios (ou ambos aleatórios) - que chamamos de efeitos aleatórios - e uma reta “geral média” que é o fixed effect - a média dos random effects. Sumarizando, modelos mistos geram dois tipos de efeitos: um efeito fixo que representa a média dos coeficientes dos efeitos aleatórios, ou a própria estimativa, caso não haja um efeito aleatório esperado.

Dessa forma, a escolha entre modelos fixos e mistos deve ser fundamentada na estrutura dos dados e na presença de autocorrelação, garantindo que a análise seja robusta e confiável.

2 Dataset de Estudo

Para ilustrar o funcionamento dos Modelos Mistos, vamos modelar dados do conjunto de dados Math do pacote flexplot. Este pacote oferece funcionalidades robustas que serão amplamente utilizadas neste projeto.

Para instalar o pacote flexplot, basta executar o seguinte código em seu ambiente R:

devtools::install_github("dustinfife/flexplot", ref="development") 
require(flexplot)

2.1 Conjunto de Dados Math do Flexplot

O conjunto de dados Math do pacote flexplot contém informações sobre o desempenho dos alunos em matemática. Ele é frequentemente utilizado para demonstrar diferentes técnicas de modelagem e visualização de dados.

2.2 Estrutura do Conjunto de Dados

A tabela abaixo apresenta as principais variáveis do conjunto de dados Math:

Variável Tipo de Dados Descrição
School Categórica Nome da escola onde o aluno está matriculado
Minority Categórica Indica se o aluno é parte de uma minoria (sim/não)
Sex Categórica Gênero do aluno (masculino, feminino)
SES Categórica Status socioeconômico do aluno
MathAch Numérica Nota do aluno em matemática (0 a 100)
MEANSES Numérica Média do status socioeconômico dos alunos da escola

2.3 Carregando o Conjunto de Dados

Para carregar o conjunto de dados Math, utilize o seguinte código:

library(flexplot)
data(Math)
  • Para nosso projeto vamos prever a nota do aluno em matemática (MathAch = variável dependente/resposta)

2.3.1 Pacote Para Regressões

Para usarmos Mixed Models Lineares (note que não se aplica à modelos mistos generalizados), usamos o pacote lme4 do R

library(lme4)

Onde a fórmula é

\[ y \sim x + (x \mid \text{cluster}). \]

O grupo \((x \mid \text{cluster})\) especifica para o modelo qual variável tem o efeito aleatório por qual cluster.

O efeito aleatório pode ser aninhado, como em

\[ (x + x_1 + \ldots + x_n \mid \text{Cluster}), \]

ou pode ser representado apenas com o intercepto

\[ (1 \mid \text{Cluster}). \]

Se algum preditor fixo for passado e você quiser manter o intercepto como fixo, use:

\[ (-1 + x_1 \mid \text{Cluster}). \]

2.4 Referências

Este trabalho é baseado no canal do youtube Quant Psych cujo conteúdo é excelente para aprofundar melhor as ideias aqui aplicadas. Recomendo fortemente os vídeos do professor https://www.youtube.com/@QuantPsych

devtools::install_github("dustinfife/flexplot", ref="development") 
require(flexplot)
set.seed(123)

data(math)

head(math, n=20)
##    School Minority    Sex    SES MathAch MEANSES
## 1    1224       No Female -1.528   5.876  -0.428
## 2    1224       No Female -0.588  19.708  -0.428
## 3    1224       No   Male -0.528  20.349  -0.428
## 4    1224       No   Male -0.668   8.781  -0.428
## 5    1224       No   Male -0.158  17.898  -0.428
## 6    1224       No   Male  0.022   4.583  -0.428
## 7    1224       No Female -0.618  -2.832  -0.428
## 8    1224       No   Male -0.998   0.523  -0.428
## 9    1224       No Female -0.888   1.527  -0.428
## 10   1224       No   Male -0.458  21.521  -0.428
## 11   1224       No Female -1.448   9.475  -0.428
## 12   1224       No Female -0.658  16.057  -0.428
## 13   1224       No   Male -0.468  21.178  -0.428
## 14   1224       No Female -0.988  20.178  -0.428
## 15   1224       No   Male  0.332  20.349  -0.428
## 16   1224       No Female -0.678  20.508  -0.428
## 17   1224       No   Male -0.298  19.338  -0.428
## 18   1224      Yes   Male -1.528   4.145  -0.428
## 19   1224       No Female  0.042   2.927  -0.428
## 20   1224       No   Male -0.078  16.405  -0.428
#Transformando School em uma variável Factor
math = math %>% 
  mutate(School = factor(School))


#Amostrando 3 Escolas de forma "aleatória" (há viés) para verificar plots a seguir
#Se não fosse feito cada plot teria 160 escolas no eixo x (impossível de visualizar)
math2 = math %>% 
  mutate(School = factor(School)) %>%
  filter(School == 2277 | School == 2458 | School == 2336)



flexplot(MathAch ~ School, data = math2) #boxplot (formula = x ~ y + quebras)

#Aqui estamos usando a fórmula (y ~ x1 + A | Minotiry)
#onde x1 é a variável de foco, A nossos clusters (vamos ver o porquê) 
#e Minotiry são a variável que esperamos que 
flexplot(MathAch ~ SES + School | Minority, data = math2, method='lm')

#Podemos notar diferentes variâncias por escolas
#os slopes (ou ângulos) de cada relação de x e Y depende de qual escola estamos analisando

#portanto a escola é nosso random effect em X - dado que os dados foram amostrados por
#escolas no estudo, portanto nossos clusters

#Minority parece também sofrer influência da escola, dado a escola Minotiry tem
#comportamentos diferentes, portanto School também pode ser um efeito aleatório para Minotiry

flexplot(MathAch ~ SES + Minority, data = math2, method='lm')

#Veja que não temos, de fato nenhum efeito interativo de Minotiry com a relação entre
#SES e MathAch

3 Porque School é um random effect e não interativo?

No meu estudo, a escola é um efeito aleatório, pois atua como um cluster de amostragem, criando dependência entre x (SES) e y (MathAch). As escolas influenciam de forma aleatória essa relação. Em outras palavras, temos grupos (clusters) nos dados que geram relações aleatórias entre nosso preditor e nossa resposta.

Por outro lado, as relações interativas são influenciadas por fatores que não causam variação nos slopes (ângulos) de maneira aleatória, mas sim de forma determinística. Portanto, esses efeitos não são causados por grupos ou clusters, mas sim por relações diretas entre as variáveis.

3.1 Identificação de Clusters e Modelagem Baseline

Ao trabalhar com modelos mistos, o primeiro passo é identificar o cluster. Um cluster é o nível mais alto de granularidade onde ele não varia em relação a nenhuma outra categoria ou valor, ou seja, as outras variáveis variam em função deste cluster. Um outra forma de pensá-lo seria entendê-lo como uma amostragem de uma população. Aqui Escolas seriam amostras de todas a população (alunos)

3.1.1 Identificação do Cluster

Nos nossos dados, utilizamos o identificador School, que é a variável que altera os dados. Portanto, School é o nosso cluster.

3.1.1.1 Dicas Importantes:

  1. Variáveis Constantes em Clusters: Variáveis que não mudam dentro de um cluster, como MEANSES, não devem ser consideradas como efeitos aleatórios. Ao aplicá-las, você pode receber avisos de singularidade. Evite essa prática.

  2. Hierarquia das Categorias: Clusters não podem ser categorias que variam dentro de uma categoria mais alta em seu conjunto de dados.

  3. Teoria dos Preditores: Decidir se um preditor terá efeito aleatório muitas vezes depende da teoria que está sendo testada. Por exemplo, acreditamos que um preditor terá seu coeficiente angular igual para todos os seus clusters?

3.1.2 Aplicação do Modelo de Baseline

Uma vez definido o cluster, podemos aplicar um modelo de baseline, cuja ideia é verificar o valor do Índice de Concordância Intraclasse (ICC). O ICC nos informa o quanto da variância da variável resposta (y) é explicada pelos clusters simples.

O modelo de baseline pode ser definido da seguinte maneira:

baseline <- lmer(y ~ 1 + (1 | School), data = data)
baseline = lmer(MathAch~1 + (1|School), data=math)


icc(baseline)
## # Intraclass Correlation Coefficient
## 
##     Adjusted ICC: 0.180
##   Unadjusted ICC: 0.180
#18% do modelo é devido aos clusters, então podemos definir como importante o uso de Mixed #Models

Portanto, podemos notar que os clusters explicam quase 28% da variância dos dados. Observamos que os dados são autocorrelacionados, uma vez que a escola influencia as notas dos alunos, ou o nível de matemática, considerando que cada escola possui programas de estudos, objetivos e professores diversos.

Nosso próximo passo é adicionar variáveis preditoras e definir se seu efeito é fixo (fixed effect), ou seja, um slope fixo, ou efeito aleatório (random effect) criado pelo cluster.

Vamos começar com a variável SES

Se o efeito for fixo, isso implica que o status socioeconômico do aluno influencia a nota de matemática independentemente da escola. Na teoria, podemos imaginar que esse efeito não diferirá entre as escolas; no entanto, essa hipótese poderá ser testada ao comparar modelos.

#SES tem efeito aleatório ou fixo?

fixed_slopes = lmer(MathAch ~ SES + (1|School), data=math)
random_slopes = lmer(MathAch ~ SES + (SES|School), data=math)

visualize(random_slopes, plot='model')

#Veja como mesmo com random effect o modelo manteve slopes muito próximos
#A linha preta é o fixed effect (a média dos efeitos aleatórios) e cada linha colorida
#representa o fit para cada School ou random effects



compare.fits(MathAch ~ SES | School, data = math, 
             fixed_slopes,
             random_slopes,
             clusters = 3) #clusters = 3 para ver apenas 3 escolas (menor poluição visual)

#Ambos os fits são muitos próximos no slope, isso significa que as visualmente efeitos fixos
#e aleatórios parecem não divergir por escola dado a variável SES, isso #significa que não
#há, possivelmente, uma explicação significativa da variância de SES dado a Escola do aluno,
#ou seja, podemos acreditar, nesse caso, que os #slopes (ou seja, o valor da nota de
#matemática está igualmente relacionada à SES quando é igual a 1 para todas as escolas) não
#sofrem influência do cluster

#Comparando modelos de fixed_slopes para random_slopes
model.comparison(fixed_slopes, random_slopes)
## $statistics
##                    aic      bic bayes.factor     p
## fixed_slopes  46653.17 46680.69      661.309 0.104
## random_slopes 46652.40 46693.68        0.002      
## 
## $predicted_differences
##    0%   25%   50%   75%  100% 
## 0.000 0.024 0.070 0.162 2.008 
## 
## $r_squared_change
##     Residual  (Intercept) 
## -0.005545779  0.012541409
#O model.comparison comprova nossa tese acima dado o visual, porém vamos explorar um pouco
#melhor

4 Pontos Importantes do model.comparison (r_squared_change)

4.1 Interpretação do R² Change

4.1.1 Residual

O valor de -0.005545779 para o Residual indica que, ao mudar de um modelo de efeitos aleatórios (onde o SES tem um slope aleatório) para um modelo com apenas intercepto aleatório, o modelo com efeito aleatório não melhora a explicação da variância de MathAch (\(y\)). Na verdade, ele está até piorando ligeiramente a capacidade de explicar essa variância.

Isso sugere que o modelo de efeitos aleatórios para SES não está capturando uma variação significativa que não já é explicada pelo modelo de intercepto fixo.

4.1.2 Intercept

O valor de 0.012541409 para o Intercept mostra que, mesmo quando se considera apenas a variação do intercepto, há uma leve melhora na explicação da variância. Isso implica que, embora o modelo de intercepto aleatório adicione alguma variância explicada, essa contribuição é pequena.

Portanto tanto visualmente como por estimativa, nosso modelo não ganha muito com SES sendo considerado como random effect, sua variância não é bem explicada por slopes aleatórios, porém op intercept no modelo continua sendo bem utilizado

4.2 2. Pontos Importantes do model.comparison (predicted_differences)

4.2.1 Interpretação do Predict Difference

Nesta seção, observamos as diferenças de 0 no primeiro quartil, 0.02 no 25º quartil, 0.07 no 50º quartil, 0.1 no 75º quartil e 2 no 100º quartil. Isso pode implicar que o modelo apresenta diferenças de previsões mais altas à medida que aumentamos os quartis. Em outras palavras, as previsões tendem a mudar significativamente em quartis altos. Porém as diferenças são pequenas entre os modelos

4.3 3. Pontos Importantes do model.comparison (statistics)

4.3.1 Interpretação das Estatisticas

Modelos com AIC e BIC menores comumente são preferíveis mas o foco é o bayer.factor que indica compatibilidade do modelo com os dados, aqui os dados estão a favor do modelo de slope para SES fixo com 661.309

#modelo com variável Minority e SES fixo
minority = lmer(MathAch ~ SES + Minority + (1 | School), data=math)

compare.fits(MathAch ~ SES | Minority + School, data = math,
             minority, fixed_slopes)

model.comparison(fixed_slopes, minority)
## $statistics
##                   aic      bic bayes.factor      p
## fixed_slopes 46653.17 46680.69 0.000000e+00 <2e-16
## minority     46459.31 46493.71 4.000398e+40       
## 
## $predicted_differences
##    0%   25%   50%   75%  100% 
## 0.000 0.200 0.458 0.962 3.007 
## 
## $r_squared_change
##    Residual (Intercept) 
##  0.02393021  0.17472191
#Veja o ganho em explicação da variância dos resíduos e interceptos
#Ou seja, nosso modelo com minority explica melhor os dados do que apenas com o nível
#socioeconomico

Agora vamos complementar nosso modelo com Minotiry sofrendo efeito aleatório de School

minority_random = lmer(MathAch ~ SES + Minority + (1 + Minority| School), data=math)

model.comparison(minority_random,minority)
## refitting model(s) with ML (instead of REML)
## $statistics
##                      aic      bic bayes.factor     p
## minority_random 46450.75 46498.90        0.074 0.002
## minority        46459.31 46493.71       13.434      
## 
## $predicted_differences
##    0%   25%   50%   75%  100% 
## 0.000 0.028 0.069 0.200 1.419 
## 
## $r_squared_change
##    Residual (Intercept) 
## 0.005602589 0.138238586
#Nosso modelo prior seria o reduzido, sem Minotiry como random effect, veja as mudanças na
#variância capturas pelo modelo sem a relação, é melhor, por mais que tenhamos mais resíduos,
#mas a mudança deste é bem pequena

summary(minority)
## Linear mixed model fit by REML ['lmerMod']
## Formula: MathAch ~ SES + Minority + (1 | School)
##    Data: math
## 
## REML criterion at convergence: 46449.3
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1721 -0.7236  0.0288  0.7565  2.9794 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  School   (Intercept)  3.935   1.984   
##  Residual             36.148   6.012   
## Number of obs: 7185, groups:  School, 160
## 
## Fixed effects:
##             Estimate Std. Error t value
## (Intercept)  13.4660     0.1822   73.90
## SES           2.1276     0.1060   20.07
## MinorityYes  -2.9382     0.2070  -14.19
## 
## Correlation of Fixed Effects:
##             (Intr) SES   
## SES         -0.057       
## MinorityYes -0.311  0.193

5 Análise do summary do Modelo

No resumo do modelo, teremos acesso a três partes importantes:

  1. R Squared: Representa a variância explicada pelos Resíduos e pelo Intercepto.

  2. Fixed Effect: Refere-se aos coeficientes médios dos efeitos aleatórios. Caso algum preditor não tenha sido definido como efeito aleatório, o coeficiente correspondente é simplesmente a estimativa para \(y \sim x\), ou seja, slopes fixos para todos os observadores.

  3. Random Effects: Apresenta o campo std. dev, que indica a variação da média dos coeficientes dos efeitos fixos. No nosso caso, permitimos apenas que o intercepto varie Observamos que o intercepto do Fixed Effect é 13.46 (média), enquanto o intercepto do Random Effect possui um valor de 3.9. Isso implica que o intercepto varia entre \(13.46 \pm 3.9\). Já os slopes para SES e Minority são fixos para todos os observadores.

5.1 Verificação do \(R^2\)

Um passo importante é verificar o \(R^2\), que não é apresentado no resumo do modelo. Podemos calcular os coeficientes de determinação utilizando o pacote r-performance com a função r2(modelo).


library(performance)

Essa função retornará duas métricas:

  1. Conditional \(R^2\): Indica a proporção da variância total que é explicada pelo modelo, incluindo tanto os efeitos fixos quanto os efeitos aleatórios.

  2. Marginal \(R^2\): Refere-se à proporção da variância total que é explicada apenas pelos efeitos fixos do modelo. Neste caso, cerca de 12% da variância é explicada.

r2_performance <- r2(minority)
print(r2_performance)
## # R2 for Mixed Models
## 
##   Conditional R2: 0.210
##      Marginal R2: 0.124
#nosso modelo explica 9% a mais que um modelo apenas de efeitos fixos, sem o intercepto
#variar, representado pelo R2 Marginal

5.2 Conditional \(R^2\) (R² Condicional)

Neste modelo, o \(R^2\) condicional (como o valor de 0.210 que você obteve) reflete a variância explicada tanto pelos efeitos fixos (como SES e Minority) quanto pela variação entre as escolas (efeito aleatório do intercepto). Assim, aproximadamente 21% da variância em MathAch é explicada por todo o modelo, considerando as diferenças entre as escolas.

5.3 Marginal \(R^2\) (R² Marginal)

O \(R^2\) marginal (como o valor de 0.124) reflete apenas a variância explicada pelos efeitos fixos. Isso indica que cerca de 12% da variância é explicada apenas pelos preditores fixos, desconsiderando a variação do intercepto entre as escolas.

hist(resid(minority), freq = FALSE, breaks=50, col='steelblue')
lines(density(resid(minority)), col='tomato')

#um pouco de assimetria à esquerda
library(e1071)
skewness(resid(minority))
## [1] -0.1432044
kurtosis(resid(minority))
## [1] -0.5250705
#Criando modelo linear para comparação 
#library(performance)


linear = lm(MathAch ~ SES + Minority, data = math)
# Comparar os modelos
compare_performance(linear, minority)
## # Comparison of Model Performance Indices
## 
## Name     |   Model |   AIC (weights) |  AICc (weights) |   BIC (weights) |  RMSE | Sigma |    R2 | R2 (adj.) | R2 (cond.) | R2 (marg.) |   ICC
## ----------------------------------------------------------------------------------------------------------------------------------------------
## linear   |      lm | 46843.2 (<.001) | 46843.2 (<.001) | 46870.7 (<.001) | 6.298 | 6.300 | 0.161 |     0.161 |            |            |      
## minority | lmerMod | 46453.6 (>.999) | 46453.6 (>.999) | 46488.0 (>.999) | 5.956 | 6.012 |       |           |      0.210 |      0.124 | 0.098

6 Conclusões

6.1 Efeito Fixo vs. Efeito Aleatório

A comparação entre modelos com diferentes efeitos (slope fixo para SES versus SES como efeito aleatório) mostrou que o efeito aleatório para SES não acrescentou muita variância explicada. Isso sugere que a relação entre SES e o resultado (MathAch) é relativamente consistente entre as escolas, ou seja, o slope de SES é semelhante entre os grupos (escolas).

O efeito aleatório que melhorou ligeiramente o ajuste foi o do intercepto, indicando que há diferenças significativas entre as escolas no desempenho médio (MathAch), mas que a relação entre SES e MathAch não varia muito entre as escolas.

6.2 2. Interpretação Geral

  • SES tem um efeito significativo, porém consistente entre diferentes escolas. Sua relação com o desempenho acadêmico é robusta e não varia muito entre os grupos (escolas).
  • As escolas variam em termos de desempenho médio (MathAch), como esperado, o que justifica o uso de interceptos aleatórios.
  • O modelo com apenas intercepto aleatório foi suficiente para capturar a variação entre as escolas, e adicionar efeitos aleatórios ao slope de SES não trouxe melhorias substanciais.

6.3 3. Análise dos Resíduos

O modelo apresentou resíduos com média 0 e assimetria à esquerda.

6.3.1 4. Próximos Passos

  • Explorar a linearidade do modelo.
  • Avaliar a heterocedasticidade.
  • Verificar a aplicabilidade de outras variáveis preditoras ao modelo.