ANOVA e Vários Preditores

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

João O. Santos

joao.santos@ispa.pt

ISPA

Dois Preditores Quantitativos

Como varia a massa com o comprimento da barbatana e do bico?

Uma linha = um espécime de pinguim; duas medidas no mesmo animal são duas colunas, não duas réplicas.

Preparar os Dados

library(car)
library(dplyr)
library(emmeans)
library(ggplot2)
library(palmerpenguins)
library(tidyr)
options(contrasts = c("contr.sum", "contr.poly"))
ds <- penguins |>
  select(species, sex, body_mass_g, flipper_length_mm, bill_length_mm) |>
  drop_na()

Os mesmos casos completos em todas as comparações.

Ver as Duas Relações

Relações separadas: ainda sem ter em conta o outro comprimento.

Somar Contribuições

\[ \widehat{\text{massa}}=b_0+b_L\,\text{barbatana}+b_B\,\text{bico} \]

Os \(b\) estimam os parâmetros \(\beta\) do modelo.

\(b_L\): g/mm de barbatana, mantendo o comprimento do bico constante.

\(b_B\): g/mm de bico, mantendo o comprimento da barbatana constante.

Sem interacção: cada declive é constante ao longo do outro preditor.

Ajustar o Modelo Aditivo

m_multi <- lm(body_mass_g ~ flipper_length_mm + bill_length_mm, data = ds)

O intercepto prevê a massa com ambos os comprimentos a zero: fora dos dados, não um animal real.

Ver o Ajuste Conjunto

As rectas são cortes do modelo; nem todas as combinações têm dados próximos.

Comparação Global

Demonstração do F global do resumo: dois modelos para perceber o teste, não para acrescentar uma etapa ao script.

m0 <- lm(body_mass_g ~ 1, data = ds)
anova(m0, m_multi)

Analysis of Variance Table

Model 1: body_mass_g ~ 1
Model 2: body_mass_g ~ flipper_length_mm + bill_length_mm
  Res.Df       RSS Df Sum of Sq      F    Pr(>F)    
1    332 215259666                                  
2    330  51071963  2 164187703 530.45 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Compacto: uma média. Aumentado: dois declives e intercepto.

Restrição conjunta: \(\beta_L=\beta_B=0\). Justifica-se acrescentar dois parâmetros?

Resumo do Modelo

summary(m_multi)

Call:
lm(formula = body_mass_g ~ flipper_length_mm + bill_length_mm, 
    data = ds)

Residuals:
     Min       1Q   Median       3Q      Max 
-1083.08  -282.65   -17.18   242.95  1293.24 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)       -5836.299    312.604 -18.670   <2e-16 ***
flipper_length_mm    48.890      2.034  24.034   <2e-16 ***
bill_length_mm        4.959      5.214   0.951    0.342    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 393.4 on 330 degrees of freedom
Multiple R-squared:  0.7627,    Adjusted R-squared:  0.7613 
F-statistic: 530.4 on 2 and 330 DF,  p-value: < 2.2e-16

O teste global tem 2 gl no numerador; não testa cada declive separadamente.

Relatar o Modelo Global

A massa aumentou com o comprimento da barbatana, ajustando para o bico: 48.9 g/mm, SE = 2.0 g/mm. O modelo aditivo melhorou o ajuste relativamente a uma média comum, F_{2, 330} = 530.45, P < .001.

O teste global não diz que ambos os preditores contribuem.

O Que Acrescenta a Barbatana?

Demonstração da linha flipper_length_mm de Anova(m_multi). anova() aqui serve para fins didácticos.

m_no_flipper <- lm(body_mass_g ~ bill_length_mm, data = ds)
anova(m_no_flipper, m_multi)

Analysis of Variance Table

Model 1: body_mass_g ~ bill_length_mm
Model 2: body_mass_g ~ flipper_length_mm + bill_length_mm
  Res.Df       RSS Df Sum of Sq      F    Pr(>F)    
1    331 140467133                                  
2    330  51071963  1  89395170 577.62 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Restrição: \(\beta_L=0\), mantendo o bico no modelo.

O Que Acrescenta o Bico?

Demonstração da linha bill_length_mm da mesma tabela.

m_no_bill <- lm(body_mass_g ~ flipper_length_mm, data = ds)
anova(m_no_bill, m_multi)

Analysis of Variance Table

