Modelos Lineares no R

Modelos Lineares

Estatística Inferencial

emmeans

Tutorial

Um fluxo prático para transformar perguntas de investigação em modelos, estimativas e comparações.

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

Este tutorial é para analisar dados. A teoria que une teste-t, regressão, ANOVA, ANCOVA e interacções está no capítulo de comparação de modelos. Para a mecânica de erro-padrão, t e valor-p, consulte Fundamentos da Inferência; para unidades, replicação e alcance causal, consulte Delineamento da Investigação. Aqui partimos dessas ideias e concentramo-nos num fluxo de trabalho mais simples:

pergunta → modelo → estimativa → comparação → gráfico → interpretação

Não precisamos de começar pelo nome do teste. Começamos por perguntar o que queremos saber e escrevemos um modelo que represente essa pergunta.

As secções seguintes reutilizam os mesmos dados dos pinguins, acrescentando uma só decisão de cada vez: primeiro comparamos grupos, depois acrescentamos um preditor quantitativo e, por fim, perguntamos se uma diferença depende de outra variável. Em cada caso percorremos o mesmo caminho: formulamos a pergunta, ajustamos o modelo, obtemos a estimativa, comparamos o que precisa de ser comparado, fazemos um gráfico e interpretamos o resultado à luz do delineamento.

Por exemplo, em Dois Grupos perguntamos pela massa nos dois sexos, ajustamos body_mass_g ~ sex, pedimos as médias e a diferença estimada com emmeans(), mostramos os pontos e a média, e só então interpretamos a comparação. Em Preditor Quantitativo, a mesma sequência termina num declive e numa recta ajustada. O modelo muda. O modo de raciocinar mantém-se.

Antes do Modelo

Confirme sempre:

qual é a unidade de observação;
se cada linha pode ser tratada como independente;
quais variáveis são quantitativas e quais são factores;
se há observações repetidas da mesma pessoa, animal, local ou item.

Se a mesma unidade aparece várias vezes, consulte medidas repetidas.

Ter o mesmo nível de um factor em várias linhas é outra coisa. Vários espécimes podem ser Adelie e continuar a ser unidades diferentes. Nesse caso, species é um preditor entre unidades. Procure medidas repetidas ou pseudorreplicação quando as linhas repetem a mesma unidade, ou quando partilham um agrupamento relevante que o modelo está a ignorar.

Confirme também se o identificador é realmente um identificador. Duas linhas com IDs diferentes não garantem duas unidades diferentes se o mesmo organismo, participante ou utilizador puder ter sido registado novamente sob outro código.

1 Preparação

Este tutorial usa palmerpenguins, ggplot2, emmeans e car. Se ainda não os instalou, veja a secção sobre pacotes. Num script de análise, carregamos as dependências no início:

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

variaveis <- c("body_mass_g", "flipper_length_mm", "species", "sex")
dados <- penguins[complete.cases(penguins[, variaveis]), variaveis]
dados$species <- droplevels(dados$species)
dados$sex <- factor(droplevels(dados$sex), levels = c("female", "male"))

# Keep factorial models centred: sum coding for species and -1/+1 for sex.
contrasts(dados$species) <- contr.sum(nlevels(dados$species))
contrasts(dados$sex) <- matrix(
  c(-1, 1),
  ncol = 1,
  dimnames = list(levels(dados$sex), "male_vs_female")
)

A sintaxe principal de lm() é:

lm(resposta ~ preditores, data = dados)

+ acrescenta termos. * inclui efeitos principais e a interacção. Para uma explicação completa das fórmulas e da matemática por trás destas comparações, volte ao capítulo teórico.

2 Uma Média

Pergunta

No ficheiro Love4Taylor.csv, a média de Love4Taylor difere de zero? Zero é uma referência com significado nesta escala.

taylor <- read.csv("../data/Love4Taylor.csv")
taylor <- taylor[complete.cases(taylor$Love4Taylor), , drop = FALSE]

Modelo

m_zero <- lm(Love4Taylor ~ 0, data = taylor)
m_media <- lm(Love4Taylor ~ 1, data = taylor)

