Medidas Repetidas e Modelos Mistos

Porque várias linhas da mesma unidade não são várias unidades independentes, e como a ANOVA e os modelos mistos representam essa dependência

Autores

Afiliações

João O. Santos

ISPA

Cristina Mendonça

ISPA; WJCR

Vitória Melita

FP-UL

Mariona Pascual Peñas

UIB

João Raposo

ISPA

Marta Barros

FP-UL

0 A Pergunta Antes do Modelo

Suponham que medimos o nível de stress das mesmas pessoas de manhã, à tarde e à noite. A pergunta de investigação é simples:

O nível médio de stress muda ao longo do dia?

Se tivermos 35 pessoas e três medições por pessoa, a base de dados tem 105 linhas. Mas não temos 105 pessoas independentes. Temos 35 unidades, cada uma observada três vezes.

Esta é a primeira ideia a fixar:

o número de linhas não nos diz, por si só, quantas unidades independentes temos.

O mesmo problema aparece quando seguimos o mesmo peixe ao longo do tempo, medimos várias vezes a mesma colónia, colocamos vários organismos no mesmo tanque ou observamos vários quadrados dentro do mesmo recife.

Antes de pensar em ANOVA, modelos mistos ou sintaxe de R, perguntem:

1. o que representa uma linha?
2. que linhas pertencem à mesma unidade?
3. em que nível varia o preditor que nos interessa?

Essas três perguntas resolvem uma parte surpreendentemente grande daquilo que costuma parecer assustador neste capítulo.

Para a implementação prática, usamos afex::aov_4(). Aqui vamos primeiro perceber o que o modelo precisa de fazer.

Ver as Ligações Entre as Linhas

O gráfico seguinte liga as três observações da mesma pessoa. Não olhem primeiro para a média geral. Olhem para as linhas individuais.

library(ggplot2)

ggplot(ds, aes(x = momento, y = resposta, group = unidade)) +
  geom_line(alpha = 0.25) +
  geom_point(alpha = 0.5) +
  theme_classic() +
  labs(x = "Momento", y = "Resposta")

Algumas pessoas podem ter valores sempre mais altos; outras, sempre mais baixos. Isso é variação entre pessoas. A pergunta sobre a hora do dia, porém, é uma pergunta sobre mudança dentro da mesma pessoa.

Se misturarmos essas duas fontes de variação, estamos a usar o erro errado para a pergunta que queremos fazer.

1 ANOVA de Medidas Repetidas

Dados = Modelo + Erro, Outra Vez

No capítulo anterior usámos

\[ \text{Dados}=\text{Modelo}+\text{Erro}. \]

A ideia não muda. O que muda é que agora sabemos uma coisa sobre o erro: observações da mesma unidade partilham história e características. Portanto, a parte que fica por explicar tem estrutura.

O Modelo Ingénuo

Podíamos ajustar:

m_ingenuo <- lm(resposta ~ momento, data = ds)

Este modelo sabe que existem manhã, tarde e noite. Não sabe que três linhas podem pertencer à mesma pessoa.

Assim, diferenças estáveis entre pessoas ficam misturadas com a variação que queremos usar para avaliar a mudança ao longo do dia.

A palavra «ingénuo» não quer dizer que os dados estejam mal recolhidos. Quer dizer apenas que o modelo esqueceu parte do delineamento.

Pôr a Unidade no Modelo

Uma forma muito concreta de ver a solução clássica é dar a cada pessoa o seu próprio nível basal.

Primeiro ajustamos um modelo que conhece apenas as pessoas:

m_unidade <- lm(resposta ~ unidade, data = ds)

Depois libertamos a possibilidade de existir uma diferença média entre momentos:

m_unidade_momento <- lm(resposta ~ unidade + momento, data = ds)

Agora comparem as duas perguntas:

• modelo compacto: as pessoas podem ter médias diferentes, mas não existe uma mudança média entre momentos;
• modelo aumentado: as pessoas podem ter médias diferentes e permitimos também uma mudança média entre momentos.

A comparação

anova(m_unidade, m_unidade_momento)

Res.Df	RSS	Df	Sum of Sq	F	Pr(>F)
70	20286.35	NA	NA	NA	NA
68	13711.37	2	6574.98	16.30393	1.6e-06

pergunta exactamente aquilo que queríamos:

Depois de retirarmos as diferenças médias entre pessoas, ainda ganhamos explicação ao permitir diferenças entre manhã, tarde e noite?

Esta é a ponte mais importante do capítulo. A ANOVA de medidas repetidas não é uma espécie diferente de estatística. Continua a ser uma forma de decidir qual é o modelo compacto, qual é o aumentado e qual é o erro relevante para essa comparação.

O F de Medidas Repetidas

Num delineamento unifactorial, completo e equilibrado, com uma observação por pessoa e momento, podemos escrever:

\[ SS_{total} = SS_{unidade} + SS_{momento} + SS_{unidade\times momento}. \]

A parte unidade × momento diz-nos quanto os perfis individuais se afastam do perfil médio. Com uma única observação por célula, ela também contém variação residual que não conseguimos separar.

Para testar momento, usamos:

\[ F= \frac{MS_{momento}} {MS_{unidade\times momento}}. \]

Não é uma fórmula nova. É o mesmo raciocínio do F do capítulo anterior: melhoria atribuível à parte do modelo que libertámos, dividida pelo erro que ainda ficou e que é apropriado para essa pergunta.

Nos nossos dados,

\[ MS_{momento} = \frac{SS_{momento}}{df_{momento}} = \frac{6574.98}{2} = 3287.49, \]

e

\[ MS_{erro} = \frac{SS_{unidade\times momento}}{df_{erro}} = \frac{1.371137\times 10^{4}}{68} = 201.64. \]

Logo,

\[ F = \frac{3287.49} {201.64} = 16.3. \]

O software transforma depois este F num valor-p. Não precisamos de redefinir aqui o valor-p: usamos exactamente a interpretação introduzida em Do t ao Valor-p.

O que mudou neste capítulo foi qual é o erro adequado ao delineamento e, portanto, como construímos a estatística F; a lógica da inferência é a mesma.

ImportanteO Que a Unidade Fez

Ao pôr unidade no modelo, não «apagámos pessoas» nem corrigimos os dados.

Separámos duas perguntas:

• quanto variam as pessoas entre si?
• quanto varia a resposta dentro da mesma pessoa entre momentos?

Para testar momento, a segunda é a fonte de variação relevante.

Com Dois Momentos: O Teste-t Emparelhado

Se cada unidade tiver apenas dois momentos, podemos calcular uma diferença para cada pessoa:

\[ D_i=Y_{i,2}-Y_{i,1}. \]

A pergunta «há uma mudança média?» passa a ser «a média das diferenças é zero?». Portanto, o teste-t emparelhado é simplesmente um teste-t de uma amostra aplicado aos valores (D_i).

Voltamos assim ao que já conhecemos:

\[ t=\frac{\bar D}{SE_{\bar D}}. \]

Com dois níveis, a ANOVA de medidas repetidas e o teste-t emparelhado testam a mesma hipótese e, outra vez,

\[ F=t^2. \]

Esta ligação é útil porque mostra que «teste-t emparelhado» e «ANOVA de medidas repetidas» não são dois truques sem relação. Um é o caso mais pequeno do outro.

A Forma Clássica no R

Em R base podemos declarar explicitamente o estrato de erro:

aov_classica <- aov(
  resposta ~ momento + Error(unidade / momento),
  data = ds
)

summary(aov_classica)

Error: unidade
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals 34 142630    4195               

Error: unidade:momento
          Df Sum Sq Mean Sq F value   Pr(>F)    
