Dependência e Modelos Mistos

Desenho Experimental e Análise Avançada de Dados Ecológicos

João O. Santos

joao.santos@ispa.pt

ISPA

Descarregar os dados de briozoários.

Voltamos ao Primeiro Dia

peixe → tanque → tratamento

As linhas nunca foram automaticamente independentes.

Hoje

agrupamento · intercepto aleatório · medidas repetidas

Declives aleatórios ficam como extensão.

Doze Tanques

Exemplo simulado

6 controlo · 6 aquecidos · 5 peixes/tanque

Desenhar cinco peixes dentro de cada um de doze recipientes. Uma linha é um peixe; o tratamento varia entre tanques. Contar primeiro as unidades que recebem o tratamento, só depois as observações. Na demonstração, os peixes são distintos e os tanques são as réplicas; confirmar num estudo real se partilham reservatório, circuito de água ou outro agrupamento.

Conhecer os Tanques

CSV simulado: uma linha por peixe, crescimento em unidades arbitrárias.

ds <- read.csv("../../data/tank_growth_demo.csv") |>
  mutate(tank = factor(tank))
head(ds)

  tank treatment    growth
1    1   Control  7.381509
2    1   Control  7.935572
3    1   Control  7.820514
4    1   Control  7.434780
5    1   Control  8.759421
6    2   Control 12.446396

Este CSV conserva a simulação usada nesta aula: semente 12, 12 desvios de tanque de uma normal com SD 1.4, crescimento basal 10, acréscimo 2 nos seis tanques aquecidos e erro normal com SD 0.7 em cada um dos 60 peixes. Ler o ficheiro permite estudar a análise sem misturar a geração dos dados com o modelo que queremos ajustar.

Ver os Tanques

Preparar os Contrastes

options(contrasts = c("contr.sum", "contr.poly"))

A codificação dos factores não resolve a dependência.

Soma-zero para factores não ordenados, contrastes polinomiais para ordenados. A opção aplica-se aos modelos ajustados a partir daí nesta sessão. O tratamento varia entre tanques; a estrutura de agrupamento continua a ser necessária. aov_4() também assegura contrastes soma-zero nos factores entre unidades. Preparação no tutorial.

Solução Simples

Agregar ao nível da unidade independente.

tank_averages <- ds |>
  group_by(tank, treatment) |>
  summarise(growth = mean(growth), .groups = "drop")

lm(growth ~ treatment, data = tank_averages)

Call:
lm(formula = growth ~ treatment, data = tank_averages)

Coefficients:
(Intercept)   treatment1  
     10.291       -1.212  

Funciona — Mas Comprime

5 peixes → 1 média do tanque

Perdemos informação dentro do tanque.

Modelo dos Tanques

Uma linha = um peixe · tratamento entre tanques · crescimento em unidades arbitrárias neste exemplo simulado.

m_tank <- mixed(growth ~ treatment + (1 | tank), data = ds,
                method = "S", progress = FALSE)
m_tank

Mixed Model Anova Table (Type 3 tests, S-method)

Model: growth ~ treatment + (1 | tank)
Data: ds
     Effect       df      F p.value
1 treatment 1, 10.00 7.74 *    .019
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1

treatment_means <- emmeans(m_tank, ~ treatment)
treatment_means

 treatment emmean    SE df lower.CL upper.CL
 Control     9.08 0.616 10     7.71     10.5
 Heated     11.50 0.616 10    10.13     12.9

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

summary(pairs(treatment_means), infer = TRUE)

 contrast         estimate    SE df lower.CL upper.CL t.ratio p.value
 Control - Heated    -2.42 0.871 10    -4.37   -0.483  -2.782  0.0194

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

Script dos Tanques

Descarreguem o CSV simulado e adaptem o caminho em read.csv() à vossa pasta de trabalho.

library(afex)
library(dplyr)
library(emmeans)
library(ggplot2)

options(contrasts = c("contr.sum", "contr.poly"))

ds <- read.csv("../../data/tank_growth_demo.csv") |>
  mutate(tank = factor(tank))
m <- mixed(growth ~ treatment + (1 | tank), data = ds, method = "S")
treatment_means <- emmeans(m, ~ treatment)
g <- ggplot(ds, aes(tank, growth, fill = treatment)) + geom_boxplot()

