Modelos Lineares no R

Modelos Lineares
Estatística Inferencial
emmeans
Tutorial
Do t-teste à ANCOVA: como fazer (quase) tudo com apenas uma função (lm).
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 Introdução

Já vimos a teoria: Dados = Modelo + Erro.

Agora vamos à prática.

A ideia deste capítulo é construir uma escada de complexidade, começando pelo modelo mais simples possível e acabando em interacções e covariáveis. O “segredo” que muitos livros de estatística clássica não deixam claro é este: são tudo regressões lineares. Na maior parte dos casos, a função lm() (Linear Model) chega.

ImportanteNota sobre medidas repetidas

Este tutorial cobre apenas designs entre-sujeitos.

Se mediu a mesma pessoa várias vezes (e.g., Antes vs. Depois), pare aqui: um lm() simples viola o pressuposto da independência.

Para isso vai precisar de ANOVAs de medições repetidas ou de modelos mistos (linear mixed models), que aparecem no capítulo seguinte.

Preparação

Vamos carregar os pacotes e definir algumas opções.

library(car)      # Anova() com SS Tipo III
library(emmeans)  # Médias marginais e contrastes
library(ggplot2)

# Em designs factorais, esta opção evita surpresas ao usar SS Tipo III
options(contrasts = c("contr.sum", "contr.poly"))
  • Em delineamentos factorais, a combinação SS Tipo III + codificação por omissão ou dummy (contr.treatment) pode produzir resultados difíceis de interpretar (e por vezes indesejados, dependendo do desenho e do desequilíbrio amostral).

  • A solução pragmática, muito comum em psicologia experimental quando se quer uma tabela ANOVA Tipo III, é usar sum-to-zero contrasts (contr.sum) que são ortogonais.

Se quiser saber mais:

Quando tiver dúvidas sobre comparações/contrastes, esta página costuma resolver muitos casos:

1 Nível 1: Teste-t para Uma Amostra

Aqui a lógica é: prever uma variável usando apenas um intercepto (a média).

\(Y_i = \beta_0 + e_i\)

Testar Contra 0

Vamos começar com um dataset que tem um zero que faz sentido: indiferença.

# Ler dados (caminho relativo para funcionar localmente e em CI/CD)
taylor <- read.csv("../data/Love4Taylor.csv")

# Espreitar
head(taylor)
ID Love4Taylor
pp01 1.50
pp02 2.15
pp03 1.32
pp04 1.17
pp05 -0.28
pp06 1.72
Nota

A ideia: 0 = indiferença, > 0 = gosto, < 0 = não gosto.

Ou seja: aqui testar contra 0 não é uma parvoíce.

# Modelo apenas com intercepto
m <- lm(Love4Taylor ~ 1, taylor)

summary(m)

Call:
lm(formula = Love4Taylor ~ 1, data = taylor)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.2552 -0.8202  0.1648  0.7548  2.4348 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.5652     0.1549   3.649 0.000638 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.095 on 49 degrees of freedom

Testar Contra Outro Valor

No entanto, há casos em que testar contra zero é bastante parvo…

library(palmerpenguins)

ds <- penguins

m <- lm(body_mass_g ~ 1, ds)

summary(m)

Call:
lm(formula = body_mass_g ~ 1, data = ds)

Residuals:
    Min      1Q  Median      3Q     Max 
-1501.8  -651.8  -151.8   548.2  2098.2 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4201.75      43.36   96.89   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 802 on 341 degrees of freedom
  (2 observations deleted due to missingness)

O que acabámos de testar?

Que a massa dos pinguins na população difere de zero…ou seja, eles existem…pesam alguma coisa…

No entanto, é fácil usar a mesma lógica para testar hipóteses mais interessantes, por exemplo, será que a massa dos pinguins difere de 4kg (4000g).

No modelo linear, testar contra um valor c é o mesmo que centrar os dados em c e voltar a testar contra 0.

m <- lm(body_mass_g - 4000 ~ 1, ds)

summary(m)

Call:
lm(formula = body_mass_g - 4000 ~ 1, data = ds)

Residuals:
    Min      1Q  Median      3Q     Max 
-1501.8  -651.8  -151.8   548.2  2098.2 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   201.75      43.36   4.652  4.7e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 802 on 341 degrees of freedom
  (2 observations deleted due to missingness)
