Modelos com Interacções

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

João O. Santos
joao.santos@ispa.pt

ISPA

2025-10-28

Two Categorical IVs

  • What if we have two categorical IVs

  • We can code each one with the same coding scheme

  • Then we can test for the main effect of each

  • But, now we’re going to test for an interaction

Modelos com Múltiplos Preditores

  • Já vimos como podemos querer comparar a média de mais que dois grupos

  • Ou seja, ver o efeito duma variável categórica numa variável dependente quantitativa

  • Na regressão múltipla vimos como podemos testar o efeito de várias variáveis independentes numa única VD

  • Vimos que ao incluir mais que uma VI no modelo vamos poder testar se cada uma tem efeito (na VD), mesmo quando controlamos estatisticamente para o efeito das restantes VIs

  • Contudo, ainda não testámos se há efeitos de interacção/moderação

Interacção/Moderação

  • Quando efeito duma VI na VD depende de outra VI

  • Exemplos:

    • O efeito duma intervenção tem efeitos opostos em jovens adultos e em idosos
    • Um grupo tem um desempenho melhor na leitura do que na escrita, enquanto outro grupo tem melhor desempenho na escrita do que na leitura

Delineamentos Fatoriais

  • Já vimos delineamentos com duas ou mais amostras independentes (e.g., placebo vs. controlo; tratamento A vs. B vs. placebo)

  • Já vimos delineamentos com duas ou mais medidas repetidas (e.g., pré vs. pós; pré vs. pós vs. follow-up)

  • Mas ainda não vimos delineamentos com mais que uma variável independente categórica

Delineamentos Fatoriais (Incluindo Mistos)

  • Imaginemos que eu quero ver o efeito da faixa etária (jovem adulto: 18 <= idade < 25; adulto de meia: 25 <= idade < 65) e do bloco da prova de Stroop (congruente; incongruente)

  • Podemos testar se a média dos tempos de resposta é diferente consoante a faixa etária (controlando para o efeito do bloco)—efeito principal da idade

  • Podemos testar se a média dos tempos de resposta é diferente consoante o bloco (controlando para o efeito da idade)—efeito principal do bloco

  • Mas podemos querer testar se a diferença das média dos tempos de respota em cada bloco são diferentes consoante a faixa etária—efeito de interacção

Delineamentos Fatoriais

  • Quando nos referimos a delineamentos fatoriais é frequente usarmos notações como:

    • 2X2 (uma VI com dois níveis e outra VI com 2 níveis)
    • 3X4 (uma VI com três níveis e outra VI com 4)

  • Falamos também em células do delineamento porque podemos construir tabelas com os níveis duma VI nas colunas e os de outra nas linhas

  • Podemos ter delineamentos com mais de duas VIs, tendo interacções de ordem superior (e.g., triplas)

    • Mas queremos evitar ter delineamentos demasiado complexos

Delineamentos Fatoriais Mistos

  • Podemos ter delineamentos fatoriais só com medições repetidas—todas as variáveis manipuladas intra-pp

  • Podemos ter delineamentos fatoriais só com amostras independentes—todas as variáveis manipuladas entre-pp

  • Podemos ter delineamentos fatoriais mistos—pelo menos uma VI intra-pp e uma entre-pp

Delineamentos Fatoriais

  • Vantagens:
    • Permitem testar efeitos de interacção
    • Permitem testar o efeito duma VI controlando para a/s outra/s
    • Podemos então responder a questões de investigações mais complexas
  • Desvantagens:
    • Acrescentam alguma complexidade
    • Tendem a aumentar o número de células do delineamento aumentando o N de participantes e/ou medidas repetidas necessário para manter o poder de teste

Delineamentos Mistos

  • É difícil pensar em vantagens e desvantagens genéricas

  • Quanto mais experiência tivermos a planear, conduzir e analisar dados de estudos, mais capazes seremos de pesar os prós e os contras de vários delineamentos e quais as variáveis que devemos manipular intra-pp e quais as que devemos manipular entre-pp

Tipo de Soma de Quadrados [Curiosidade]

  • Existem várias formas de fazer as somas de quadrados nas tabelas da ANOVA
    • Cada tipo de soma de quadrados compara modelos diferentes e/ou de forma diferente
  • Geralmente falando para delineamentos fatoriais (ortogonais) queremos o Tipo II ou III (ver Hector & collaborators, 2010)
    • Tradicionalmente nas ciências sociais recorre-se ao Tipo III
  • O JASP e o SPSS usam o Tipo III
    • O R usa o Tipo I