# Print results and plot
print(m)
print(treatment_means)
pairs(treatment_means) |> summary(infer = TRUE) |> print()
print(g)

Resultados dos Tanques

Mixed Model Anova Table (Type 3 tests, S-method)

Model: growth ~ treatment + (1 | tank)
Data: ds
     Effect       df      F p.value
1 treatment 1, 10.00 7.74 *    .019
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1

 treatment emmean    SE df lower.CL upper.CL
 Control     9.08 0.616 10     7.71     10.5
 Heated     11.50 0.616 10    10.13     12.9

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

 contrast         estimate    SE df lower.CL upper.CL t.ratio p.value
 Control - Heated    -2.42 0.871 10    -4.37   -0.483  -2.782  0.0194

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

Relatar os Tanques

O crescimento estimado foi menor no controlo do que no aquecimento: Control − Heated = −2.42 unidades, 95% CI [−4.37, −0.48], t₁₀ = −2.78, P = .019. São dados simulados: há seis tanques por tratamento, não 30 réplicas. Os cinco peixes por tanque são subamostras. O modelo supõe tanques independentes e resíduos independentes dentro de cada tanque depois de considerar o seu intercepto; não representa um eventual circuito de água partilhado.

Com contrastes soma-zero, o coeficiente de tratamento é metade da diferença entre médias; pairs() devolve a diferença inteira, Control − Heated.

O Que Varia?

intercepto do tanque

Cada tanque pode ter um nível basal diferente.

Uma Recta Por Grupo

\[ \beta_{0j}=\beta_0+b_{0j} \]

intercepto global + desvio do grupo

Declives Também Variam

\[ \beta_{1j}=\beta_1+b_{1j} \]

declive global + desvio da unidade

Extensão em R

Sintaxe esquemática: exige temperatura variável dentro de cada colónia; não se aplica ao CSV dos tanques.

Neste exemplo esquemático, a temperatura tem de variar dentro de cada colónia para estimarmos declives por colónia. Comparar com o tratamento nos tanques: cada tanque recebe uma só condição, por isso não temos um declive de tratamento dentro de cada tanque. Que quantidade e distribuição de medições permitiriam estimar a variação dos declives?

mixed(
  metabolism ~ temperature +
    (temperature | colony),
  data = ds
)

Não É Magia

O delineamento tem de fornecer informação para estimar essa variação.

Medidas Repetidas

ds <- read.csv("../../data/bryozoan_wrangled.csv") |>
  filter(Stage %in% c("EARLY", "LATE")) |>
  mutate(Unit = interaction(Species, Id, drop = TRUE))

A mesma unidade aparece em mais de uma linha.

Às Vezes Não Sabemos

mesma unidade · IDs diferentes

Identificação falhada pode esconder a repetição.

Dois casos hipotéticos: um peixe com marca ou chip reconhecido nas duas visitas; um peixe parecido, sem marca, que recebe um novo ID em cada captura. No segundo caso, duas linhas podem ser o mesmo animal. Que registo permitiria confirmar a identidade antes de construir Unit?

Linhas da Mesma Unidade

Species	Id	Stage	Unit
Burgula	ID372	EARLY	Burgula.ID372
Burgula	ID372	LATE	Burgula.ID372
Watersipora	ID173	EARLY	Watersipora.ID173
Watersipora	ID173	LATE	Watersipora.ID173

Unit: chave construída ≠ identidade biológica verificada

Uma linha é uma medição; duas linhas partilham a unidade em estádios distintos. Mostramos a primeira unidade de cada espécie: quatro linhas reais. A chave construída é interaction(Species, Id, drop = TRUE). Neste ficheiro os Id dos dois grupos não colidem; incluir Species torna o agrupamento explícito, não demonstra que os organismos foram reconhecidos. Essa confirmação depende dos registos de recolha e dos métodos do estudo. Há 255 Unit (197 Burgula, 58 Watersipora), com dois estádios por unidade.

Briozoários

De Onde Vieram os Briozoários?

Pettersen et al. (2015)

Why Does Offspring Size Affect Performance? Integrating Metabolic Scaling With Life-History Theory