Estimativa e Comparação

A estimativa que procuramos é a média. Comparamos o modelo dessa média com o modelo que fixa a resposta em zero e mostramos a estimativa com o seu intervalo:

anova(m_zero, m_media)

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

coef(m_media)

(Intercept) 
     0.5652

confint(m_media)

                2.5 %    97.5 %
(Intercept) 0.2539361 0.8764639

Não nos interessa apenas saber se o valor-p atravessou um limiar. A estimativa é a média. O intervalo mostra a incerteza associada a essa estimativa sob o modelo.

Outra Referência

Para testar a média contra uma referência mu0, recentre a resposta:

mu0 <- 3
m_ref <- lm(I(Love4Taylor - mu0) ~ 1, data = taylor)
summary(m_ref)
confint(m_ref)

O intercepto passa a ser a diferença média relativamente a mu0.

3 Dois Grupos

Pergunta

Nos dados dos pinguins, qual é a diferença estimada de massa corporal entre os dois níveis de sex?

Modelo. Ajustamos uma média para cada nível:

m_sexo <- lm(body_mass_g ~ sex, data = dados)

Estimativa e comparação. Em vez de interpretar a codificação do coeficiente à mão, pedimos ao modelo as médias previstas e a diferença entre elas:

emmeans(m_sexo, list(pairwise ~ sex))

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

$`pairwise differences of sex`
 1             estimate SE  df t.ratio p.value
 female - male     -683 80 331  -8.542 <0.0001

Gráfico. Verificamos se a diferença resumida pelo modelo faz sentido face aos dados:

ggplot(dados, aes(x = sex, y = body_mass_g)) +
  geom_jitter(width = 0.12, alpha = 0.35) +
  stat_summary(fun = mean, geom = "point", size = 3) +
  theme_classic() +
  labs(x = "Sexo", y = "Massa corporal (g)")

Interpretação

A diferença estimada compara as massas previstas para os dois sexos no modelo sem outros preditores. O gráfico permite ver se essa diferença representa a separação geral dos dados ou se esconde grande sobreposição. Como sex não foi atribuído experimentalmente, a comparação descreve estes pinguins. Não estima o efeito causal de alterar o sexo.

4 Preditor Quantitativo

Pergunta

Como varia a massa corporal prevista com o comprimento da barbatana?

Modelo. Ajustamos uma recta aos dados:

m_barbatana <- lm(body_mass_g ~ flipper_length_mm, data = dados)

Estimativa e comparação. O declive responde directamente à pergunta: mudança prevista na massa por unidade adicional de flipper_length_mm. O intervalo quantifica a incerteza dessa estimativa.

summary(m_barbatana)


Call:
lm(formula = body_mass_g ~ flipper_length_mm, data = dados)

Residuals:
     Min       1Q   Median       3Q      Max 
-1057.33  -259.79   -12.24   242.97  1293.89 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)       -5872.09     310.29  -18.93   <2e-16 ***
flipper_length_mm    50.15       1.54   32.56   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 393.3 on 331 degrees of freedom
Multiple R-squared:  0.7621,    Adjusted R-squared:  0.7614 
F-statistic:  1060 on 1 and 331 DF,  p-value: < 2.2e-16

confint(m_barbatana)

                        2.5 %      97.5 %
(Intercept)       -6482.47224 -5261.71313
flipper_length_mm    47.12339    53.18314

`emmeans()` Sem Grupos: `~ 1`

Também podemos pedir uma única estimativa prevista, sem a dividir por grupos:

emmeans(m_barbatana, ~ 1)

 1       emmean   SE  df lower.CL upper.CL
 overall   4207 21.6 331     4165     4249

Confidence level used: 0.95

~ 1 Não Significa «Intercepto»

Em emmeans, ~ 1 significa apenas uma estimativa sem grupos. Se o modelo tiver preditores quantitativos, estes são por defeito avaliados nos valores da grelha de referência — normalmente as suas médias — e não necessariamente em zero. Portanto, não chame automaticamente a este resultado «o intercepto».

