Interacções

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

João O. Santos

joao.santos@ispa.pt

ISPA

Uma Pergunta

Os nutrientes aumentam o crescimento dos corais da mesma forma em condições de baixa e alta luminosidade?

Exemplo simulado · crescimento em mm ao fim de oito semanas

Antes do Modelo

Uma linha = crescimento de um fragmento de coral em oito semanas.

• Fragmentos provenientes de colónias diferentes
• Um aquário por fragmento; luminosidade e circuito de água independentes
• Condições atribuídas aleatoriamente

64 observações independentes · 16 por combinação.

Dois Factores, Quatro Células

Condições de Luminosidade	Controlo	Adição de Nutrientes
Baixa	16 observações	16 observações
Alta	16 observações	16 observações

2 × 2: duas condições de luminosidade × duas condições de nutrientes.

Factores entre observações: uma combinação por fragmento; sem repetir o fragmento noutras condições.

Preparar o Script

library(car)
library(dplyr)
library(emmeans)
library(ggplot2)
library(palmerpenguins)
library(tidyr)
options(contrasts = c("contr.sum", "contr.poly"))

Contrastes soma-zero · diferenças médias com peso igual para os níveis do outro factor.

Interacção

O efeito de uma variável depende do nível de outra.

Nos corais, perguntamos se o efeito dos nutrientes depende das condições de luminosidade.

As linhas não precisam de se cruzar. Basta que a distância entre elas mude.

Perguntas de Interacção

• O aquecimento altera o efeito da acidificação no crescimento de larvas?
• A adição de nutrientes muda a cobertura de macroalgas de forma diferente com e sem herbívoros?
• A relação entre salinidade e sobrevivência difere entre espécies?

Em todas: «depende de quê?» já faz parte da pergunta biológica.

Oito Combinações

Dois efeitos principais e uma interacção: \(2^3=8\) combinações. Simulámos um cenário para cada combinação.

Os títulos indicam os efeitos presentes na construção. Nestes exemplos, os testes também correspondem aos títulos: p < .05 para os efeitos presentes e p ≥ .05 para os restantes.

Cada cenário aparece em duas versões dos mesmos dados: as mesmas quatro médias e os mesmos intervalos de confiança. Só trocamos os factores de lugar.

Dados construídos · São exemplos para ler padrões, não resultados empíricos.

Ler os Centros

Barra contínua: diferença entre os centros das linhas. Círculos vazios e barra tracejada: médias à esquerda e à direita; diferença do factor no eixo x.

Pesos iguais nas médias; nenhuma destas barras é um intervalo de confiança.

Sem Efeitos

Crescimento semelhante nas quatro condições.

Só Luminosidade

Mais luz → mais crescimento, com ou sem adição de nutrientes.

Só Nutrientes

Adição de nutrientes → mais crescimento, com baixa ou alta luminosidade.

Dois Efeitos Principais

Mais luz e mais nutrientes → mais crescimento; efeitos aditivos.

Só Interacção

Nutrientes prejudiciais com pouca luz, benéficos com muita luz; efeitos opostos que se anulam em média.

Luminosidade e Interacção

Mais crescimento com mais luz, em média; nutrientes prejudiciais com pouca luz, benéficos com muita luz.

Nutrientes e Interacção

Adição de nutrientes → mais crescimento, sobretudo com muita luz.

Três Efeitos

Mais luz e mais nutrientes → mais crescimento; ganho com nutrientes maior com muita luz.

O Que Mudou?

Vimos as oito combinações de resultados dos testes, não todas as formas possíveis de interacção.

A mesma tabela de «significativo/não significativo» pode acompanhar padrões biológicos muito diferentes. Voltem sempre ao gráfico e à pergunta.

Os oito conjuntos reutilizam o mesmo ruído para facilitar a comparação entre padrões. Não são experiências independentes nem resultados empíricos. O script de geração guarda as sementes, as condições e as verificações dos oito padrões. Uma diferença não significativa não prova ausência de efeito, e a interacção não se define por um teste simples significativo e outro não significativo.

Conhecer os Corais

CSV simulado · um registo por fragmento.

ds <- read.csv("../../data/coral_growth_2x2_demo.csv")

head(ds)

  microcosm light nutrients growth_mm
1       C01  High   Control  3.971323
2       C02   Low   Control  3.470618
3       C03   Low  Enriched  5.171945
4       C04   Low   Control  4.183034
5       C05  High  Enriched  6.999499
6       C06  High   Control  3.579312

growth_mm: crescimento nas oito semanas; Low/High: baixa/alta luminosidade; Control/Enriched: controlo/adição de nutrientes.

Um Modelo com Interacção