Model 1: body_mass_g ~ flipper_length_mm
Model 2: body_mass_g ~ flipper_length_mm + bill_length_mm
  Res.Df      RSS Df Sum of Sq      F Pr(>F)
1    331 51211963                           
2    330 51071963  1    140000 0.9046 0.3422

Restrição: \(\beta_B=0\), mantendo a barbatana no modelo.

Uma Tabela, Testes Parciais

Anova(m_multi, type = "III")

Anova Table (Type III tests)

Response: body_mass_g
                    Sum Sq  Df  F value Pr(>F)    
(Intercept)       53945591   1 348.5679 <2e-16 ***
flipper_length_mm 89395170   1 577.6243 <2e-16 ***
bill_length_mm      140000   1   0.9046 0.3422    
Residuals         51071963 330                    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

• Cada linha: um preditor, mantendo o outro no modelo
• Estes Fs: testes das duas demonstrações, não o F global
• Na análise: o resumo já dá os testes dos dois declives
• Tabela ANOVA: sem hipótese nova; sem sinal ou magnitude dos coeficientes
• Aqui: ligar as duas saídas, não juntar testes ao script

Anova() vem do pacote car, carregado no início. Não confundir com anova(model) com um só modelo, que dá testes sequenciais. Neste modelo aditivo, tipos II e III testam as mesmas contribuições dos preditores. A linha do intercepto não responde à pergunta biológica. Somas de quadrados.

Quando F É t²

Compare os valores de t no resumo com os valores de F na tabela anterior.

Testes parciais: cada restrição tem 1 gl. O F repete o teste-t do coeficiente, com \(F=t^2\) e o mesmo valor-P bilateral.

Teste F global: 2 gl para \(\beta_L=\beta_B=0\) em conjunto. Nenhum dos ts dos coeficientes responde sozinho a essa pergunta.

É aqui que o F começa a fazer mais do que repetir um t ao quadrado.

A identidade usa o mesmo modelo linear e a mesma estimativa da variância residual. Não estender a identidade a testes diferentes só porque se chamam t, F ou qui-quadrado. Em particular, o teste de um factor com k > 2 níveis terá mais de um grau de liberdade; nenhum t isolado substitui esse teste.

Relatar os Declives

A comprimento do bico constante, a massa aumentou 48.9 g por mm de barbatana, SE = 2.0 g/mm, t₃₃₀ = 24.03, P < .001. Para o bico, a barbatana constante, o declive estimado foi 5.0 g/mm, SE = 5.2 g/mm, t₃₃₀ = 0.95, P = .342.

Ainda não distinguimos espécies. A amostragem limita a generalização; ajustar não torna a relação causal.

Script Quantitativo

library(ggplot2)
library(palmerpenguins)
library(tidyr)

ds <- penguins |>
  drop_na(species, sex, body_mass_g, flipper_length_mm, bill_length_mm)
model <- lm(body_mass_g ~ flipper_length_mm + bill_length_mm, data = ds)
g <- ggplot(ds, aes(flipper_length_mm, body_mass_g)) + geom_point()

# Print results and plot
summary(model) |> print()
print(g)

Um modelo e um resumo: os ts testam cada declive; o F testa os dois em conjunto. A tabela ANOVA e os modelos reduzidos ficaram nas demonstrações.

Resultados dos Declives

Call:
lm(formula = body_mass_g ~ flipper_length_mm + bill_length_mm, 
    data = ds)

Residuals:
     Min       1Q   Median       3Q      Max 
-1083.08  -282.65   -17.18   242.95  1293.24 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)       -5836.299    312.604 -18.670   <2e-16 ***
flipper_length_mm    48.890      2.034  24.034   <2e-16 ***
bill_length_mm        4.959      5.214   0.951    0.342    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 393.4 on 330 degrees of freedom
Multiple R-squared:  0.7627,    Adjusted R-squared:  0.7613 
F-statistic: 530.4 on 2 and 330 DF,  p-value: < 2.2e-16

Agora: Três Grupos

A massa corporal difere entre espécies?

Várias Linhas por Espécie

linha = espécime · species = factor entre unidades

Repetir o rótulo da espécie não é repetir o espécime.

Espécimes diferentes podem partilhar uma espécie sem pseudorreplicação. Duas medições do mesmo espécime repetem a unidade. Perguntar ainda por local, ocasião de amostragem e outros agrupamentos relevantes: espécies iguais ou IDs diferentes não garantem independência.

