A Abordagem de Comparação de Modelos

Modelos Lineares

Abordagem de Comparação de Modelos

Como descobrir que o teste-t, a ANOVA e a regressão são o mesmo modelo com nomes diferentes

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

Evitem “Cargo-Cult Statistics”

Se já alguma vez sentiram que a estatística é uma montra de testes (um para cada ocasião), não estão sozinhos. A tentação é decorar “qual teste usar quando”, mas isso leva ao que Stark & Saltelli (2018) chamam cargo-cult statistics: “o mimar ritualístico da estatística em vez de prática consciente”.

Este capítulo oferece uma alternativa: pensar nos testes estatísticos equanto comparações entre modelos (Correll et al., 2025). Esta abordagem unificadora não é necessariamente superior ou inferior a outras perspectivas, mas é particularmente útil para quem trabalha com R, pois permite usar lm() (ou outras funções que suportam a sintaxe de fórmulas do R) para uma vasta gama de análises, em vez de andar à procura de funções com o nome específico do “teste” que queremos usar.

Nota Histórica: Fórmulas Específicas para Cada Teste?

No passado, ter fórmulas específicas para cada “teste” facilitava cálculos manuais. Calcular um t “à mão” era mais simples do que estimar o declive da regressão manualmente.

Contudo, nos dias de hoje, a maioria dos computadores corre tudo como regressões—apenas mudam a forma como apresentam os resultados conforme o menu/função que o utilizador usou para correr a análise.

Curiosidade: o manual da função aov() do R diz abertamente que aov() apenas chama lm().

A Velocidade da Álgebra Matricial

Os computadores modernos são incrivelmente rápidos a realizar operações matriciais—a mesma matemática que alimenta jogos 3D, inteligência artificial, e regressões lineares.

Uma regressão linear pode ser resolvida por: \(\hat{\beta} = (X^TX)^{-1}X^TY\)

Esta operação matricial, mesmo com milhares de observações e dezenas de preditores, demora milissegundos num portátil moderno. É por isso que não há penalização computacional em usar lm() para tudo.

O objectivo deste capítulo não é decorar output nem coleccionar funções; é perceber o que cada modelo/teste está a testar.

1 Fundamentos da Inferência Estatística

Antes de começarmos a comparar modelos, precisamos de compreender o que significa “fazer inferência” e porque é que os testes estatísticos funcionam.

Dedução vs Inferência

Dedução (lógica formal):

Todo o A é B
C é A
∴ C é B

A dedução é verdadeira sempre que as suas premissas forem verdadeiras e a conclusão derivar logicamente das premissas.

Inferência (testes de hipóteses):

O João é professor de estatística
O João é parvo
∴ Todos os professores de estatística são parvos

As inferências nunca são garantidas mesmo que as suas premissas sejam verdadeiras e o raciocínio seja lógico.

O Que Isto Significa Para Nós

A inferência estatística tenta testar hipóteses sobre parâmetros populacionais (μ, β, σ²) a partir de estatísticas amostrais (\(\bar{x}\), b, s²).

Nunca temos a certeza absoluta. O melhor que podemos fazer é quantificar a incerteza e controlar a taxa de erros a longo prazo.

O Teorema do Limite Central

Antes de testarmos hipóteses sobre modelos, precisamos de compreender porque é que esses testes funcionam mesmo quando os dados não são perfeitos.

O Problema

Imaginemos que recolhemos uma amostra de 30 participantes e calculamos a média do seu tempo de reacção. Essa média é apenas uma possível média que poderíamos ter obtido—se recolhêssemos outra amostra diferente, obteríamos uma média ligeiramente diferente.

A questão fundamental: Como é que sabemos se a média que observámos é suficientemente diferente de zero (ou de qualquer outro valor de H₀) para concluirmos que a média populacional não é zero (ou outro valor)?

A Distribuição Amostral

O Teorema do Limite Central (CLT) diz-nos algo extraordinário:

Teorema do Limite Central

Se recolhermos muitas amostras da mesma população e calcularmos a média de cada amostra, essas médias amostrais seguirão uma distribuição normal, mesmo que a população original não seja normal.

Propriedades da distribuição de médias amostrais:

Centrada na verdadeira média populacional (μ)
Desvio padrão = \(\frac{\sigma}{\sqrt{N}}\) (chamado erro padrão, SE)
Forma cada vez mais normal à medida que N aumenta

O erro padrão (SE) é um conceito fundamental: Representa o erro de estimação da média. Quando calculamos \(t = \frac{\bar{x}}{SE}\), estamos literalmente a comparar o sinal (a média) com o ruído (o erro dessa estimação).

Demonstração Interativa (opcional)

Experimentem:

Ir medindo vários peixes
Clickar em completar a amostra
Clickar em calcular a média
Clickar em calcular as médias de múltiplas amostras
Alterar os parâmetros (média, desvio-padrão, dimensão da amostra) e repetir o passo 4. O que se altera? O que se mantém?

Link directo para a aplicação

Demonstração Visual em R

Aqui vamos ver, lado a lado, (1) a distribuição populacional e (2) a distribuição amostral da média.

Nota importante: Quando mantemos N fixo e apenas aumentamos o número de amostras simuladas, a distribuição “não fica mais normal”—nós é que a vemos com mais nitidez (a aproximação Monte Carlo fica mais suave/estável).

Código

library(ggplot2)

set.seed(42)

simulate_sampling_means <- function(population, sample_size, n_samples) {
    replicate(n_samples, mean(sample(population, sample_size, replace = TRUE)))
}

make_sampling_df <- function(population, sample_size, n_samples_vec) {
    do.call(rbind, lapply(n_samples_vec, function(ns) {
        data.frame(
            mean = simulate_sampling_means(population, sample_size, ns),
            n_samples = factor(ns, levels = n_samples_vec)
        )
    }))
}

plot_population <- function(population, title, bins = 40, x_limits = NULL) {
    p <- ggplot(data.frame(x = population), aes(x = x)) +
        geom_histogram(aes(y = after_stat(density)), bins = bins,
                       fill = "grey70", color = "white") +
        geom_density(linewidth = 0.8) +
        theme_classic() +
        labs(title = title, x = "Valor", y = "Densidade")

    if (!is.null(x_limits)) {
        p <- p + coord_cartesian(xlim = x_limits)
    }

    return(p)
}

plot_sampling <- function(df_means, title, bins = 30, x_limits = NULL) {
    facet_labels <- paste0(unique(df_means$n_samples), " amostras")
    p <- ggplot(df_means, aes(x = mean)) +
        geom_histogram(aes(y = after_stat(density)), bins = bins,
                       fill = "steelblue", color = "white", alpha = 0.8) +
        geom_density(linewidth = 0.8) +
        facet_wrap(~ n_samples, nrow = 1,
                   labeller = labeller(n_samples = facet_labels)) +
        theme_classic() +
        labs(title = title, x = "Média amostral", y = "Densidade")

    if (!is.null(x_limits)) {
        p <- p + coord_cartesian(xlim = x_limits)
    }

    return(p)
}

n_samples_vec <- c(10, 50, 200, 2000)

População Normal (N pequeno)

Código

pop_normal <- rnorm(50000, mean = 0, sd = 1)
plot_population(pop_normal, "População: Normal(0, 1)", x_limits = c(-4, 4))

Código

N <- 5

df_means_normal <- make_sampling_df(pop_normal, sample_size = N, n_samples_vec)
plot_sampling(df_means_normal,
              paste0("Distribuição amostral da média (Normal; N = ", N, ")"),
              x_limits = c(-2, 2))

População Uniforme (N pequeno)

Código

pop_unif <- runif(50000, min = -1.5, max = 1.5)
plot_population(pop_unif, "População: Uniforme(-1.5, 1.5)", x_limits = c(-2, 2))

Código

N <- 5

df_means_unif <- make_sampling_df(pop_unif, sample_size = N, n_samples_vec)
plot_sampling(df_means_unif,
              paste0("Distribuição amostral da média (Uniforme; N = ", N, ")"),
              x_limits = c(-1.5, 1.5))

População Exponencial (N pequeno)

Código

pop_exp <- rexp(50000, rate = 1)
plot_population(pop_exp, "População: Exponencial(rate = 1)", x_limits = c(0, 8))

Código

N <- 5

df_means_exp_small <- make_sampling_df(pop_exp, sample_size = N, n_samples_vec)
plot_sampling(df_means_exp_small,
              paste0("Distribuição amostral da média (Exponencial; N = ", N, ")"),
              x_limits = c(0, 4))

População Exponencial (N maior)

Código

N <- 30

df_means_exp_large <- make_sampling_df(pop_exp, sample_size = N, n_samples_vec)
plot_sampling(df_means_exp_large,
              paste0("Distribuição amostral da média (Exponencial; N = ", N, ")"),
              x_limits = c(0, 3))

Perspectiva Crítica: Limites e Abusos do CLT

N ≥ 30 É Uma Regra Arbitrária

A famosa “regra” de que precisamos de N ≥ 30 para “invocar o CLT” é excessivamente simplista e pode ser enganadora.

Realidade mais complexa:

A convergência sempre acontece (mesmo com amostras pequenas), mas a velocidade depende da forma da distribuição populacional
Distribuições simétricas (ex: uniforme): Convergência rápida, já com N = 5-10
Distribuições moderadamente assimétricas: N = 20-30 pode ser suficiente
Distribuições extremas (outliers, caudas pesadas, bimodais): Podem requerer N > 100 ou nunca convergir adequadamente