Ler o Delineamento

15 minutos · grupos de 3 ou 5

Procurem nos Métodos como foram acompanhados os organismos.

1. Que medições pertencem ao mesmo indivíduo?
2. Que factor é intra-unidades?
3. O que aconteceria à incerteza se tratássemos essas linhas como independentes?

Intra- e Inter-Unidades

Species → inter-unidades

Stage → intra-unidades

O Que Se Repete?

muitas Unit na mesma Species → normal

a mesma Unit em vários Stage → dependência

O rótulo do grupo pode repetir-se. A unidade é que não pode fingir ser nova.

A Pergunta Não Mudou

A mudança entre estádios difere entre espécies?

ANOVA de Medidas Repetidas

m_rm <- aov_4(
  Log10MR ~ Species * Stage + (Stage | Unit),
  data = ds
)
m_rm

Anova Table (Type 3 tests)

Response: Log10MR
         Effect     df  MSE          F  ges p.value
1       Species 1, 253 0.03 288.52 *** .414   <.001
2         Stage 1, 253 0.02 757.70 *** .533   <.001
3 Species:Stage 1, 253 0.02  16.84 *** .025   <.001
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1

(Stage | Unit) especifica as medidas repetidas dentro de Unit em aov_4. Relaciona-se com a ideia de variação dentro de unidades, mas não deve ser lido como uma troca automática pela sintaxe de declives aleatórios de lmer.

Mudança por Espécie

stage_means <- emmeans(m_rm, ~ Stage | Species)
stage_changes <- pairs(stage_means)
summary(stage_means)

Species = Burgula:
 Stage  emmean     SE  df lower.CL upper.CL
 EARLY -0.7906 0.0125 253   -0.815   -0.766
 LATE  -0.4443 0.0102 253   -0.464   -0.424

Species = Watersipora:
 Stage  emmean     SE  df lower.CL upper.CL
 EARLY -0.5315 0.0231 253   -0.577   -0.486
 LATE  -0.0639 0.0187 253   -0.101   -0.027

Confidence level used: 0.95 

summary(stage_changes, infer = TRUE)

Species = Burgula:
 contrast     estimate     SE  df lower.CL upper.CL t.ratio p.value
 EARLY - LATE   -0.346 0.0141 253   -0.374   -0.318 -24.557 <0.0001

Species = Watersipora:
 contrast     estimate     SE  df lower.CL upper.CL t.ratio p.value
 EARLY - LATE   -0.468 0.0260 253   -0.519   -0.416 -17.993 <0.0001

Confidence level used: 0.95 

summary(pairs(stage_changes, by = NULL), infer = TRUE)

 contrast                                            estimate     SE  df
 (EARLY - LATE Burgula) - (EARLY - LATE Watersipora)    0.121 0.0296 253
 lower.CL upper.CL t.ratio p.value
   0.0631     0.18   4.103 <0.0001

Confidence level used: 0.95 

• EARLY − LATE negativo: valores maiores em LATE
• Última comparação: comparar as mudanças, não dois valores-P
• Intervalos pontuais de 95%, na escala log10

Script Completo

Descarreguem o CSV e adaptem o caminho em read.csv() à vossa pasta de trabalho. A chave junta linhas; os Métodos é que confirmam a identidade.

library(afex)
library(dplyr)
library(emmeans)
library(ggplot2)

options(contrasts = c("contr.sum", "contr.poly"))

ds <- read.csv("../../data/bryozoan_wrangled.csv") |>
  filter(Stage %in% c("EARLY", "LATE")) |>
  mutate(Unit = interaction(Species, Id, drop = TRUE))
m_rm <- aov_4(Log10MR ~ Species * Stage + (Stage | Unit), data = ds)
stage_means <- emmeans(m_rm, ~ Stage | Species)
stage_changes <- pairs(stage_means)
g <- ggplot(ds, aes(Stage, Log10MR, group = Unit)) +
  geom_line(alpha = .2) + facet_wrap(~ Species)

# Print results and plot
print(m_rm)
print(stage_means)
summary(stage_changes, infer = TRUE) |> print()
pairs(stage_changes, by = NULL) |> summary(infer = TRUE) |> print()
print(g)