m_coral <- lm(growth_mm ~ light * nutrients, data = ds)

light * nutrients inclui luminosidade + nutrientes + luminosidade × nutrientes.

Não temos de provar primeiro que há efeitos principais para estudar uma interacção: foi essa a pergunta que motivou a experiência.

As Médias das Quatro Células

emm <- emmeans(m_coral, ~ light * nutrients)
print(emm)

 light nutrients emmean    SE df lower.CL upper.CL
 High  Control     4.94 0.165 60     4.61     5.27
 Low   Control     4.16 0.165 60     3.83     4.49
 High  Enriched    8.07 0.165 60     7.74     8.40
 Low   Enriched    4.51 0.165 60     4.18     4.84

Confidence level used: 0.95 

coral_plot <- as.data.frame(emm)

• emmeans(): médias estimadas e sua incerteza
• as.data.frame(): tabela para o ggplot()

O Nosso Padrão

O que nos diz este padrão?

Descrevam a distância entre as linhas em cada condição de luminosidade.

Nutrientes por Luminosidade

• Barras pretas: diferenças entre nutrientes na mesma condição de luminosidade
• Barras coloridas: intervalos de confiança das médias

Trocar a Pergunta

• Comparação: luminosidade, mantendo os nutrientes fixos
• Barras pretas: mesmas quatro médias, noutra direcção

Diferença de Diferenças

\[ (\widehat\mu_{\mathrm{High,Enriched}}-\widehat\mu_{\mathrm{High,Control}})- (\widehat\mu_{\mathrm{Low,Enriched}}-\widehat\mu_{\mathrm{Low,Control}}) \]

3.13 − 0.34 ≈ 2.79 mm: quanto muda a diferença entre nutrientes ao mudar as condições de luminosidade.

Esta é a interacção no 2 × 2, não uma comparação entre valores-p.

Garcia-Marques et al. (2014): interpretem o padrão completo das quatro células.

Efeito Principal: Fazer a Média

• Contínua: entre centros das linhas · tracejada: entre médias no eixo x
• Pesos iguais ao fazer a média sobre o outro factor
• Diferenças entre médias, não intervalos de confiança

Três Perguntas Diferentes

Pergunta	O Que Comparamos?
Luminosidade, em média?	Médias por luminosidade após fazer a média sobre nutrientes
Nutrientes, em média?	Médias por nutrientes após fazer a média sobre luminosidade
Interacção?	A diferença entre nutrientes nas duas condições de luminosidade

Um efeito principal resume o padrão. Não substitui olhar para as quatro células.

Ler a Tabela

Anova(m_coral, type = "III") |> print()

Anova Table (Type III tests)

Response: growth_mm
                 Sum Sq Df  F value    Pr(>F)    
(Intercept)     1880.51  1 4332.454 < 2.2e-16 ***
light             75.48  1  173.885 < 2.2e-16 ***
nutrients         48.25  1  111.162 2.776e-15 ***
light:nutrients   31.01  1   71.434 8.338e-12 ***
Residuals         26.04 60                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

light:nutrients: a diferença entre nutrientes é igual nas duas condições de luminosidade?

Cada termo deste 2 × 2 tem 1 gl: o seu F ainda coincide com o \(t^2\) do coeficiente correspondente. O teste global do resumo tem 3 gl e é outra pergunta.

Comparar Nutrientes

emm <- emmeans(m_coral, pairwise ~ nutrients | light)
print(emm)

$emmeans
light = High:
 nutrients emmean    SE df lower.CL upper.CL
 Control     4.94 0.165 60     4.61     5.27
 Enriched    8.07 0.165 60     7.74     8.40

light = Low:
 nutrients emmean    SE df lower.CL upper.CL
 Control     4.16 0.165 60     3.83     4.49
 Enriched    4.51 0.165 60     4.18     4.84

Confidence level used: 0.95 

$contrasts
light = High:
 contrast           estimate    SE df t.ratio p.value
 Control - Enriched   -3.129 0.233 60 -13.432 <0.0001

light = Low:
 contrast           estimate    SE df t.ratio p.value
 Control - Enriched   -0.344 0.233 60  -1.479  0.1444

pairwise pede médias e comparações no mesmo resultado. Em a | b, comparamos a dentro de cada nível de b. Trocar a ordem muda a pergunta.

Comparar Luminosidade

emm <- emmeans(m_coral, pairwise ~ light | nutrients)
print(emm)

$emmeans
nutrients = Control:
 light emmean    SE df lower.CL upper.CL
 High    4.94 0.165 60     4.61     5.27
 Low     4.16 0.165 60     3.83     4.49

nutrients = Enriched:
 light emmean    SE df lower.CL upper.CL
 High    8.07 0.165 60     7.74     8.40
 Low     4.51 0.165 60     4.18     4.84