Padrões de Interacção

  • Há vários padrões possíveis para interacções com variáveis quantitativas

  • A VI1 pode ter efeito na VD só num nível da VI2, mas não no outro (e.g., o stress só afectar o desempenho dos homens, mas não das mulheres)

  • A VI1 pode ter um efeito na VD numa direcção num nível da VI2 e ter o efeito oposto no outro nível (e.g., o stress afectar negativamente os homens e positivamente as mulheres)

  • A VI1 pode ter um efeito na VD na mesma direcção em cada nível da VI2, mas o efeito ter uma magnitude diferente em cada nível (e.g., o stress piorar o desempenho de tanto homens como mulheres, mas ser particularmente nocivo nos homens).

Interpretação de Padrões

  • Vejamos agora a representação de vários padrões de resultados possíveis, com e sem interações e efeitos principais…

  • Vamos ver todas as combinações de efeitos significativos e não signficativos na tabela duma anova fatorial com dois fatores:

    • Apenas efeito da VI1
    • Apenas efeito da VI2
    • Apenas efeito de interacção
    • Efeito da VI1 e de interacção
    • Efeito da VI2 e de interacção
    • Todos os efeitos significativos

Apenas Efeito Principal de Condição

library(ggplot2)

mk_graph <- function(ds, x = NULL, line = NULL) {
    ggplot(ds, aes(y = EMM, x = .data[[x]], group = .data[[line]],
                   shape = .data[[line]], color = .data[[line]])) +
    geom_line(linewidth = 1.5, alpha = 0.85, show.legend = FALSE) +
    geom_point(size = 5, alpha = 0.85) +
    ylim(0, 10) +
    ggtitle(paste0(x, " by ", line)) +
    theme_classic() +
    theme(text = element_text(family = "Times New Roman", size = 16),
          plot.title = element_text(hjust = 0.5))
}
emms <- data.frame(Condition = c(rep("Control", 2), rep("Experimental", 2)),
                   Task = rep(c("Reading", "Writing"), 2),
                   EMM = c(4, 4.1, 8, 8.1))
mk_graph(emms, "Condition", "Task")

mk_graph(emms, "Task", "Condition")

  • Sem efeito principal de Tarefa
  • Sem efeito de Interacção
  • Grupo de controlo mais rápido que o experimental independentemente da tarefa

Apenas Efeito Principal de Tarefa

emms <- data.frame(Condition = c(rep("Control", 2), rep("Experimental", 2)),
                   Task = rep(c("Reading", "Writing"), 2),
                   EMM = c(4, 8, 4.1, 8.1))
mk_graph(emms, "Condition", "Task")

mk_graph(emms, "Task", "Condition")

  • Sem efeito principal de Condição
  • Sem efeito de Interacção
  • Leitura mais rápida que escrita independentemente do grupo

Apenas Efeito de Interacção

emms <- data.frame(Condition = c(rep("Control", 2), rep("Experimental", 2)),
                   Task = rep(c("Reading", "Writing"), 2),
                   EMM = c(8, 4, 4, 8))
mk_graph(emms, "Condition", "Task")

mk_graph(emms, "Task", "Condition")

  • Sem efeito principal de Condição
  • Sem efeito principal de Tarefa
  • Há diferenças entre tarefas em cada grupo e entre grupos em cada tarefa
  • Controlo rápido na escrita e lento na leitura. Experimental rápido na leitura e lento na escrita. Experimental demorou na leitura o que o de controlo demorou na escrita e na escrita o que o de controlo demorou na leitura.

Efeito Principal de Condição e Efeito de Interacção

emms <- data.frame(Condition = c(rep("Control", 2), rep("Experimental", 2)),
                   Task = rep(c("Reading", "Writing"), 2),
                   EMM = c(2, 5, 8, 5))
mk_graph(emms, "Condition", "Task")

mk_graph(emms, "Task", "Condition")

  • Sem efeito principal de Tarefa
  • Há diferenças entre condiçoes só na leitura, não há na escrita
  • O grupo de controlo é tão rápido a escrever como o experimental, mas é mais rápido a ler do que experimental

Efeito Principal de Tarefa e Efeito de Interacção

emms <- data.frame(Condition = c(rep("Control", 2), rep("Experimental", 2)),
                   Task = rep(c("Reading", "Writing"), 2),
                   EMM = c(2, 8, 5, 5))
mk_graph(emms, "Condition", "Task")