momento    2   6575    3287    16.3 1.64e-06 ***
Residuals 68  13711     202                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

No trabalho corrente do site preferimos a interface do afex:

m_afex <- afex::aov_4(
  resposta ~ momento + (momento | unidade),
  data = ds
)

afex::nice(m_afex)

Effect	df	MSE	F	ges	p.value
momento	1.99, 67.58	202.90	16.30 ***	.040	<.001

A parte (momento | unidade) diz ao aov_4() quais são as medidas repetidas dentro da unidade. A função trata depois da contabilidade da ANOVA.

O objectivo desta página não é decorar as duas sintaxes. É poder olhar para uma delas e reconhecer a pergunta estatística que está a ser feita.

Porque o Truque do lm() Tem Limites

O modelo

y ~ unidade + momento

é uma excelente ponte pedagógica para um único factor within. Mas não é uma solução geral para todos os delineamentos repetidos.

Com dois factores repetidos, por exemplo, A e B, os testes clássicos usam fontes de erro diferentes para A, B e A × B. Aparecem componentes como

\[ SS_{unidade\times A},\qquad SS_{unidade\times B},\qquad SS_{unidade\times A\times B}. \]

Um único resíduo do modelo y ~ unidade + A * B não reproduz automaticamente todos esses denominadores.

Também há outro problema. Suponham que cada colónia pertence a uma espécie e é medida várias vezes. species não varia dentro da colónia. Se dermos um coeficiente fixo separado a cada colónia, podemos absorver precisamente as diferenças entre colónias que seriam necessárias para estudar species.

Portanto, a construção com lm() serve para ver a lógica, não para substituir todas as ANOVAs repetidas ou todos os modelos mistos.

Esfericidade

Com dois momentos há apenas uma diferença possível. Com três momentos já temos, por exemplo:

• manhã − tarde;
• manhã − noite;
• tarde − noite.

A ANOVA univariada clássica com três ou mais níveis precisa de uma condição sobre a forma como estas diferenças variam. Essa condição chama-se esfericidade.

Uma forma intuitiva de a ler é esta: o teste clássico não quer que uma diferença entre dois níveis seja tremendamente mais variável do que outra de uma forma que o seu denominador não representa.

Quando a esfericidade não é plausível, correcções como Greenhouse–Geisser e Huynh–Feldt ajustam os graus de liberdade do teste (Greenhouse & Geisser, 1959; Huynh & Feldt, 1976). Não alteram os dados e não tornam as observações independentes.

NotaDetalhe Matemático: O Que É Esfericidade?

Para um factor within com vários níveis, a esfericidade pode ser expressa em termos da matriz de covariância das medidas repetidas. Uma formulação equivalente exige que todos os contrastes normalizados do efeito tenham a mesma variância.

No caso simples de um único factor repetido, isto corresponde à conhecida condição de igualdade das variâncias das diferenças entre pares de níveis.

As correcções Greenhouse–Geisser e Huynh–Feldt estimam quanto a estrutura se afasta da condição ideal e reduzem os graus de liberdade usados no teste F.

Quando a ANOVA Começa a Ficar Apertada

A ANOVA de medidas repetidas é especialmente cómoda quando temos poucas condições bem definidas, um delineamento quase completo e uma observação ou um resumo por unidade × condição.

Começa a ficar menos natural quando aparecem:

• observações em falta;
• números diferentes de medições por unidade;
• tempos irregulares;
• várias unidades de agrupamento;
• ensaios individuais dentro de condições;
• efeitos que podem variar de unidade para unidade.

É aqui que os modelos mistos deixam de parecer uma técnica nova e passam a parecer a continuação natural da mesma pergunta: como representamos a estrutura que sabemos que existe nos dados?

2 Modelos Mistos

Um modelo misto mantém os efeitos que representam a nossa pergunta científica, mas acrescenta componentes que representam variação entre unidades.

Interceptos Aleatórios