Confidence level used: 0.95 

$contrasts
nutrients = Control:
 contrast   estimate    SE df t.ratio p.value
 High - Low     0.78 0.233 60   3.348  0.0014

nutrients = Enriched:
 contrast   estimate    SE df t.ratio p.value
 High - Low     3.56 0.233 60  15.301 <0.0001

As comparações seguem a ordem dos níveis: Control − Enriched e Low − High. Tukey aplica-se dentro de cada grupo, não entre todos os pares das quatro células.

O ajuste de Tukey é o predefinido para as comparações pairwise e aplica-se separadamente em cada grupo definido por |. Com apenas duas médias por grupo, há um único contraste: o ajuste não altera o seu valor-p. Não reunimos os contrastes de grupos diferentes numa única família.

Relatar o Crescimento

O efeito da adição de nutrientes no crescimento dos corais varia com a luminosidade, F(1, 60) = 71.43, p < .001. A adição de nutrientes favorece mais o crescimento sob alta do que sob baixa luminosidade: alta luminosidade, M_controlo = 4.94 mm, SE = 0.16; M_nutrientes = 8.07 mm, SE = 0.16; t(60) = −13.43, p < .001; baixa luminosidade, M_controlo = 4.16 mm, SE = 0.16; M_nutrientes = 4.51 mm, SE = 0.16; t(60) = −1.48, p = .144.

• Comentário: este exemplo usa dados simulados, não resultados empíricos.

Não Façam Isto

significativo aqui + não significativo ali ≠ interacção

O gráfico ajuda a perceber o padrão. O teste da interacção compara directamente as diferenças, em vez de comparar os seus valores-p.

Script dos Corais

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

library(car)
library(dplyr)
library(emmeans)
library(ggplot2)

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

ds <- read.csv("../../data/coral_growth_2x2_demo.csv")
m <- lm(growth_mm ~ light * nutrients, data = ds)
emm <- list(~ light, ~ nutrients, pairwise ~ nutrients | light, pairwise ~ light | nutrients) |>
  emmeans(object = m, specs = _)
g <- emmeans(m, ~ light * nutrients) |> data.frame() |> # means for the 4 cells
     ggplot(aes(x = light, y = emmean, group = nutrients, colour = nutrients, shape = nutrients)) +
     geom_line() + geom_point(size = 3) +
     geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = .08) +
     ylim(0, 10) +
     theme_classic()

# Print out results and plot
Anova(m, type = "III") |> print()
print(emm)
print(g)

Um modelo, o gráfico das quatro médias e as comparações da nossa pergunta.

Resultados dos Corais

Anova Table (Type III tests)

Response: growth_mm
                 Sum Sq Df  F value    Pr(>F)    
(Intercept)     1880.51  1 4332.454 < 2.2e-16 ***
light             75.48  1  173.885 < 2.2e-16 ***
nutrients         48.25  1  111.162 2.776e-15 ***
light:nutrients   31.01  1   71.434 8.338e-12 ***
Residuals         26.04 60                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

$`emmeans of light`
 light emmean    SE df lower.CL upper.CL
 High    6.51 0.116 60     6.27     6.74
 Low     4.33 0.116 60     4.10     4.57

Results are averaged over the levels of: nutrients 
Confidence level used: 0.95 

$`emmeans of nutrients`
 nutrients emmean    SE df lower.CL upper.CL
 Control     4.55 0.116 60     4.32     4.79
 Enriched    6.29 0.116 60     6.06     6.52

Results are averaged over the levels of: light 
Confidence level used: 0.95 

$`emmeans of nutrients | light`
light = High:
 nutrients emmean    SE df lower.CL upper.CL
 Control     4.94 0.165 60     4.61     5.27
 Enriched    8.07 0.165 60     7.74     8.40

light = Low:
 nutrients emmean    SE df lower.CL upper.CL
 Control     4.16 0.165 60     3.83     4.49
 Enriched    4.51 0.165 60     4.18     4.84

Confidence level used: 0.95 

$`pairwise differences of nutrients | light`
light = High:
 2                  estimate    SE df t.ratio p.value
 Control - Enriched   -3.129 0.233 60 -13.432 <0.0001

light = Low:
 2                  estimate    SE df t.ratio p.value
 Control - Enriched   -0.344 0.233 60  -1.479  0.1444


$`emmeans of light | nutrients`
nutrients = Control:
 light emmean    SE df lower.CL upper.CL
 High    4.94 0.165 60     4.61     5.27
 Low     4.16 0.165 60     3.83     4.49