Porquê \(k-1\)?

Para descrever três médias, precisamos de três parâmetros. Um já é o intercepto: restam dois para distinguir as espécies.

Níveis do Factor, \(k\)	Intercepto	Coeficientes de Espécie
2	1	1
3	1	2
4	1	3

Com intercepto, um factor de \(k\) níveis usa \(k-1\) coeficientes. Não acrescentamos um coeficiente por espécime!

A contagem pressupõe níveis observados e um modelo de posto completo. Com intercepto, incluir todos os k indicadores duplicaria informação: a soma desses indicadores já é a coluna de uns. Sem intercepto, response ~ 0 + species pode estimar directamente k médias; não é um modelo com uma média a menos. Aqui mantemos o intercepto e a codificação soma-zero.

Três Níveis, Dois Contrastes

\[ \widehat{Y}=b_0+b_A C_A+b_C C_C \]

m_species <- lm(body_mass_g ~ species, data = ds)
contrasts(ds$species)

          [,1] [,2]
Adelie       1    0
Chinstrap    0    1
Gentoo      -1   -1

Soma-zero não tem grupo de referência. \(b_0\) é a média das três médias, com pesos iguais; \(b_A\) e \(b_C\) são os desvios de Adelie e Chinstrap.

Ler Cada Linha

Espécie	\(C_A\)	\(C_C\)	Massa Prevista
Adelie	1	0	\(b_0+b_A\)
Chinstrap	0	1	\(b_0+b_C\)
Gentoo	−1	−1	\(b_0-b_A-b_C\)

Gentoo não é uma referência escondida. O seu desvio é \(-b_A-b_C\): os três desvios somam zero.

Das Contas às Médias

coef(m_species)

(Intercept)    species1    species2 
  4177.2299   -471.0655   -444.1416 

Espécie	Substituir na Equação	Massa Prevista
Adelie	\(4177.23-471.07\)	≈ 3706 g
Chinstrap	\(4177.23-444.14\)	≈ 3733 g
Gentoo	\(4177.23+471.07+444.14\)	≈ 5092 g

Os coeficientes são interpretáveis. Fazer estas contas à mão é que se torna pouco prático — sobretudo quando quisermos também a incerteza.

O R Faz as Contas

emm_species <- emmeans(m_species, ~ species)
emm_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 

As mesmas três médias, agora com erros-padrão e IC 95%. emmeans() combina os coeficientes e a sua incerteza.

Que Testa Cada Coeficiente?

Se \(\bar\mu=(\mu_A+\mu_C+\mu_G)/3\), os testes-t de soma-zero avaliam:

Coeficiente	Hipótese
species1, \(b_A\)	Adelie tem a média \(\bar\mu\)
species2, \(b_C\)	Chinstrap tem a média \(\bar\mu\)

Não são testes de Adelie contra Chinstrap ou Gentoo. Para perguntar se há diferenças entre espécies, testamos os dois coeficientes em conjunto. Depois pedimos as médias e as comparações que nos interessam.

Curiosidade: Onde Está o Zero?

Codificação	Coeficiente	O Que Testa o Seu t?
Soma-zero	species1	Adelie = média das três médias
Soma-zero	species2	Chinstrap = média das três médias
0/1, referência Adelie	speciesChinstrap	Chinstrap = Adelie
0/1, referência Adelie	speciesGentoo	Gentoo = Adelie

O intercepto passa da média das três médias para a média de Adelie.

Mudam os coeficientes e as perguntas dos seus ts, não as médias previstas nem o teste F conjunto de espécie.

A codificação de tratamento, o padrão inicial do R, pode ser pedida só para um modelo com contrasts = list(species = “contr.treatment”). Não precisamos de ensinar duas receitas paralelas: a tabela mostra por que um t individual não é uma pergunta invariável sobre a espécie inteira. Demonstração e detalhes de codificação.

Duas Restrições em Conjunto

Demonstração do teste de espécie: que hipótese representa a linha com 2 gl?

anova(m0, m_species)

Analysis of Variance Table

Model 1: body_mass_g ~ 1
Model 2: body_mass_g ~ species
  Res.Df       RSS Df Sum of Sq      F    Pr(>F)    