# Ou críamos uma nova coluna e usamos essa coluna...
ds$mass_minus_4000 <- ds$body_mass_g - 4000

m <- lm(mass_minus_4000 ~ 1, ds)

summary(m)

Call:
lm(formula = mass_minus_4000 ~ 1, data = ds)

Residuals:
    Min      1Q  Median      3Q     Max 
-1501.8  -651.8  -151.8   548.2  2098.2 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   201.75      43.36   4.652  4.7e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 802 on 341 degrees of freedom
  (2 observations deleted due to missingness)

2 Nível 2: amostras emparelhadas

Um teste-t para amostras emparelhado nada mais é que um teste-t a uma amostra aplicado à diferença entre momentos de medida.

\(\Delta = (\text{pós} - \text{pré}),\qquad \Delta_i = \beta_0 + e_i\)

Exemplo: reduzir risco de “luta”

Um mini-programa de intervenção para reduzir o risco de confrontos (em humanos…ou pinguins com mau feitio) foi avaliado ao medir-se o risco de confronto antes e depois do tratamento.

fight <- read.csv("../data/fight_risk.csv")
head(fight)
pp pre post diff
pp01 3.155373 1.659283 1.4960897
pp02 5.260346 5.163075 0.0972709
pp03 0.000000 0.000000 0.0000000
pp04 0.000000 0.000000 0.0000000
pp05 6.098303 5.454262 0.6440405
pp06 3.482212 3.150342 0.3318697
m <- lm(diff ~ 1, fight)

summary(m)

Call:
lm(formula = diff ~ 1, data = fight)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.4281 -0.4303  0.0094  0.3242  1.1333 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.48053    0.05703   8.427 1.28e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.5101 on 79 degrees of freedom

O gráfico abaixo ajuda a visualizar o emparelhamento.

Código
library(tidyr)

fight_long <- pivot_longer(fight, c("pre", "post"),
                           names_to = "momento", values_to = "risco")

ggplot(fight_long, aes(x = momento, y = risco, group = pp)) +
geom_line(alpha = 0.35) +
geom_point() +
theme_classic()

3 Nível 3: duas amostras independentes (2 grupos)

Aqui já temos um preditor categórico com 2 níveis.

\(Y_i = \beta_0 + \beta_1 X + e_i\)

Vamos usar palmerpenguins para um exemplo rápido (é muito visual).

library(palmerpenguins)

ds <- penguins

# Garantir que sex tem apenas níveis válidos (o dataset original tem NAs)
ds <- ds[!is.na(ds$sex), ]
m <- lm(body_mass_g ~ sex, ds)

Anova(m, type = 3)
Sum Sq Df F value Pr(>F)
(Intercept) 5884776622 1 11043.50022 0
sex 38878897 1 72.96099 0
Residuals 176380769 331 NA NA

Médias marginais

Como há apenas 2 níveis, pedir pairwise costuma é redundante. Vamos pedir só as médias.

emmeans(m, ~ sex)
 sex    emmean   SE  df lower.CL upper.CL
 female   3862 56.8 331     3750     3974
 male     4546 56.3 331     4435     4656

Confidence level used: 0.95 

4 Nível 4: One-way ANOVA (k > 2)

Agora temos um factor com mais de 2 níveis (e portanto faz sentido pedir comparações múltiplas).

m <- lm(body_mass_g ~ species, ds)
Anova(m, type = 3)
Sum Sq Df F value Pr(>F)
(Intercept) 5242015989 1 24687.8683 0
species 145190219 2 341.8949 0
Residuals 70069447 330 NA NA

Médias + post-hoc numa chamada (eficiente)

# Numa lista só: médias e comparações
emm <- emmeans(m, list(pairwise ~ species))

print(emm)
$`emmeans of species`
 species   emmean   SE  df lower.CL upper.CL
 Adelie      3706 38.1 330     3631     3781
 Chinstrap   3733 55.9 330     3623     3843
 Gentoo      5092 42.2 330     5009     5176

Confidence level used: 0.95 

$`pairwise differences of species`
 1                  estimate   SE  df t.ratio p.value
 Adelie - Chinstrap    -26.9 67.7 330  -0.398  0.9164
 Adelie - Gentoo     -1386.3 56.9 330 -24.359 <0.0001
 Chinstrap - Gentoo  -1359.3 70.0 330 -19.406 <0.0001

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