nutrients = Enriched:
 light emmean    SE df lower.CL upper.CL
 High    8.07 0.165 60     7.74     8.40
 Low     4.51 0.165 60     4.18     4.84

Confidence level used: 0.95 

$`pairwise differences of light | nutrients`
nutrients = Control:
 2          estimate    SE df t.ratio p.value
 High - Low     0.78 0.233 60   3.348  0.0014

nutrients = Enriched:
 2          estimate    SE df t.ratio p.value
 High - Low     3.56 0.233 60  15.301 <0.0001

Dos Corais aos Pinguins

A diferença de massa entre sexos é igual em todas as espécies?

• 2 × 2 → 2 × 3: dois coeficientes para espécie e dois para sexo × espécie
• Dados observacionais: sem atribuição aleatória das condições

Preparar os Pinguins

library(dplyr)
library(emmeans)
library(palmerpenguins)
library(tidyr)

ds <- penguins |>
  select(sex, species, body_mass_g, flipper_length_mm) |>
  drop_na()
m_int <- lm(body_mass_g ~ sex * species, data = ds)
emm <- emmeans(m_int, ~ sex * species)
penguin_plot <- as.data.frame(emm)

• Uma linha = um espécime
• Repetir espécie ou sexo ≠ repetir o indivíduo
• Mesmos casos completos ao longo do exemplo

Ver Primeiro

Antes de resumir por sexo ou por espécie, descrevam as seis médias.

Codificação e Interpretação

coef(m_int)

  (Intercept)          sex1      species1      species2 sex1:species1 
   4173.84711    -315.25282    -467.68273    -440.75888     -22.07595 
sex1:species2 
    109.37047 

sex1: metade de female − male, com igual peso para cada espécie.

Os dois coeficientes de espécie descrevem desvios das médias marginais; as duas interacções descrevem como a diferença entre sexos varia por espécie.

Continuamos com soma-zero, definida no início. Com níveis female, male, sex1 não é a diferença inteira nem apenas a diferença em Adelie. No modelo barbatana × espécie mais adiante, o coeficiente de barbatana é a média dos três declives, com pesos iguais, e as interacções são desvios desse declive. Interpretação no tutorial.

A Interacção na Tabela

Anova(m_int, type = "III") |> print()

Anova Table (Type III tests)

Response: body_mass_g
                Sum Sq  Df   F value    Pr(>F)    
(Intercept) 5232595969   1 54661.828 < 2.2e-16 ***
sex           29851220   1   311.838 < 2.2e-16 ***
species      143001222   2   746.924 < 2.2e-16 ***
sex:species    1676557   2     8.757 0.0001973 ***
Residuals     31302628 327                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

sex:species: a diferença entre sexos é igual nas três espécies?

Aqui o F tem 2 gl no numerador: testa dois coeficientes em conjunto. Nenhum dos ts da tabela de coeficientes o substitui.

Comparações de Modelos

Efeito	Modelos Comparados
sex	Completo vs. interacção + espécie, sem sexo
species	Completo vs. interacção + sexo, sem espécie
sex:species	Completo vs. sexo + espécie, sem interacção
(Intercept)	Completo vs. média das seis células = 0

Esta tabela explica as hipóteses. Não precisamos de ajustar um modelo reduzido por linha no script.

Mantêm-se pesos iguais para os níveis do outro factor. São restrições no mesmo modelo completo, não passos de uma sequência. Com interacção, apagar sex ou species da fórmula não constrói necessariamente o teste Tipo III: o R pode voltar a representar o mesmo espaço através das colunas da interacção. A tabela substitui essa implementação, que não pertence ao percurso prático desta aula. O intercepto não responde à pergunta biológica. Tipo III.

Quatro Perguntas

Queremos…	Pedido a emmeans()
Sexo, em média	~ sex
Espécie, em média	pairwise ~ species
Sexo dentro de cada espécie	pairwise ~ sex | species
Espécie dentro de cada sexo	pairwise ~ species | sex

Nas duas primeiras, fazemos a média sobre o outro factor. Nas duas últimas, comparamos os níveis de um factor dentro de cada nível do outro.

Sexo Dentro de Espécie

emm <- emmeans(m_int, pairwise ~ sex | species)
print(emm)

$emmeans
species = Adelie:
 sex    emmean   SE  df lower.CL upper.CL
 female   3369 36.2 327     3298     3440
 male     4043 36.2 327     3972     4115

species = Chinstrap:
 sex    emmean   SE  df lower.CL upper.CL
 female   3527 53.1 327     3423     3632
 male     3939 53.1 327     3835     4043

species = Gentoo:
 sex    emmean   SE  df lower.CL upper.CL
 female   4680 40.6 327     4600     4760
 male     5485 39.6 327     5407     5563