O Abuso do “Invocar o CLT”

Um problema comum em metodologia é usar o CLT como desculpa para ignorar violações graves de pressupostos:

“Os meus dados são terríveis, mas tenho N > 30, logo invoco o CLT!”
Esta atitude transformou o CLT numa espécie de meme metodológico

Limites reais do CLT:

Outliers extremos: Podem dominar a média mesmo com N = 1000
Dependências nos dados: O CLT assume observações independentes. Medidas repetidas, clustering, autocorrelação violam isto
Distribuições com variância infinita: Cauchy, algumas power-law distributions—o CLT não se aplica

Soluções:

Estatística robusta: M-estimadores, medianas aparadas
Modelos apropriados: Mixed models para medidas repetidas
Bootstrap: Simular a distribuição amostral empiricamente

CLT e Modelos Lineares

O CLT aplica-se não só às médias, mas também a outros parâmetros estimados, como declives em regressão e diferenças entre médias.

Para mais detalhes técnicos:

Por outras palavras: Os coeficientes (interceptos, declives) dos nossos modelos lineares também seguem distribuições normais, permitindo teste de hipóteses.

Estatísticas de Teste: Sinal vs Ruído

Agora que compreendemos que as médias amostrais seguem uma distribuição conhecida (graças ao CLT), podemos perguntar: como decidir se um efeito observado é real ou apenas ruído?

A Ideia Central

Genericamente, uma estatística de teste diz-nos o quão distantes estamos do que seria de esperar se não houvesse efeito (hipótese nula).

As estatísticas de teste tendem a ser calculadas a partir da discrepância dum dado sumário dos dados face ao que seria de esperar se não houvesse efeito, ponderando o erro/ruído dos dados.

Informalmente, as estatísticas de teste comparam a magnitude do sinal com a magnitude do ruído, sendo tão maiores (em módulo) quanto mais sinal detectarmos e menos ruído tivermos:

\[\text{Estatística de Teste} = \frac{\text{Sinal (efeito observado)}}{\text{Ruído (erro padrão)}}\]

Quanto maior o valor da estatística de teste, mais evidência temos contra a hipótese nula (a ausência de efeito).

2 A Equação Fundamental

Agora que estabelecemos os fundamentos da inferência estatística, podemos introduzir a equação que unifica todos os modelos lineares.

Quase tudo o que vamos fazer cabe nesta frase:

\(\text{Dados} = \text{Modelo} + \text{Erro}\)

Quando temos uma variável dependente (VD) quantitativa, os dados são a VD. O modelo são as previsões \(\hat{Y}\) e o erro são os resíduos \((e_i = Y_i - \hat{Y}_i)\).

A regra do jogo é:

Queremos menos erro.
Mas “menos erro” quase sempre vem com mais parâmetros (mais complexidade).
Uma forma de pensar no teste F é pensar que pergunta: vale a pena pagar essa complexidade, para comprar o aumento no desempenho?

O Que é o Erro?

O erro que vamos usar aqui é o erro quadrático total:

\(SSE = \sum (Y_i - \hat{Y}_i)^2\)

Porquê o quadrado? Porque assim erros positivos e negativos não se anulam, e penalizamos os erros maiores. A transformação mantém a ordem (é monotónica) mas amplifica as distâncias (\(1^2 = 1; 2^2 = 4, 5^2 = 25, 10^2 = 100\)). Os estimadores da regressão linear são definidos exactamente para minimizar este SSE.

3 O Modelo Mais Simples

O modelo mais simples que podemos propor para prever a VD é o modelo com apenas um intercepto:

\(\hat{Y}_i = \beta_0\)

Aqui \(\beta_0\) é, na prática, a média (ou a média de uma versão centrada).

Vamos usar um exemplo em que “zero” tem significado: o índice Love4Taylor vai de -5 a 5 e o 0 representa indiferença.

library(ggplot2)

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

head(ds)

ID	Love4Taylor
pp01	1.50
pp02	2.15
pp03	1.32
pp04	1.17
pp05	-0.28
pp06	1.72

Código

M <- mean(ds$Love4Taylor)
SD <- sd(ds$Love4Taylor)

ggplot(ds, aes(x = Love4Taylor)) +
geom_histogram(bins = 20) +
geom_vline(xintercept = 0, linetype = "dashed") +
geom_vline(xintercept = M, color = "#F8766D", linetype = "dashed") +
xlim(-5, 5) +
theme_classic() +
labs(x = "Love4Taylor", y = "Frequência",
     subtitle = paste0("M = ", round(M, 2), "; DP = ", round(SD, 2),
                       ". Tracejado preto = 0."))

Compacto vs Aumentado

Modelo compacto (m0): “as pessoas são indiferentes”, ou seja, a previsão é 0.
Modelo aumentado (m1): “vamos estimar um intercepto a partir dos dados”, ou seja, a previsão será a nossa média.

# Modelo aumentado: estima intercepto (média)
m1 <- lm(Love4Taylor ~ 1, data = ds)

summary(m1)


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

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

Se querem ver a lógica de comparação de modelos de forma explícita, comparem:

# `~ 0` define que o modelo não tem intercepto
m0 <- lm(Love4Taylor ~ 0, data = ds)

anova(m0, m1)

Res.Df	RSS	Df	Sum of Sq	F	Pr(>F)
50	74.75060	NA	NA	NA	NA
49	58.77805	1	15.97255	13.31543	0.0006379

O t de Student

Genericamente, um teste-t é um rácio entre uma dada estimativa (e.g., uma média, um declive) e o seu erro de estimação (erro padrão).

\(t = \frac{M}{SE_M}; t = \frac{\beta}{SE_{\beta}}\)

O t de ~~Student~~ William Gosset

O teste-t foi desenvolvido pelo William Gosset que trabalhava para a famosa cervejaria Guiness. Por trabalhar para a Guiness, publicou sobre o pseudónimo de Student e desde então o teste acabou por ficar conhecido com t de Student.

A próxima vez que beberem um pint de Guiness podem pensar no William Gosset e a próxima vez que computarem um teste-t podem beber uma pint de Guiness (ou talvez não…se beberem não computem?…sejam responsáveis, bebam com moderação…).

O F de Snedcor

O teste F pondera a redução do erro ao passarmos do modelo compacto para o aumentado (SSE), penalizando a complexidade adicional do aumentado (\(MSR = \frac{SSE_{m0} - SSE_{m1}}{np_{m1} - np_{m0}}\)), assim como o erro que fica por explicar, considerando os graus de liberdade da estimativa/do erro (\(MSE = \frac{SSE_{m1}}{N - np_{m1}}\)).

Assim:

\(F = \frac{\frac{SSE_{m0} - SSE_{m1}}{np_{m1} - np_{m0}}}{\frac{SSE_{m1}}{N - np_{m1}}} = \frac{MSR}{MSE}\)

Onde:

\(SSE_{m0}\) = soma do erro quadrado do modelo compacto
\(SSE_{m1}\) = soma do erro quadrado do modelo aumentado
\(np_{m0}\) = número de parâmetros do modelo compacto
\(np_{m1}\) = número de parâmetros do modelo aumentado

Código

SSE_m0 <- sum(resid(m0)^2)
SSE_m1 <- sum(resid(m1)^2)
np_m1 <- length(coef(m1))
np_m0 <- length(coef(m0))
MSR <- (SSE_m0 - SSE_m1) / (np_m1 - np_m0)
MSE <- SSE_m1 / (nrow(ds) - np_m1)

c(SSE_m0 = SSE_m0, SSE_m1 = SSE_m1, np_m0 = np_m0, np_m1 = np_m1,
  MSR = MSR, MSE = MSE, F = MSR / MSE)

   SSE_m0    SSE_m1     np_m0     np_m1       MSR       MSE         F 
74.750600 58.778048  0.000000  1.000000 15.972552  1.199552 13.315431

Se quiserem ver a ligação com o teste-t: quando o modelo difere do compacto por um parâmetro, acontece a equivalência clássica:

F com 1 grau de liberdade no numerador/factor é o quadrado de t.

Código

t_val <- summary(m1)$coefficients[1, "t value"]
F_val <- anova(m0, m1)$F[2]

c(t = as.numeric(t_val), F = as.numeric(F_val), t2 = as.numeric(t_val^2))

        t         F        t2 
 3.649032 13.315431 13.315431

4 Regressão e Correlação

Agora passamos ao caso “clássico”: prever a VD com uma VI quantitativa. Vamos usar o exemplo Manatees: mortes de peixes-boi (ManateeDeaths) em função do número de barcos registados (Powerboats).

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

head(ds)

Year	ManateeDeaths	Powerboats
1982	13	447
1983	21	460
1984	24	481
1985	16	498
1986	24	513
1987	20	512

Código

ggplot(ds, aes(x = Powerboats, y = ManateeDeaths)) +
geom_point() +
geom_smooth(method = "lm", se = FALSE) +
theme_classic() +
labs(x = "Powerboats (milhares)", y = "ManateeDeaths")

O modelo aumentado é:

\(\hat{Y}_i = \beta_0 + \beta_1 X_i\)

E o modelo compacto é o modelo da média:

\(\hat{Y}_i = \beta_0\)

m0 <- lm(ManateeDeaths ~ 1, data = ds)
m1 <- lm(ManateeDeaths ~ Powerboats, ds)

