Um Preditor

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

João O. Santos

joao.santos@ispa.pt

ISPA

Da Média à Recta

\[ \underbrace{\widehat{Y}_i=b_0}_{\text{mesma previsão para todos}} \quad\longrightarrow\quad \underbrace{\widehat{Y}_i=b_0+b_1X_i}_{\text{previsão varia com }X} \]

O que acrescenta o declive?

Preparar os Dados

Usamos os mesmos casos completos para comparar os dois preditores.

library(car)
library(dplyr)
library(emmeans)
library(ggplot2)
library(palmerpenguins)
library(parameters)

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

ds <- penguins |>
  filter(!is.na(sex), !is.na(body_mass_g), !is.na(flipper_length_mm))

Um Preditor Contínuo

Pinguins com barbatanas maiores têm maior massa?

Ler a Recta

\[ \widehat{\text{massa}}_i=b_0+b_1\,\text{barbatana}_i \]

\(b_0\): massa prevista quando o comprimento é zero.

\(b_1\): mudança na massa prevista por mais 1 mm de barbatana, em g/mm.

Os \(b\) são estimativas dos parâmetros \(\beta\) do modelo.

O zero está fora dos comprimentos observados: não interpretar o intercepto como a massa de um animal real sem barbatana. Centrar o comprimento, por exemplo subtraindo 200 mm, permitiria ler o intercepto nesse comprimento, sem mudar o declive nem as previsões.

A letra b representa uma estimativa; a recta prevê a média condicional, não a massa exacta de cada pinguim. Não estamos ainda a distinguir espécies: a associação conjunta não tem de ser igual à associação dentro de cada espécie. Essa pergunta reaparece na aula de vários preditores.

Ajustar e Estimar

m_flipper <- lm(body_mass_g ~ flipper_length_mm, data = ds)
summary(m_flipper)

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

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

model_parameters(m_flipper)

Parameter         | Coefficient |     SE |               95% CI | t(331) |      p
---------------------------------------------------------------------------------
(Intercept)       |    -5872.09 | 310.29 | [-6482.47, -5261.71] | -18.92 | < .001
flipper length mm |       50.15 |   1.54 | [   47.12,    53.18] |  32.56 | < .001

O Que a Recta Prevê

Mais 20 mm → diferença prevista de \(20b_1\) gramas.

O Que a Média Prevê

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

(Intercept) 
   4207.057 

A mesma massa prevista para todos os comprimentos.

Duas Previsões

Comparar com Uma Média

Demonstração do teste do declive, não um passo extra da análise. anova() serve aqui para fins didácticos.

anova(m0, m_flipper)

Analysis of Variance Table

Model 1: body_mass_g ~ 1
Model 2: body_mass_g ~ flipper_length_mm
  Res.Df       RSS Df Sum of Sq      F    Pr(>F)    
1    332 215259666                                  
2    331  51211963  1 164047703 1060.3 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Restrição: \(\beta_1=0\). Testada pelo t do declive e, neste modelo, pelo F global. O código mostra o significado desses testes.

Recuperar a comparação compacto/aumentado, sem repetir a derivação de F nem as tabelas de RSS da aula de comparação de modelos. São as mesmas observações em ambos os modelos; só a regra de previsão muda. O modelo compacto fixa o parâmetro beta1 em zero. O teste não pergunta se o declive estimado b1 é literalmente zero; avalia a hipótese sobre beta1.

Associação ou Causalidade?

Randall Munroe, xkcd 552: Correlation, CC BY-NC 2.5.

O comprimento foi medido ou manipulado?

Agora: Declive

20 minutos · grupos de 3 ou 5

Modelos Lineares no R — Exercícios

Façam exercício 3. Identifiquem a unidade do declive e o que a comparação com o modelo só com intercepto permite concluir.

Script Quantitativo

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

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

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

Para esta pergunta, o resumo basta: declive, erro-padrão e teste. O modelo nulo e a tabela ANOVA ficaram na demonstração; não precisamos de repetir o mesmo teste de várias maneiras no script.

Resultados do Declive

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

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

Ler o Primeiro Resumo

summary(m_flipper)

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

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

model_parameters(m_flipper)

Parameter         | Coefficient |     SE |               95% CI | t(331) |      p
---------------------------------------------------------------------------------
(Intercept)       |    -5872.09 | 310.29 | [-6482.47, -5261.71] | -18.92 | < .001
flipper length mm |       50.15 |   1.54 | [   47.12,    53.18] |  32.56 | < .001

Anova(m_flipper)

Anova Table (Type II tests)

Response: body_mass_g
                     Sum Sq  Df F value    Pr(>F)    
flipper_length_mm 164047703   1  1060.3 < 2.2e-16 ***
Residuals          51211963 331                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Localizem declive e EP, depois \(R^2\) e F. O teste do intercepto não responde à pergunta sobre barbatanas maiores.

Outra Tabela Para Quê?