1    332 215259666                                  
2    330  70069447  2 145190219 341.89 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Compacto: um intercepto. Aumentado: intercepto + dois contrastes.

\(H_0: \beta_A=\beta_C=0\) equivale a três médias iguais. É uma ANOVA com espécie como factor: o mesmo raciocínio omnibus, agora com categorias.

Ver as Previsões

Tracejado: média comum. Verde: uma média por espécie.

Dois ts Não Chegam

summary(m_species)

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

Residuals:
     Min       1Q   Median       3Q      Max 
-1142.44  -342.44   -33.09   307.56  1207.56 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4177.23      26.59  157.12   <2e-16 ***
species1     -471.07      34.52  -13.65   <2e-16 ***
species2     -444.14      41.80  -10.62   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 460.8 on 330 degrees of freedom
Multiple R-squared:  0.6745,    Adjusted R-squared:  0.6725 
F-statistic: 341.9 on 2 and 330 DF,  p-value: < 2.2e-16

Os ts de species1 e species2 testam \(\beta_A=0\) e \(\beta_C=0\) separadamente. Não são o teste conjunto de igualdade das três médias.

Uma Linha, Dois Coeficientes

Anova(m_species, type = "III")

Anova Table (Type III tests)

Response: body_mass_g
                Sum Sq  Df  F value    Pr(>F)    
(Intercept) 5242015989   1 24687.87 < 2.2e-16 ***
species      145190219   2   341.89 < 2.2e-16 ***
Residuals     70069447 330                       
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Agora a tabela ANOVA agrupa o que a tabela de coeficientes separa: a linha species testa \(\beta_A=\beta_C=0\), com 2 gl no numerador.

Com \(k\) níveis, o factor tem \(k-1\) coeficientes. Para \(k>2\), nenhum t isolado substitui o teste conjunto do factor.

Como só há um preditor, este é também o F global do resumo: não é outro teste.

Cumprimos a promessa da aula de um preditor: o termo já não corresponde a um único coeficiente. Não estamos a obter informação nova sobre os dados, mas a responder a uma pergunta conjunta que nenhum t individual responde. Como este modelo só tem um preditor, o F global no fim de summary() já testa a mesma hipótese que a linha species. Não apresentar essa repetição como mais uma confirmação; no modelo com factor e comprimento, a tabela vai separar os dois termos, ao contrário do F global.

Relatar a ANOVA

Comparámos a massa corporal das três espécies através de uma ANOVA com espécie como factor. Os resultados indicaram diferenças entre espécies, F(2, 330) = 341.89, P < .001.

Não identifica quais as médias que diferem. Nenhum t isolado dos dois contrastes substitui este F com 2 gl.

Que Espécies Diferem?

emm_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 

species_pairs <- pairs(emm_species)
species_pairs

 contrast           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 

confint(species_pairs)

 contrast           estimate   SE  df lower.CL upper.CL
 Adelie - Chinstrap    -26.9 67.7 330     -186      132
 Adelie - Gentoo     -1386.3 56.9 330    -1520    -1252
 Chinstrap - Gentoo  -1359.3 70.0 330    -1524    -1194

Confidence level used: 0.95 
Conf-level adjustment: tukey method for comparing a family of 3 estimates 

Primeiro médias em gramas e IC 95%; depois diferenças na ordem impressa. Tukey ajusta a família das três comparações.

Relatar as Comparações

Mais concretamente, Gentoo apresentou maior massa (M = 5092 g, SE = 42.2 g) do que Adelie (M = 3706 g, SE = 38.1 g), t(330) = 24.36, P < .001, e Chinstrap (M = 3733 g, SE = 55.9 g), t(330) = 19.41, P < .001.

Não encontrámos evidência de diferença entre Adelie e Chinstrap, t(330) = −0.40, P = .916. Os valores-P têm ajuste de Tukey. Não encontrar evidência de diferença não estabelece igualdade.

As médias e os erros-padrão vêm do modelo com variância residual comum; os t referem-se a Gentoo menos Adelie, Gentoo menos Chinstrap e Adelie menos Chinstrap, respectivamente. A tabela imprime os dois primeiros no sentido oposto; inverter o contraste muda o sinal, não o valor-P bilateral. Os dados são observacionais e a generalização depende da amostragem.

Script Entre Espécies

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

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

ds <- penguins |>
  drop_na(species, sex, body_mass_g, flipper_length_mm, bill_length_mm)