anova(m0, m1)

Res.Df	RSS	Df	Sum of Sq	F	Pr(>F)
34	29273.886	NA	NA	NA	NA
33	3886.267	1	25387.62	215.5774	0

summary(m1)


Call:
lm(formula = ManateeDeaths ~ Powerboats, data = ds)

Residuals:
    Min      1Q  Median      3Q     Max 
-21.023  -5.645  -0.885   6.522  28.252 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -57.12918    8.05671  -7.091 4.05e-08 ***
Powerboats    0.15237    0.01038  14.683 5.00e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 10.85 on 33 degrees of freedom
Multiple R-squared:  0.8672,    Adjusted R-squared:  0.8632 
F-statistic: 215.6 on 1 and 33 DF,  p-value: 4.998e-16

\(R^2\): Proporção de Redução do Erro

O \(R^2\) quantifica “quanto” do erro total conseguimos reduzir ao passar do modelo compacto para o aumentado:

\(R^2 = 1 - \frac{SSE_{m1}}{SSE_{m0}} = \frac{SSE_{m0} - SSE_{m1}}{SSE_{m0}}\)

Código

SSE_m0 <- sum(resid(m0)^2)
SSE_m1 <- sum(resid(m1)^2)
R2 <- 1 - (SSE_m1 / SSE_m0)

cat("R² =", round(R2, 3), "\n")

R² = 0.867

Código

cat("Reduzimos", round(R2 * 100, 1),
    "% do erro ao incluir Powerboats no modelo.")

Reduzimos 86.7 % do erro ao incluir Powerboats no modelo.

r de Pearson

Muitas vezes o Pearson \(r\) é visto apenas como a raiz de \(R^2\). Mas, fundamentalmente, a correlação de Pearson é o declive da regressão quando as variáveis estão padronizadas (scaled).

Isto significa que se o preditor subir 1 desvio-padrão, a VD sobe exactamente \(r\) desvios-padrão.

# 1. Correlação de Pearson clássica
r_val <- cor(ds$ManateeDeaths, ds$Powerboats)

# 2. Regressão com variáveis padronizadas (Z-scores)
m_scaled <- lm(scale(ManateeDeaths) ~ scale(Powerboats), ds)
beta_std <- coef(m_scaled)[2]

cat("Pearson r =", round(r_val, 3), "\n")

Pearson r = 0.931

cat("Declive Padronizado (Beta) =", round(beta_std, 3), "\n")

Declive Padronizado (Beta) = 0.931

cat("São idênticos!")

São idênticos!

Notem que os valores são idênticos. A correlação é simplesmente uma regressão linear “sem unidades”.

5 Teste-t para Amostras Independentes

Se a VI tem dois níveis, isto é equivalente ao “teste-t para amostras independentes”. Usamos o dataset penguins para comparar o peso entre sexos.

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

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

Esquemas de Codificação

Quando temos uma variável categórica (factor), o R precisa de a transformar em números para ajustar o modelo linear. Existem várias formas de fazer isto, cada uma com interpretações diferentes dos coeficientes.

Link directo para UCLA página sobre codificação de contrastes

`contr.treatment` (Dummy Coding)

Esta é a codificação por defeito do R. Um grupo é a “referência” (codificado como 0) e o(s) outro(s) grupo(s) são codificados como 1.

Interpretação dos coeficientes:

\(\beta_0\) = média do grupo de referência
\(\beta_1\) = diferença entre o grupo codificado como 1 e o grupo de referência

# Usar dummy coding
options(contrasts = c("contr.treatment", "contr.poly"))

m_dummy <- lm(body_mass_g ~ sex, data = ds)
summary(m_dummy)$coefficients

             Estimate Std. Error  t value      Pr(>|t|)
(Intercept) 3862.2727   56.82895 67.96312 1.703018e-196
sexmale      683.4118   80.00868  8.54172  4.897247e-16

Calcular EMMs com `contr.treatment`

b0 <- coef(m_dummy)[1]  # Intercepto
b1 <- coef(m_dummy)[2]  # Declive

# female é a referência (0), male é 1
EMM_female <- b0
EMM_male <- b0 + b1

cat("EMM female (grupo de referência):", round(EMM_female, 2), "\n")

EMM female (grupo de referência): 3862.27

cat("EMM male (referência + declive):", round(EMM_male, 2), "\n")

EMM male (referência + declive): 4545.68

`contr.sum` (Effect/Contrast Coding)

Nesta codificação, os grupos são codificados como -1, 0, +1 (ou fracções quando k > 2), de forma que a soma dos códigos seja zero.

Interpretação dos coeficientes:

\(\beta_0\) = média global (grande média, não condicionada aos grupos)
\(\beta_1\) = desvio do primeiro grupo em relação à média global

# Usar contrast coding (soma)
options(contrasts = c("contr.sum", "contr.poly"))

m_sum <- lm(body_mass_g ~ sex, data = ds)
summary(m_sum)$coefficients

             Estimate Std. Error   t value      Pr(>|t|)
(Intercept) 4203.9786   40.00434 105.08806 2.647312e-256
sex1        -341.7059   40.00434  -8.54172  4.897247e-16

Calcular EMMs com `contr.sum`

b0 <- coef(m_sum)[1]  # Intercepto (média global)
b1 <- coef(m_sum)[2]  # Desvio do primeiro grupo

# Com contr.sum, o primeiro grupo (female) é -1, o segundo (male) é +1
# Mas R só estima o desvio do primeiro grupo, o do segundo é -b1
EMM_female <- b0 + b1
EMM_male <- b0 - b1

cat("Média global (intercepto):", round(b0, 2), "\n")

Média global (intercepto): 4203.98

cat("EMM female (média + desvio):", round(EMM_female, 2), "\n")

EMM female (média + desvio): 3862.27

cat("EMM male (média - desvio):", round(EMM_male, 2), "\n")

EMM male (média - desvio): 4545.68

Porquê Usar `contr.sum`?

Para modelos com interacções e ANOVAs com múltiplos factores, contr.sum é essencial:

Ortogonalidade: Os contrastes são ortogonais, o que significa que os efeitos principais e interacções são testados de forma independente.
Tipo III Somas de Quadrados: As ANOVAs Tipo III (que testam cada efeito controlando para todos os outros) requerem contrastes ortogonais para produzir os resultados esperados.
Interpretação dos efeitos principais: Em modelos com interacções, os coeficientes dos efeitos principais representam médias marginais (médias colapsando sobre os níveis do outro factor).

Regra Prática

Para a grande maioria das análises em psicologia/ciências sociais, usem contr.sum:

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

Coloquem isto no início do vosso script, antes de ajustar qualquer modelo com factores.

Verificar com `emmeans`

A boa notícia é que não precisamos de fazer os cálculos manualmente. O pacote emmeans faz isto por nós, independentemente da codificação usada:

# Médias marginais estimadas (funciona com qualquer codificação)
emmeans(m_sum, ~ 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

Porquê reportar EMMs em vez de coeficientes?

Evita erros: Os cálculos manuais são tediosos e propensos a erros, especialmente com k > 2 ou interacções.
Interpretação directa: EMMs são as médias (ajustadas) de cada grupo, que é o que queremos comunicar.
Consistência: EMMs são as mesmas independentemente da codificação usada, facilitando a comparação entre estudos.

Comparação de Modelos

Agora que entendemos as codificações, vejamos a comparação de modelos:

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

anova(m0, m_sum)

Res.Df	RSS	Df	Sum of Sq	F	Pr(>F)
332	215259666	NA	NA	NA	NA
331	176380769	1	38878897	72.96099	0

E a equivalência \(t^2 = F\):

Código

t_sex <- summary(m_sum)$coefficients[2, "t value"]
F_sex <- anova(m0, m_sum)$F[2]

c(t = round(t_sex, 3), F = round(F_sex, 3), t_squared = round(t_sex^2, 3))

        t         F t_squared 
   -8.542    72.961    72.961

6 Regressão Múltipla

Antes de passarmos para factores com k > 2, vejamos como a comparação de modelos funciona quando temos múltiplos preditores quantitativos.

Isto é importante porque:

Mostra como o teste F pode comparar modelos que diferem em mais de um parâmetro (teste omnibus).
Introduz a ideia de controlar estatisticamente para outros preditores (testes parciais).
Prepara-nos conceptualmente para ANOVAs com k > 2, onde um factor categórico corresponde a múltiplos parâmetros (k - 1 contrastes).

Vamos usar o exemplo WaterQualityTesting:

ds <- read.csv("../data/WaterQualityTesting.csv")
colnames(ds) <- c("Sample", "pH", "Temp", "NTU", "Oxygen", "Conductivity")

head(ds)

Sample	pH	Temp	NTU	Oxygen	Conductivity
1	7.25	23.1	4.5	7.8	342
2	7.11	22.3	5.1	6.2	335
3	7.03	21.5	3.9	8.3	356
4	7.38	22.9	3.2	9.5	327
5	7.45	20.7	3.8	8.1	352
6	6.89	23.6	4.6	7.2	320

Teste Omnibus

O teste omnibus pergunta: “Será que o conjunto de todos os preditores reduz significativamente o erro?”

m0 <- lm(Oxygen ~ 1, data = ds)
m_full <- lm(Oxygen ~ pH + Temp + NTU + Conductivity, ds)

anova(m0, m_full)

Res.Df	RSS	Df	Sum of Sq	F	Pr(>F)
499	337.4916	NA	NA	NA	NA
495	104.8176	4	232.6739	274.6999	0

Interpretação: O modelo com os 4 preditores reduz significativamente o erro em comparação com o modelo que só tem a média (\(F(4, 49) = 60.66, p < .001\)).

Finalmente: \(df_1 > 1\)

Repararam? Este é o primeiro exemplo onde o numerador do F tem mais de 1 grau de liberdade (\(df_1 = 4\)).

Até aqui, todos os modelos diferiam por um parâmetro, logo \(df_1 = 1\) e \(F = t^2\).

Agora, com 4 preditores adicionados simultaneamente, \(df_1 = 4\) e já não há equivalência com um teste-t. O F omnibus testa se “algum” dos 4 preditores contribui, mas não diz “qual” ou “quantos”.

ANOVA Tipo III

Agora queremos saber: “Será que cada preditor contribui significativamente, mesmo quando controlamos para os outros?”

Para isso, usamos a ANOVA Tipo III, que compara o modelo completo com modelos que têm todos os preditores excepto um:

library(parameters)

Anova(m_full, type = "III")

	Sum Sq	Df	F value	Pr(>F)
(Intercept)	32.2456898	1	152.2798585	0.000000
pH	21.9911133	1	103.8527517	0.000000
Temp	0.0894279	1	0.4223222	0.516082
NTU	12.9663764	1	61.2335471	0.000000
Conductivity	46.6706267	1	220.4014393	0.000000
Residuals	104.8176469	495	NA	NA

Cada linha corresponde a uma comparação de modelos:

# Modelo sem o Intercepto (mas com todos os declives)
m_no_intercept <- lm(Oxygen ~ 0 + Temp + NTU + Conductivity, ds)
# Modelo sem pH (mas com todos os outros declives)
m_no_pH <- lm(Oxygen ~ Temp + NTU + Conductivity, ds)
m_no_Temp <- lm(Oxygen ~ pH + NTU + Conductivity, ds)
m_no_NTU <- lm(Oxygen ~ pH + Temp + Conductivity, ds)
m_no_Conductivity <- lm(Oxygen ~ pH + Temp + NTU, ds)

# Comparar com modelo completo
anova(m_no_intercept, m_full)

Res.Df	RSS	Df	Sum of Sq	F	Pr(>F)
497	137.6328	NA	NA	NA	NA
495	104.8176	2	32.81513	77.48452	0

anova(m_no_pH, m_full)

Res.Df	RSS	Df	Sum of Sq	F	Pr(>F)
496	126.8088	NA	NA	NA	NA
495	104.8176	1	21.99111	103.8528	0

anova(m_no_Temp, m_full)

Res.Df	RSS	Df	Sum of Sq	F	Pr(>F)
496	104.9071	NA	NA	NA	NA
495	104.8176	1	0.0894279	0.4223222	0.516082

anova(m_no_NTU, m_full)

Res.Df	RSS	Df	Sum of Sq	F	Pr(>F)
496	117.7840	NA	NA	NA	NA
495	104.8176	1	12.96638	61.23355	0

anova(m_no_Conductivity, m_full)

Res.Df	RSS	Df	Sum of Sq	F	Pr(>F)
496	151.4883	NA	NA	NA	NA
495	104.8176	1	46.67063	220.4014	0

Notem: Agora todos os \(df_1 = 1\) (cada preditor adiciona 1 parâmetro), logo cada F é o quadrado do respectivo t.

Podemos ver isto nos coeficientes. O pacote parameters apresenta-os de forma mais legível:

parameters(m_full)

Parameter	Coefficient	SE	CI	CI_low	CI_high	t	df_error	p
(Intercept)	-19.6351392	1.5911560	0.95	-22.7613915	-16.5088869	-12.3401725	495	0.000000
pH	2.6264210	0.2577243	0.95	2.1200526	3.1327894	10.1908170	495	0.000000
Temp	-0.0160866	0.0247539	0.95	-0.0647223	0.0325490	-0.6498632	495	0.516082
NTU	-0.4179501	0.0534109	0.95	-0.5228901	-0.3130101	-7.8251867	495	0.000000
Conductivity	0.0328333	0.0022116	0.95	0.0284880	0.0371786	14.8459233	495	0.000000

Reparem: as tabelas estão a contar a mesma história! A ANOVA Tipo III testa cada preditor controlando para os outros (comparação de modelos), e os coeficientes mostram a direcção e magnitude desses efeitos (com sinal e em unidades originais).

Os valores-t ao quadrado correspondem exactamente aos F da tabela ANOVA Tipo III. Por exemplo, para NTU:

Código

t_NTU <- parameters(m_full)$t[parameters(m_full)$Parameter == "NTU"]
F_NTU <- Anova(m_full, type = "III")["NTU", "F value"]

c(t = round(t_NTU, 3), F = round(F_NTU, 3), t_squared = round(t_NTU^2, 3))

        t         F t_squared 
   -7.825    61.234    61.234

O Que Significa “Controlar Para”?

Quando dizemos “controlar estatisticamente para X”, estamos simplesmente a perguntar: “Será que Y contribui para reduzir o erro MESMO quando X já está no modelo?”

Não há magia aqui—é apenas comparação de modelos:

Modelo sem Y (mas com X): Oxygen ~ pH + Temp + Conductivity
Modelo com Y e X: Oxygen ~ pH + Temp + Conductivity + NTU

Se o modelo com Y reduz significativamente o erro, dizemos que “Y tem um efeito significativo, controlando para X”.

Tipo I vs Tipo II vs Tipo III

Existem três formas principais de calcular as somas de quadrados em ANOVAs:

Tipo I (Sequential):

Testa cada preditor na ordem em que foi adicionado ao modelo.
O efeito de A é testado contra o modelo nulo.
O efeito de B é testado contra o modelo com apenas A.
Problema: O resultado depende da ordem dos preditores na fórmula.

Tipo II:

Testa cada efeito principal controlando para os outros efeitos principais, mas não para interacções que envolvam esse efeito.
Apropriado quando não há interacções significativas.

Tipo III:

Testa cada efeito controlando para todos os outros efeitos E interacções.
Apropriado quando há (ou pode haver) interacções.
Requer contrastes ortogonais (e.g., contr.sum) para funcionar correctamente.

Recomendação: Para modelos com interacções (ou que possam ter interacções), usem Tipo III com contr.sum. Caso contrário, Tipo II é mais parcimonioso.

7 One-Way ANOVA (k > 2)

Com mais de dois níveis, aparece a grande diferença prática:

O teste F é omnibus—diz: “é altamente improvável que NÃO haja diferenças algures”—MAS não diz “onde”.

Vamos usar as espécies de pinguins (k = 3):

ds <- penguins
ds <- subset(ds, complete.cases(ds))

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

m0 <- lm(body_mass_g ~ 1, data = ds)
m1 <- lm(body_mass_g ~ species, ds)

anova(m0, m1)

Res.Df	RSS	Df	Sum of Sq	F	Pr(>F)
332	215259666	NA	NA	NA	NA
330	70069447	2	145190219	341.8949	0

Novamente: \(df_1 > 1\), logo \(F \neq t^2\)

Com k = 3 grupos, precisamos de k - 1 = 2 contrastes para representar o factor species.

Portanto, \(df_1 = 2\) e não há equivalência simples com um teste-t.

O F omnibus testa se “há diferenças entre os 3 grupos”, mas não especifica quais pares diferem.

Coeficientes: Contrastes, Não Efeitos

Se olharmos para os coeficientes estimados:

parameters(m1)

Parameter	Coefficient	SE	CI	CI_low	CI_high	t	df_error
(Intercept)	4177.2299	26.58561	0.95	4124.9312	4229.5285	157.12374	330
species1	-471.0655	34.51912	0.95	-538.9708	-403.1602	-13.64651	330
species2	-444.1416	41.80474	0.95	-526.3790	-361.9042	-10.62419	330

Vemos 2 declives (para k = 3 grupos). Estes representam os contrastes codificados (variáveis numéricas “falsas” criadas pelo esquema de codificação), não o “efeito” substantivo da variável species.

Estes declives são:

species1: desvio de Adelie em relação à média global
species2: desvio de Chinstrap em relação à média global

(O desvio de Gentoo é inferido: \(-\text{species1} - \text{species2}\))

Estes coeficientes são interpretáveis, mas não são triviais e raramente correspondem ao que queremos comunicar (“Adelie é X gramas mais pesado que Chinstrap”).

Solução: Se o F omnibus for significativo, passamos para comparações múltiplas (pairwise) usando emmeans, que dá as diferenças directas entre pares de grupos—muito mais útil para comunicar resultados!

Como Funcionam os Graus de Liberdade?

Quando temos um factor com k níveis, precisamos de k - 1 contrastes (declives) para o representar no modelo:

Modelo aumentado (\(m_1\)): \(\beta_0\) (intercepto) + \(\beta_1, \beta_2, ..., \beta_{k-1}\) (contrastes) = k parâmetros
Modelo compacto (\(m_0\)): \(\beta_0\) (intercepto) = 1 parâmetro

Portanto:

\(df_1\) (numerador do F) = \(np_{m1} - np_{m0} = k - 1\)
\(df_2\) (denominador do F) = \(N - np_{m1} = N - k\)

No nosso exemplo com espécies (k = 3):

Código

k <- nlevels(ds$species)
N <- nrow(ds)

cat("k (número de níveis) =", k, "\n")

k (número de níveis) = 3

Código

cat("N (número de observações) =", N, "\n")

N (número de observações) = 333

Código

cat("df1 = k - 1 =", k - 1, "\n")

df1 = k - 1 = 2

Código

cat("df2 = N - k =", N - k, "\n")

df2 = N - k = 330

Porquê k - 1 contrastes?

Com k = 3 grupos, se soubermos as médias de 2 grupos E a média global, podemos deduzir a média do 3º grupo. Por isso, só precisamos de estimar k - 1 parâmetros para representar um factor com k níveis.

Comparações Múltiplas

Como o teste omnibus não diz onde estão as diferenças, precisamos de comparações múltiplas:

\(\text{número de pares} = \frac{k(k-1)}{2}\)

emmeans(m1, pairwise ~ species)

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

8 Múltiplos Factores: Interacções

Em modelos complexos, cada linha da tabela da ANOVA (Tipo III) corresponde a uma comparação de modelos: o modelo completo vs o modelo sem aquele termo específico.

Quando temos dois (ou mais) factores, podemos testar:

Efeito principal de A: A contribui, controlando para B e A:B?
Efeito principal de B: B contribui, controlando para A e A:B?
Interacção A:B: A interacção contribui, controlando para A e B?

Sintaxe de Interacções no R

Na sintaxe de fórmulas do R:

A * B expande para A + B + A:B (efeitos principais + interacção)
A + B são apenas efeitos aditivos (sem interacção)
A:B é apenas a interacção (raramente usado sozinho)

Para modelos com múltiplas interacções:

A * B * C = A + B + C + A:B + A:C + B:C + A:B:C
(A + B + C)^2 = todas as interacções duplas, mas não a tripla

Vamos usar um exemplo 2 × 3 com sex * species:

ds <- penguins
ds <- subset(ds, complete.cases(ds))

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

m_add <- lm(body_mass_g ~ sex + species, data = ds)
m_int <- lm(body_mass_g ~ sex * species, ds)

# Teste da interacção
anova(m_add, m_int)

Res.Df	RSS	Df	Sum of Sq	F	Pr(>F)
329	32979185	NA	NA	NA	NA
327	31302628	2	1676557	8.756997	0.0001973

A tabela ANOVA Tipo III mostra todos os testes simultaneamente:

Anova(m_int, type = "III")

	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

parameters(m_int)

Parameter	Coefficient	SE	CI	CI_low	CI_high	t	df_error	p
(Intercept)	4173.84711	17.85231	0.95	4138.72724	4208.96699	233.7986911	327	0.0000000
sex1	-315.25282	17.85231	0.95	-350.37269	-280.13295	-17.6589355	327	0.0000000
species1	-467.68273	23.17885	0.95	-513.28121	-422.08425	-20.1771319	327	0.0000000
species2	-440.75888	28.07051	0.95	-495.98045	-385.53731	-15.7018480	327	0.0000000
sex1:species1	-22.07595	23.17885	0.95	-67.67443	23.52253	-0.9524176	327	0.3415887
sex1:species2	109.37047	28.07051	0.95	54.14890	164.59204	3.8962765	327	0.0001185

O Que É uma Interacção?

Uma interacção ocorre quando o efeito de A na VD depende do nível de B (e vice-versa).

Por exemplo:

Machos são mais pesados que fêmeas em todas as espécies → sem interacção
Machos são mais pesados que fêmeas em Gentoos, mas não há diferença em Adelies → interacção!

Visualizar a Interacção

Vamos criar uma função auxiliar para fazer os gráficos de interacção:

# Função auxiliar para gráficos de interacção
# Se 'data' estiver vazio ou sem linhas, não adiciona camada de jitter
plot_interaction <- function(data, emms, x_var, group_var, y_var, title) {
    p <- ggplot(emms, aes(x = .data[[x_var]], y = emmean,
                          color = .data[[group_var]], shape = .data[[group_var]]))

    # Adicionar jitter apenas se houver dados brutos
    if (!is.null(data) && nrow(data) > 0) {
        p <- p + geom_jitter(data = data, aes(y = .data[[y_var]]),
                             alpha = 0.15, width = 0.15)
    }

    p <- p +
        geom_point(size = 3) +
        geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.2) +
        geom_line(aes(group = .data[[group_var]])) +
        theme_classic() +
        labs(x = x_var, y = y_var, title = title)

    return(p)
}