mk_graph(emms, "Task", "Condition")

  • Sem efeito principal de Condição
  • Há diferenças entre condições no grupo de controlo mas não no experimental
  • O grupo de controlo é rápido a ler e lento a escrever, o grupo experimental demora o mesmo tempo em ambas as tarefas, sendo mais lento a ler que o do controlo, mas mais rápido a escrever que eles

Todos os Efeitos Significativos

emms <- data.frame(Condition = c(rep("Control", 2), rep("Experimental", 2)),
                   Task = rep(c("Reading", "Writing"), 2),
                   EMM = c(2, 8, 2, 2))
mk_graph(emms, "Condition", "Task")

mk_graph(emms, "Task", "Condition")

  • Há diferenças significativas entre condições apenas na escrita, há diferenças entre tarefas apenas no controlo
  • O grupo de controlo é mais rápido a ler que a escrever, enquanto o grupo experimental demora o mesmo tempo em ambos, sendo o tempo que demora em ambas as tarefas igual ao tempo que o grupo de controlo demora na sua tarefa mais rápida (a leitura)

Infinitas Possibilidades

  • Vimos todas as combinações possíveis de efeitos numa tabela da ANOVA com dois fatores

  • Acabámos de ver exemplos com um delineamento 2X2—ambos os fatores só com dois níveis—o delineamento fatorial mais simples

    • Mas com os pinguins já temos um exemplo dum delineamento 3X2

  • Há um número infinito de padrões possíveis, mas já ficam a conhecer uns quantos…

  • Pratiquem agora usando as bases de dados da “Data Library”->“ANOVA” do JASP que tenham mais que uma VI e/ou as do “Ecampus”->“Bases de Dados”->“Delineamentos Fatoriais”

  • Releiam artigos que conheçam com delineamentos fatoriais, olhando com atenção para a secção de resultados para:

    • Tentarem compreender o padrão
    • Tentarem perceber como os autores reportaram os resultados
      • Conseguem pensar em alterações que fariam ao texto para o tornar mais claro e/ou conciso?

Será que o Peso Depende do Sexo e da Espécie?

Proposed Model

library(car)
library(palmerpenguins)

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

ds <- penguins
ds <- subset(ds, !is.na(sex))

# Long form
model <- lm(body_mass_g ~ sex + species + species:sex, ds)
# Short form
model <- lm(body_mass_g ~ sex * species, ds)

ANOVA Tables

  • There are different ways we can compute the model comparisons

  • We can compare the models sequentially (Type I)

  • We can omit the interaction in models testing for main effects (Type

  • We can include all effects (even the interaction) in all comparsions, except the one where we test for the interaction (Type III)

ANOVA Tables

  • Type I is usually not what we want, because the order of the variables in the model matters

  • There is some debate between Type II and Type III

  • The default in R is type I

    • so watch out!

  • SPSS, JASP, etc, use Type III

ANOVA Tables

  • I prefer type III, I’m ok with type II, I would avoid type I unless you know what you’re doing and how to explain it

  • Type II vs III should only matter for really unbalanced designs with interactions

  • For more information see Hector and collaborators (2010)

ANOVA Tables

The car package allows us to get an type II and III sums of squares ANOVA tables:

library(car)

Anova(model, type = 3)
Sum Sq Df F value Pr(>F)
(Intercept) 5232595969 1 54661.827958 0.0000000
sex 29851220 1 311.838003 0.0000000
species 143001222 2 746.924492 0.0000000
sex:species 1676557 2 8.756997 0.0001973
Residuals 31302628 327 NA NA

Complete Script

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

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

# Get data and tidy it
ds <- penguins
ds <- subset(ds, !is.na(sex))

# Fit model
model <- lm(body_mass_g ~ sex * species, ds)

# Type III sums of squares ANOVA table
results <- Anova(model, type = 3)
# Estimated marginal means for decomposing the main effect of sex
emm_sex <- emmeans(model, ~ sex)
# Estimated marginal means for decomposing the main effect of species
emm_species <- emmeans(model, ~ species)
# Estimated marginal means for decomposing the interaction (2 X 3 = 6)
emm_int <- emmeans(model, ~ sex * species)
# Pairwise comparisons for the main effect of species
pc_species <- pairs(emm_species)
# Pairwise comparisons for one IV by the levels of the other IV
pc_sex_by_species <- pairs(emmeans(model, ~ sex | species))
pc_species_by_sex <- pairs(emmeans(model, ~ species | sex))

# Plot the emmeans
g_data <- data.frame(emm_int)
g_species_sex <- ggplot(g_data, aes(x = sex, y = emmean, color = species)) +
                 geom_point() + geom_line(aes(group = species)) +
                 geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
                               width = 0.25) +
                 theme_classic()
