Modelos Lineares no R
lm).
0 Introdução
Já vimos a teoria: Dados = Modelo + Erro.
Agora vamos à prática.
A ideia deste capítulo é construir uma escada de complexidade, começando pelo modelo mais simples possível e acabando em interacções e covariáveis. O “segredo” que muitos livros de estatística clássica não deixam claro é este: são tudo regressões lineares. Na maior parte dos casos, a função lm() (Linear Model) chega.
Este tutorial cobre apenas designs entre-sujeitos.
Se mediu a mesma pessoa várias vezes (e.g., Antes vs. Depois), pare aqui: um lm() simples viola o pressuposto da independência.
Para isso vai precisar de ANOVAs de medições repetidas ou de modelos mistos (linear mixed models), que aparecem no capítulo seguinte.
Preparação
Vamos carregar os pacotes e definir algumas opções.
contr.sum?
Em delineamentos factorais, a combinação SS Tipo III + codificação por omissão ou dummy (
contr.treatment) pode produzir resultados difíceis de interpretar (e por vezes indesejados, dependendo do desenho e do desequilíbrio amostral).A solução pragmática, muito comum em psicologia experimental quando se quer uma tabela ANOVA Tipo III, é usar sum-to-zero contrasts (
contr.sum) que são ortogonais.
Se quiser saber mais:
emmeans (a documentação que vale ouro)
Quando tiver dúvidas sobre comparações/contrastes, esta página costuma resolver muitos casos:
1 Nível 1: Teste-t para Uma Amostra
Aqui a lógica é: prever uma variável usando apenas um intercepto (a média).
\(Y_i = \beta_0 + e_i\)
Testar Contra 0
Vamos começar com um dataset que tem um zero que faz sentido: indiferença.
| ID | Love4Taylor |
|---|---|
| pp01 | 1.50 |
| pp02 | 2.15 |
| pp03 | 1.32 |
| pp04 | 1.17 |
| pp05 | -0.28 |
| pp06 | 1.72 |
A ideia: 0 = indiferença, > 0 = gosto, < 0 = não gosto.
Ou seja: aqui testar contra 0 não é uma parvoíce.
Call:
lm(formula = Love4Taylor ~ 1, data = taylor)
Residuals:
Min 1Q Median 3Q Max
-2.2552 -0.8202 0.1648 0.7548 2.4348
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.5652 0.1549 3.649 0.000638 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.095 on 49 degrees of freedom
Testar Contra Outro Valor
No entanto, há casos em que testar contra zero é bastante parvo…
library(palmerpenguins)
ds <- penguins
m <- lm(body_mass_g ~ 1, ds)
summary(m)
Call:
lm(formula = body_mass_g ~ 1, data = ds)
Residuals:
Min 1Q Median 3Q Max
-1501.8 -651.8 -151.8 548.2 2098.2
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4201.75 43.36 96.89 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 802 on 341 degrees of freedom
(2 observations deleted due to missingness)
O que acabámos de testar?
Que a massa dos pinguins na população difere de zero…ou seja, eles existem…pesam alguma coisa…
No entanto, é fácil usar a mesma lógica para testar hipóteses mais interessantes, por exemplo, será que a massa dos pinguins difere de 4kg (4000g).
No modelo linear, testar contra um valor c é o mesmo que centrar os dados em c e voltar a testar contra 0.
Call:
lm(formula = body_mass_g - 4000 ~ 1, data = ds)
Residuals:
Min 1Q Median 3Q Max
-1501.8 -651.8 -151.8 548.2 2098.2
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 201.75 43.36 4.652 4.7e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 802 on 341 degrees of freedom
(2 observations deleted due to missingness)
Call:
lm(formula = mass_minus_4000 ~ 1, data = ds)
Residuals:
Min 1Q Median 3Q Max
-1501.8 -651.8 -151.8 548.2 2098.2
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 201.75 43.36 4.652 4.7e-06 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 802 on 341 degrees of freedom
(2 observations deleted due to missingness)
2 Nível 2: amostras emparelhadas
Um teste-t para amostras emparelhado nada mais é que um teste-t a uma amostra aplicado à diferença entre momentos de medida.
\(\Delta = (\text{pós} - \text{pré}),\qquad \Delta_i = \beta_0 + e_i\)
Exemplo: reduzir risco de “luta”
Um mini-programa de intervenção para reduzir o risco de confrontos (em humanos…ou pinguins com mau feitio) foi avaliado ao medir-se o risco de confronto antes e depois do tratamento.
| pp | pre | post | diff |
|---|---|---|---|
| pp01 | 3.155373 | 1.659283 | 1.4960897 |
| pp02 | 5.260346 | 5.163075 | 0.0972709 |
| pp03 | 0.000000 | 0.000000 | 0.0000000 |
| pp04 | 0.000000 | 0.000000 | 0.0000000 |
| pp05 | 6.098303 | 5.454262 | 0.6440405 |
| pp06 | 3.482212 | 3.150342 | 0.3318697 |
Call:
lm(formula = diff ~ 1, data = fight)
Residuals:
Min 1Q Median 3Q Max
-1.4281 -0.4303 0.0094 0.3242 1.1333
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.48053 0.05703 8.427 1.28e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.5101 on 79 degrees of freedom
O gráfico abaixo ajuda a visualizar o emparelhamento.
Código
library(tidyr)
fight_long <- pivot_longer(fight, c("pre", "post"),
names_to = "momento", values_to = "risco")
ggplot(fight_long, aes(x = momento, y = risco, group = pp)) +
geom_line(alpha = 0.35) +
geom_point() +
theme_classic()3 Nível 3: duas amostras independentes (2 grupos)
Aqui já temos um preditor categórico com 2 níveis.
\(Y_i = \beta_0 + \beta_1 X + e_i\)
Vamos usar palmerpenguins para um exemplo rápido (é muito visual).
library(palmerpenguins)
ds <- penguins
# Garantir que sex tem apenas níveis válidos (o dataset original tem NAs)
ds <- ds[!is.na(ds$sex), ]
m <- lm(body_mass_g ~ sex, ds)
Anova(m, type = 3)| Sum Sq | Df | F value | Pr(>F) | |
|---|---|---|---|---|
| (Intercept) | 5884776622 | 1 | 11043.50022 | 0 |
| sex | 38878897 | 1 | 72.96099 | 0 |
| Residuals | 176380769 | 331 | NA | NA |
Médias marginais
Como há apenas 2 níveis, pedir pairwise costuma é redundante. Vamos pedir só as médias.
emmeans(m, ~ sex) sex emmean SE df lower.CL upper.CL
female 3862 56.8 331 3750 3974
male 4546 56.3 331 4435 4656
Confidence level used: 0.95
4 Nível 4: One-way ANOVA (k > 2)
Agora temos um factor com mais de 2 níveis (e portanto faz sentido pedir comparações múltiplas).
| Sum Sq | Df | F value | Pr(>F) | |
|---|---|---|---|---|
| (Intercept) | 5242015989 | 1 | 24687.8683 | 0 |
| species | 145190219 | 2 | 341.8949 | 0 |
| Residuals | 70069447 | 330 | NA | NA |
Médias + post-hoc numa chamada (eficiente)
$`emmeans of species`
species emmean SE df lower.CL upper.CL
Adelie 3706 38.1 330 3631 3781
Chinstrap 3733 55.9 330 3623 3843
Gentoo 5092 42.2 330 5009 5176
Confidence level used: 0.95
$`pairwise differences of species`
1 estimate SE df t.ratio p.value
Adelie - Chinstrap -26.9 67.7 330 -0.398 0.9164
Adelie - Gentoo -1386.3 56.9 330 -24.359 <0.0001
Chinstrap - Gentoo -1359.3 70.0 330 -19.406 <0.0001
P value adjustment: tukey method for comparing a family of 3 estimates
Ajustes de Valores-p
Em emmeans, as comparações (e os seus ajustes) aparecem de forma mais transparente quando usamos summary().
# Comparações com FDR (BH)
emmeans(m, pairwise ~ species, adjust = "fdr")$emmeans
species emmean SE df lower.CL upper.CL
Adelie 3706 38.1 330 3631 3781
Chinstrap 3733 55.9 330 3623 3843
Gentoo 5092 42.2 330 5009 5176
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Adelie - Chinstrap -26.9 67.7 330 -0.398 0.6909
Adelie - Gentoo -1386.3 56.9 330 -24.359 <0.0001
Chinstrap - Gentoo -1359.3 70.0 330 -19.406 <0.0001
P value adjustment: fdr method for 3 tests
Não existe um ajuste “mágico”.
Tukey: muito comum em ANOVA clássica (família de comparações entre todos os pares).
Holm / BH: alternativas frequentes quando há muitas comparações.
O importante é ser explícito e coerente com o objectivo da análise.
5 Nível 5: ANOVA factorial (interacções)
A pergunta: o efeito de species depende de sex?
| Sum Sq | Df | F value | Pr(>F) | |
|---|---|---|---|---|
| (Intercept) | 5232595969 | 1 | 54661.827958 | 0.0000000 |
| species | 143001222 | 2 | 746.924492 | 0.0000000 |
| sex | 29851220 | 1 | 311.838003 | 0.0000000 |
| species:sex | 1676557 | 2 | 8.756997 | 0.0001973 |
| Residuals | 31302628 | 327 | NA | NA |
Na presença de uma interacção, não interprete efeitos principais como se estivessem sozinhos.
Ignorar a interacção é não ouvir o que os dados nos querem dizer e contar apenas uma parte da história.
O que normalmente queremos são efeitos simples / comparações condicionais.
Efeitos Simples/Comparações Condicionais
Queremos:
Diferenças entre espécies dentro de cada sexo
Diferenças entre sexos dentro de cada espécie
$`emmeans of species | sex`
sex = female:
species emmean SE df lower.CL upper.CL
Adelie 3369 36.2 327 3298 3440
Chinstrap 3527 53.1 327 3423 3632
Gentoo 4680 40.6 327 4600 4760
sex = male:
species emmean SE df lower.CL upper.CL
Adelie 4043 36.2 327 3972 4115
Chinstrap 3939 53.1 327 3835 4043
Gentoo 5485 39.6 327 5407 5563
Confidence level used: 0.95
$`pairwise differences of species | sex`
sex = female:
2 estimate SE df t.ratio p.value
Adelie - Chinstrap -158 64.2 327 -2.465 0.0377
Adelie - Gentoo -1311 54.4 327 -24.088 <0.0001
Chinstrap - Gentoo -1153 66.8 327 -17.246 <0.0001
sex = male:
2 estimate SE df t.ratio p.value
Adelie - Chinstrap 105 64.2 327 1.627 0.2357
Adelie - Gentoo -1441 53.7 327 -26.855 <0.0001
Chinstrap - Gentoo -1546 66.2 327 -23.345 <0.0001
P value adjustment: tukey method for comparing a family of 3 estimates
$`emmeans of sex | species`
species = Adelie:
sex emmean SE df lower.CL upper.CL
female 3369 36.2 327 3298 3440
male 4043 36.2 327 3972 4115
species = Chinstrap:
sex emmean SE df lower.CL upper.CL
female 3527 53.1 327 3423 3632
male 3939 53.1 327 3835 4043
species = Gentoo:
sex emmean SE df lower.CL upper.CL
female 4680 40.6 327 4600 4760
male 5485 39.6 327 5407 5563
Confidence level used: 0.95
$`pairwise differences of sex | species`
species = Adelie:
2 estimate SE df t.ratio p.value
female - male -675 51.2 327 -13.174 <0.0001
species = Chinstrap:
2 estimate SE df t.ratio p.value
female - male -412 75.0 327 -5.487 <0.0001
species = Gentoo:
2 estimate SE df t.ratio p.value
female - male -805 56.7 327 -14.188 <0.0001
Visualizar a interacção (obrigatório)
Abaixo estão duas visões complementares:
- Dados brutos (o mais “honesto” possível)
# Dados brutos (exemplo)
ggplot(ds, aes(x = species, y = body_mass_g, fill = sex, color = sex)) +
geom_violin(alpha = .5, color = "black") +
geom_jitter(position = position_jitterdodge(dodge.width = .9, jitter.width = .15),
show.legend = FALSE) +
theme_classic()- Médias marginais (o que o modelo está realmente a prever)
# Médias marginais estimadas
for_g <- as.data.frame(emmeans(m, ~ species * sex))
ggplot(for_g, aes(x = species, y = emmean, ymin = lower.CL, ymax = upper.CL,
color = sex, group = sex)) +
geom_line() + geom_point(size = 2.5) + geom_errorbar(width = 0.15) +
theme_classic()6 Nível 6: ANCOVA (factor + covariável)
Agora juntamos uma covariável contínua (ex.: flipper_length_mm).
| Sum Sq | Df | F value | Pr(>F) | |
|---|---|---|---|---|
| (Intercept) | 8.052763e+02 | 1 | 0.0097542 | 0.9213869 |
| species | 1.262607e+07 | 2 | 76.4688266 | 0.0000000 |
| sex | 1.274795e+07 | 1 | 154.4138757 | 0.0000000 |
| flipper_length_mm | 4.389049e+06 | 1 | 53.1638663 | 0.0000000 |
| species:sex | 1.739989e+06 | 2 | 10.5381054 | 0.0000367 |
| Residuals | 2.691358e+07 | 326 | NA | NA |
Em ANOVA factorial experimental, é comum testar interacções por defeito (e usar SS Tipo III para avaliar efeitos “controlando” o resto).
Em regressão múltipla, muita gente coloca apenas efeitos principais (por vezes por hábito, por vezes para reduzir complexidade).
Nenhuma destas tradições é automaticamente “errada”.
Evitar interacções pode ser prudente (menos parâmetros; menos espaço para pescarias).
Mas ignorar interacções quando elas existem pode esconder o fenómeno.
Interacções com covariáveis (quando considerar)
Um caso clássico é permitir que a relação entre flipper_length_mm e body_mass_g varie por espécie/sexo.
| Sum Sq | Df | F value | Pr(>F) | |
|---|---|---|---|---|
| (Intercept) | 8.052763e+02 | 1 | 0.0097542 | 0.9213869 |
| species | 1.262607e+07 | 2 | 76.4688266 | 0.0000000 |
| sex | 1.274795e+07 | 1 | 154.4138757 | 0.0000000 |
| flipper_length_mm | 4.389049e+06 | 1 | 53.1638663 | 0.0000000 |
| species:sex | 1.739989e+06 | 2 | 10.5381054 | 0.0000367 |
| Residuals | 2.691358e+07 | 326 | NA | NA |
7 Cábula do emmeans
| Situação | Pedido típico | Nota |
|---|---|---|
| 1 amostra | summary(lm(y ~ 1)) |
Intercepto = média |
| 2 grupos | emmeans(m, ~ A) |
Post-hoc é redundante |
| k > 2 | emmeans(m, list(~ A, pairwise ~ A)) |
Post-hoc faz sentido |
| Interacção |
emmeans(m, ~ A | B) + pairs()/contrast()
|
Efeitos simples/condicionais |
8 emmeans: Alguns Truques
Nesta secção não vou entrar em teoria (isso fica para o capítulo teórico). A ideia é apenas dar-lhe botões que resolvem a maioria das dúvidas práticas quando está a decompor efeitos.
Gráfico Rápido de Interacções
Quando tem uma interacção, é muito fácil interpretar mal uma tabela. O emmip() é uma forma rápida de ver o que o modelo está a prever.
O emmip() é óptimo para ver padrões, mas para um gráfico de artigo normalmente vai querer um ggplot2 mais personalizado.
Contrastes Planeados
Nem sempre a pergunta científica é “quem difere de quem?”. Às vezes a pergunta é um contraste específico, por exemplo:
- Gentoo vs a média de (Adelie e Chinstrap).
levels(ds$species)[1] "Adelie" "Chinstrap" "Gentoo"
$emmeans
species emmean SE df lower.CL upper.CL
Adelie 3706 25.6 327 3656 3757
Chinstrap 3733 37.5 327 3659 3807
Gentoo 5082 28.4 327 5026 5138
Results are averaged over the levels of: sex
Confidence level used: 0.95
$contrasts
contrast estimate SE df t.ratio p.value
Gentoo_vs_Outros 1363 36.3 327 37.495 <0.0001
Results are averaged over the levels of: sex
Confirme sempre a ordem dos níveis (a ordem do vector tem de bater certo)
Ajustes de Valores-p
contrast estimate SE df t.ratio p.value
female - male -631 35.7 327 -17.659 <0.0001
Results are averaged over the levels of: species
pairs(emm_sp, adjust = "fdr") contrast estimate SE df t.ratio p.value
female - male -631 35.7 327 -17.659 <0.0001
Results are averaged over the levels of: species
pairs(emm_sp, adjust = "bonferroni") contrast estimate SE df t.ratio p.value
female - male -631 35.7 327 -17.659 <0.0001
Results are averaged over the levels of: species
O método “certo” depende do objectivo.
Se está a fazer todas as comparações possíveis porque sim, ajustes mais conservadores são uma defesa contra falsos positivos.
Se tem contrastes planeados (pré-especificados), faz sentido reportar esses contrastes com muito mais foco.
9 Exercícios
- No dataset
Love4Taylor.csv, teste se a média é diferente de 0 e de um outro valor (idealmente um que tenha algum significado).
Clique para ver as soluções
Call:
lm(formula = Love4Taylor ~ 1, data = taylor)
Residuals:
Min 1Q Median 3Q Max
-2.2552 -0.8202 0.1648 0.7548 2.4348
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.5652 0.1549 3.649 0.000638 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.095 on 49 degrees of freedom
Call:
lm(formula = love_vs_neg1 ~ 1, data = taylor)
Residuals:
Min 1Q Median 3Q Max
-2.2552 -0.8202 0.1648 0.7548 2.4348
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 1.5652 0.1549 10.11 1.42e-13 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 1.095 on 49 degrees of freedom
- No dataset
fight_risk.csv, faça:
o modelo
lm(diff ~ 1)as descritivas do pré e pós
uma frase que reporte o resultado (uma ou duas linhas).
Clique para ver as soluções
Call:
lm(formula = diff ~ 1, data = fight)
Residuals:
Min 1Q Median 3Q Max
-1.4281 -0.4303 0.0094 0.3242 1.1333
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 0.48053 0.05703 8.427 1.28e-12 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.5101 on 79 degrees of freedom
mean(fight$pre)[1] 3.043889
sd(fight$pre)[1] 1.560794
mean(fight$post)[1] 2.563358
sd(fight$post)[1] 1.574395
Exemplo de como reportar os resultados:
Encontrámos evidência significativa de que o tratamento reduziu o risco de comportamentos violentos, t(29) = 7.45, p < .001; Mpre = 50.51, SD = 7.80, Mpós = 45.86, SD = 7.65.
- Com
palmerpenguins, ajustelm(body_mass_g ~ species * sex)e faça:
comparações entre espécies dentro de cada sexo
comparação entre sexos dentro de cada espécie (uma por espécie)
dois gráficos: dados brutos e emmeans.
Clique para ver as soluções
| Sum Sq | Df | F value | Pr(>F) | |
|---|---|---|---|---|
| (Intercept) | 5232595969 | 1 | 54661.827958 | 0.0000000 |
| species | 143001222 | 2 | 746.924492 | 0.0000000 |
| sex | 29851220 | 1 | 311.838003 | 0.0000000 |
| species:sex | 1676557 | 2 | 8.756997 | 0.0001973 |
| Residuals | 31302628 | 327 | NA | NA |
# (a) species dentro de sex
sp_within_sex <- emmeans(m, list(pairwise ~ species | sex, pairwise ~ sex | species))
# gráficos (exemplos/sugestões)
ggplot(ds, aes(x = species, y = body_mass_g, fill = sex, color = sex)) +
geom_violin(alpha = .5, color = "black") +
geom_jitter(position = position_jitterdodge(dodge.width = .9, jitter.width = .15),
show.legend = FALSE) +
theme_classic()for_g <- as.data.frame(emmeans(m, ~ species * sex))
ggplot(for_g, aes(x = species, y = emmean, ymin = lower.CL, ymax = upper.CL,
color = sex, group = sex)) +
geom_line() + geom_point(size = 2.5) + geom_errorbar(width = 0.15) +
theme_classic()- Faça uma ANCOVA com
flipper_length_mme compare um modelo sem e com interacções com a covariável.
Clique para ver as soluções
| Sum Sq | Df | F value | Pr(>F) | |
|---|---|---|---|---|
| (Intercept) | 47510.60 | 1 | 0.5858908 | 0.4445749 |
| species | 216230.99 | 2 | 1.3332578 | 0.2650728 |
| sex | 16183.84 | 1 | 0.1995757 | 0.6553648 |
| flipper_length_mm | 4492633.05 | 1 | 55.4022165 | 0.0000000 |
| species:sex | 634560.21 | 2 | 3.9126323 | 0.0209482 |
| species:flipper_length_mm | 279932.75 | 2 | 1.7260362 | 0.1796363 |
| sex:flipper_length_mm | 73438.41 | 1 | 0.9056272 | 0.3419933 |
| species:sex:flipper_length_mm | 596743.03 | 2 | 3.6794555 | 0.0263072 |
| Residuals | 26030280.01 | 321 | NA | NA |
anova(m_main, m_int)| Res.Df | RSS | Df | Sum of Sq | F | Pr(>F) |
|---|---|---|---|---|---|
| 326 | 26913579 | NA | NA | NA | NA |
| 321 | 26030280 | 5 | 883299.4 | 2.178533 | 0.0563138 |