Ajustes de Valores-p

Em emmeans, as comparações (e os seus ajustes) aparecem de forma mais transparente quando usamos summary().

# Comparações com FDR (BH)
emmeans(m, pairwise ~ species, adjust = "fdr")
$emmeans
 species   emmean   SE  df lower.CL upper.CL
 Adelie      3706 38.1 330     3631     3781
 Chinstrap   3733 55.9 330     3623     3843
 Gentoo      5092 42.2 330     5009     5176

Confidence level used: 0.95 

$contrasts
 contrast           estimate   SE  df t.ratio p.value
 Adelie - Chinstrap    -26.9 67.7 330  -0.398  0.6909
 Adelie - Gentoo     -1386.3 56.9 330 -24.359 <0.0001
 Chinstrap - Gentoo  -1359.3 70.0 330 -19.406 <0.0001

P value adjustment: fdr method for 3 tests 
Nota

Não existe um ajuste “mágico”.

  • Tukey: muito comum em ANOVA clássica (família de comparações entre todos os pares).

  • Holm / BH: alternativas frequentes quando há muitas comparações.

O importante é ser explícito e coerente com o objectivo da análise.

5 Nível 5: ANOVA factorial (interacções)

A pergunta: o efeito de species depende de sex?

m <- lm(body_mass_g ~ species * sex, ds)

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

Na presença de uma interacção, não interprete efeitos principais como se estivessem sozinhos.

Ignorar a interacção é não ouvir o que os dados nos querem dizer e contar apenas uma parte da história.

O que normalmente queremos são efeitos simples / comparações condicionais.

Efeitos Simples/Comparações Condicionais

Queremos:

  • Diferenças entre espécies dentro de cada sexo

  • Diferenças entre sexos dentro de cada espécie

emmeans(m, list(pairwise ~ species | sex, pairwise ~ sex | species))
$`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 

$`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

Visualizar a interacção (obrigatório)

Abaixo estão duas visões complementares:

  1. Dados brutos (o mais “honesto” possível)
# Dados brutos (exemplo)
ggplot(ds, aes(x = species, y = body_mass_g, fill = sex, color = sex)) +
geom_violin(alpha = .5, color = "black") +
geom_jitter(position = position_jitterdodge(dodge.width = .9, jitter.width = .15),
            show.legend = FALSE) +
theme_classic()

  1. Médias marginais (o que o modelo está realmente a prever)
# Médias marginais estimadas
for_g <- as.data.frame(emmeans(m, ~ species * sex))

ggplot(for_g, aes(x = species, y = emmean, ymin = lower.CL, ymax = upper.CL,
                  color = sex, group = sex)) +
geom_line() + geom_point(size = 2.5) + geom_errorbar(width = 0.15) +
theme_classic()

6 Nível 6: ANCOVA (factor + covariável)

Agora juntamos uma covariável contínua (ex.: flipper_length_mm).

m_main <- lm(body_mass_g ~ species * sex + flipper_length_mm, ds)

Anova(m_main, type = 3)
Sum Sq Df F value Pr(>F)
(Intercept) 8.052763e+02 1 0.0097542 0.9213869
species 1.262607e+07 2 76.4688266 0.0000000
sex 1.274795e+07 1 154.4138757 0.0000000
flipper_length_mm 4.389049e+06 1 53.1638663 0.0000000
species:sex 1.739989e+06 2 10.5381054 0.0000367
Residuals 2.691358e+07 326 NA NA
NotaTradição experimental vs regressão
  • Em ANOVA factorial experimental, é comum testar interacções por defeito (e usar SS Tipo III para avaliar efeitos “controlando” o resto).

  • Em regressão múltipla, muita gente coloca apenas efeitos principais (por vezes por hábito, por vezes para reduzir complexidade).

Nenhuma destas tradições é automaticamente “errada”.

  • Evitar interacções pode ser prudente (menos parâmetros; menos espaço para pescarias).

  • Mas ignorar interacções quando elas existem pode esconder o fenómeno.

Interacções com covariáveis (quando considerar)