model <- lm(body_mass_g ~ species, data = ds)
emm <- emmeans(model, ~ species)
g <- ggplot(ds, aes(species, body_mass_g)) + geom_boxplot()

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

Resultados Entre Espécies

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

Residuals:
     Min       1Q   Median       3Q      Max 
-1142.44  -342.44   -33.09   307.56  1207.56 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4177.23      26.59  157.12   <2e-16 ***
species1     -471.07      34.52  -13.65   <2e-16 ***
species2     -444.14      41.80  -10.62   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 460.8 on 330 degrees of freedom
Multiple R-squared:  0.6745,    Adjusted R-squared:  0.6725 
F-statistic: 341.9 on 2 and 330 DF,  p-value: < 2.2e-16

 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 

 contrast           estimate   SE  df lower.CL upper.CL t.ratio p.value
 Adelie - Chinstrap    -26.9 67.7 330     -186      132  -0.398  0.9164
 Adelie - Gentoo     -1386.3 56.9 330    -1520    -1252 -24.359 <0.0001
 Chinstrap - Gentoo  -1359.3 70.0 330    -1524    -1194 -19.406 <0.0001

Confidence level used: 0.95 
Conf-level adjustment: tukey method for comparing a family of 3 estimates 
P value adjustment: tukey method for comparing a family of 3 estimates 

Combinar Factor e Comprimento

E se compararmos espécies ao mesmo comprimento da barbatana?

m_add <- lm(body_mass_g ~ species + flipper_length_mm, data = ds)
summary(m_add)

Call:
lm(formula = body_mass_g ~ species + flipper_length_mm, data = ds)

Residuals:
   Min     1Q Median     3Q    Max 
-898.8 -252.0  -24.8  229.8 1191.6 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)       -3986.80     619.53  -6.435 4.36e-10 ***
species1            -26.38      43.81  -0.602    0.547    
species2           -231.76      37.50  -6.180 1.89e-09 ***
flipper_length_mm    40.61       3.08  13.186  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 373.3 on 329 degrees of freedom
Multiple R-squared:  0.787, Adjusted R-squared:  0.7851 
F-statistic: 405.3 on 3 and 329 DF,  p-value: < 2.2e-16

Anova(m_add, type = "III")

Anova Table (Type III tests)

Response: body_mass_g
                    Sum Sq  Df F value    Pr(>F)    
(Intercept)        5770339   1  41.412 4.355e-10 ***
species            5368818   2  19.265 1.225e-08 ***
flipper_length_mm 24226302   1 173.864 < 2.2e-16 ***
Residuals         45843144 329                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

ANCOVA aditiva: dois contrastes + um declive comum.

Na mesma tabela, species reúne dois coeficientes num F com 2 gl; flipper_length_mm tem um só coeficiente e o seu F repete o \(t^2\).

O F global do resumo não separa estas duas perguntas.

Ver as Rectas Paralelas

Declive comum; diferenças entre espécies constantes.

Relatar o Ajuste Aditivo

Considerando a espécie, a massa aumentou 40.6 g por mm de barbatana, SE = 3.1 g/mm, t₃₂₉ = 13.19, P < .001. A espécie contribuiu depois de considerar o comprimento, F_{2, 329} = 19.27, P < .001.

Comparar a comprimento constante não é comparar as médias brutas.

Comparações Ajustadas

adjusted_species <- emmeans(m_add, ~ species,
                            at = list(flipper_length_mm = 205))
adjusted_species

 species   emmean   SE  df lower.CL upper.CL
 Adelie      4311 55.3 329     4202     4420
 Chinstrap   4106 53.4 329     4001     4211
 Gentoo      4596 50.9 329     4495     4696

Confidence level used: 0.95 

adjusted_pairs <- pairs(adjusted_species)
adjusted_pairs

 contrast           estimate   SE  df t.ratio p.value
 Adelie - Chinstrap      205 57.6 329   3.568  0.0012
 Adelie - Gentoo        -285 95.4 329  -2.981  0.0086
 Chinstrap - Gentoo     -490 87.0 329  -5.631 <0.0001

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

confint(adjusted_pairs)

 contrast           estimate   SE  df lower.CL upper.CL
 Adelie - Chinstrap      205 57.6 329     69.8    340.9
 Adelie - Gentoo        -285 95.4 329   -509.2    -59.8
 Chinstrap - Gentoo     -490 87.0 329   -694.7   -285.1