g_sex_species <- ggplot(g_data, aes(x = species, y = emmean, color = sex)) +
                 geom_point() + geom_line(aes(group = sex)) +
                 geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
                               width = 0.25) +
                 theme_classic()

Exercise

  1. Replicate the analysis on your machine

  2. Interpret the results

  3. Report the results

Solutions: 1

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

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

# Get data and tidy it
ds <- penguins
ds <- subset(ds, !is.na(sex))

# Fit model
model <- lm(body_mass_g ~ sex * species, ds)

# Type III sums of squares ANOVA table
results <- Anova(model, type = 3)
# Estimated marginal means for decomposing the main effect of sex
emm_sex <- emmeans(model, ~ sex)
# Estimated marginal means for decomposing the main effect of species
emm_species <- emmeans(model, ~ species)
# Estimated marginal means for decomposing the interaction (2 X 3 = 6)
emm_int <- emmeans(model, ~ sex * species)
# Pairwise comparisons for the main effect of species
pc_species <- pairs(emm_species)
# Pairwise comparisons for one IV by the levels of the other IV
pc_sex_by_species <- pairs(emmeans(model, ~ sex | species))
pc_species_by_sex <- pairs(emmeans(model, ~ species | sex))

# Plot the emmeans
g_data <- data.frame(emm_int)
g_species_sex <- ggplot(g_data, aes(x = sex, y = emmean, color = species)) +
                 geom_point() + geom_line(aes(group = species)) +
                 geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
                               width = 0.25) +
                 theme_classic()
g_sex_spcies <- ggplot(g_data, aes(x = species, y = emmean, color = sex)) +
                geom_point() + geom_line(aes(group = sex)) +
                geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
                              width = 0.25) +
                theme_classic()

Solutions: 2

  • Males are heavier than females across species

  • When we compare males with males, and females with females, there is significant evidence that Gentoos are heavier than Adelies or Chinstraps.

  • There is significant evidence that female Chinstraps are heavier than female Adelies.

  • There is no significant difference between the body mass of male Chinstraps and male Adelies (if anything male Adelies are heavier than male Chinstraps)

Going One Step Further

Looks like the difference between sexes is a smaller for Chinstraps than the other species

print(g_sex_species)

Going One Step Further

emmeans allows us to do a clever test for this

# Test for differences in the size of differences (advanced and complex)
pairs(pc_sex_by_species, by = NULL)
 contrast                                           estimate   SE  df t.ratio
 (female - male Adelie) - (female - male Chinstrap)     -263 90.8 327  -2.894
 (female - male Adelie) - (female - male Gentoo)         130 76.4 327   1.706
 (female - male Chinstrap) - (female - male Gentoo)      393 94.1 327   4.181
 p.value
  0.0113
  0.2043
  0.0001

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

Going One Step Further

  • Reading the results is tricky because the difference between females and males is negative…

  • Still, as we suspected, we see the size of the difference between males and females is smaller for Chinstraps than for Adelies and Gentoos

  • There is no significant evidence that the size of the differences between males and females differs between Chinstraps and Gentoos

Going One Step Further

  • The above test is not something most people think of, nor something I expect you to do

  • I just wanted to show you how you can test interesting hypothesis when you apply statistics in a clever way, and use the full power of R and emmeans to do it

Solutions: 3

  • Brace yourselves it’s going to feel like a lot…

  • Remember that the way to eat an elephant is one piece at a time…

Solutions: 3

To test for differences in body mass between sexes and species, we computed a 2X3 ANOVA (sex: male and female; species: Adelie, Chinstrap, Gentoo).

Solutions: 3

Results showed a significant interaction between sex and species, F2, 327 = 8.76, P < .001. Planned comparisons, with Tukey’s correction for multiple comparisons applied, show significant evidence that female Gentoos (M = 4679.74, SE = 40.63) are heavier than female Chinstraps (M = 3527.21, SE = 53.06; t = 17.25, P < .001), and than female Adelies (M = 3368.84, SE = 36.21; t = 24.09, P < .001). There is also significant evidence that female Chinstraps are heavier than female Adelies, t = 2.47, P = .038.

Solutions: 3

