Medidas repetidas no R

Modelos lineares
Medidas repetidas
afex
modelos mistos
emmeans
Quando a mesma pessoa aparece várias vezes nos teus dados…ANOVAs de medidas repetidas e modelos mistos no R
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

Setup

library(ggplot2)
library(dplyr)
library(tidyr)

library(afex)
library(emmeans)

# Recommended defaults
options(contrasts = c("contr.sum", "contr.poly"))
afex_options(es_aov = "pes")

Um delineamento de medidas repetidas significa que mediste a mesma unidade (participante, animal, turma…) mais do que uma vez. Se ignorares isto e correres um lm()/uma ANOVA entre-sujeitos “normal”, podes subestimar a incerteza e acabar com valores-p que parecem “bons demais para ser verdade”.

Os Dados

Podes descarregar os ficheiros CSV que acompanham este tutorial.

paired    <- read.csv("../data/rm_paired.csv")
oneway    <- read.csv("../data/rm_oneway.csv")
mixed     <- read.csv("../data/rm_mixed.csv")
within2_2 <- read.csv("../data/rm_within2_2.csv")

aov() + Error(...)

O R base consegue fazer ANOVAs de medidas repetidas com aov(... + Error(pp / fator_within)). O ponto importante é que tens de dizer explicitamente ao R o que pertence à mesma pessoa/unidade de observação.

Teste-t emparelhado como ANOVA MR

Um teste-t emparelhado é simplesmente uma ANOVA one-way de medidas repetidas com 2 níveis.

paired <- pivot_longer(paired, c("pre", "post"),
                       names_to = "time", values_to = "score")

m_aov_paired <- aov(score ~ time + Error(pp / time), paired)

summary(m_aov_paired)

Error: pp
          Df Sum Sq Mean Sq F value Pr(>F)
Residuals 29   3293   113.5               

Error: pp:time
          Df Sum Sq Mean Sq F value   Pr(>F)    
time       1  323.8   323.8    55.5 3.29e-08 ***
Residuals 29  169.2     5.8                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

One-way RM (k = 3)

m_aov_oneway <- aov(stress ~ time + Error(pp / time), oneway)
summary(m_aov_oneway)

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

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

Delineamento Misto (2b × 2w)

Aqui group é entre-participantes (between-subjects) e season é intra-participantes (within-subjects).

m_aov_mixed <- aov(score ~ group * season + Error(pp / season), mixed)
summary(m_aov_mixed)

Error: pp
          Df Sum Sq Mean Sq F value Pr(>F)
group      1    108   107.8    0.21  0.649
Residuals 46  23607   513.2               

Error: pp:season
             Df Sum Sq Mean Sq F value   Pr(>F)    
season        1 1771.4  1771.4   76.92 2.22e-11 ***
group:season  1  392.9   392.9   17.06 0.000151 ***
Residuals    46 1059.3    23.0                     
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
NotaPorque é que o aov() base não é “o fim da história”

O aov() base é ótimo para aprender a ideia central: “põe o fator repetido dentro do Error(pp / ...)”.

Mas não te ajuda muito com problemas práticos recorrentes: correções de esfericidade, múltiplas observações por participante × célula (dados por ensaio), dados em falta, estruturas complexas de efeitos aleatórios, etc.

afex::aov_4()

O aov_4() dá-te resultados de ANOVA de medidas repetidas, mas:

  • Exige uma coluna de ID (neste caso, pp), porque se tiveres múltiplas observações por participante e por célula do delineamneto, o afex pode agregar (normalmente, fazendo a média) em vez de fingir que são observações independentes. Mesmo que não haja medições repetidas o afex::aov_4() exige esta coluna, porque ele próprio confirma se há medidas repetidas quando te esqueceste das considerar, ou se não há medidas repetidas (ou se nem todos os participantes estão em todas as condições intra-pp) quando disseste que havia.

  • Reporta tabelas da ANOVA com soma de quadrados tipo III e suporta reportar graus de liberdade corrigidos por esfericidade para efeitos within-subjects.

One-way RM com aov_4() 

m_aov4_oneway <- aov_4(stress ~ time + (time | pp), oneway)
nice(m_aov4_oneway)
Effect df MSE F pes p.value
time 1.99, 67.58 202.90 16.30 *** .324 <.001

Delineamento Misto com aov_4() 