Um caso clássico é permitir que a relação entre flipper_length_mm e body_mass_g varie por espécie/sexo.

m_int <- lm(body_mass_g ~ species * sex * flipper_length_mm, ds)

Anova(m_main, type = 3)
Sum Sq Df F value Pr(>F)
(Intercept) 8.052763e+02 1 0.0097542 0.9213869
species 1.262607e+07 2 76.4688266 0.0000000
sex 1.274795e+07 1 154.4138757 0.0000000
flipper_length_mm 4.389049e+06 1 53.1638663 0.0000000
species:sex 1.739989e+06 2 10.5381054 0.0000367
Residuals 2.691358e+07 326 NA NA

7 Cábula do emmeans

Situação Pedido típico Nota
1 amostra summary(lm(y ~ 1)) Intercepto = média
2 grupos emmeans(m, ~ A) Post-hoc é redundante
k > 2 emmeans(m, list(~ A, pairwise ~ A)) Post-hoc faz sentido
Interacção emmeans(m, ~ A | B) + pairs()/contrast() Efeitos simples/condicionais

8 emmeans: Alguns Truques

Nesta secção não vou entrar em teoria (isso fica para o capítulo teórico). A ideia é apenas dar-lhe botões que resolvem a maioria das dúvidas práticas quando está a decompor efeitos.

Gráfico Rápido de Interacções

Quando tem uma interacção, é muito fácil interpretar mal uma tabela. O emmip() é uma forma rápida de ver o que o modelo está a prever.

m <- lm(body_mass_g ~ species * sex, ds)

emmip(m, sex ~ species, CIs = TRUE)

Dica

O emmip() é óptimo para ver padrões, mas para um gráfico de artigo normalmente vai querer um ggplot2 mais personalizado.

Contrastes Planeados

Nem sempre a pergunta científica é “quem difere de quem?”. Às vezes a pergunta é um contraste específico, por exemplo:

  • Gentoo vs a média de (Adelie e Chinstrap).
levels(ds$species)
[1] "Adelie"    "Chinstrap" "Gentoo"   
# Contraste planeado (vector de comprimento 3)
# Ordem assumida: Adelie, Chinstrap, Gentoo
emmeans(m, ~ species, contr = list(Gentoo_vs_Outros = c(-0.5, -0.5, 1)))
$emmeans
 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 

$contrasts
 contrast         estimate   SE  df t.ratio p.value
 Gentoo_vs_Outros     1363 36.3 327  37.495 <0.0001

Results are averaged over the levels of: sex 
Aviso

Confirme sempre a ordem dos níveis (a ordem do vector tem de bater certo)

Ajustes de Valores-p

emm_sp <- emmeans(m, ~ sex)

# O mesmo conjunto de comparações, com ajustes diferentes
pairs(emm_sp, adjust = "tukey")
 contrast      estimate   SE  df t.ratio p.value
 female - male     -631 35.7 327 -17.659 <0.0001

Results are averaged over the levels of: species 
pairs(emm_sp, adjust = "fdr")
 contrast      estimate   SE  df t.ratio p.value
 female - male     -631 35.7 327 -17.659 <0.0001

Results are averaged over the levels of: species 
pairs(emm_sp, adjust = "bonferroni")
 contrast      estimate   SE  df t.ratio p.value
 female - male     -631 35.7 327 -17.659 <0.0001

Results are averaged over the levels of: species 
Nota

O método “certo” depende do objectivo.

  • Se está a fazer todas as comparações possíveis porque sim, ajustes mais conservadores são uma defesa contra falsos positivos.

  • Se tem contrastes planeados (pré-especificados), faz sentido reportar esses contrastes com muito mais foco.

9 Exercícios

  1. No dataset Love4Taylor.csv, teste se a média é diferente de 0 e de um outro valor (idealmente um que tenha algum significado).
Clique para ver as soluções
# Repetir o que fizemos no Nível 1

taylor <- read.csv("../data/Love4Taylor.csv")

# Contra 0
m <- lm(Love4Taylor ~ 1, taylor)

summary(m)

Call:
lm(formula = Love4Taylor ~ 1, data = taylor)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.2552 -0.8202  0.1648  0.7548  2.4348 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   0.5652     0.1549   3.649 0.000638 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.095 on 49 degrees of freedom
# Contra outro valor (e.g., -1)
taylor$love_vs_neg1 <- taylor$Love4Taylor + 1