Confidence level used: 0.95 
Conf-level adjustment: tukey method for comparing a family of 3 estimates 

Fixámos o comprimento em 205 mm, dentro do intervalo partilhado pelas três espécies. at indica o valor em que pedimos as médias previstas.

Relatar as Médias Ajustadas

A 205 mm de barbatana, Gentoo apresentou maior massa prevista (M = 4596 g, SE = 50.9 g) do que Adelie (M = 4311 g, SE = 55.3 g), t(329) = 2.98, P = .009, e Chinstrap (M = 4106 g, SE = 53.4 g), t(329) = 5.63, P < .001. Adelie também excedeu Chinstrap, t(329) = 3.57, P = .001. Os valores-P têm ajuste de Tukey.

Ao mesmo comprimento, sob a hipótese de declive comum — não as médias brutas de cada espécie.

Script Aditivo

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

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

ds <- penguins |>
  drop_na(species, sex, body_mass_g, flipper_length_mm, bill_length_mm)
m <- lm(body_mass_g ~ species + flipper_length_mm, data = ds)
adjusted_species <- emmeans(m, ~ species, at = list(flipper_length_mm = 205))
g <- ggplot(ds, aes(flipper_length_mm, body_mass_g, colour = species)) + geom_point()

# Print results and plot
summary(m) |> print()
Anova(m, type = "III") |> print()
print(adjusted_species)
pairs(adjusted_species) |> summary(infer = TRUE) |> print()
print(g)

Um declive comum; diferenças entre espécies ao mesmo comprimento.

Resultados Aditivos

Call:
lm(formula = body_mass_g ~ species + flipper_length_mm, data = ds)

Residuals:
   Min     1Q Median     3Q    Max 
-898.8 -252.0  -24.8  229.8 1191.6 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)       -3986.80     619.53  -6.435 4.36e-10 ***
species1            -26.38      43.81  -0.602    0.547    
species2           -231.76      37.50  -6.180 1.89e-09 ***
flipper_length_mm    40.61       3.08  13.186  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 373.3 on 329 degrees of freedom
Multiple R-squared:  0.787, Adjusted R-squared:  0.7851 
F-statistic: 405.3 on 3 and 329 DF,  p-value: < 2.2e-16

Anova Table (Type III tests)

Response: body_mass_g
                    Sum Sq  Df F value    Pr(>F)    
(Intercept)        5770339   1  41.412 4.355e-10 ***
species            5368818   2  19.265 1.225e-08 ***
flipper_length_mm 24226302   1 173.864 < 2.2e-16 ***
Residuals         45843144 329                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

 species   emmean   SE  df lower.CL upper.CL
 Adelie      4311 55.3 329     4202     4420
 Chinstrap   4106 53.4 329     4001     4211
 Gentoo      4596 50.9 329     4495     4696

Confidence level used: 0.95 

 contrast           estimate   SE  df lower.CL upper.CL t.ratio p.value
 Adelie - Chinstrap      205 57.6 329     69.8    340.9   3.568  0.0012
 Adelie - Gentoo        -285 95.4 329   -509.2    -59.8  -2.981  0.0086
 Chinstrap - Gentoo     -490 87.0 329   -694.7   -285.1  -5.631 <0.0001

Confidence level used: 0.95 
Conf-level adjustment: tukey method for comparing a family of 3 estimates 
P value adjustment: tukey method for comparing a family of 3 estimates 

Agora: Vários Grupos

20 minutos · grupos de 3 ou 5

Modelos Lineares no R — Exercícios

Façam exercício 2.

Agora: Pensar o Modelo

15 minutos · grupos de 3 ou 5

Comparação de Modelos — Exercícios

Façam exercícios 3 e 4.

Desafio no R

15 minutos · grupos de 3 ou 5

Exercício de interpretação, não um passo extra do script.

Comparem body_mass_g ~ flipper_length_mm com body_mass_g ~ species + flipper_length_mm e identifiquem a linha de Anova() que corresponde a essa comparação.

1. Quantos parâmetros foram acrescentados?
2. Qual é a restrição do modelo compacto?
3. Que conclusão continua a não ser causal?

Se Precisarem de Teoria

Comparação de Modelos — Regressão Múltipla

A Seguir

E se o declive não for igual em todas as espécies?