Confidence level used: 0.95 

$contrasts
species = Adelie:
 contrast      estimate   SE  df t.ratio p.value
 female - male     -675 51.2 327 -13.174 <0.0001

species = Chinstrap:
 contrast      estimate   SE  df t.ratio p.value
 female - male     -412 75.0 327  -5.487 <0.0001

species = Gentoo:
 contrast      estimate   SE  df t.ratio p.value
 female - male     -805 56.7 327 -14.188 <0.0001

Uma comparação por espécie: fêmeas − machos. O teste global de interacção já avaliou se estas diferenças variam.

Outras Distâncias no Mesmo Gráfico

Agora a pergunta é espécie dentro de sexo. Não surgiram dados novos; mudou a comparação que queremos ler.

Espécie Dentro de Sexo

emm <- emmeans(m_int, pairwise ~ species | sex)
print(emm)

$emmeans
sex = female:
 species   emmean   SE  df lower.CL upper.CL
 Adelie      3369 36.2 327     3298     3440
 Chinstrap   3527 53.1 327     3423     3632
 Gentoo      4680 40.6 327     4600     4760

sex = male:
 species   emmean   SE  df lower.CL upper.CL
 Adelie      4043 36.2 327     3972     4115
 Chinstrap   3939 53.1 327     3835     4043
 Gentoo      5485 39.6 327     5407     5563

Confidence level used: 0.95 

$contrasts
sex = female:
 contrast           estimate   SE  df t.ratio p.value
 Adelie - Chinstrap     -158 64.2 327  -2.465  0.0377
 Adelie - Gentoo       -1311 54.4 327 -24.088 <0.0001
 Chinstrap - Gentoo    -1153 66.8 327 -17.246 <0.0001

sex = male:
 contrast           estimate   SE  df t.ratio p.value
 Adelie - Chinstrap      105 64.2 327   1.627  0.2357
 Adelie - Gentoo       -1441 53.7 327 -26.855 <0.0001
 Chinstrap - Gentoo    -1546 66.2 327 -23.345 <0.0001

P value adjustment: tukey method for comparing a family of 3 estimates 

São três comparações entre espécies em cada sexo. Tukey ajusta cada família de três; não a família de todos os pares das seis células.

E as Médias Marginais?

emm <- list(~ sex, pairwise ~ species,
            pairwise ~ sex | species, pairwise ~ species | sex) |>
  emmeans(object = m_int, specs = _)
print(emm)

$`emmeans of sex`
 sex    emmean   SE  df lower.CL upper.CL
 female   3859 25.3 327     3809     3908
 male     4489 25.2 327     4440     4539

Results are averaged over the levels of: species 
Confidence level used: 0.95 

$`emmeans of species`
 species   emmean   SE  df lower.CL upper.CL
 Adelie      3706 25.6 327     3656     3757
 Chinstrap   3733 37.5 327     3659     3807
 Gentoo      5082 28.4 327     5026     5138

Results are averaged over the levels of: sex 
Confidence level used: 0.95 

$`pairwise differences of species`
 1                  estimate   SE  df t.ratio p.value
 Adelie - Chinstrap    -26.9 45.4 327  -0.593  0.8241
 Adelie - Gentoo     -1376.1 38.2 327 -36.007 <0.0001
 Chinstrap - Gentoo  -1349.2 47.0 327 -28.682 <0.0001

Results are averaged over the levels of: sex 
P value adjustment: tukey method for comparing a family of 3 estimates 

$`emmeans of sex | species`
species = Adelie:
 sex    emmean   SE  df lower.CL upper.CL
 female   3369 36.2 327     3298     3440
 male     4043 36.2 327     3972     4115

species = Chinstrap:
 sex    emmean   SE  df lower.CL upper.CL
 female   3527 53.1 327     3423     3632
 male     3939 53.1 327     3835     4043

species = Gentoo:
 sex    emmean   SE  df lower.CL upper.CL
 female   4680 40.6 327     4600     4760
 male     5485 39.6 327     5407     5563

Confidence level used: 0.95 

$`pairwise differences of sex | species`
species = Adelie:
 2             estimate   SE  df t.ratio p.value
 female - male     -675 51.2 327 -13.174 <0.0001

species = Chinstrap:
 2             estimate   SE  df t.ratio p.value
 female - male     -412 75.0 327  -5.487 <0.0001

species = Gentoo:
 2             estimate   SE  df t.ratio p.value
 female - male     -805 56.7 327 -14.188 <0.0001


$`emmeans of species | sex`
sex = female:
 species   emmean   SE  df lower.CL upper.CL
 Adelie      3369 36.2 327     3298     3440
 Chinstrap   3527 53.1 327     3423     3632
 Gentoo      4680 40.6 327     4600     4760