Para testar este declive, a informação repete-se.

Saída	Hipótese
t de flipper_length_mm no resumo	Declive zero
F de flipper_length_mm na ANOVA	Declive zero

Neste modelo, \(F=t^2\), com o mesmo valor-P bilateral. Não precisamos de relatar ambos.

Porquê aprender a ler a tabela ANOVA? Com três espécies, um único factor precisará de dois coeficientes. Aí, um t de cada vez já não responde à pergunta sobre o factor inteiro.

Relatar a Associação

Nos pinguins observados, a massa prevista aumentou 50.2 g por mm de barbatana, 95% CI [47.1, 53.2] g/mm; F_{1, 331} = 1060.3, P < .001. É uma associação conjunta das três espécies, não causal nem ajustada para espécie; a generalização depende da amostragem.

E se X For um Nome?

A massa corporal difere entre os dois sexos observados?

Mantemos os mesmos pinguins e a mesma resposta. Cada linha representa um espécime; repetir o rótulo do sexo não significa repetir a unidade. A comparação é entre grupos observados, não entre tratamentos atribuídos.

Como Entram os Nomes?

\[ \widehat{Y}=b_0+b_1X \]

Grupo	\(X\)	Previsão
female	0	\(b_0\)
male	1	\(b_0+b_1\)

Codificação por contrastes: os nomes ficam nos dados; o modelo usa números.

O Zero Define a Referência

\[ b_0=\widehat{\mu}_{\mathrm{female}},\qquad b_1=\widehat{\mu}_{\mathrm{male}}-\widehat{\mu}_{\mathrm{female}} \]

# Treatment coding only for this comparison.
m_sex01 <- lm(body_mass_g ~ sex, data = ds,
              contrasts = list(sex = "contr.treatment"))
summary(m_sex01)

Call:
lm(formula = body_mass_g ~ sex, data = ds, contrasts = list(sex = "contr.treatment"))

Residuals:
    Min      1Q  Median      3Q     Max 
-1295.7  -595.7  -237.3   737.7  1754.3 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  3862.27      56.83  67.963  < 2e-16 ***
sexmale       683.41      80.01   8.542  4.9e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 730 on 331 degrees of freedom
Multiple R-squared:  0.1806,    Adjusted R-squared:  0.1781 
F-statistic: 72.96 on 1 and 331 DF,  p-value: 4.897e-16

Esta é a codificação de tratamento, o padrão inicial do R para factores não ordenados. Pedimo-la explicitamente só neste modelo: a opção definida no início da aula usa soma-zero. Neste conjunto, female é o primeiro nível. A coluna sexmale é o X numérico da tabela anterior. Substituir 0 e 1 na equação antes de interpretar a saída.

E se o Zero Ficar a Meio?

Grupo	\(X\)	Previsão
female	\(+1\)	\(b_0+b_1\)
male	\(-1\)	\(b_0-b_1\)

\[ b_0=\frac{\widehat{\mu}_{\mathrm{female}}+\widehat{\mu}_{\mathrm{male}}}{2}, \qquad 2b_1=\widehat{\mu}_{\mathrm{female}}-\widehat{\mu}_{\mathrm{male}} \]

Intercepto: ponto médio das duas médias. Coeficiente: metade da diferença.

Os símbolos b0 e b1 designam os coeficientes da codificação agora mostrada: os seus valores mudam relativamente a 0/1. Este intercepto é a média não ponderada das duas médias de grupo, não necessariamente a média de todos os indivíduos quando os grupos têm tamanhos diferentes.

Escolher a Codificação no R

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

       [,1]
female    1
male     -1

m_sex <- lm(body_mass_g ~ sex, data = ds)
coef(m_sex)

(Intercept)        sex1 
  4203.9786   -341.7059 

Mudam os coeficientes, não as médias previstas nem os resíduos.

Ler as linhas female e male da matriz antes de ler os coeficientes. contr.sum codifica factores não ordenados; contr.poly aplica contrastes polinomiais aos ordenados. A opção vale para os modelos ajustados a partir daí nesta sessão, salvo codificações explícitas. Aqui, contr.sum(2) dá +1 a female e −1 a male: sex1 é metade de female − male. Não atribuir uma direcção à diferença sem olhar para a matriz. Preparação e interpretação.

Estimar as Médias

emm_sex <- emmeans(m_sex, ~ sex)
emm_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 

A massa prevista de cada sexo, em gramas, com erro-padrão e IC 95%.

Com dois grupos há um único par. O teste da diferença já está no resumo; não precisamos de pedir comparações múltiplas.

Dois Grupos, a Mesma Comparação

Demonstração da linha sex na tabela ANOVA. O script não precisa de ajustar novamente o modelo só com intercepto.

anova(m0, m_sex)

Analysis of Variance Table

Model 1: body_mass_g ~ 1
Model 2: body_mass_g ~ sex
  Res.Df       RSS Df Sum of Sq      F    Pr(>F)    