Num modelo com uma covariável quantitativa, isto não pede automaticamente o intercepto. Por defeito, emmeans constrói a grelha de referência colocando a covariável no seu valor médio. Assim, esta chamada devolve a massa prevista para um pinguim com comprimento da barbatana igual à média observada.

Se quisermos pedir exactamente a previsão associada ao intercepto do modelo, temos de avaliar a recta em zero:

coef(m_barbatana)["(Intercept)"]

(Intercept) 
  -5872.093

emmeans(
  m_barbatana,
  ~ 1,
  at = list(flipper_length_mm = 0)
)

 1       emmean  SE  df lower.CL upper.CL
 overall  -5872 310 331    -6482    -5262

Confidence level used: 0.95

As duas estimativas coincidem: num modelo y = beta0 + beta1 * x, quando x = 0 a previsão é beta0. Neste exemplo, contudo, uma barbatana com 0 mm está muito fora dos dados e o intercepto original não tem uma interpretação biológica útil.

Podemos recentrar o preditor para dar ao zero um significado melhor:

dados$flipper_200 <- dados$flipper_length_mm - 200

m_barbatana_200 <- lm(
  body_mass_g ~ flipper_200,
  data = dados
)

coef(m_barbatana_200)["(Intercept)"]

(Intercept) 
   4158.561

emmeans(
  m_barbatana_200,
  ~ 1,
  at = list(flipper_200 = 0)
)

 1       emmean   SE  df lower.CL upper.CL
 overall   4159 21.6 331     4116     4201

Confidence level used: 0.95

Agora zero significa 200 mm de barbatana. O intercepto e a estimativa do emmeans voltam a ser a mesma quantidade, mas desta vez descrevem uma previsão num ponto plausível e fácil de interpretar. Em modelos Gaussianos com ligação identidade, não precisamos de type = "response": a escala do modelo já é a escala da resposta.

Gráfico. Mostramos os dados, a recta estimada e a respectiva faixa de incerteza:

ggplot(dados, aes(x = flipper_length_mm, y = body_mass_g)) +
  geom_point(alpha = 0.45) +
  geom_smooth(method = "lm", se = TRUE) +
  theme_classic() +
  labs(x = "Comprimento da barbatana (mm)", y = "Massa corporal (g)")

Interpretação

O declive é a mudança prevista na massa corporal por cada milímetro adicional de barbatana, dentro do intervalo observado. A recta ajuda a verificar se uma relação linear é uma descrição razoável e se há zonas com poucos dados. Não interprete automaticamente a associação como causal: estes dados não são uma experiência em que se tenha manipulado o comprimento da barbatana.

5 Vários Preditores

Pergunta

A associação entre comprimento da barbatana e massa corporal mantém-se quando também incluímos species?

m_aditivo <- lm(body_mass_g ~ species + flipper_length_mm, data = dados)

m_sem_especie <- lm(body_mass_g ~ flipper_length_mm, data = dados)
anova(m_sem_especie, m_aditivo)

Res.Df	RSS	Df	Sum of Sq	F	Pr(>F)
331	51211963	NA	NA	NA	NA
329	45843144	2	5368818	19.26505	0

Para uma tabela dos termos do modelo, use uma função que corresponda à hipótese que pretende testar. Se estiver a seguir uma análise Tipo III com factores, veja a explicação de contrastes e somas de quadrados antes de escolher a parametrização.

Para comunicar os resultados, costuma ser mais útil pedir previsões ou médias ajustadas que correspondam à pergunta científica:

emmeans(m_aditivo, list(pairwise ~ species))

$`emmeans of species`
 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 

$`pairwise differences of species`
 1                  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

emtrends(m_aditivo, ~ species, var = "flipper_length_mm")

 species   flipper_length_mm.trend   SE  df lower.CL upper.CL
 Adelie                       40.6 3.08 329     34.5     46.7
 Chinstrap                    40.6 3.08 329     34.5     46.7
 Gentoo                       40.6 3.08 329     34.5     46.7

Confidence level used: 0.95