Código

# Médias marginais estimadas para todas as células
emm_cells <- as.data.frame(emmeans(m_int, ~ sex * species))

# Gráfico 1: x = species, linhas = sex
plot_interaction(ds, emm_cells, "species", "sex", "body_mass_g",
                 "Espécie por Sexo")

Código

# Gráfico 2: x = sex, linhas = species
plot_interaction(ds, emm_cells, "sex", "species", "body_mass_g",
                 "Sexo por Espécie")

Visualizem SEMPRE as Interacções

Mesmo quando “só querem” saber os números, façam um gráfico antes de interpretar uma interacção. Uma interacção mal visualizada dá origem a textos muito confiantes… e muito errados.

Decompor a Interacção

Quando a interacção é significativa, a pergunta útil não é “há efeito principal de sexo?” mas sim:

Diferenças entre espécies dentro de cada sexo (species | sex)
Diferenças entre sexos dentro de cada espécie (sex | species)

emmeans(m_int, 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

Padrões de Interacção (2 × 2)

Vejamos agora todos os padrões possíveis em delineamentos factoriais 2 × 2, para percebermos como interacções e efeitos principais se relacionam.

Vamos usar exemplos hipotéticos de tempo de resposta (RT) em duas condições (Control vs Experimental) e duas tarefas (Reading vs Writing).

Apenas Efeito Principal de Condição

Código

emms <- data.frame(Condition = c("Control", "Experimental", "Control", "Experimental"),
                   Task = c("Reading", "Reading", "Writing", "Writing"),
                   emmean = c(4, 4.1, 8, 8.1),
                   lower.CL = c(3.5, 3.6, 7.5, 7.6),
                   upper.CL = c(4.5, 4.6, 8.5, 8.6))

plot_interaction(NULL, emms, "Condition", "Task", "emmean",
                 "Condição por Tarefa") + ylim(0, 10)

Código

plot_interaction(NULL, emms, "Task", "Condition", "emmean",
                 "Tarefa por Condição") + ylim(0, 10)

Padrão:

Efeito principal de Condição
Efeito principal de Tarefa
Interacção

Interpretação: O grupo de controlo é mais rápido que o experimental, independentemente da tarefa.

Apenas Efeito Principal de Tarefa

Código

emms$emmean <- c(4, 8, 4.1, 8.1)
emms$lower.CL <- c(3.5, 7.5, 3.6, 7.6)
emms$upper.CL <- c(4.5, 8.5, 4.6, 8.6)

plot_interaction(NULL, emms, "Condition", "Task", "emmean",
                 "Condição por Tarefa") + ylim(0, 10)

Código

plot_interaction(NULL, emms, "Task", "Condition", "emmean",
                 "Tarefa por Condição") + ylim(0, 10)

Padrão:

Efeito principal de Condição
Efeito principal de Tarefa
Interacção

Interpretação: Leitura é mais rápida que escrita, independentemente do grupo.

Apenas Interacção (Crossover)

Código

emms$emmean <- c(8, 4, 4, 8)
emms$lower.CL <- c(7.5, 3.5, 3.5, 7.5)
emms$upper.CL <- c(8.5, 4.5, 4.5, 8.5)

plot_interaction(NULL, emms, "Condition", "Task", "emmean",
                 "Condição por Tarefa") + ylim(0, 10)

Código

plot_interaction(NULL, emms, "Task", "Condition", "emmean",
                 "Tarefa por Condição") + ylim(0, 10)

Padrão:

Efeito principal de Condição
Efeito principal de Tarefa
Interacção

Interpretação: Controlo é rápido na escrita e lento na leitura; Experimental é rápido na leitura e lento na escrita. As linhas cruzam-se (crossover interaction).

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

Código

emms$emmean <- c(2, 5, 8, 5)
emms$lower.CL <- c(1.5, 4.5, 7.5, 4.5)
emms$upper.CL <- c(2.5, 5.5, 8.5, 5.5)

plot_interaction(NULL, emms, "Condition", "Task", "emmean",
                 "Condição por Tarefa") + ylim(0, 10)

Código

plot_interaction(NULL, emms, "Task", "Condition", "emmean",
                 "Tarefa por Condição") + ylim(0, 10)

Padrão:

Efeito principal de Condição
Efeito principal de Tarefa
Interacção

Interpretação: Controlo é mais rápido que Experimental apenas na leitura, não há diferença na escrita. As linhas não são paralelas.

Efeito Principal de Tarefa + Interacção

Código

emms$emmean <- c(2, 8, 5, 5)
emms$lower.CL <- c(1.5, 7.5, 4.5, 4.5)
emms$upper.CL <- c(2.5, 8.5, 5.5, 5.5)

plot_interaction(NULL, emms, "Condition", "Task", "emmean",
                 "Condição por Tarefa") + ylim(0, 10)

Código

plot_interaction(NULL, emms, "Task", "Condition", "emmean",
                 "Tarefa por Condição") + ylim(0, 10)

Padrão:

Efeito principal de Condição
Efeito principal de Tarefa
Interacção

Interpretação: Há diferenças entre tarefas apenas no grupo de controlo, não no experimental.

Todos os Efeitos Significativos

Código

emms$emmean <- c(2, 8, 2, 2)
emms$lower.CL <- c(1.5, 7.5, 1.5, 1.5)
emms$upper.CL <- c(2.5, 8.5, 2.5, 2.5)

plot_interaction(NULL, emms, "Condition", "Task", "emmean",
                 "Condição por Tarefa") + ylim(0, 10)

Código

plot_interaction(NULL, emms, "Task", "Condition", "emmean",
                 "Tarefa por Condição") + ylim(0, 10)

Padrão:

Efeito principal de Condição
Efeito principal de Tarefa
Interacção

Interpretação: Há diferenças entre condições apenas na escrita. Há diferenças entre tarefas apenas no controlo. Experimental demora o mesmo em ambas as tarefas (igual ao tempo que controlo demora na leitura).

Lição Importante

Estes padrões mostram que efeitos principais podem ser enganadores quando há interacção.

Sempre que a interacção for significativa, interpretem-na primeiro e contextualizem os efeitos principais à luz da interacção.

9 ANCOVA: Factores + Covariáveis

ANCOVA (Analysis of Covariance) é simplesmente um modelo que mistura preditores categóricos (factores) e contínuos (covariáveis). Na prática, é uma regressão múltipla com tipos mistos de preditores.

Vamos usar o exemplo dos pinguins: será que o efeito de species no peso se mantém quando controlamos para flipper_length_mm (comprimento da barbatana)?

Modelo Sem Covariável

Primeiro, o modelo só com o factor:

ds <- penguins
ds <- subset(ds, complete.cases(ds))

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

m_factor <- lm(body_mass_g ~ species, data = ds)

Anova(m_factor, type = "III")

	Sum Sq	Df	F value	Pr(>F)
(Intercept)	5242015989	1	24687.8683	0
species	145190219	2	341.8949	0
Residuals	70069447	330	NA	NA

parameters(m_factor)

Parameter	Coefficient	SE	CI	CI_low	CI_high	t	df_error
(Intercept)	4177.2299	26.58561	0.95	4124.9312	4229.5285	157.12374	330
species1	-471.0655	34.51912	0.95	-538.9708	-403.1602	-13.64651	330
species2	-444.1416	41.80474	0.95	-526.3790	-361.9042	-10.62419	330

emmeans(m_factor, pairwise ~ species)

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

ANCOVA: Factor + Covariável (Sem Interacção)

A tradição em ANCOVA é incluir a covariável sem interacção com o factor:

m_ancova <- lm(body_mass_g ~ species + flipper_length_mm, ds)

Anova(m_ancova, type = "III")

	Sum Sq	Df	F value	Pr(>F)
(Intercept)	5770339	1	41.41167	0
species	5368818	2	19.26505	0
flipper_length_mm	24226302	1	173.86359	0
Residuals	45843144	329	NA	NA

parameters(m_ancova)

Parameter	Coefficient	SE	CI	CI_low	CI_high	t	df_error	p
(Intercept)	-3986.79618	619.530421	0.95	-5205.53685	-2768.05551	-6.4351903	329	0.0000000
species1	-26.38271	43.809840	0.95	-112.56546	59.80004	-0.6022096	329	0.5474494
species2	-231.75819	37.500785	0.95	-305.52976	-157.98662	-6.1800890	329	0.0000000
flipper_length_mm	40.60617	3.079553	0.95	34.54807	46.66426	13.1857342	329	0.0000000

emmeans(m_ancova, pairwise ~ species)

$emmeans
 species   emmean   SE  df lower.CL upper.CL
 Adelie      4147 45.5 329     4058     4237
 Chinstrap   3942 48.0 329     3848     4036
 Gentoo      4432 60.7 329     4312     4551

Confidence level used: 0.95 

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

Interpretação:

O efeito de species mantém-se significativo, mesmo controlando para flipper_length_mm.
As EMMs são agora ajustadas para o comprimento médio da barbatana (é como se comparássemos pinguins de todas as espécies “com barbatanas do mesmo tamanho”).
Por cada mm adicional de barbatana, o peso aumenta ~49g, independentemente da espécie (porque não há interacção).

E Se Houver Interacção?

Mas será que a relação entre comprimento da barbatana e peso é a mesma em todas as espécies? Talvez não! Podemos testar incluindo a interacção:

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

Anova(m_ancova_int, type = "III")

	Sum Sq	Df	F value	Pr(>F)
(Intercept)	5663696	1	41.720187	0.0000000
species	1229661	2	4.528993	0.0114771
flipper_length_mm	22649648	1	166.842898	0.0000000
species:flipper_length_mm	1451475	2	5.345963	0.0051931
Residuals	44391669	327	NA	NA

parameters(m_ancova_int)

Parameter	Coefficient	SE	CI	CI_low	CI_high	t	df_error	p
(Intercept)	-4073.162597	630.606760	0.95	-5313.72066	-2832.6045354	-6.459117	327	0.0000000
species1	1565.074854	814.329012	0.95	-36.90991	3167.0596176	1.921920	327	0.0554846
species2	1035.966823	952.637242	0.95	-838.10413	2910.0377753	1.087473	327	0.2776289
flipper_length_mm	40.475910	3.133594	0.95	34.31136	46.6404577	12.916768	327	0.0000000
species1:flipper_length_mm	-7.786996	4.142046	0.95	-15.93541	0.3614234	-1.879988	327	0.0609986
species2:flipper_length_mm	-5.902516	4.806025	0.95	-15.35715	3.5521135	-1.228149	327	0.2202742

A interacção é significativa! Isto significa que o “declive” (efeito de flipper_length_mm no peso) difere entre espécies.

Podemos visualizar:

Código

ggplot(ds, aes(x = flipper_length_mm, y = body_mass_g,
               color = species, shape = species)) +
geom_point(alpha = 0.3) +
geom_smooth(method = "lm", se = FALSE) +
theme_classic() +
    labs(x = "Comprimento da Barbatana (mm)", y = "Peso (g)",
         title = "Declives Diferentes por Espécie (Interacção Significativa)")

Quando a interacção é significativa, as EMMs dependem do valor da covariável. Uma escolha comum é olhar para valores “típicos” (média ± 1 DP):

fl_mean <- mean(ds$flipper_length_mm)
fl_sd <- sd(ds$flipper_length_mm)
fl_vals <- fl_mean + c(-1, 0, 1) * fl_sd

# Comparações entre espécies para flipper baixo / médio / alto
emmeans(m_ancova_int, pairwise ~ species,
        at = list(flipper_length_mm = fl_vals))

$emmeans
 species   emmean   SE  df lower.CL upper.CL
 Adelie      4061 59.4 327     3944     4178
 Chinstrap   3911 55.2 327     3802     4020
 Gentoo      4211 90.3 327     4034     4389

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

$contrasts
 contrast           estimate    SE  df t.ratio p.value
 Adelie - Chinstrap      150  81.1 327   1.854  0.1539
 Adelie - Gentoo        -150 108.0 327  -1.387  0.3489
 Chinstrap - Gentoo     -300 106.0 327  -2.836  0.0134

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

# Diferenças de declives (slopes) entre espécies
emtrends(m_ancova_int, pairwise ~ species, var = "flipper_length_mm")

$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

Mostrar as EMMs (média ± DP)

emm_typical <- as.data.frame(
    emmeans(m_ancova_int, ~ species | flipper_length_mm,
            at = list(flipper_length_mm = fl_vals))
)

ggplot(emm_typical,
       aes(x = flipper_length_mm, y = emmean,
           color = species, shape = species)) +
    geom_point(size = 2.5) +
    geom_line(aes(group = species)) +
    geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0) +
    theme_classic() +
    labs(x = "flipper_length_mm (valores típicos)",
         y = "body_mass_g (EMM)",
         title = "EMMs por espécie em flipper baixo/médio/alto")