1    332 215259666                                  
2    331 176380769  1  38878897 72.961 4.897e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Restrição: \(\beta_1=0\) → médias iguais.

É o teste-t de Student: \(F=t^2\). Ainda temos só um coeficiente para o factor: a redundância mantém-se com dois grupos.

Não voltar a apresentar um teste-t como método independente antes do modelo: a novidade é reconhecer a comparação de grupos dentro da equação já conhecida. Dois Grupos.

Agora: Contraste

10 minutos · grupos de 3 ou 5

Comparação de Modelos — Exercícios

Façam exercício 2: substituam \(X=-1\) e \(X=+1\) na equação e subtraiam as previsões.

E se trocarmos os sinais dos dois grupos?

Fixar A = −1 e B = +1: a diferença B − A é 2b1. Se trocarmos os sinais, o coeficiente muda de sinal, mas as previsões para cada grupo mantêm-se quando voltamos a ajustar o modelo.

Script de Dois Grupos

Outra pergunta, não uma extensão do modelo anterior. O resumo dá o teste da diferença; emmeans() dá as duas médias. A comparação explícita com m0 e a tabela ANOVA ficaram na demonstração.

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

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

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

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

Resultados de Dois Grupos

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

Residuals:
    Min      1Q  Median      3Q     Max 
-1295.7  -595.7  -237.3   737.7  1754.3 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   4204.0       40.0 105.088  < 2e-16 ***
sex1          -341.7       40.0  -8.542  4.9e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 730 on 331 degrees of freedom
Multiple R-squared:  0.1806,    Adjusted R-squared:  0.1781 
F-statistic: 72.96 on 1 and 331 DF,  p-value: 4.897e-16

 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 

Ler o Segundo Resumo

summary(m_sex)

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

Residuals:
    Min      1Q  Median      3Q     Max 
-1295.7  -595.7  -237.3   737.7  1754.3 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   4204.0       40.0 105.088  < 2e-16 ***
sex1          -341.7       40.0  -8.542  4.9e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 730 on 331 degrees of freedom
Multiple R-squared:  0.1806,    Adjusted R-squared:  0.1781 
F-statistic: 72.96 on 1 and 331 DF,  p-value: 4.897e-16

model_parameters(m_sex)

Parameter   | Coefficient |    SE |             95% CI | t(331) |      p
------------------------------------------------------------------------
(Intercept) |     4203.98 | 40.00 | [4125.28, 4282.67] | 105.09 | < .001
sex [1]     |     -341.71 | 40.00 | [-420.40, -263.01] |  -8.54 | < .001

Anova(m_sex)

Anova Table (Type II tests)

Response: body_mass_g
             Sum Sq  Df F value    Pr(>F)    
sex        38878897   1  72.961 4.897e-16 ***
Residuals 176380769 331                      
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

sex1 é metade da diferença female − male. O sinal depende da codificação; o teste global pergunta se as duas médias são iguais.

Ler as Duas Médias

emm_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 

Que massa estimamos para cada sexo? Com que incerteza?

Podemos subtrair as médias para perceber a diferença. Não precisamos de interpretar o coeficiente soma-zero para recuperar cada média.

Relatar os Dois Grupos

As fêmeas apresentaram menor massa corporal (M = 3862 g, SE = 56.8 g) do que os machos (M = 4546 g, SE = 56.3 g), t(331) = −8.54, P < .001.

O teste refere-se a female − male. O modelo não ajusta para espécie; estes dados observacionais não estabelecem causalidade.

Estes exemplos ligam o método à pergunta e relatam direcção, magnitude, unidade e estatística de teste. O coeficiente soma-zero é metade da diferença; o seu t coincide com o t da diferença inteira.

Entrega Rápida

10 minutos · grupos de 3 ou 5

Para cada um dos dois preditores, escrevam:

unidade observada · significado de \(b_0\) e \(b_1\) · estimativa e unidade

Na comparação de grupos, indiquem a codificação usada.

Que informação do delineamento falta para interpretar causalmente a associação?

Uma Possível Resposta

Preditor	Intercepto	Coeficiente do Preditor
Barbatana	Massa prevista a 0 mm, fora dos dados	Mudança de massa por mm, em g/mm
Sexo, soma-zero	Média das duas médias de sexo	Metade de female − male, em g

Na recta, \(b_0\approx-5872\) g e \(b_1\approx50.2\) g/mm. Com sexo, \(b_0\approx4204\) g e \(b_1\approx-341.7\) g.

Uma linha = um espécime. No primeiro caso, relatamos o declive e a sua incerteza; no segundo, as médias de cada sexo e o teste da diferença.

Sem conhecer a amostragem e o delineamento, uma associação não identifica uma causa. A espécie também pode ajudar a explicar ambas as relações.

A Seguir

E se X tiver três grupos? E se tivermos mais de um X?