Suponham várias medições de cada colónia. Algumas colónias podem ter valores sistematicamente altos e outras sistematicamente baixos.

Um intercepto aleatório permite essa variação de nível basal:

lme4::lmer(y ~ x + (1 | colonia), data = dados)

A parte (1 | colonia) diz, em linguagem informal:

«As colónias podem começar em níveis diferentes; quero estimar quanto esses níveis variam.»

As observações da mesma colónia passam a partilhar essa componente. É assim que o modelo representa parte da correlação entre linhas da mesma unidade.

NotaDetalhe Matemático: Intercepto Aleatório

Um modelo simples pode ser escrito como

\[ Y_{ij} = \beta_0+\beta_1X_{ij}+u_{0j}+\epsilon_{ij}, \]

com

\[ u_{0j}\sim N(0,\tau_0^2). \]

O índice (j) identifica a unidade. Duas observações da mesma unidade partilham (u_{0j}), e é essa componente partilhada que produz covariância entre elas no modelo marginal.

Partial Pooling

Há três ideias úteis para comparar.

Se ignorarmos completamente a colónia, tratamos todas as observações como se a unidade não importasse.

Se dermos a cada colónia um coeficiente fixo completamente separado, cada uma é estimada quase como um pequeno mundo próprio.

O modelo misto fica entre esses extremos. Estima a variação entre colónias e usa essa informação para estimar os desvios de cada uma. A isto chama-se agrupamento parcial (partial pooling).

Unidades com pouca informação tendem a ser mais puxadas para o padrão geral; unidades com muita informação conseguem afastar-se mais quando os dados o justificam.

NotaDetalhe Matemático: A Forma Matricial

Um modelo linear misto pode ser escrito como

\[ \mathbf y = \mathbf X\boldsymbol\beta + \mathbf Z\mathbf u + \boldsymbol\epsilon. \]

(X) contém os efeitos fixos; (Zu) distribui os efeitos aleatórios pelas observações. A inferência exige estimar também as variâncias e covariâncias associadas a (u), normalmente por máxima verosimilhança ou REML.

Os algoritmos usados para o fazer são importantes para quem estuda computação estatística, mas não são necessários para interpretar (1 | colonia).

Declives Aleatórios

Uma unidade pode diferir não só no ponto de partida, mas também na resposta ao preditor.

Talvez todas as colónias tenham, em média, maior metabolismo a temperaturas mais altas, mas algumas reajam muito e outras pouco. Nesse caso faz sentido perguntar se o declive varia entre colónias:

lme4::lmer(y ~ temperatura + (temperatura | colonia), data = dados)

O efeito fixo de temperatura descreve a tendência média. O declive aleatório permite que cada colónia se afaste dessa tendência.

Há uma regra conceptual muito simples:

só pode variar entre unidades um efeito que tenha sido observado dentro dessas unidades.

Se cada colónia só aparece a uma única temperatura, os dados não nos mostram um declive de temperatura para cada colónia. Um modelo não consegue inventar essa informação.

NotaDetalhe Matemático: Intercepto e Declive Aleatórios

Podemos escrever

\[ Y_{ij} = (\beta_0+u_{0j}) + (\beta_1+u_{1j})X_{ij} + \epsilon_{ij}. \]

O modelo pode ainda estimar a covariância entre (u_{0j}) e (u_{1j}): unidades com níveis basais mais altos podem ter, por exemplo, declives sistematicamente maiores ou menores. Essa estrutura adicional precisa de dados suficientes para ser estimada.

Efeitos Fixos e Aleatórios

Não usem a regra «poucos níveis = fixo; muitos níveis = aleatório». O número de níveis afecta aquilo que conseguimos estimar, mas não define sozinho a pergunta.

Perguntem antes:

Quero estimar diferenças específicas entre estes níveis, ou quero representar a variação entre unidades de uma população mais ampla?