ggplot(dados, aes(flipper_length_mm, body_mass_g, colour = species)) +
  geom_point(alpha = 0.4) +
  geom_smooth(method = "lm", se = TRUE) +
  theme_classic() +
  labs(x = "Comprimento da barbatana (mm)", y = "Massa corporal (g)",
       colour = "Espécie")

Interpretação

Aqui a pergunta mudou: o declive de flipper_length_mm é estimado depois de incluir species, e as médias de espécie são comparadas no mesmo valor de referência do preditor quantitativo. O gráfico mostra se rectas aditivas são uma descrição plausível. Se as rectas diferirem na inclinação, a pergunta já não é apenas «qual é a diferença entre espécies?». Precisamos da interacção apresentada mais à frente.

Ajustamento e Confundimento

Incluir uma variável num modelo não resolve, por si só, o confundimento. Num modelo múltiplo, a interpretação de cada coeficiente depende dos restantes termos que também estão no modelo. Isso é uma propriedade estatística, não uma afirmação causal sobre o que aconteceria se manipulássemos uma das variáveis.

6 Vários Grupos

Pergunta

A massa corporal prevista varia entre espécies?

m_especie <- lm(body_mass_g ~ species, data = dados)

O teste omnibus pode ser obtido comparando com o modelo de uma só média:

m_especie_0 <- lm(body_mass_g ~ 1, data = dados)
anova(m_especie_0, m_especie)

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

Depois, escolha as comparações que respondem à pergunta. Se quer comparar todos os pares:

emmeans(m_especie, list(pairwise ~ species))

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

Se tinha planeado uma comparação mais específica, escreva-a directamente em vez de gerar todos os pares só porque estão disponíveis.

ggplot(dados, aes(x = species, y = body_mass_g, colour = sex)) +
  geom_jitter(
    position = position_jitterdodge(dodge.width = 0.5, jitter.width = 0.12),
    alpha = 0.3
  ) +
  stat_summary(
    fun = mean,
    geom = "point",
    position = position_dodge(width = 0.5),
    size = 3
  ) +
  theme_classic() +
  labs(x = "Espécie", y = "Massa corporal (g)", colour = "Sexo")

Interpretação

O omnibus responde apenas se o modelo com médias específicas por espécie melhora o modelo de uma só média. As comparações e o gráfico dizem quais espécies diferem e quanta sobreposição existe. Se a pergunta tiver sido planeada para um par de espécies, esse contraste é preferível a escolher uma comparação depois de ver todos os resultados.

7 Interacção Entre Factores

Pergunta

A diferença entre níveis de sex é igual em todas as espécies?

Primeiro ajustamos o modelo que obriga as diferenças a serem aditivas e depois o modelo que permite uma interacção.

m_sem_interaccao <- lm(body_mass_g ~ species + sex, data = dados)
m_interaccao <- lm(body_mass_g ~ species * sex, data = dados)

# Direct test of the interaction.
anova(m_sem_interaccao, m_interaccao)

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

# Default omnibus table for the factorial model used in this material.
Anova(m_interaccao, type = "III")

	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 comparação dos dois modelos isola a interacção. A tabela Tipo III mantém species, sex e species:sex no mesmo modelo, usando os contrastes centrados definidos no início.

É essa a razão da nossa preferência aqui: quando perguntamos pelo efeito médio de species ou de sex, não obrigamos primeiro a interacção a ser zero. Os restantes termos continuam no modelo e funcionam como nuisance parameters para a pergunta que estamos a testar.

Se a interacção for relevante, não pare na linha da tabela. Veja as células e as comparações condicionais:

emmeans(
  m_interaccao,
  list(
    pairwise ~ species,
    ~ sex,
    pairwise ~ species | sex,
    pairwise ~ sex | species
  )
)

$`emmeans of species`
 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 

$`pairwise differences of species`
 1                  estimate   SE  df t.ratio p.value
 Adelie - Chinstrap    -26.9 45.4 327  -0.593  0.8241
 Adelie - Gentoo     -1376.1 38.2 327 -36.007 <0.0001
 Chinstrap - Gentoo  -1349.2 47.0 327 -28.682 <0.0001

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