m <- lm(love_vs_neg1 ~ 1, taylor)

summary(m)

Call:
lm(formula = love_vs_neg1 ~ 1, data = taylor)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.2552 -0.8202  0.1648  0.7548  2.4348 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   1.5652     0.1549   10.11 1.42e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.095 on 49 degrees of freedom
  1. No dataset fight_risk.csv, faça:
  • o modelo lm(diff ~ 1)

  • as descritivas do pré e pós

  • uma frase que reporte o resultado (uma ou duas linhas).

Clique para ver as soluções
library(tidyr)

fight <- read.csv("../data/fight_risk.csv")

m <- lm(diff ~ 1, fight)

summary(m)

Call:
lm(formula = diff ~ 1, data = fight)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.4281 -0.4303  0.0094  0.3242  1.1333 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  0.48053    0.05703   8.427 1.28e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.5101 on 79 degrees of freedom
mean(fight$pre)
[1] 3.043889
sd(fight$pre)
[1] 1.560794
mean(fight$post)
[1] 2.563358
sd(fight$post)
[1] 1.574395

Exemplo de como reportar os resultados:

Encontrámos evidência significativa de que o tratamento reduziu o risco de comportamentos violentos, t(29) = 7.45, p < .001; Mpre = 50.51, SD = 7.80, Mpós = 45.86, SD = 7.65.

  1. Com palmerpenguins, ajuste lm(body_mass_g ~ species * sex) e faça:
  • comparações entre espécies dentro de cada sexo

  • comparação entre sexos dentro de cada espécie (uma por espécie)

  • dois gráficos: dados brutos e emmeans.

Clique para ver as soluções
library(car)
library(palmerpenguins)

ds <- penguins

ds <- ds[!is.na(ds$sex), ]
m <- lm(body_mass_g ~ species * sex, ds)

Anova(m, type = 3)
Sum Sq Df F value Pr(>F)
(Intercept) 5232595969 1 54661.827958 0.0000000
species 143001222 2 746.924492 0.0000000
sex 29851220 1 311.838003 0.0000000
species:sex 1676557 2 8.756997 0.0001973
Residuals 31302628 327 NA NA
# (a) species dentro de sex
sp_within_sex <- emmeans(m, list(pairwise ~ species | sex, pairwise ~ sex | species))

# gráficos (exemplos/sugestões)
ggplot(ds, aes(x = species, y = body_mass_g, fill = sex, color = sex)) +
geom_violin(alpha = .5, color = "black") +
geom_jitter(position = position_jitterdodge(dodge.width = .9, jitter.width = .15),
            show.legend = FALSE) +
theme_classic()

for_g <- as.data.frame(emmeans(m, ~ species * sex))

ggplot(for_g, aes(x = species, y = emmean, ymin = lower.CL, ymax = upper.CL,
                  color = sex, group = sex)) +
geom_line() + geom_point(size = 2.5) + geom_errorbar(width = 0.15) +
theme_classic()

  1. Faça uma ANCOVA com flipper_length_mm e compare um modelo sem e com interacções com a covariável.
Clique para ver as soluções
library(car)
library(palmerpenguins)

m_main  <- lm(body_mass_g ~ species * sex + flipper_length_mm, ds)
m_int <- lm(body_mass_g ~ species * sex * flipper_length_mm, ds)

Anova(m_int, type = 3)
Sum Sq Df F value Pr(>F)
(Intercept) 47510.60 1 0.5858908 0.4445749
species 216230.99 2 1.3332578 0.2650728
sex 16183.84 1 0.1995757 0.6553648
flipper_length_mm 4492633.05 1 55.4022165 0.0000000
species:sex 634560.21 2 3.9126323 0.0209482
species:flipper_length_mm 279932.75 2 1.7260362 0.1796363
sex:flipper_length_mm 73438.41 1 0.9056272 0.3419933
species:sex:flipper_length_mm 596743.03 2 3.6794555 0.0263072
Residuals 26030280.01 321 NA NA
anova(m_main, m_int)
Res.Df RSS Df Sum of Sq F Pr(>F)
326 26913579 NA NA NA NA
321 26030280 5 883299.4 2.178533 0.0563138