Resultados das Mudanças

Anova Table (Type 3 tests)

Response: Log10MR
         Effect     df  MSE          F  ges p.value
1       Species 1, 253 0.03 288.52 *** .414   <.001
2         Stage 1, 253 0.02 757.70 *** .533   <.001
3 Species:Stage 1, 253 0.02  16.84 *** .025   <.001
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1

Species = Burgula:
 Stage  emmean     SE  df lower.CL upper.CL
 EARLY -0.7906 0.0125 253   -0.815   -0.766
 LATE  -0.4443 0.0102 253   -0.464   -0.424

Species = Watersipora:
 Stage  emmean     SE  df lower.CL upper.CL
 EARLY -0.5315 0.0231 253   -0.577   -0.486
 LATE  -0.0639 0.0187 253   -0.101   -0.027

Confidence level used: 0.95 

Species = Burgula:
 contrast     estimate     SE  df lower.CL upper.CL t.ratio p.value
 EARLY - LATE   -0.346 0.0141 253   -0.374   -0.318 -24.557 <0.0001

Species = Watersipora:
 contrast     estimate     SE  df lower.CL upper.CL t.ratio p.value
 EARLY - LATE   -0.468 0.0260 253   -0.519   -0.416 -17.993 <0.0001

Confidence level used: 0.95 

 contrast                                            estimate     SE  df
 (EARLY - LATE Burgula) - (EARLY - LATE Watersipora)    0.121 0.0296 253
 lower.CL upper.CL t.ratio p.value
   0.0631     0.18   4.103 <0.0001

Confidence level used: 0.95 

Ler a Tabela

m_rm

Anova Table (Type 3 tests)

Response: Log10MR
         Effect     df  MSE          F  ges p.value
1       Species 1, 253 0.03 288.52 *** .414   <.001
2         Stage 1, 253 0.02 757.70 *** .533   <.001
3 Species:Stage 1, 253 0.02  16.84 *** .025   <.001
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '+' 0.1 ' ' 1

Species:Stage: a mudança EARLY → LATE é igual nas duas espécies? A tabela respeita a repetição dentro de Unit; não é a ANOVA de 510 unidades independentes.

O modelo representa a dependência dentro de Unit e supõe unidades independentes. A ideia de testar restrições continua a aplicar-se, mas não copiar para aqui as comparações lm() da aula anterior ignorando a covariância intra-unidade. O teste de interacção usa a variação das mudanças dentro das unidades. Com dois estádios não há uma condição adicional de esfericidade a verificar.

Relatar em Biologia

A diferença EARLY − LATE na taxa metabólica log10 foi -0.346 em Burgula, 95% CI [-0.374, -0.318], e -0.468 em Watersipora, 95% CI [-0.519, -0.416]. A diferença entre estas mudanças foi 0.121 unidades log10, 95% CI [0.063, 0.180]; a interacção na ANOVA de medidas repetidas foi F_{1, 253} = 16.84, P < .001. As medições pertencem a 255 unidades; a generalização depende da amostragem e da confirmação da identidade nos registos do estudo.

Um Tanque por Tratamento?

Nenhum modelo inventa réplicas que não foram recolhidas.

Três Perguntas Outra Vez

1. O que representa uma linha?
2. Que linhas pertencem ao mesmo grupo?
3. A que nível varia o preditor?

Agora: Delineamento Primeiro

15 minutos · grupos de 3 ou 5

Medidas Repetidas no R — Exercícios

Façam exercícios 5, 6 e 7.

Depois: Medidas Repetidas

20 minutos · grupos de 3 ou 5

Medidas Repetidas no R — Exercícios

Façam exercícios 3 e 4.

Explicar o Modelo

10 minutos · grupos de 3 ou 5

Escolham tanques ou briozoários e expliquem, sem usar a palavra «aleatório»:

1. que linhas andam juntas;
2. o que pode variar entre grupos;
3. porque tratar todas as linhas como independentes mudaria a incerteza.

Se Precisarem de Teoria

Medidas Repetidas e Modelos Mistos

O Que Mudou?

Agora conseguimos manter as observações no nível mais fino sem fingir que todas as linhas são independentes.