$`emmeans of sex`
 sex    emmean   SE  df lower.CL upper.CL
 female   3859 25.3 327     3809     3908
 male     4489 25.2 327     4440     4539

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

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

ggplot(dados, aes(x = species, y = body_mass_g, colour = sex)) +
  geom_jitter(
    position = position_jitterdodge(dodge.width = 0.5, jitter.width = 0.12),
    alpha = 0.3
  ) +
  stat_summary(
    fun = mean,
    geom = "point",
    position = position_dodge(width = 0.5),
    size = 3
  ) +
  theme_classic() +
  labs(x = "Espécie", y = "Massa corporal (g)", colour = "Sexo")

Interpretação

A diferença entre sexos não é necessariamente a mesma em todas as espécies. A comparação entre os modelos pergunta se permitir essa diferença entre espécies melhora a explicação dos dados. As médias por célula, os contrastes condicionais e o gráfico mostram em que espécies a diferença aparece e evitam reduzir uma interacção a dois valores-p separados.

8 ANCOVA

Pergunta

As espécies diferem na massa corporal quando comparamos previsões para o mesmo comprimento de barbatana?

O modelo tradicional de ANCOVA começa por impor o mesmo declive a todas as espécies:

m_paralelo <- lm(
  body_mass_g ~ species + flipper_length_mm,
  data = dados
)

Antes de interpretar essa ANCOVA, precisamos de perguntar se a hipótese de declives paralelos é razoável. Permitimos então que o declive varie por espécie:

m_declives <- lm(
  body_mass_g ~ species * flipper_length_mm,
  data = dados
)

anova(m_paralelo, m_declives)

Res.Df	RSS	Df	Sum of Sq	F	Pr(>F)
329	45843144	NA	NA	NA	NA
327	44391669	2	1451475	5.345963	0.0051931

Se o modelo com interacção for necessário para responder à pergunta, examine os declives directamente:

trends <- emtrends(m_declives, ~ species, var = "flipper_length_mm")
summary(trends, infer = c(TRUE, TRUE))

species	flipper_length_mm.trend	SE	df	lower.CL	upper.CL	t.ratio	p.value
Adelie	32.68891	4.691630	327	23.45933	41.91850	6.967497	0e+00
Chinstrap	34.57339	6.311529	327	22.15707	46.98972	5.477816	1e-07
Gentoo	54.16542	5.150527	327	44.03307	64.29777	10.516481	0e+00

summary(pairs(trends), infer = c(TRUE, TRUE))

contrast	estimate	SE	df	lower.CL	upper.CL	t.ratio	p.value
Adelie - Chinstrap	-1.88448	7.864273	327	-20.40041	16.6314458	-0.2396254	0.9688451
Adelie - Gentoo	-21.47651	6.967016	327	-37.87990	-5.0731138	-3.0825974	0.0062871
Chinstrap - Gentoo	-19.59203	8.146369	327	-38.77213	-0.4119247	-2.4050013	0.0440144

E veja o padrão no gráfico:

ggplot(dados, aes(
  x = flipper_length_mm,
  y = body_mass_g,
  colour = species
)) +
  geom_point(alpha = 0.4) +
  geom_smooth(method = "lm", se = TRUE) +
  theme_classic() +
  labs(
    x = "Comprimento da barbatana (mm)",
    y = "Massa corporal (g)",
    colour = "Espécie"
  )

Interpretação

A ANCOVA responde à comparação entre espécies para o mesmo comprimento de barbatana apenas enquanto os declives forem suficientemente compatíveis. A comparação entre m_paralelo e m_declives, as estimativas dos declives e o gráfico verificam essa condição. Se os declives variarem, não há uma única diferença entre espécies que resuma adequadamente todos os comprimentos observados. Devemos interpretar os declives ou previsões em valores substantivamente úteis, sem extrapolar.

9 Interacção Quantitativa

Para uma interacção entre duas variáveis quantitativas, o fluxo é o mesmo. Ajuste primeiro o modelo aditivo, depois permita a interacção:

# Exemplo genérico: substitua y, x e z por variáveis do seu conjunto de dados.
m_add <- lm(y ~ x + z, data = dados)
m_int <- lm(y ~ x * z, data = dados)

anova(m_add, m_int)

Se a interacção for necessária, não há um único «efeito de x» que possamos apresentar sem qualificação. Estime o declive de x em valores de z que sejam úteis para a pergunta e estejam dentro do intervalo coberto pelos dados.

A página de comparação de modelos explica por que este declive muda com z. emmeans::emtrends() permite pedir esses declives sem fazer a álgebra à mão.

Estimativa, Gráfico e Interpretação

Depois de comparar os modelos, peça os declives de x em valores de z que façam sentido para a investigação e represente as previsões. Assim, a conclusão pode dizer em que valores de z a associação muda, em vez de apresentar um efeito médio que não responde à pergunta. Se z não variar suficientemente na amostra, essa leitura será incerta. O gráfico deve tornar essa limitação visível.

10 Diagnóstico

Um modelo executado sem erro não é automaticamente um bom modelo. Depois de especificar a pergunta e ajustar o modelo, examine:

linearidade quando há preditores quantitativos;
padrão e variância dos resíduos;
observações muito influentes;
independência de acordo com o delineamento;
se as previsões estão a extrapolar para regiões sem dados.

O tutorial de pressupostos e diagnóstico desenvolve este passo.

11 Da Pergunta ao Modelo

Pergunta	Modelo típico	O que pedir depois
A média difere de uma referência?	`y ~ 1` depois de recentrar	intercepto + IC
Dois grupos diferem?	`y ~ grupo`	`emmeans()` + `pairs()`
Y varia com X?	`y ~ x`	declive + IC + gráfico
X acrescenta informação depois de incluir Z?	`y ~ x + z`	coeficiente/contraste condicional
Vários grupos diferem?	`y ~ grupo`	omnibus + contrastes planeados
Uma diferença depende de outro factor?	`y ~ A * B`	células + efeitos simples
Um declive depende do grupo?	`y ~ grupo * x`	`emtrends()` + gráfico
Um declive depende de outra variável quantitativa?	`y ~ x * z`	declives simples

12 Exercícios

Em Love4Taylor.csv, escolha uma referência com significado diferente de zero. Recentre a resposta e interprete o intercepto e o intervalo de confiança.
Em penguins, ajuste body_mass_g ~ species. Antes de consultar as médias, escolha uma diferença entre espécies que lhe interessaria. Use emmeans() e peça directamente esse contraste; explique por que o escolheu antes de ver os resultados.
Ajuste body_mass_g ~ flipper_length_mm. Interprete o declive em unidades reais e faça um gráfico com os dados observados e a recta ajustada.
Compare body_mass_g ~ species + sex com body_mass_g ~ species * sex. Depois peça, numa só chamada, emmeans(m, list(pairwise ~ species, ~ sex, pairwise ~ species | sex, pairwise ~ sex | species)). Distinga efeitos médios de comparações condicionais. Se olhar apenas para uma diferença média entre sexos, que informação pode perder?
Compare a ANCOVA com declive comum e o modelo species * flipper_length_mm. Use emtrends() para explicar o que muda entre os modelos.

Solução 1: Uma Referência Com Significado

Por exemplo, se 3 for uma referência substantivamente útil na escala, recentre a resposta nesse valor. O intercepto passa a ser a diferença média relativamente a 3, e o intervalo é o intervalo dessa diferença:

mu0 <- 3
m_ref <- lm(I(Love4Taylor - mu0) ~ 1, data = taylor)

coef(m_ref)
confint(m_ref)

Uma estimativa positiva indica uma média acima de 3; uma negativa indica uma média abaixo de 3. O intervalo mostra quais diferenças relativamente a 3 são compatíveis com o modelo.

Solução 2: Comparação Planeada Entre Espécies

Antes de consultar as médias, podemos escolher comparar Adelie com Gentoo. Ajustamos o modelo e pedimos esse contraste, em vez de calcular todas as comparações por pares:

m_especie <- lm(body_mass_g ~ species, data = dados)
emm_especie <- emmeans(m_especie, ~ species)
levels(emm_especie)