Se queremos comparar especificamente duas espécies, a diferença entre essas espécies costuma ser um efeito fixo.

Se amostrámos vários recifes para representar variação espacial e não queremos uma conclusão separada para cada recife, recife pode funcionar como efeito aleatório.

A mesma variável pode desempenhar papéis diferentes em estudos diferentes.

3 Estrutura e Replicação

Within e Between Dependem da Unidade

Uma variável é within ou between relativamente a uma unidade.

Num estudo com peixes dentro de tanques:

• o tempo pode variar dentro de cada peixe;
• o tratamento pode variar entre tanques;
• vários peixes podem partilhar o mesmo tanque.

Portanto, dizer apenas «tratamento é between» é incompleto. Between o quê?

Esta pergunta é especialmente importante em dados biológicos, onde podemos ter organismos, tanques, parcelas, locais e campanhas todos na mesma base de dados.

Tanques, Peixes e Pseudorreplicação

Imaginem um tratamento térmico atribuído a tanques. Cada tanque contém vários peixes.

A estrutura é:

peixe → tanque → tratamento

Uma linha pode representar um peixe, mas o tratamento varia ao nível do tanque.

Se tivermos vários tanques por tratamento, um modelo possível é:

lme4::lmer(crescimento ~ tratamento + (1 | tanque), data = dados)

Os peixes do mesmo tanque partilham ambiente e tratamento, e o efeito aleatório representa essa semelhança.

Mas suponham agora que existe apenas um tanque de controlo e um tanque aquecido. Medir 30 peixes em cada tanque não cria 30 réplicas independentes do tratamento. Tratamento e tanque estão perfeitamente confundidos.

Este é o problema clássico da pseudorreplicação (Hurlbert, 1984).

ImportanteO Modelo Não Cria Réplicas

Um modelo misto pode representar dependência entre unidades que existem.

Não pode criar tanques que nunca foram replicados.

Antes de escolher a fórmula, identifiquem sempre a que nível o preditor foi atribuído ou varia independentemente.

Estruturas Aninhadas e Cruzadas

Nem todo o agrupamento tem a mesma geometria.

Se vários quadrados pertencem a cada recife, podemos ter uma estrutura aninhada:

lme4::lmer(cobertura ~ impacto + (1 | recife/quadrado), data = dados)

Se várias fotografias são classificadas por vários observadores, fotografias e observadores podem estar cruzados:

lme4::lmer(
  cobertura ~ tratamento + (1 | fotografia) + (1 | observador),
  data = dados
)

A sintaxe vem depois da pergunta:

que observações partilham a mesma fonte de variação?

Quanto Deve Variar?

A escolha da estrutura aleatória é uma área em que é fácil transformar princípios úteis em slogans.

Uma estrutura demasiado pobre pode ignorar variação importante. Uma estrutura demasiado rica pode pedir aos dados variâncias e covariâncias que eles não conseguem estimar.

Convergência e singularidade são sinais a investigar, não oráculos que decidem sozinhos o modelo científico.

NotaDetalhe Técnico: Maximal, Parcimonioso e Delineamento

Uma linha de trabalho, associada sobretudo a Barr e colaboradores, enfatizou estruturas aleatórias tão ricas quanto o delineamento justifica para proteger a inferência sobre efeitos fixos. Outra linha, associada ao trabalho de Bates e colaboradores, mostrou os problemas de sobreparametrização e defendeu simplificação informada pelos dados e pela identificabilidade.

Na prática, uma boa análise começa pelo delineamento: que interceptos e declives podem realmente variar entre as unidades amostradas? Depois verifica se há informação suficiente para estimar essa estrutura e documenta quaisquer simplificações. Harrison et al. (2018) discutem estes compromissos no contexto ecológico.

4 Um Roteiro Para Ler Estes Modelos

Quando virem uma análise de medidas repetidas ou um modelo misto, não comecem pela fórmula. Façam estas perguntas:

1. O que representa uma linha?
2. Qual é a unidade que aparece repetidamente?
3. Qual é a pergunta científica?
4. O preditor varia dentro ou entre que unidades?
5. Que diferenças estáveis entre unidades precisam de ser separadas do erro?
6. O efeito do preditor pode plausivelmente variar entre unidades?
7. Há outros agrupamentos, aninhados ou cruzados?
8. O delineamento contém réplicas independentes no nível em que queremos fazer inferência?

Depois disso, a fórmula deixa de parecer uma língua secreta. Passa a ser uma descrição abreviada das respostas.

5 Síntese

A passagem de modelos lineares para medidas repetidas não exige abandonar a ideia central do capítulo anterior.

Continuamos a perguntar:

Dados = Modelo + Erro. Que parte do modelo representa a pergunta e que variabilidade é o erro certo para a avaliar?

Quando uma unidade aparece várias vezes, o erro deixa de poder ser tratado como se todas as linhas fossem independentes.

A ANOVA de medidas repetidas separa diferenças entre unidades da variabilidade dentro das unidades. Num caso simples, conseguimos ver isso directamente pondo a unidade no modelo e comparando um modelo compacto com um aumentado.

Os modelos mistos levam a mesma ideia mais longe. Permitem representar interceptos, declives e várias fontes de agrupamento sem transformar cada unidade num coeficiente fixo separado.

E há um limite que nenhum método ultrapassa: um modelo pode reconhecer a estrutura do delineamento; não pode criar replicação que o delineamento não teve.

No tutorial prático usamos estas ideias para ajustar ANOVAs com aov_4() e primeiros modelos mistos.

6 Exercícios

1. Trinta pessoas são medidas antes e depois de uma intervenção. Explique por que há 60 linhas mas apenas 30 pares independentes. Mostre como o problema pode ser reduzido a uma média de diferenças.
2. No exemplo manhã/tarde/noite, explique por palavras o que fica no modelo compacto resposta ~ unidade e o que é libertado quando acrescentamos momento.
3. Explique por que o denominador do F para momento deve reflectir variabilidade dentro das unidades e não diferenças estáveis entre unidades.
4. Um tratamento é atribuído a quatro tanques e há vinte peixes em cada tanque. Qual é a unidade de replicação do tratamento? O que representam os peixes?
5. Explique a diferença substantiva entre (1 | colonia) e (temperatura | colonia).
6. Dê um exemplo de agrupamento aninhado e outro de agrupamento cruzado.
7. Um modelo misto é singular. Diga três perguntas que faria antes de remover termos automaticamente.

7 Leituras

Greenhouse, S. W., & Geisser, S. (1959). On methods in the analysis of profile data. Psychometrika, 24(2), 95–112. https://doi.org/10.1007/BF02289823

Harrison, X. A., Donaldson, L., Correa-Cano, M. E., Evans, J., Fisher, D. N., Goodwin, C. E. D., Robinson, B. S., Hodgson, D. J., & Inger, R. (2018). A brief introduction to mixed effects modelling and multi-model inference in ecology. PeerJ, 6, e4794. https://doi.org/10.7717/peerj.4794

Bolker, B. M., Brooks, M. E., Clark, C. J., Geange, S. W., Poulsen, J. R., Stevens, M. H. H., & White, J.-S. S. (2009). Generalized linear mixed models: A practical guide for ecology and evolution. Trends in Ecology & Evolution, 24(3), 127–135. https://doi.org/10.1016/j.tree.2008.10.008

Hurlbert, S. H. (1984). Pseudoreplication and the design of ecological field experiments. Ecological Monographs, 54(2), 187–211. https://doi.org/10.2307/1942661

Huynh, H., & Feldt, L. S. (1976). Estimation of the Box correction for degrees of freedom from sample data in randomized block and split-plot designs. Journal of Educational Statistics, 1(1), 69–82. https://doi.org/10.3102/10769986001001069