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?

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.

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 = peng)
summary(m_flipper)$coefficients

                     Estimate Std. Error   t value      Pr(>|t|)
(Intercept)       -5872.09268 310.285155 -18.92483  1.183941e-54
flipper_length_mm    50.15327   1.540231  32.56217 3.132836e-105

confint(m_flipper)

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

O Que a Recta Prevê

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

Comparar com Uma Média

m0 <- lm(body_mass_g ~ 1, data = peng)
anova(m0, m_flipper)

Res.Df	RSS	Df	Sum of Sq	F	Pr(>F)
332	215259666	NA	NA	NA	NA
331	51211963	1	164047703	1060.295	0

Restrição: \(b_1=0\). O declive reduz erro suficiente para justificar o parâmetro?

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.

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.

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}} \]

m_sex01 <- lm(body_mass_g ~ sex, data = peng)
summary(m_sex01)$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

Esta é a codificação de tratamento, por defeito no R para factores não ordenados. A preparação fixa female como primeiro nível e explicita esta codificação, para a demonstração não depender de opções de outra sessão. 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{male}}-\widehat{\mu}_{\mathrm{female}} \]

O intercepto fica a meio das duas médias; o 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

contrasts(peng$sex) <- cbind(half_difference = c(-1, 1))
contrasts(peng$sex)

       half_difference
female              -1
male                 1

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

       (Intercept) sexhalf_difference 
         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. O contraste é soma-zero; o sinal foi escolhido explicitamente para que 2b1 seja male − female. contr.sum(2) usa a ordem de sinais inversa: não atribuir uma direcção à diferença sem olhar para a matriz.

Estimar a Diferença

sex_means <- emmeans(m_sex, ~ sex)
contrast(sex_means, list("male - female" = c(-1, 1))) |>
  confint()

contrast	estimate	SE	df	lower.CL	upper.CL
male - female	683.4118	80.00868	331	526.0222	840.8014

Uma diferença em gramas e o seu IC, qualquer que seja a codificação.

Dois Grupos, a Mesma Comparação

anova(m0, m_sex)

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

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

Com variância residual comum, é o teste-t de Student: \(F=t^2\).

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. A equivalência exacta é com Student, não com o teste de Welch que t.test() usa por defeito. Welch estima a incerteza de outra forma; pode dar outro valor-p e outros graus de liberdade, embora estime a mesma diferença de médias. 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.

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?

A Seguir

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