Mostrar os declives (o que `emtrends()` testa)

tr_species <- as.data.frame(emtrends(m_ancova_int, ~ species, var = "flipper_length_mm"))

slope_col <- grep("\\.trend$", names(tr_species), value = TRUE)[1]

ggplot(tr_species, aes(x = species, y = .data[[slope_col]], color = species)) +
    geom_point(size = 2.5) +
    geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.15) +
    theme_classic() +
    labs(x = "Espécie", y = "Declive (g por mm)",
         title = "Declives estimados por espécie (emtrends)") +
    guides(color = "none")

Notas Importantes sobre ANCOVA

Não Adicionem Covariáveis Levianamente

“Controlar para” não é mágico—é apenas adicionar ao modelo.

Se um efeito deixa de ser significativo quando adicionamos uma covariável, isso não significa que o efeito original era “falso” ou “inútil”. Apenas significa que:

O “tamanho” do efeito é menor quando ajustamos para a covariável.
Na prática, se a covariável estiver disponível, usá-la pode ser uma forma mais eficiente de prever a VD.
Mas teoricamente, o facto de o efeito ser significativo por si só ainda pode ser relevante para a teoria em teste.

Usem o cérebro. Evitem cargo-cult statistics.

Quando Incluir Interacções?

A tradição em ANCOVA é não incluir interacções entre factores e covariáveis. Mas:

Incluímos interacções entre factores por defeito (e.g., sex * species).
Então porquê não incluir species * flipper_length_mm? Apenas por tradição?

Pensem criticamente: Se há razão teórica para crer que a relação entre a covariável e a VD pode diferir entre grupos, testem a interacção.

Nota prática:

Cada variável quantitativa adiciona 1 parâmetro ao modelo (não contando interacções).
Cada variável categórica com k níveis adiciona k - 1 parâmetros.

Por isso, adicionar muitos factores (especialmente com muitos níveis) e todas as suas interacções pode tornar o modelo demasiado complexo. Usem o cérebro: incluíam o que faz sentido teoricamente e testável com os vossos dados.

10 Regressão Múltipla Quantitativa com Interacções

Quando ambos os preditores são quantitativos, uma interacção significa que o “efeito” (declive) de uma variável depende do nível da outra.

Vamos voltar ao exemplo WaterQuality e testar se o efeito da turbidez (NTU) no oxigénio depende da temperatura (Temp).

Ajustar e Comparar Modelos

library(car)
library(emmeans)

# Recarregar dataset

ds <- read.csv("../data/WaterQualityTesting.csv")
colnames(ds) <- c("Sample", "pH", "Temp", "NTU", "Oxygen", "Conductivity")

# Modelo aditivo (sem interacção)
m_add <- lm(Oxygen ~ pH + Temp + NTU + Conductivity, data = ds)

# Modelo com interacção (Temp × NTU)
m_int <- lm(Oxygen ~ pH + Temp * NTU + Conductivity, data = ds)

anova(m_add, m_int)

Res.Df	RSS	Df	Sum of Sq	F	Pr(>F)
495	104.8176	NA	NA	NA	NA
494	104.5080	1	0.3096057	1.463478	0.2269561

Anova(m_int, type = "III")

	Sum Sq	Df	F value	Pr(>F)
(Intercept)	3.8389747	1	18.146484	0.0000245
pH	22.1839190	1	104.861366	0.0000000
Temp	0.2770492	1	1.309586	0.2530237
NTU	0.1753988	1	0.829094	0.3629787
Conductivity	44.1001261	1	208.457282	0.0000000
Temp:NTU	0.3096057	1	1.463478	0.2269561
Residuals	104.5080412	494	NA	NA

Decompor a Interacção com `emtrends()`

Se a interacção for significativa, uma forma útil de interpretar é olhar para declives simples: qual é o declive de NTU em valores baixos/médios/altos de Temp (média ± 1 DP)?

temp_mean <- mean(ds$Temp)
temp_sd <- sd(ds$Temp)
temp_vals <- temp_mean + c(-1, 0, 1) * temp_sd

# Estimar o declive de NTU em Temp baixo / médio / alto
emtrends(m_int, ~ Temp, var = "NTU",
         at = list(Temp = temp_vals))

 Temp NTU.trend     SE  df lower.CL upper.CL
 21.2    -0.334 0.0878 494   -0.506   -0.161
 22.1    -0.402 0.0549 494   -0.510   -0.294
 23.0    -0.471 0.0691 494   -0.607   -0.335

Confidence level used: 0.95

# Comparar estes declives entre si (diferenças de declives)
emtrends(m_int, pairwise ~ Temp, var = "NTU",
         at = list(Temp = temp_vals))

$emtrends
 Temp NTU.trend     SE  df lower.CL upper.CL
 21.2    -0.334 0.0878 494   -0.506   -0.161
 22.1    -0.402 0.0549 494   -0.510   -0.294
 23.0    -0.471 0.0691 494   -0.607   -0.335

Confidence level used: 0.95 

$contrasts
 contrast                                    estimate     SE  df t.ratio
 Temp21.1512765843986 - Temp22.0544            0.0687 0.0568 494   1.210
 Temp21.1512765843986 - Temp22.9575234156014   0.1374 0.1140 494   1.210
 Temp22.0544 - Temp22.9575234156014            0.0687 0.0568 494   1.210
 p.value
  0.4479
  0.4479
  0.4479

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

Visualizar o que os declives simples representam

temp_lab <- c("Temp baixa (M-1DP)", "Temp média (M)", "Temp alta (M+1DP)")

# Curvas preditas: Oxygen ~ NTU para Temp baixa/média/alta (mantendo os restantes preditores em valores típicos)
newgrid <- expand.grid(
    Temp = temp_vals,
    NTU = seq(quantile(ds$NTU, 0.05), quantile(ds$NTU, 0.95), length.out = 60),
    pH = mean(ds$pH),
    Conductivity = mean(ds$Conductivity)
)

newgrid$Temp_level <- factor(newgrid$Temp, levels = temp_vals, labels = temp_lab)
newgrid$Oxygen_hat <- predict(m_int, newdata = newgrid)

ggplot(newgrid, aes(x = NTU, y = Oxygen_hat, color = Temp_level)) +
    geom_line(linewidth = 0.9) +
    theme_classic() +
    labs(x = "NTU", y = "Oxygen (predito)",
         color = "Temperatura",
         title = "Declives simples: efeito de NTU em diferentes níveis de Temp")

Visualizar os declives estimados por `emtrends()`

sl_ntu_by_temp <- as.data.frame(emtrends(m_int, ~ Temp, var = "NTU", at = list(Temp = temp_vals)))

sl_ntu_by_temp$Temp_level <- factor(sl_ntu_by_temp$Temp, levels = temp_vals, labels = temp_lab)

slope_col <- grep("\\.trend$", names(sl_ntu_by_temp), value = TRUE)[1]

ggplot(sl_ntu_by_temp,
       aes(x = Temp_level, y = .data[[slope_col]], color = Temp_level)) +
    geom_point(size = 2.5) +
    geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.12) +
    theme_classic() +
    labs(x = "Temp (valores típicos)", y = "Declive de NTU (Oxygen por NTU)",
         title = "Declives simples estimados (emtrends)") +
    guides(color = "none")

Podemos inverter a pergunta: qual é o declive de Temp em valores baixos/médios/altos de NTU?

ntu_mean <- mean(ds$NTU)
ntu_sd <- sd(ds$NTU)
ntu_vals <- ntu_mean + c(-1, 0, 1) * ntu_sd

emtrends(m_int, ~ NTU, var = "Temp",
         at = list(NTU = ntu_vals))

  NTU Temp.trend     SE  df lower.CL upper.CL
 3.77     0.0189 0.0381 494  -0.0559   0.0937
 4.17    -0.0113 0.0251 494  -0.0606   0.0379
 4.57    -0.0416 0.0325 494  -0.1054   0.0223

Confidence level used: 0.95

ntu_lab <- c("NTU baixa (M-1DP)", "NTU média (M)", "NTU alta (M+1DP)")

sl_temp_by_ntu <- as.data.frame(emtrends(m_int, ~ NTU, var = "Temp", at = list(NTU = ntu_vals)))
sl_temp_by_ntu$NTU_level <- factor(sl_temp_by_ntu$NTU, levels = ntu_vals, labels = ntu_lab)

slope_col <- grep("\\.trend$", names(sl_temp_by_ntu), value = TRUE)[1]

ggplot(sl_temp_by_ntu,
       aes(x = NTU_level, y = .data[[slope_col]], color = NTU_level)) +
    geom_point(size = 2.5) +
    geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL), width = 0.12) +
    theme_classic() +
    labs(x = "NTU (valores típicos)", y = "Declive de Temp (Oxygen por unidade de Temp)",
         title = "Declives simples estimados (Temp) em diferentes NTU") +
    guides(color = "none")

11 Medidas Repetidas e Modelos Mistos

Até aqui, todos os modelos assumiam que cada linha da base de dados representa uma unidade independente. Em delineamentos de medidas repetidas, a mesma pessoa aparece várias vezes (e.g., pré e pós), o que quebra a independência.

A boa notícia é que isto não requer um “tipo diferente” de estatística. A lógica de comparação de modelos mantém-se, mas precisamos de modelar explicitamente a estrutura de agrupamento (“quem é quem”) através de efeitos aleatórios.

Capítulo Dedicado

Dado que medidas repetidas e modelos mistos envolvem decisões conceptuais importantes (estrutura de efeitos aleatórios, debate maximal vs parsimonious, crossed random effects), criámos um capítulo teórico dedicado:

👉 Capítulo 4: Medidas Repetidas e Modelos Mistos

Esse capítulo explica:

Porquê o erro tem estrutura (decomposição em componentes)
Random intercepts e random slopes
Crossed random effects (participantes E itens)
O debate maximal vs parsimonious (Barr et al., Bates et al., Singmann & Kellen)
Como medidas repetidas se encaixam na lógica de comparação de modelos

Para aspectos práticos (sintaxe do R, comparação de pacotes aov() vs afex::aov_4() vs lme4::lmer(), interpretação de output), consultem o tutorial prático sobre medidas repetidas.

12 Exercícios

1) Tradução: “que modelo é este?”

Completai a tabela (podeis usar rascunho em papel):

Nome comum	Fórmula no `lm()`
teste-t (1 amostra)
correlação/regressão simples
teste-t (2 grupos)
one-way ANOVA
regressão múltipla
ANCOVA
moderação/interacção (2 variáveis)

Clique para ver as soluções