In males, we also find significant evidence of Gentoos (M = 5484.84, SE = 39.61) being heavier than both Chinstraps (M = 3938.97, SE = 53.06;t = 23.35,P < .001), and Adelies (M = 4043.49, SE = 36.21; t = 26.85, P < .001). However, we find no significant differences between the body mass of male Adelies and male Chinstraps, t = 1.63, P = .236. Unsurprisingly, multiple comparisons of sex across species show that, in all species, males are heavier than females, all Ps < .001 (with Tukey’s correction for multiple comparisons applied).

Solutions: 3

We also ran planned comparisons (with Tukey’s correction) testing for differences across species of the size of the difference in body mass between males and females. Results showed, that the difference between males and females is smaller for Chinstraps than for Adelies (t = -2.89, P = .011) and Gentoos (t = -4.180, P < .001). Yet, there is no evidence of difference in the magnitude of differences between sexes when we compare Adelies with Gentoos, t = 1.710, P = .204.

  • This is the analysis I didn’t expect you to include

Solutions: 3

The ANOVA also yielded a significant main effect of species, F1, 327 = 746.92, P < .001. When we look at the pairwise comparisons of the estimated marginal means for each species (with both sexes aggregated) we see significant evidence of Gentoos being heavier (M = 5082.29, SE = 28.37) than both Chinstraps (M = 3733.09, SE = 37.52; t = 28.68, P < .001) and Adelies (M = 3706.16, SE = 25.61; t = 36.01 P < .001). Yet, there is no significant evidence of differences between Chinstraps and Adelies, t < 1.

Solutions: 3

Finally, we see an unsurprising main effect of sex, F1, 327 = 311.84, P < .001. This merely tells us that, when we aggregate across species, males (M = 4489.10, SE = 25.16) are heavier than females (M = 3858.59, SE = 25.34)

Solutions: 3

In summary, results show that male penguins are heavier than females, regardless of species; with this difference being smaller in Chinstraps than the other species. Further, Gentoos are heavier than both Adelies and Chinstraps, regardless of sex. However, female Chinstraps are heavier than female Adelies, but male Chinstraps are not heavier than female Adelies (if anything they are lighter, but the difference is not significant).

  • The part in italic refers to the analysis I didn’t expect you to include

Two Quantitative IVs

  • We already looked at multiple quantitative IVs
    • But, we didn’t test for interactions then
  • Now we’re gonna test for interactions between two quantitative variables!

  • Link if slide is bronken

Lobsters

  • Do the number of traps and license holders predict the Tonnes of lobsters fished in Maine?

  • Is there an interaction between the two?

Lobsters

ds <- read.delim("../../data/lobsters-2016.tsv")

model <- lm(Tonnes ~ Traps.M. * license.holders, ds)

summary(model)

Call:
lm(formula = Tonnes ~ Traps.M. * license.holders, data = ds)

Residuals:
     Min       1Q   Median       3Q      Max 
-17746.5  -4741.8    -20.3   3122.0  20247.4 

Coefficients:
                           Estimate Std. Error t value Pr(>|t|)    
(Intercept)              -31459.246  20033.487  -1.570 0.121347    
Traps.M.                  44161.052  11385.640   3.879 0.000254 ***
license.holders               5.088      3.334   1.526 0.131982    
Traps.M.:license.holders     -5.256      1.846  -2.848 0.005942 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 8529 on 63 degrees of freedom
Multiple R-squared:  0.6745,    Adjusted R-squared:  0.659 
F-statistic: 43.52 on 3 and 63 DF,  p-value: 2.337e-15

Scaling

  • To understand interactions with quantitative variables is better to scale them (i.e., M = 0, SD = 1)

  • Then, we can refit the model with the scaled variables to understand what’s happening

ds$s_Tonnes <- scale(ds$Tonnes)
ds$s_Traps <- scale(ds$Traps.M.)
ds$s_License <- scale(ds$license.holders)

model <- lm(s_Tonnes ~ s_Traps * s_License, ds)

Visualizing the Interaction

# Dictomize variables
ds$d_T <- ifelse(ds$Traps.M. > mean(ds$Traps.M.), "> M", "<= M")
ds$d_L <- ifelse(ds$license.holders > mean(ds$license.holders), "> M", "<= M")

g_TL <- ggplot(ds, aes(x = Traps.M., y = Tonnes, color = d_L)) +
        geom_point() + geom_smooth(method = "lm") +
        theme_classic()

g_LT <- ggplot(ds, aes(x = license.holders, y = Tonnes, color = d_T)) +
        geom_point() + geom_smooth(method = "lm") +
        theme_classic()

Visualizing the Interaction

print(g_TL)

Visualizing the Interaction

print(g_LT)