A estatística está em tudo: Generalized Estimating Equations e Mixed Models
- Fernanda Kelly
- há 19 horas
- 20 min de leitura
Vamos conversar sobre alguns métodos estatísticos MUITOOOO utilizados na área da economia, saúde, agronômica e também em estudos longitudinais?
Generalized Estimating Equations e Mixed Models são modelos que geralmente são empregados em situações em que devemos levar em consideração nas coletas e em casos longitudinais a estrutura de correlação entre indivíduos ou objetos de estudo, ou seja, sob brandas suposições de dependência. No universo de Machine Learning, esse tipo de análise também é conhecida como Modelagens Multiníveis ou Hierárquicos e, essas nomenclaturas trazem exatamente a ideia de se tratar de dados agrupados e estruturas hierárquicas complexas ou não.
Mas que tipo de hierarquia/multinível estamos falando?
Como esse assunto é um pouco mais complexo, para responder os questionamentos desse post, vou utilizar um exemplo da área agronômica, a qual fiz um trabalho na graduação em 2017.
Diggle, Liang, e Zeger (1994), apresentaram uma base de dados que apresenta a coleta de dados relacionados a proteína contida no leite de vacas. A ideia é caracterizar a mudança no nível dessas proteínas e fatores que influenciam essas mudanças, como por exemplo, a dieta a qual as vacas estão submetidas. Estes dados fazem parte da pasta do pacote nlme do software R, descrevendo o nível de proteína presente no leite de 79 vacas australianas. A coleta de dados foi efetuada após o nascimento das vacas e, em seguida, administrou-se três tipos de dietas para três grupos diferentes de vacas. Repetiu-se a coleta de observações do nível de proteína durante 19 semanas.
Lendo o experimento acima, já é possível identificar níveis hierárquicos de possíveis causas do aumento ou a diminuição de proteína no leite dessas vacas, por exemplo, a dieta. Há três níveis de dieta, três grupos de vacas e a coleta é ao longo do tempo. E claro, considerando que os dados são experimentais, quando as medidas são tomadas na mesma unidade experimental ou parcela, que no nosso caso são os grupos de vacas, diz-se que tem-se dados provenientes de um planejamento com medidas repetidas. Quando essas medidas repetidas são tomadas ao longo (ou de profundidades, ou distâncias), diz-se que essas medidas repetidas constituem dados longitudinais (LIMA, 1996).
Conseguiu pegar as características de modelagem?
Podemos dizer que essa modelagem é um tanto sofisticada e, talvez, um pouco complicada, mas é importante você conseguir entender que esse tipo de modelagem possui características marcantes e que não devem ser esquecidas, como as medidas repetidas ao longo do tempo.
Conforme Eisenhart (1947) apud Rocha (2015) os modelos de regressão podem ser classificados em três categorias:
Modelos de efeitos fixos: são modelos em que todos os componentes, exceto o erro experimental, são considerados fixos;
Modelos de efeitos aleatórios: são modelos em que todos os componentes, exceto a média geral, são considerados aleatórios;
Modelos de efeitos aleatórios e fixos: são formados pela combinação dos dois modelos anteriores, isto é, são constituídos por termos de efeitos fixos e aleatórios, além da média geral e do erro experimental, respectivamente. São também chamados de modelos mistos.
Pera aí Fer, modelos mistos? Nera longitudinal?
Então...
Se você é desse universo, talvez seja um pouco mais fácil, mas para quem não é, faz sentido falar sobre correlação entre as medidas tomadas na mesma unidade experimental.
Senhorrrrrrrrrrrrrrrrrrrrrrrr...
Pois é, na teoria é mais difícil mesmo. E, para exemplificar essa correlação, vou trazer o que vivenciamos em 2020, a pandemia da COVID-19. A COVID-19 é uma infecção respiratória, certo? E ela tem suas características de expansão, seja pelo espirro ou contato direto ou indireto com a pessoa infectada.
Vamos pensar o hospital como uma unidade experimental?
Nesse hospital haviam três tipos de Unidade de Terapia Intensiva. Em cada unidade tínhamos perfis de gravidade diferente entre os pacientes e, ainda tínhamos àqueles que já estavam internados no hospital antes da COVID-19 e não estavam infectados. Agora vamos pensar que a infecção respiratória (COVID-19) é a correlação. Se a gente juntar uma UTI em que os pacientes estivessem com sintomas leves da infecção com os não infectados, provavelmente essa infecção iria se expandir aos não infectados, pois naquela ala de UTI os pacientes já estavam auto correlacionados pela infecção. Ou seja, a ala de UTI possui sua própria correlação e os seus pacientes também, e ainda assim, uma ala de UTI pode ter pacientes mais graves do que outra ala.
Mas, como o hospital fica nessa loucurada de correlação?
O hospital é uma unidade experimental em que os pacientes internados em suas dependências, independente de ala, se tornam correlacionados. E agora vem a dúvida...
Fê, é muita correlação, como modelar isso?
É aí que entra os Modelos Lineares Mistos, pois estes modelos são M-U-I-T-O úteis para a análise de dados longitudinais, visto que a parte aleatória do modelo e a correlação da unidade experimental devem ser acomodadas corretamente na modelagem.
E de qual modelagem estamos falando?
A ideia central desse tipo de modelagem é separar o efeito médio da população e a variação individual. Esses modelos são amplamente usados porque conseguem modelar trajetórias individuais ao longo do tempo.
Os pacotes mais utilizados para ajustar modelos nesse cenário são os descritos abaixo.
library(nlme)
library(lme4)Vamos de exemplo?
Com os pacotes devidamente carregados, vamos iniciar com a percepção dos dados.
dados_milk <- nlme::MilkNa base de dados Milk, há 4 variáveis, sendo elas:
protein: Vetor contínuo com observações da proteína contida no leite;
Time: Identificador das semanas de 1 a19;
Cow: Fator ordenado dado para identificar cada vaca analisada, sendo ele individual para as 79 vacas;
Diet: Fator com níveis identificando a dieta para cada vaca. As dietas são: barley, lupins e barley+lupins.
É importante ressaltar que não foram todas as vacas que completaram o ciclo de 19 semanas observadas, mas isso não significa que teremos valores ausentes.
Fê, tá faltando uma análise exploratória, né?
Sabe que eu também acho? Vamos de gráfico.
ggplot(dados_milk, aes(x = factor(Time),
y = protein,
fill = Diet)) +
geom_boxplot() +
labs(
x = "Tempo",
y = "Nível de proteína",
fill = "Dieta"
) +
scale_fill_manual(values = c("lightgray","lightblue","lightgreen")) + theme_minimal()
ggplot(dados_milk, aes(x = Diet,
y = protein,
fill = Diet)) +
geom_boxplot() +
labs(
x = "Tempo",
y = "Nível de proteína",
fill = "Dieta"
) +
scale_fill_manual(values = c("lightgray","lightblue","lightgreen")) + theme_minimal()
O que podemos interpretar com esses dois gráficos?
A parte boa do gráfico Box-Plot é a fácil visualização das médias de cada uma das dietas ao longo do tempo ou individualmente. Ao longo das 19 semanas de observação, a dieta baseada em barley apresentou os maiores níveis médios de proteína e maior estabilidade temporal. A dieta barley+lupins apresentou valores intermediários, enquanto a dieta composta apenas por lupins apresentou os menores níveis de proteína e uma tendência de redução ao longo do tempo.
Abaixo é possível observar o caminho percorrido por cada uma das 79 vacas alocadas nas dietas Barley, Barley+Lupins e Lupins.
dadosGD <- nlme::groupedData(protein ~ Time|Cow,
outer = ~Diet,
data = dados_milk)
plot(dadosGD,
outer = TRUE)
Mas não podemos inferir apenas com a visualização...
... e é agora que a brincadeira começa a ficar boa.
Essa modelagem nos dá a opção de trabalhar com efeitos fixos, aleatórios ou ambos na mesma modelagem. Expliquei lá encima o que era cada um, porém não expliquei como aplicá-los e quando faz sentido.
O Efeito Fixo (EF) é um parâmetro que representa um valor ou característica específica da população que queremos estimar. Ou seja, ele é constante, desconhecido, mas não é considerado uma variável aleatória, por exemplo: gênero, tratamento médico, dieta nutritiva. Já o Efeito Aleatório (EA) representa variação entre unidades, assumida como proveniente de uma distribuição probabilística. Ou seja, cada unidade tem um valor próprio e esses valores vêm de uma distribuição, como citado em exemplos, escolas, hospitais, pacientes/indivíduos. Porém tudo depende da sua questão de pesquisa e, claro,
por que isso importa em dados longitudinais?
Os hospitais podem ser consideradas de efeito fixo ou aleatório. Se é aleatório, não queremos estimar cada hospital especificamente, queremos estimar a variância entre hospitais.
Entendeu?
Se não, ao longo da modelagem essa ideia pode ficar um pouco mais palpável e lhe ajudar a entender melhor esse conceito de EF e EA, mas é de extrema importância ter entendimento de negócio ou do experimento que estamos avaliando.
Para as modelagens longitudinais, utilizaremos a função lmer do pacote lme4. Esse pacote é um dos clássicos da modelagem. Sua expressão é a seguinte:
resp ~ EFexpr + (EAexpr1 | factor1) + (EAexpr2 | factor2) + ...É igual a uma regressão, mas com um plus (rs). O EA é identificado quando sua expressão está entre parênteses "()" é separado por "|" (barra vertical). Uma maneira de pensar sobre o operador de barra vertical é como um tipo especial de interação entre a matriz do modelo e o fator de agrupamento. Essa interação garante que as colunas da matriz do modelo tenham efeitos diferentes para cada nível do fator de agrupamento.
Mas e a expressão (1 | factor)?
Essa é a expressão mais simples possível, onde cada nível do grupo factor é o seu próprio intercepto aleatório. Abaixo é possível observar alguns exemplos de expressões. É importante ressaltar que o pacote lme4 assume que todos os coeficientes associados com termos de EA são correlacionados, caso seja de interesse especificar que estes não são correlacionados, usa-se "||".
No nosso exemplo, o primeiro modelo a ser testado é:
modA <- lme4::lmer(formula = protein ~ Time*Diet + (1| Cow),
data = dados_milk,
REML = FALSE)O que são esses parâmetros, Fê?
No modelo A (modA), estamos considerando que a nossa variável resposta é o nível de proteína. Em relação as variáveis explicativas, como Efeito Fixo temos a interação entre Time e Diet e Efeito Aleatório a Cow como identidade. Esse modelo analisa como o teor de proteína do leite varia ao longo do tempo e de acordo com a dieta das vacas, considerando que cada vaca é uma unidade repetida. Já os parâmetros da função,
REML: É um escalar lógico, em que as estimativas devem ser escolhidas para otimizar o critério REML (em vez da log-verossimilhança). Ele recebe TRUE ou FALSE.
O REML é a abreviação de Restricted Maximum Likelihood (REML), que é um método para estimar variâncias em modelos mistos que remove o efeito dos parâmetros fixos antes de estimar as variâncias.
O REML estima as variâncias com base apenas na parte “aleatória” dos dados.
E quando decidimos não utilizar a REML, consequentemente, estamos utilizando a Maximum Likelihood (ML), que é, geralmente, utilizada em regressões e métodos lineares.
Com o intuito de comparar os métodos, abaixo se encontra o resultado do modelo A.
anova(modA)Analysis of Variance Table npar Sum Sq Mean Sq F value Time 1 1.30825 1.30825 17.5882 Diet 2 1.31209 0.65605 8.8200 Time:Diet 2 0.51855 0.25928 3.4857 |
modB <- lme4::lmer(formula = protein ~ Time*Diet + (1| Cow),
data = dados_milk,
REML = TRUE)
anova(modB)Analysis of Variance Table npar Sum Sq Mean Sq F value Time 1 1.30804 1.30804 17.5441 Diet 2 1.26337 0.63168 8.4725 Time:Diet 2 0.51576 0.25788 3.4588 |
Viu que não tem muita diferença nesse experimento?
Os efeitos de tempo e dieta foram significativos sob ambos os métodos de estimação (ML e REML), indicando robustez dos resultados à forma de estimação dos componentes de variância. Mas, aqui tem um pulo do gato, pois quando a ideia é comparar os EF da modelagem, utilizamos REML igual FALSE e, para aleatórios, usamos igual a TRUE. E, apesar dos resultados serem próximos, eles não são iguais.
Para você entender melhor essa comparação, vamos ajustar um novo modelo com EF diferente do modA e modB com REML igual a FALSE. O modelo abaixo possui a mesma estrutura fixa, mas a aleatória passa a ser (Time | Cow), ou seja, intercepto aleatório E inclinação aleatória para o tempo.
cada vaca tem seu próprio nível inicial de proteína;
cada vaca evolui ao longo do tempo de forma diferente.
Explicando de uma forma mais tranquila é dizendo que essa sintaxe equivale a:
(1 + Time | Cow)
E isso demonstra que na modelagem há três componentes aleatórios (EA):
Variância do intercepto: diferença inicial entre vacas.
Variância da inclinação: diferença na resposta ao tempo.
Covariância intercepto–inclinação:
Indica se vacas com proteína inicial alta:
aumentam mais rápido?
mais devagar?
não há relação?
modC <- lme4::lmer(formula = protein ~ Time*Diet + (Time | Cow),
data = dados_milk,
REML = FALSE)
anova(modC)Analysis of Variance Table npar Sum Sq Mean Sq F value Time 1 0.96010 0.96010 15.887 Diet 2 0.99617 0.49808 8.242 Time:Diet 2 0.09378 0.04689 0.776 |
O que podemos fazer agora é a comparação do modelo A e C através do AIC e BIC, vamos lá?
anova(modA, modC)Data: dados_milk Models: modA: protein ~ Time * Diet + (1 | Cow) modC: protein ~ Time * Diet + (Time | Cow) npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq) modA 8 489.15 530.74 -236.58 473.15 modC 10 360.58 412.56 -170.29 340.58 132.57 2 < 2.2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 |
De acordo com o Critério de Informação de Akaike (AIC) que é uma métrica estatística usada para selecionar o melhor modelo entre vários candidatos, equilibrando a qualidade do ajuste e a complexidade, o modelo escolhido é o modC, pois este possui os menores valores de AIC e BIC e é o estatisticamente significativo (Pr(>Chisq)).
Escolhemos a modelagem aleatória, certo?
Agora falta a fixa.
modC1 <- lme4::lmer(formula = protein ~ Time + (Time| Cow),
data = dados_milk,
REML = F)
modC2 <- lme4::lmer(formula = protein ~ Time + Diet + (Time| Cow),
data = dados_milk,
REML = F)
modC3 <- lme4::lmer(formula = protein ~ Time * Diet + (Time| Cow),
data = dados_milk,
REML = F)Comparando os modelos,
anova(modC1, modC2, modC3)Data: dados_milk Models: modC1: protein ~ Time + (Time | Cow) modC2: protein ~ Time + Diet + (Time | Cow) modC3: protein ~ Time * Diet + (Time | Cow) npar AIC BIC logLik -2*log(L) Chisq Df Pr(>Chisq) modC1 6 368.72 399.91 -178.36 356.72 modC2 8 358.12 399.70 -171.06 342.12 14.6082 2 0.0006728 *** modC3 10 360.58 412.56 -170.29 340.58 1.5343 2 0.4643393 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 |
Assim escolhemos o modC2, pois é estatisticamente significante, tem os menores valores de AIC e BIC. O modC3 não foi escolhido, mas é interessante reparar que não há interação entre Time e Diet, indicando que o efeito do tempo é assumido igual para todas as dietas.
Agora vamos analisar o nosso modelo?
É um tabelão, mas essa tabela traz muita informação e é incrível para a análise ser completa e robusta.
summary(modC2)Linear mixed model fit by maximum likelihood ['lmerMod'] Formula: protein ~ Time + Diet + (Time | Cow) Data: dados_milk AIC BIC logLik -2*log(L) df.resid 358.1 399.7 -171.1 342.1 1329 Scaled residuals: Min 1Q Median 3Q Max -3.1431 -0.5930 -0.0248 0.5467 3.6001 Random effects: Groups Name Variance Std.Dev. Corr Cow (Intercept) 0.0700823 0.26473 Time 0.0006278 0.02506 -0.78 Residual 0.0604166 0.24580 Number of obs: 1337, groups: Cow, 79 Fixed effects: Estimate Std. Error t value (Intercept) 3.616574 0.043924 82.337 Time -0.012491 0.003153 -3.961 Dietbarley+lupins -0.094164 0.048706 -1.933 Dietlupins -0.197188 0.048729 -4.047 Correlation of Fixed Effects: (Intr) Time Dtbrl+ Time -0.601 Dtbrly+lpns -0.576 0.001 Dietlupins -0.576 0.001 0.519 optimizer (nloptwrap) convergence code: 0 (OK). |
Fer, mas não tem pressupostos?
Há pressuposições (e muitas), visto que temos Efeitos Fixos e Aleatórios na mesma modelagem, mas veja os resíduos em Scaled residuals, a distribuição me parece ser aproximadamente simétrica, pois há poucos outliers extremos e o ajuste parece razoável.
O que não podemos esquecer é que a variância/covariância e correlações são de E-X-T-R-E-M-A importância nas modelagens que fazemos, visto que na maioria dos casos a ideia é conseguir modelar a variação entre variáveis/indivíduos de interesse ao longo do tempo ou não.A saída do modelo ajustado já nos conta muita coisa, por exemplo, no caso dos EF:
Intercepto da vaca: Há diferenças relevantes no nível inicial de proteína entre vacas.
Variância = 0.0701
Desvio padrão ≈ 0.265
Inclinação do tempo por vaca: As vacas diferem na taxa de mudança ao longo do tempo, mas é pouco, trazendo a questão que existe heterogeneidade dinâmica, porém pequena.
Variância = 0.000628
Desvio padrão ≈ 0.025
Correlação intercepto-inclinação: Esse é um resultado MUITO interessante, pois indica que vacas com proteína inicial alta tendem a diminuir mais rápido ao longo do tempo ou o inverso, vacas com nível baixo tendem a cair menos (ou até melhorar)
Corr = −0.78
Variância residual: Representa variação dentro da mesma vaca ao longo do tempo.
Variância = 0.0604
Desvio padrão ≈ 0.246
Mas, ainda temos que verificar outras coisinhas, tá?
Geralmente, não tem como correr da linearidade e especificação do modelo, normalidade dos resíduos e, nesse caso, normalidade dos efeitos aleatórios. A homoscedasticidade residual nem preciso falar, né? E claro, nós temos indivíduos (vacas) ao longo do tempo, identificar a autocorrelação residual é crucial. Só que você não precisa se desesperar, muitos desses diagnósticos são visuais.
Podemos iniciar com a correlação intraclasse, que no nosso caso temos o intercepto aleatório, slope aleatório e a correlação entre eles.
vc <- as.data.frame(nlme::VarCorr(modC))
var_intercept <- vc$vcov[vc$grp == "Cow" & vc$var1 == "(Intercept)"]
var_resid <- vc$vcov[vc$grp == "Residual"]
ICC <- var_intercept / (var_intercept + var_resid)Resultado:
0.53360979 -0.09175861O valor 0.53360979 é o ICC associado ao intercepto (no tempo zero ou médio), que traz para nós a interpretação de que cerca de 53% da variabilidade da proteína é explicada pelas diferenças entre vacas. Essa porcentagem é alta, implicando que as vacas são muito diferentes entre si e há forte dependência intra-vaca, trazendo, novamente, que o uso da Modelagem Mista é absolutamente necessária. Já o segundo valor -0.09175861 reflete o efeito do slope aleatório mais a correlação.
E o que isso indica?
Esse resultado pode indicar:
mudança da correlação ao longo do tempo
dependência diminuindo com o tempo
estimativa instável
correlação intercepto–slope forte
Isso é um tanto RUIM e, para tentar entender um pouco mais, vamos estudar a dependência temporal, ou seja, a autocorrelação residual.
acf(resid(modC))
Diante desse gráfico, podemos observar que há vacas com autocorrelação positiva alta e, mesmo que o modelo misto explique parte da dependência intra-vaca, ainda assim sobra dependência temporal para algumas vacas.
Você achou mesmo que trabalhar com Modelos Mistos ia te livrar de aprender um pouco sobre Séries Temporais?
Pois é, não foi dessa vez. O nosso modelo demonstra que há dependência temporal residual, mesmo que Modelos Mistos captem diferenças entre vacas, ele pode não capturar totalmente a dinâmica temporal.
E agora, Fê?
Vamos modelar com correlação residual. Veja a fórmula abaixo.
modC1 <- nlme::lme(protein ~ Time + Diet,
random = ~ Time | Cow,
correlation = corAR1(form = ~ Time | Cow),
data = dados_milk)Agora o modelo assume AR(1) dentro de cada vaca. O resultado se encontra a seguir.
summary(modC1)Random effects: Formula: ~Time | Cow Structure: General positive-definite, Log-Cholesky parametrization StdDev Corr (Intercept) 0.18420039 (Intr) Time 0.01685216 -0.721 Residual 0.29470326 Correlation Structure: ARMA(1,0) Formula: ~Time | Cow Parameter estimate(s): Phi1 0.5803862 Fixed effects: protein ~ Time + Diet Correlation: (Intr) Time Dtbrl+ Time -0.596 Dietbarley+lupins -0.579 0.002 Dietlupins -0.580 0.004 0.518 Standardized Within-Group Residuals: Min Q1 Med Q3 Max -2.84166066 -0.67057009 -0.08829411 0.52443728 3.12557000 Number of Observations: 1337 Number of Groups: 79 AIC BIC logLik 142.3767 189.1333 -62.18833 |
Value Std.Error DF t-value p-value (Intercept) 3.699696 0.04545448 1257 81.39343 0.0000 Time -0.018369 0.00315239 1257 -5.82713 0.0000 Dietbarley+lupins -0.095369 0.05067196 76 -1.88209 0.0637 Dietlupins -0.200787 0.05071025 76 -3.95949 0.0002 |
AIC BIC logLik 142.3767 189.1333 -62.18833 |
E a matriz de variâncias do modelo?
Podemos observar que as vacas começam com níveis moderadamente diferentes de proteína (0.03393 (SD = 0.1842), mas as vacas evoluem de forma relativamente parecida (0.000284 (SD = 0.01685)), ou seja, a principal diferença entre vacas está no nível inicial, não na dinâmica temporal.
Fê, e o AR(1)? Foi pra onde?
A dependência temporal está separa, pois ela não entra na variância residual e sim na estrutura da correlação. O Phi igual a 0.58, descreve uma dependência moderada-alta, ou seja, se o resíduo foi alto hoje, tende a continuar alto na próxima medição.
Existe persistência temporal dentro de cada vaca.
Trazendo o coeficiente de correlação intraclasse,
nlme::VarCorr(modC1)ICC = variância entre vacas ( 0.0339) / variância entre vacas ( 0.08685) + variância dentro da vaca (0.0339) = 0.28
Note que ele é menor comparado ao melo sem o ajuste de AR1, pois agora o AR1 captura parte da dependência.
Voltando alguns passos...
O que podemos responder com esse modelo?
Como o teor de proteína do leite muda ao longo do tempo
Se a dieta influencia esse teor
Considerando que cada vaca é diferente
E que medições próximas no tempo são parecidas
Ou seja,
“Como tempo e dieta afetam a proteína do leite, levando em conta as diferenças entre vacas e a evolução ao longo do tempo.”
Se este estudo fosse uma história, cada vaca seria uma personagem com personalidade própria, algumas começam a lactação produzindo leite riquíssimo em proteína, enquanto outras são mais modestas desde o início. Em média, porém, a “vaca típica” inicia com teor proteico estimado em 3,6997, nosso ponto de partida da narrativa. À medida que o tempo passa, ocorre um pequeno declínio natural na proteína (β = −0,01837), como se a lactação fosse lentamente “consumindo energia” ao longo dos capítulos.
A dieta também entra como coadjuvante importante: em comparação com a dieta de referência (barley), a combinação barley + lupins tende a reduzir levemente a proteína (β = −0,09537; p = 0,0637), um efeito tímido, quase sussurrado pela estatística, enquanto a dieta composta apenas por lupins provoca uma queda mais evidente e consistente (β = −0,20079; p = 0,0002), como um plot twist nutricional nada favorável ao teor proteico.
Mas o grande charme do modelo está em reconhecer que as vacas não seguem um roteiro único. Há diferenças claras no nível inicial entre elas (desvio-padrão do intercepto = 0,1842), mostrando que algumas já entram na história com “superpoderes proteicos”, enquanto outras começam em modo mais econômico. Também existem pequenas diferenças na velocidade com que cada vaca muda ao longo do tempo (desvio-padrão da inclinação = 0,01685), embora a maioria caminhe em ritmos semelhantes. Um detalhe particularmente intrigante é a forte correlação negativa entre nível inicial e taxa de mudança (−0,721), ou seja, as vacas que começam no topo tendem a descer mais rápido, quase como atletas que largam em disparada, mas não sustentam o sprint até o final. Dentro de cada vaca, ainda há variações naturais de uma medição para outra (desvio-padrão residual = 0,2947), lembrando que sistemas biológicos nunca são perfeitamente previsíveis.
Além disso, o tempo tem memória: a estrutura AR(1) estimou φ₁ = 0,5804, indicando que medições próximas são moderadamente parecidas, indicando que se a proteína está alta hoje, há boa chance de continuar relativamente alta amanhã, como se cada vaca carregasse consigo um “estado recente” que influencia o próximo momento. O estudo acompanhou 1.337 observações de 79 vacas, e os resíduos permaneceram dentro de limites aceitáveis (aproximadamente entre −2,84 e 3,13), sugerindo que o modelo capturou bem a trama principal sem grandes incoerências.
No fim das contas, a história que emerge é bastante humana, ou melhor, bovina: a proteína do leite tende a diminuir com o tempo, dietas com lupins podem acelerar essa queda, cada vaca possui sua própria identidade produtiva e o passado recente ajuda a prever o presente. Em outras palavras, não existe “uma vaca média” perfeita, mas sim um elenco diverso, cada uma com seu próprio arco narrativo ao longo da lactação .
Fê, mas não respondeu muita coisa, né? Por exemplo, qual a melhor dieta?
Então...
Nesse caso nós modelamos as vaquinhas, observando que cada uma possui sua própria curva em relação a dieta e o Modelo Misto não vai nos responder com precisão esse questionamento e, por isso, esse post é sobre GEE também.
No GEE, podemos responder a seguinte pergunta:
Em média, a dieta aumenta a proteína das vaquinhas?
E claro, essa decisão é dependente da estrutura dos nossos dados, o que até aqui já foi B-A-S-T-A-N-T-E discutido. No GEE, dados correlacionados são modelados utilizando o mesmo preditor linear e função de ligação que seriam utilizados no caso independente. Entretanto, a estrutura de covariância das medidas repetidas também é modelada. As associações intra-indivíduos, são consideradas através das matrizes auxiliares de correlação intra-indivíduo, que incorporam a estrutura de dependência diretamente na estimação dos betas, essas matrizes auxiliares de correlação são, comumente (ou estatisticamente), denominadas por estruturas de correlações de trabalho.
Ai Fê, de novo essa história de correlação?
Sim! Dessa história não tem como correr se a estrutura dos dados é essa: indivíduo ao longo do tempo. Se em Modelos Mistos nós tínhamos a opção de modelagem com REML para estimar os componentes de variância (variâncias dos efeitos aleatórios, variância residual, estrutura de covariância implícita, sem precisar "chutar" uma correlação), já na GEE utilizaremos correlações de trabalho.
Há sete formas de se trabalhar com as correlações de trabalho quando se trata do pacote gee.
install.packages("gee")
install.packages("geepack")Sua estrutura em uma modelagem segue abaixo.
gee(formula, id, data, subset, na.action, R = NULL, b = NULL, tol = 0.001, maxiter = 25, family = gaussian, corstr = "independence", Mv = 1, silent = TRUE, contrasts = NULL, scale.fix = FALSE, scale.value = 1, v4.4compat = FALSE) É notório que há uma quantidade enorme de parâmetros que devem ser declarados nessa modelagem, mas, evidenciando a correlação de trabalho, o nosso parâmetro de interesse é o corstr que se baseia em uma cadeia de caracteres especificando a estrutura de correlação. Veja que o parâmetro corstr já se encontra com a declaração "independence" e essa será a primeira estrutura de correlação de trabalho que vamos estudar.
A estrutura de correlação de trabalho denominada por independence é uma estrutura que utiliza uma matriz identidade de ordem t.
O que essa declaração quer dizer?
O independence significa que, a partir dessa estrutura de matriz identidade, não estamos permitindo que, por exemplo, indivíduos em um hospital não sejam correlacionados devido a mesma exposição de tratamento. Ou seja, estamos modelando o modelo S-E-M a prerrogativa de correlação entre os indivíduos.
Fê, é contraditório?
É sim, mas existe essa declaração.
A segunda estrutura de correlação de trabalho é a exchangeable. A estrutura exchangeable atua com a prerrogativa de que as correlações são iguais independentemente do intervalo de tempo. Sua declaração no R é por "exchangeable" ou simplesmente "exc".
A estrutura stationary ou m-dependent structure é declarada no R por "stat_M_dep" e, quando se tem ideia do método, a nomeação dela é bem intuitiva. Essa estrutura assume que as correlações t medidas separadas são iguais, as correlações +1 medidas separadas são consideradas iguais, e assim por diante para t = 1 até t = m. As correlações com mais de m medições de distância são assumidas como zero. Quando, por exemplo, uma 'estrutura de correlação dependente de 2' é assumida, todas as correlações separadas por uma medição são consideradas as mesmas, todas as correlações separadas por duas medições são consideradas iguais, e as correlações com mais de duas medições separadas são consideradas zero. Em contrapartida, há a estrutura denominada por non stationary e sua declaração no R é dada por "non_stat_M_dep".
Uma quarta possibilidade é uma estrutura de correlação denominada por autoregressive, ou seja, as correlações separadas por uma medição são consideradas ρ; correlações separadas por duas medições são assumidas como ρ2; as correlações t medições separadas são consideradas ρt. Sua declaração no R é dada por "AR-M" e para a declaração do "M" deve-se utilizar o parâmetro "Mv" com a especificação.
Há também a estrutura de correlação denominada por "unstructured". Sua abordagem é considerar que todas as correlações são diferentes. É a opção mais flexível, porém seu número de parâmetros cresce rapidamente e a estimação pode se tornar instável. Portanto, é de extrema importância perceber qual estrutura de correlação é mais apropriada para a análise de interesse e dados de trabalho.
Fer, como escolher qual estrutura de correlação usar?
A boa notícia é que as estimativas GEE são válidas mesmo se você especificar incorretamente a estrutura de correlação.
É claro que isso pressupõe que o modelo esteja correto, mas, novamente, nenhum modelo está exatamente correto. Verifique SEMPRE como as estimativas de coeficiente e os erros padrão mudam com outras estruturas de correlação e, se as alterações forem mínimas, use a estrutura de correlação mais simples.
O GEE NÃO tenta modelar a variabilidade entre unidades.
Ele só quer:
Estimar o efeito médio populacional corretamente, mesmo que haja correlação nos dados.
Vamos testar bem rapidinho?
gee_exch <- geeglm(protein ~ Time + Diet,
id = Cow,
data = dados_milk,
corstr = "exchangeable")summary(gee_exch)Call: geeglm(formula = protein ~ Time + Diet, data = dados_milk, id = Cow, corstr = "exchangeable") Coefficients: Estimate Std.err Wald Pr(>|W|) (Intercept) 3.582748 0.042551 7089.408 < 2e-16 *** Time -0.006190 0.002737 5.116 0.023711 * Dietbarley+lupins -0.096383 0.046500 4.296 0.038195 * Dietlupins -0.204533 0.054840 13.910 0.000192 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Correlation structure = exchangeable Estimated Scale Parameters: Link = identity Estimated Correlation Parameters: Number of clusters: 79 Maximum cluster size: 19 |
O modelo GEE ajustado com estrutura de correlação exchangeable indica que o teor de proteína do leite diminui ligeiramente ao longo do tempo e é influenciado pela dieta, considerando a dependência entre medições da mesma vaca.
A proteína média inicial estimada para a dieta de referência foi de 3,5827, e observou-se uma pequena redução ao longo do tempo (β = −0,00619; p = 0,0237). Em comparação com a dieta base (barley), a dieta barley + lupins apresentou redução moderada na proteína (β = −0,09638; p = 0,0382), enquanto a dieta composta apenas por lupins mostrou queda mais acentuada e altamente significativa (β = −0,20453; p = 0,000192), indicando efeito negativo consistente sobre o teor proteico.
A estrutura exchangeable assume correlação constante entre todas as medições dentro de cada vaca, refletindo dependência intra-animal, mas sem modelar diferenças individuais específicas.
Tá, mas qual dieta usar?
A dieta de referência (baseline) é a que produz maior teor de proteína e, claro, deve-se evitar a dieta que possui apenas lupins se o objetivo é maximizar proteína.
Fê, e se o problema for custo?
Então a mistura barley + lupins pode ser aceitável nesse caso. O nosso ranking é o seguinte:
🥇 Dieta base > 🥈 mistura > 🥉 lupins
É muita informação, né? Eu também acho e fico bem zonza, mas essas modelagens são mais robustas e precisam de uma atenção maior no quesito tempo/temporabilidade dos dados e isso dá uma dificultada no entendimento do que usar na hora do aperto.
Espero que você tenha gosto e que aprofunde ainda mais nessas modelagens.
Até o próximo post.
Fernanda Kelly R. Silva | Estatística
Anexos
library(lme4)
mod_lmer <- lmer(protein ~ Time + Diet + (Time | Cow),
data = dados_milk)
library(DHARMa)
sim <- simulateResiduals(mod_lmer)
plot(sim)
# Mas aqui perdemos a estrutura AR(1).
Referências
BATES, Douglas et al. Fitting linear mixed-effects models using lme4. Journal of statistical software, v. 67, p. 1-48, 2015.
Diggle, P.J., Liang, K.-Y. and Zeger, S.L. (1994) Analysis of Longitudinal Data. Clarendon Press, Oxford.
LIMA, C.G. Análise de dados longitudinais provenientes de experimentos em blocos casualizados. 1996. 126p. Tese (Doutorado em Estatística e Experimentação Agronômica) - Escola Superior de Agricultura “Luiz de Queiroz", Universidade de São Paulo, Piracicaba, 1996.
Comentários