sex = male:
 species   emmean   SE  df lower.CL upper.CL
 Adelie      4043 36.2 327     3972     4115
 Chinstrap   3939 53.1 327     3835     4043
 Gentoo      5485 39.6 327     5407     5563

Confidence level used: 0.95 

$`pairwise differences of species | sex`
sex = female:
 2                  estimate   SE  df t.ratio p.value
 Adelie - Chinstrap     -158 64.2 327  -2.465  0.0377
 Adelie - Gentoo       -1311 54.4 327 -24.088 <0.0001
 Chinstrap - Gentoo    -1153 66.8 327 -17.246 <0.0001

sex = male:
 2                  estimate   SE  df t.ratio p.value
 Adelie - Chinstrap      105 64.2 327   1.627  0.2357
 Adelie - Gentoo       -1441 53.7 327 -26.855 <0.0001
 Chinstrap - Gentoo    -1546 66.2 327 -23.345 <0.0001

P value adjustment: tukey method for comparing a family of 3 estimates 

Nas médias marginais, cada nível do outro factor recebe o mesmo peso. Para espécie, com três níveis, pedimos também as comparações de Tukey; para sexo, descrevemos as duas médias.

Estas médias não substituem o padrão dentro de cada espécie.

Relatar o Teste

A diferença média de massa corporal entre sexos varia em função da espécie, F(2, 327) = 8.76, p < .001.

• Comentário: o nome da análise e o valor-p não substituem a descrição do padrão.

Relatar Sexo em Cada Espécie

Os machos são mais pesados do que as fêmeas nas três espécies: Adelie, M_fêmeas = 3369 g, SE = 36; M_machos = 4043 g, SE = 36; t(327) = −13.17, p < .001; Chinstrap, M_fêmeas = 3527 g, SE = 53; M_machos = 3939 g, SE = 53; t(327) = −5.49, p < .001; Gentoo, M_fêmeas = 4680 g, SE = 41; M_machos = 5485 g, SE = 40; t(327) = −14.19, p < .001.

O Padrão Entre Espécies

As fêmeas Adelie são mais leves do que as Chinstrap, M_Adelie = 3369 g, SE = 36; M_Chinstrap = 3527 g, SE = 53; t(327) = −2.47, p = .038. Já nos machos, não encontramos evidência estatisticamente significativa de que os Adelie sejam mais pesados do que os Chinstrap, M_Adelie = 4043 g, SE = 36; M_Chinstrap = 3939 g, SE = 53; t(327) = 1.63, p = .236.

• Comentário: a interacção é testada directamente, não pela diferença entre estes valores-p. Os dados são observacionais; o alcance da generalização depende da amostragem.

Os valores-p têm ajuste de Tukey entre espécies em cada sexo.

Script Factorial

library(car)
library(emmeans)
library(ggplot2)
library(palmerpenguins)
library(tidyr)

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

ds <- penguins |> drop_na(sex, species, body_mass_g, flipper_length_mm)
model <- lm(body_mass_g ~ sex * species, data = ds)
emm <- list(~ sex, pairwise ~ species,
            pairwise ~ sex | species, pairwise ~ species | sex) |>
  emmeans(object = model, specs = _)
g <- emmeans(model, ~ sex * species) |> data.frame() |>
  ggplot(aes(species, emmean, colour = sex, shape = sex, group = sex)) +
  geom_line() + geom_point() +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = .08) +
  scale_shape_manual(values = c(female = 15, male = 17))

# Print results and plot
Anova(model, type = "III") |> print()
print(emm)
print(g)

Um modelo, o gráfico das seis médias e as quatro perguntas no mesmo emm.

Resultados Factoriais

Anova Table (Type III tests)

Response: body_mass_g
                Sum Sq  Df   F value    Pr(>F)    
(Intercept) 5232595969   1 54661.828 < 2.2e-16 ***
sex           29851220   1   311.838 < 2.2e-16 ***
species      143001222   2   746.924 < 2.2e-16 ***
sex:species    1676557   2     8.757 0.0001973 ***
Residuals     31302628 327                        
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

$`emmeans of sex`
 sex    emmean   SE  df lower.CL upper.CL
 female   3859 25.3 327     3809     3908
 male     4489 25.2 327     4440     4539

Results are averaged over the levels of: species 
Confidence level used: 0.95 

$`emmeans of species`
 species   emmean   SE  df lower.CL upper.CL
 Adelie      3706 25.6 327     3656     3757
 Chinstrap   3733 37.5 327     3659     3807
 Gentoo      5082 28.4 327     5026     5138

Results are averaged over the levels of: sex 
Confidence level used: 0.95 