m_aov4_mixed <- aov_4(score ~ group * season + (season | pp), mixed)
nice(m_aov4_mixed)
Effect df MSE F pes p.value
group 1, 46 513.20 0.21 .005 .649
season 1, 46 23.03 76.92 *** .626 <.001
group:season 1, 46 23.03 17.06 *** .271 <.001

emmeans Contínua a Ser o Herói

# Simple effects: winter vs summer within each group
emmeans(m_aov4_mixed, list(pairwise ~ season | group, pairwise ~ group | season))
$`emmeans of season | group`
group = control:
 season emmean   SE df lower.CL upper.CL
 winter   98.8 3.24 46     92.3    105.3
 summer   94.2 3.44 46     87.3    101.2

group = treatment:
 season emmean   SE df lower.CL upper.CL
 winter  104.9 3.24 46     98.4    111.5
 summer   92.3 3.44 46     85.4     99.2

Confidence level used: 0.95 

$`pairwise differences of season | group`
group = control:
 2               estimate   SE df t.ratio p.value
 winter - summer     4.55 1.39 46   3.281  0.0020

group = treatment:
 2               estimate   SE df t.ratio p.value
 winter - summer    12.64 1.39 46   9.122 <0.0001


$`emmeans of group | season`
season = winter:
 group     emmean   SE df lower.CL upper.CL
 control     98.8 3.24 46     92.3    105.3
 treatment  104.9 3.24 46     98.4    111.5

season = summer:
 group     emmean   SE df lower.CL upper.CL
 control     94.2 3.44 46     87.3    101.2
 treatment   92.3 3.44 46     85.4     99.2

Confidence level used: 0.95 

$`pairwise differences of group | season`
season = winter:
 2                   estimate   SE df t.ratio p.value
 control - treatment    -6.17 4.59 46  -1.344  0.1856

season = summer:
 2                   estimate   SE df t.ratio p.value
 control - treatment     1.93 4.86 46   0.396  0.6937
DicaUma observação por célula vs muitos ensaios

Se os teus dados “brutos” têm vários ensaios por participante × condição, tens duas opções sensatas:

  • Agregar primeiro (ex.: média do RT por participante × condição).

  • Manter os dados ao nível do ensaio e passar para modelos mistos (secção seguinte), porque aí consegues modelar o ruído por ensaio e efeitos aleatórios por item/estímulo.

O afex::aov_4()/aov_car() também agrega automaticamente nos seus helpers de ANOVA quando deteta múltiplas observações por célula.

afex::mixed(): a mesma fórmula que lme4::lmer()

Já tens vindo a usar sintaxe do lme4 (o “gold standard” de modelos mistos no R) sem saber!

O afex::mixed() basicamente deixa-te manter essa sintaxe, mas dá-te testes estilo ANOVA para efeitos fixos (com defaults razoáveis para graus de liberdade / valores-p), além de alguma conveniência extra para o ajuste do modelo.

Delienemaneto 2w × 2w

# Anova de Medições Repetidas
m <- aov_4(rt ~ difficulty * duration + (difficulty * duration | pp), within2_2)


# Modelo Misto
# Com várias observações por célula do delineamento, substituir aov_4 por mixed
#m <- mixed(rt ~ difficulty * duration + (difficulty * duration | pp), within2_2)
# Com apenas uma observação por célula usar apenas interceptos aleatórios
m <- mixed(rt ~ difficulty * duration + (1 | pp), within2_2)
nice(m)
Effect df F p.value
difficulty 1, 117.00 102.01 *** <.001
duration 1, 117.00 0.01 .928
difficulty:duration 1, 117.00 6.63 * .011

dat_prop <- data.frame(grupo = factor(c(“A”, “B”, “C”)), sucessos = c(12, 20, 28), total = c(40, 40, 40)) ### Follow-ups planeados

emmeans(m, list(pairwise ~ difficulty | duration, pairwise ~ duration | difficulty))
$`emmeans of difficulty | duration`
duration = long:
 difficulty emmean SE   df lower.CL upper.CL
 easy          622 11 63.6      600      644
 hard          667 11 63.6      646      689

duration = short:
 difficulty emmean SE   df lower.CL upper.CL
 easy          607 11 63.6      585      629
 hard          684 11 63.6      662      705

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

$`pairwise differences of difficulty | duration`
duration = long:
 2           estimate   SE  df t.ratio p.value
 easy - hard    -45.5 8.55 117  -5.321 <0.0001

duration = short:
 2           estimate   SE  df t.ratio p.value
 easy - hard    -76.7 8.55 117  -8.962 <0.0001