Nome comum	Fórmula no `lm()`
teste-t (1 amostra)	`VD ~ 1` (ou `VD - mu0 ~ 1` para testar contra mu0)
correlação/regressão simples	`VD ~ VIN`
teste-t (2 grupos)	`VD ~ VIF` (factor com 2 níveis)
one-way ANOVA	`VD ~ VIF` (factor com k>2)
regressão múltipla	`VD ~ VIN1 + VIN2 + ...`
ANCOVA	`VD ~ VIF + VIN` (ou `VD ~ VIF * VIN` com interacção)
moderação/interacção	`VD ~ A * B`

2) Um teste “parcial” por comparação de modelos

No exemplo WaterQuality, calculai o teste ao preditor NTU por comparação directa:

modelo completo vs modelo sem NTU

Clique para ver as soluções

Código

ds <- read.csv("../data/WaterQualityTesting.csv")
colnames(ds) <- c("Sample", "pH", "Temp", "NTU", "Oxygen", "Conductivity")

m_full <- lm(Oxygen ~ pH + Temp + NTU + Conductivity, data = ds)
m_no_NTU <- lm(Oxygen ~ pH + Temp + Conductivity, data = ds)

anova(m_no_NTU, m_full)

Res.Df	RSS	Df	Sum of Sq	F	Pr(>F)
496	117.7840	NA	NA	NA	NA
495	104.8176	1	12.96638	61.23355	0

3) F vs t²: quando são equivalentes?

Para cada cenário abaixo, indicai se \(F = t^2\) (verdadeiro/falso) e porquê:

Comparar médias de 2 grupos (teste-t independente)
Testar se um declive de regressão é diferente de 0
Testar o efeito de um factor com 4 níveis (one-way ANOVA)
Testar uma interacção 2 × 3
Testar se um preditor quantitativo contribui numa regressão múltipla (controlando para os outros)

Clique para ver as soluções

a) Verdadeiro. Comparar 2 grupos = adicionar 1 parâmetro (1 contraste), logo \(df_1 = 1\) e \(F = t^2\).

b) Verdadeiro. Testar um declive = adicionar 1 parâmetro, logo \(df_1 = 1\) e \(F = t^2\).

c) Falso. Factor com 4 níveis requer k - 1 = 3 contrastes, logo \(df_1 = 3\) e \(F \neq t^2\). O F omnibus testa se “há diferenças algures entre os 4 grupos”, não há um único t correspondente.

d) Falso. Interacção 2 × 3 tem \((2-1) \times (3-1) = 2\) parâmetros, logo \(df_1 = 2\) e \(F \neq t^2\).

e) Verdadeiro. Testar um preditor (controlando para outros) = adicionar 1 parâmetro, logo \(df_1 = 1\) e \(F = t^2\). Isto é exactamente o que a ANOVA Tipo III faz para cada linha (cada preditor).

Regra geral: \(F = t^2\) apenas quando \(df_1 = 1\), ou seja, quando o modelo aumentado adiciona exactamente 1 parâmetro em relação ao compacto.

4) Interpretar contrastes vs EMMs

Usai o exemplo das espécies de pinguins (k = 3) do capítulo.

Ajustai o modelo body_mass_g ~ species com contr.sum.
Olhai para os coeficientes (parameters()). O que representam species1 e species2?
Usai emmeans() para obter as médias de cada espécie. São mais fáceis de interpretar que os coeficientes?
Por que razão raramente reportamos os coeficientes species1 e species2 num artigo?

Clique para ver as soluções

Código

library(palmerpenguins)
ds <- subset(penguins, complete.cases(penguins))

options(contrasts = c("contr.sum", "contr.poly"))
m <- lm(body_mass_g ~ species, data = ds)

b) species1 e species2 representam os desvios das primeiras 2 espécies (Adelie e Chinstrap) em relação à média global de todas as 3 espécies. O desvio de Gentoo é implicitamente \(-(\text{species1} + \text{species2})\).

c) Sim! EMMs dão directamente as médias de cada espécie (em gramas), que é o que queremos comunicar.

Código

library(emmeans)
emmeans(m, ~ 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

d) Porque os coeficientes species1 e species2 não correspondem a comparações teoricamente significativas (“Adelie vs média de todas as espécies” raramente é a questão de investigação). Queremos reportar:

As médias de cada grupo (EMMs)
Comparações directas entre pares de grupos (pairwise)

Ambas são facilmente obtidas com emmeans() e são muito mais interpretáveis.

5) ANCOVA com e sem interacção

Usai o dataset penguins.

Ajustai um ANCOVA para prever body_mass_g com species (factor) e bill_length_mm (covariável), sem interacção.
Agora adicionai a interacção: species * bill_length_mm.
A interacção é significativa? O que isso significa substantivamente?
Se a interacção for significativa, qual modelo deveríeis reportar?

Clique para ver as soluções

Código

library(palmerpenguins)
ds <- subset(penguins, complete.cases(penguins))

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

m_ancova <- lm(body_mass_g ~ species + bill_length_mm, data = ds)
Anova(m_ancova, type = "III")

	Sum Sq	Df	F value	Pr(>F)
(Intercept)	16351.66	1	0.116158	0.733457
species	94153501.57	2	334.421002	0.000000
bill_length_mm	23755815.48	1	168.755139	0.000000
Residuals	46313631.32	329	NA	NA

Código

m_ancova_int <- lm(body_mass_g ~ species * bill_length_mm, data = ds)
Anova(m_ancova_int, type = "III")

	Sum Sq	Df	F value	Pr(>F)
(Intercept)	117242.4	1	0.8475808	0.3579163
species	168316.0	2	0.6084036	0.5448338
bill_length_mm	21276307.8	1	153.8128523	0.0000000
species:bill_length_mm	1081048.1	2	3.9076115	0.0210333
Residuals	45232583.2	327	NA	NA

c) Sim, a interacção é significativa (\(p < .001\)). Isto significa que o declive (relação entre comprimento do bico e peso) difere entre espécies. Ou seja: um mm adicional de bico não tem o mesmo “efeito” no peso em todas as espécies.

d) Se a interacção é significativa, deveríeis reportar o modelo com interacção, pois ele representa melhor os dados. O modelo sem interacção assume declives paralelos, o que é violado quando a interacção existe.

13 Recursos

Para relembrar interceptos e declives de forma interactiva: https://joao-o-santos.shinyapps.io/intercepts_and_slopes
Uma boa visualização interactiva de regressão: https://seeing-theory.brown.edu/regression-analysis/index.html
Obrigatório: Lindeløv (2019) - Common statistical tests are linear models

Sobre cargo-cult statistics: Stark & Saltelli (2018)

0 Introdução

1 Fundamentos da Inferência Estatística

Dedução vs Inferência

O Teorema do Limite Central

O Problema

A Distribuição Amostral

Demonstração Interativa (opcional)

Demonstração Visual em R

População Normal (N pequeno)

População Uniforme (N pequeno)

População Exponencial (N pequeno)

População Exponencial (N maior)

Perspectiva Crítica: Limites e Abusos do CLT

CLT e Modelos Lineares

Estatísticas de Teste: Sinal vs Ruído

2 A Equação Fundamental

O Que é o Erro?

3 O Modelo Mais Simples

Compacto vs Aumentado

O t de Student

O F de Snedcor

4 Regressão e Correlação

\(R^2\): Proporção de Redução do Erro

r de Pearson

5 Teste-t para Amostras Independentes

Esquemas de Codificação

contr.treatment (Dummy Coding)

Calcular EMMs com contr.treatment

contr.sum (Effect/Contrast Coding)

Calcular EMMs com contr.sum

Porquê Usar contr.sum?

Verificar com emmeans

Comparação de Modelos

6 Regressão Múltipla

Teste Omnibus

ANOVA Tipo III

O Que Significa “Controlar Para”?

7 One-Way ANOVA (k > 2)

Coeficientes: Contrastes, Não Efeitos

Como Funcionam os Graus de Liberdade?

Comparações Múltiplas

8 Múltiplos Factores: Interacções

O Que É uma Interacção?

Visualizar a Interacção

Decompor a Interacção

Padrões de Interacção (2 × 2)

Apenas Efeito Principal de Condição

Apenas Efeito Principal de Tarefa

Apenas Interacção (Crossover)

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

Efeito Principal de Tarefa + Interacção

Todos os Efeitos Significativos

9 ANCOVA: Factores + Covariáveis

Modelo Sem Covariável

ANCOVA: Factor + Covariável (Sem Interacção)

E Se Houver Interacção?

Mostrar as EMMs (média ± DP)

Mostrar os declives (o que emtrends() testa)

Notas Importantes sobre ANCOVA

10 Regressão Múltipla Quantitativa com Interacções

Ajustar e Comparar Modelos

Decompor a Interacção com emtrends()

Visualizar o que os declives simples representam

Visualizar os declives estimados por emtrends()

11 Medidas Repetidas e Modelos Mistos

12 Exercícios

1) Tradução: “que modelo é este?”

2) Um teste “parcial” por comparação de modelos

3) F vs t²: quando são equivalentes?

4) Interpretar contrastes vs EMMs

5) ANCOVA com e sem interacção

13 Recursos

`contr.treatment` (Dummy Coding)

Calcular EMMs com `contr.treatment`

`contr.sum` (Effect/Contrast Coding)

Calcular EMMs com `contr.sum`

Porquê Usar `contr.sum`?

Verificar com `emmeans`

Mostrar os declives (o que `emtrends()` testa)

Decompor a Interacção com `emtrends()`

Visualizar os declives estimados por `emtrends()`