$`pairwise differences of species`
 1                  estimate   SE  df t.ratio p.value
 Adelie - Chinstrap    -26.9 45.4 327  -0.593  0.8241
 Adelie - Gentoo     -1376.1 38.2 327 -36.007 <0.0001
 Chinstrap - Gentoo  -1349.2 47.0 327 -28.682 <0.0001

Results are averaged over the levels of: sex 
P value adjustment: tukey method for comparing a family of 3 estimates 

$`emmeans of sex | species`
species = Adelie:
 sex    emmean   SE  df lower.CL upper.CL
 female   3369 36.2 327     3298     3440
 male     4043 36.2 327     3972     4115

species = Chinstrap:
 sex    emmean   SE  df lower.CL upper.CL
 female   3527 53.1 327     3423     3632
 male     3939 53.1 327     3835     4043

species = Gentoo:
 sex    emmean   SE  df lower.CL upper.CL
 female   4680 40.6 327     4600     4760
 male     5485 39.6 327     5407     5563

Confidence level used: 0.95 

$`pairwise differences of sex | species`
species = Adelie:
 2             estimate   SE  df t.ratio p.value
 female - male     -675 51.2 327 -13.174 <0.0001

species = Chinstrap:
 2             estimate   SE  df t.ratio p.value
 female - male     -412 75.0 327  -5.487 <0.0001

species = Gentoo:
 2             estimate   SE  df t.ratio p.value
 female - male     -805 56.7 327 -14.188 <0.0001


$`emmeans of species | sex`
sex = female:
 species   emmean   SE  df lower.CL upper.CL
 Adelie      3369 36.2 327     3298     3440
 Chinstrap   3527 53.1 327     3423     3632
 Gentoo      4680 40.6 327     4600     4760

sex = male:
 species   emmean   SE  df lower.CL upper.CL
 Adelie      4043 36.2 327     3972     4115
 Chinstrap   3939 53.1 327     3835     4043
 Gentoo      5485 39.6 327     5407     5563

Confidence level used: 0.95 

$`pairwise differences of species | sex`
sex = female:
 2                  estimate   SE  df t.ratio p.value
 Adelie - Chinstrap     -158 64.2 327  -2.465  0.0377
 Adelie - Gentoo       -1311 54.4 327 -24.088 <0.0001
 Chinstrap - Gentoo    -1153 66.8 327 -17.246 <0.0001

sex = male:
 2                  estimate   SE  df t.ratio p.value
 Adelie - Chinstrap      105 64.2 327   1.627  0.2357
 Adelie - Gentoo       -1441 53.7 327 -26.855 <0.0001
 Chinstrap - Gentoo    -1546 66.2 327 -23.345 <0.0001

P value adjustment: tukey method for comparing a family of 3 estimates 

Agora um Declive

A relação massa–barbatana difere entre espécies?

Antes do código, desenhem o padrão que a interacção significaria.

Ver os Declives

A interacção pergunta se os declives diferem, não se um deles tem p < .05.

Ajustar os Declives

m_slope <- lm(body_mass_g ~ flipper_length_mm * species, data = ds)

O modelo estima um declive por espécie.

Comparar Declives

Demonstração: rectas paralelas → declives iguais. Teste explicado: flipper_length_mm:species.

m_parallel <- lm(body_mass_g ~ flipper_length_mm + species, data = ds)
anova(m_parallel, m_slope)

Analysis of Variance Table

Model 1: body_mass_g ~ flipper_length_mm + species
Model 2: body_mass_g ~ flipper_length_mm * species
  Res.Df      RSS Df Sum of Sq     F   Pr(>F)   
1    329 45843144                               
2    327 44391669  2   1451475 5.346 0.005193 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Para perceber o teste, não uma etapa extra do script final.

O Teste na Tabela

Anova(m_slope, type = "III") |> print()

Anova Table (Type III tests)

Response: body_mass_g
                            Sum Sq  Df F value    Pr(>F)    
(Intercept)                5663696   1  41.720 3.814e-10 ***
flipper_length_mm         22649648   1 166.843 < 2.2e-16 ***
species                    1229661   2   4.529  0.011477 *  
flipper_length_mm:species  1451475   2   5.346  0.005193 ** 
Residuals                 44391669 327                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

flipper_length_mm:species: todos os declives são iguais? É o mesmo teste da demonstração, com 2 gl; no script basta esta tabela.

A linha species compara espécies a comprimento zero, fora dos dados; não responde à pergunta sobre diferenças entre declives. Centrar a covariável muda essa referência dos efeitos de ordem inferior, não o teste da interacção.

Estimar os Declives

emm <- emtrends(m_slope, pairwise ~ species, var = "flipper_length_mm")
print(emm)