Degrees-of-freedom method: kenward-roger 

$`emmeans of duration | difficulty`
difficulty = easy:
 duration emmean SE   df lower.CL upper.CL
 long        622 11 63.6      600      644
 short       607 11 63.6      585      629

difficulty = hard:
 duration emmean SE   df lower.CL upper.CL
 long        667 11 63.6      646      689
 short       684 11 63.6      662      705

Degrees-of-freedom method: kenward-roger 
Confidence level used: 0.95 

$`pairwise differences of duration | difficulty`
difficulty = easy:
 2            estimate   SE  df t.ratio p.value
 long - short     15.0 8.55 117   1.757  0.0815

difficulty = hard:
 2            estimate   SE  df t.ratio p.value
 long - short    -16.1 8.55 117  -1.884  0.0620

Degrees-of-freedom method: kenward-roger

Efeitos Aleatórios Maximais

Eu sigo a recomendação de Singmann e Kellen (2019): começa com a estrutura maximal de efeitos aleatórios justificada pelo design, e depois simplifica apenas se for preciso (i.e., sempre que há erros ou avisos no ajuste).

No entanto, para podermos ter declives aleatórios para uma variável temos de ter mais que uma observação por célula do delineamento. Mesmo quando temos várias medições por célula o modelo máximo pode não convergir (i.e., dar erros) ou ter problemas a convergir (i.e., dar avisos). Nessa situação, os dados podem não ter informação suficiente para estimar todos os declives aleatórios e correlações, e o modelo pode ficar singular (aviso) ou nem sequer convergir (erro).

Receita prática para simplificar

Quando o mixed() se “queixar”, simplifica por esta ordem (Singmann e Kellen, 2019):

  1. Remove correlações entre efeitos aleatórios: usa || em vez de | e define expand_re = TRUE.

  2. Remove primeiro os declives aleatórios das interações.

  3. Se ainda for preciso, remove os declives mais “pequenos” (mantém os que correspondem aos teus preditores intra-particiapantes importantes)

  4. Se não conseguires ajustar o modelo com declives aleatórios, experimenta ter apenas interceptos aleatórios (e considera se não será melhor/mais simples uma ANOVA de medições repetidas clássica).

  5. Não faças batota nos efeitos fixos só para “dar significativo”.

Importante

MAS não te esqueças que só podes ter declives aleatórios para uma dada variável/fator se tiveres mais que uma observação por célula do delineamento experimental. Se cada participante só completou uma medição em cada combinação de dificuldade X duração, não será possível ter declives aleatórios para essas variáveis no modelo.

# 0. Máximo (se tivessemos mais que uma medição por célula)
#m0 <- mixed(rt ~ difficulty * duration + (difficulty * duration | pp), within2_2)

# 1.
# Also won't unless you have more than one observation per design cell
#m1 <- mixed(rt ~ difficulty * duration + (difficulty * duration || pp),
#           within2_2, expand_re = TRUE)

# 2.
# Also won't unless you have more than one observation per design cell
#m2 <- mixed(rt ~ difficulty * duration + (difficulty + duration || pp),
#           within2_2, expand_re = TRUE)

# 3.
#m3 <- mixed(rt ~ difficulty * duration + (difficulty || pp),
#           within2_2, expand_re = TRUE)

m <- mixed(rt ~ difficulty * duration + (1 | pp), within2_2, expand_re = TRUE)

nice(m)
Effect df F p.value
difficulty 1, 117.00 102.01 *** <.001
duration 1, 117.00 0.01 .928
difficulty:duration 1, 117.00 6.63 * .011
NotaOnde está a “teoria” de modelos mistos?

Este tutorial é um tutorial focado no código e o workflow de análise.

O capítudolo teórico vai cobrir o significado de interceptos/declives aleatórios, porque é que (algumas) correlações são difíceis de estimar, o que é “partial pooling”, e como é que estes modelos se relacionam com ANOVA de medidas repetidas.

Cábula

  • ANOVA MR (R base): aov(y ~ A*B + Error(pp / (fatores within...)), data=...)
  • ANOVA MR (afex): aov_4(y ~ A*B + (fatores within | pp), data=...)
  • Modelo misto: mixed(y ~ fixos + (aleatórios | pp) + (aleatórios | item), data=...)

Se só memorizares uma coisa: tem sempre uma variável de ID (aqui, pp), e diz sempre ao R o que se repete dentro desse pp.