contraste_planeado <- contrast(
  emm_especie,
  list("Adelie - Gentoo" = c(1, 0, -1))
)
contraste_planeado

O vector de pesos corresponde à ordem dos níveis mostrada por levels(emm_especie), que neste conjunto é Adelie, Chinstrap, Gentoo. Esta comparação é planeada porque a pergunta foi escolhida antes de consultar os resultados; outras comparações podem ser feitas para exploração, mas não transformam retrospectivamente essa exploração num plano anterior.

Solução 3: Declive e Recta Ajustada

m_barbatana <- lm(body_mass_g ~ flipper_length_mm, data = dados)

coef(m_barbatana)["flipper_length_mm"]
confint(m_barbatana)["flipper_length_mm", ]

ggplot(dados, aes(flipper_length_mm, body_mass_g)) +
  geom_point(alpha = 0.45) +
  geom_smooth(method = "lm", se = TRUE) +
  theme_classic() +
  labs(
    x = "Comprimento da barbatana (mm)",
    y = "Massa corporal (g)"
  )

O declive está em gramas por milímetro: o seu sinal indica a direcção da mudança prevista e o intervalo quantifica a incerteza. Interprete-o dentro do intervalo de comprimentos observado, não como uma previsão para qualquer comprimento.

Solução 4: Interacção Entre Espécie e Sexo

m_aditivo <- lm(body_mass_g ~ species + sex, data = dados)
m_interaccao <- lm(body_mass_g ~ species * sex, data = dados)

anova(m_aditivo, m_interaccao)

emmeans(
  m_interaccao,
  list(
    pairwise ~ species,
    ~ sex,
    pairwise ~ species | sex,
    pairwise ~ sex | species
  )
)

As comparações sem | resumem diferenças médias sobre o outro factor. As comparações condicionais depois de | perguntam pela diferença dentro de cada nível. Uma única diferença média entre sexos pode esconder que a sua magnitude, ou mesmo a sua direcção, varia entre espécies.

Solução 5: Declives Comuns ou Específicos por Espécie

m_paralelo <- lm(
  body_mass_g ~ species + flipper_length_mm,
  data = dados
)
m_declives <- lm(
  body_mass_g ~ species * flipper_length_mm,
  data = dados
)

anova(m_paralelo, m_declives)

trends <- emtrends(
  m_declives,
  ~ species,
  var = "flipper_length_mm"
)
summary(trends, infer = c(TRUE, TRUE))
summary(pairs(trends), infer = c(TRUE, TRUE))

No modelo paralelo, todas as espécies partilham um declive. No modelo com interacção, emtrends() estima um declive por espécie e compara-os; se estes forem importantes para a pergunta, uma única diferença ajustada entre espécies não resume todos os comprimentos de barbatana.

13 Recursos

14 Síntese

No trabalho prático, não precisamos de reconstruir a teoria de cada teste. Precisamos de escrever claramente a pergunta, ajustar um modelo que a represente, pedir a estimativa ou contraste que responde a essa pergunta e mostrar a incerteza e o padrão nos dados.

Quando surgir uma dúvida sobre por que duas análises são equivalentes, por que um factor precisa de contrastes, o que uma interacção muda ou o que distingue Tipo I, II e III, essa explicação pertence ao capítulo teórico.

0 Introdução

1 Preparação

2 Uma Média

Pergunta

Modelo

Estimativa e Comparação

3 Dois Grupos

Pergunta

Interpretação

4 Preditor Quantitativo

Pergunta

emmeans() Sem Grupos: ~ 1

Interpretação

5 Vários Preditores

Pergunta

Interpretação

6 Vários Grupos

Pergunta

Interpretação

7 Interacção Entre Factores

Pergunta

Interpretação

8 ANCOVA

Pergunta

Interpretação

9 Interacção Quantitativa

Estimativa, Gráfico e Interpretação

10 Diagnóstico

11 Da Pergunta ao Modelo

12 Exercícios

13 Recursos

14 Síntese

`emmeans()` Sem Grupos: `~ 1`