$emtrends
 species   flipper_length_mm.trend   SE  df lower.CL upper.CL
 Adelie                       32.7 4.69 327     23.5     41.9
 Chinstrap                    34.6 6.31 327     22.2     47.0
 Gentoo                       54.2 5.15 327     44.0     64.3

Confidence level used: 0.95 

$contrasts
 contrast           estimate   SE  df t.ratio p.value
 Adelie - Chinstrap    -1.88 7.86 327  -0.240  0.9688
 Adelie - Gentoo      -21.48 6.97 327  -3.083  0.0063
 Chinstrap - Gentoo   -19.59 8.15 327  -2.405  0.0440

P value adjustment: tukey method for comparing a family of 3 estimates 

emtrends() estima os declives em g/mm; pairwise pede também as suas comparações, com ajuste de Tukey. Não comparar valores-p separados para decidir se os declives diferem.

Relatar os Declives

O declive da relação entre comprimento da barbatana e massa corporal varia entre espécies, F(2, 327) = 5.35, p = .005. O declive é maior em Gentoo do que em Adelie, B_Adelie = 32.69 g/mm, SE = 4.69; B_Chinstrap = 34.57 g/mm, SE = 6.31; B_Gentoo = 54.17 g/mm, SE = 5.15; t_{Adelie vs. Gentoo}(327) = −3.08, p = .006.

• Comentário: interpretamos os declives como associações observacionais, não como efeitos causais; não extrapolamos para além dos comprimentos amostrados.

As comparações entre declives têm ajuste de Tukey.

Script dos Declives

Um modelo com interacção: a tabela testa os termos; os contrastes estimam as diferenças entre declives. Sem repetir a demonstração com m_parallel.

library(car)
library(emmeans)
library(ggplot2)
library(palmerpenguins)
library(tidyr)

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

ds <- penguins |> drop_na(sex, species, body_mass_g, flipper_length_mm)
m_slope <- lm(body_mass_g ~ flipper_length_mm * species, data = ds)
emm <- emtrends(m_slope, pairwise ~ species, var = "flipper_length_mm")
g <- ggplot(ds, aes(flipper_length_mm, body_mass_g,
                    colour = species, shape = species)) +
  geom_point() + geom_smooth(method = "lm") +
  scale_shape_manual(values = c(Adelie = 15, Chinstrap = 17, Gentoo = 16))

# Print results and plot
Anova(m_slope, type = "III") |> print()
print(emm)
print(g)

Resultados dos Declives

Anova Table (Type III tests)

Response: body_mass_g
                            Sum Sq  Df F value    Pr(>F)    
(Intercept)                5663696   1  41.720 3.814e-10 ***
flipper_length_mm         22649648   1 166.843 < 2.2e-16 ***
species                    1229661   2   4.529  0.011477 *  
flipper_length_mm:species  1451475   2   5.346  0.005193 ** 
Residuals                 44391669 327                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

$emtrends
 species   flipper_length_mm.trend   SE  df lower.CL upper.CL
 Adelie                       32.7 4.69 327     23.5     41.9
 Chinstrap                    34.6 6.31 327     22.2     47.0
 Gentoo                       54.2 5.15 327     44.0     64.3

Confidence level used: 0.95 

$contrasts
 contrast           estimate   SE  df t.ratio p.value
 Adelie - Chinstrap    -1.88 7.86 327  -0.240  0.9688
 Adelie - Gentoo      -21.48 6.97 327  -3.083  0.0063
 Chinstrap - Gentoo   -19.59 8.15 327  -2.405  0.0440

P value adjustment: tukey method for comparing a family of 3 estimates 

Agora: Modelo

20 minutos · grupos de 3 ou 5

Modelos Lineares no R — Exercícios

Façam exercício 4.

Agora: Restrição

10 minutos · grupos de 3 ou 5

Comparação de Modelos — Exercícios

Façam exercício 6 para explicar a hipótese por trás do teste, não como mais uma etapa do script de análise.

Comparar Antes de Calcular

10 minutos · grupos de 3 ou 5

• Dois gráficos de corais: padrão aditivo de luminosidade e nutrientes vs. interacção
• Assinalem as distâncias comparadas pelo teste de interacção
• Troquem os factores nos eixos: mudou o padrão ou só a maneira de o ler?

Se Sobrar Tempo

Modelos Lineares no R — Exercícios

Façam exercício 5.

Se Precisarem de Teoria

Comparação de Modelos — Interacções e Contrastes

O capítulo explica a diferença de diferenças, os testes conjuntos e as restrições por trás de cada linha da tabela.

A Seguir

Antes de generalizar o modelo, precisamos de ver onde ele falha.