A Abordagem de Comparação de Modelos

Modelos Lineares
Abordagem de Comparação de Modelos
Como teste-t, regressão, ANOVA, ANCOVA e interacções cabem no mesmo modelo linear
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

Em Fundamentos da Inferência começámos com uma média, o seu erro-padrão, uma estatística t e um valor-p. Agora voltamos aos mesmos dados para dar o passo que organiza o resto da estatística neste site: escrever a pergunta como uma comparação entre modelos.

A pergunta continua a ser se, em média, a população é indiferente à música da Taylor Swift. Na escala Love4Taylor, zero representa indiferença, por isso a hipótese de referência é

\[ H_0:\mu=0. \]

Mais tarde mudaremos os preditores e o número de parâmetros, mas a estrutura mantém-se: que restrição representa a ausência do efeito que queremos estudar e o que ganhamos ao libertá-la?

1 O Modelo Linear

Dados, Modelo e Erro

O mesmo problema pode ser escrito sem começar pelo nome de um teste:

\[ \text{Dados}=\text{Modelo}+\text{Erro}. \]

Para cada observação,

\[ Y_i=\hat Y_i+e_i. \]

O modelo dá-nos uma previsão; o resíduo é o que ficou por explicar. Para comparar modelos lineares, resumimos esses resíduos através da soma dos quadrados:

\[ SSE=\sum_i(Y_i-\hat Y_i)^2. \]

Quanto menor o (SSE), melhor o modelo reproduz os dados. Mas um modelo com mais parâmetros tem mais liberdade para reduzir o erro. Portanto, a pergunta não pode ser apenas «qual tem menos SSE?». Precisamos de comparar a melhoria com a complexidade que acrescentámos e com o erro que ainda ficou por explicar.

No modelo linear podemos escrever todas as observações de uma vez como

$ y=X+. $

Os mínimos quadrados escolhem () para minimizar

$ (y-X)^T (y-X). $

Quando a matriz do modelo tem posto completo, este problema é quadrático e convexo: existe uma solução única. A expressão clássica é

$ =(X^TX)^{-1}X^Ty. $

Esta fórmula é óptima para perceber a álgebra, mas software estatístico sério não precisa de formar explicitamente essa inversa. Métodos como decomposição QR são numericamente mais estáveis e encontram a mesma solução de mínimos quadrados.

Esta maquinaria é útil para quem quiser perceber como lm() resolve o problema. Não é necessária para interpretar uma média, um contraste, um declive ou uma interacção.

Dois Modelos Para a Pergunta da Taylor

A nossa hipótese nula fornece imediatamente um modelo compacto.

Modelo Compacto

Se as pessoas forem, em média, indiferentes à Taylor Swift, a previsão é sempre zero:

\[ m_0:\quad \hat Y=0. \]

Este modelo não estima a média a partir dos dados. Impõe a restrição que queremos testar.

Modelo Aumentado

O modelo aumentado liberta essa restrição e estima um intercepto:

\[ m_1:\quad \hat Y=\beta_0. \]

Num modelo apenas com intercepto, a estimativa de (_0) é a média amostral. Portanto,

\[ \hat Y=0.57. \]

A pergunta estatística passou a ter uma forma muito concreta:

Vale a pena abandonar a previsão fixa de zero e estimar a média a partir dos dados?

Em R, os dois modelos são:

m0 <- lm(Love4Taylor ~ 0, data = love)
m1 <- lm(Love4Taylor ~ 1, data = love)

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

Do t ao F

O t olhou para um coeficiente dividido pelo seu erro-padrão. A comparação de modelos F faz a mesma pergunta numa forma que se generaliza quando libertamos mais do que um parâmetro.

Para dois modelos lineares aninhados,

\[ F= \frac{(SSE_0-SSE_1)/(p_1-p_0)} {SSE_1/(N-p_1)}. \]

Vamos resolver a equação com os nossos dados.

O modelo compacto tem

\[ SSE_0=74.75, \]

e o modelo aumentado tem

\[ SSE_1=58.78. \]

Libertámos um parâmetro, por isso

\[ MSR = \frac{SSE_0-SSE_1}{1} = \frac{74.75-58.78}{1} = 15.97. \]

O erro que ainda ficou no modelo aumentado é

\[ MSE = \frac{SSE_1}{N-p_1} = \frac{58.78}{49} = 1.2. \]

Logo,

\[ F = \frac{MSR}{MSE} = \frac{15.97}{1.2} = 13.32. \]

Neste exemplo os modelos diferem por um único parâmetro. Por isso,

\[ F=t^2. \]

Podemos verificar:

\[ 13.32 \approx (3.65)^2. \]

Esta equivalência é importante, mas o F é a ideia mais geral. Quando um factor com quatro níveis exige três parâmetros, ou quando comparamos dois modelos que diferem em vários termos, já não existe um único t que faça todo o trabalho. A comparação F continua exactamente com a mesma lógica: quanto erro desapareceu por termos libertado a restrição, relativamente ao erro que ficou?

O F Usa a Mesma Lógica de Inferência

Depois de calcularmos o F, não precisamos de aprender uma nova definição de valor-p. O software pergunta quão incompatível é um F pelo menos tão grande como o observado com o modelo compacto, usando a distribuição F adequada aos graus de liberdade.

A interpretação é a que já construímos em Do t ao Valor-p. O que mudou foi a estatística cuja distribuição nula usamos, não o significado do valor-p.

Quando a comparação acrescenta apenas um parâmetro, (F=t^2), e os testes t e F dão o mesmo valor-p.

ImportanteO Que Deve Ficar Deste Primeiro Exemplo

Não precisamos de decorar três técnicas independentes.

Primeiro perguntámos quanto a nossa estimativa varia entre amostras. Isso levou ao erro-padrão e ao t.

Depois escrevemos a hipótese como uma restrição num modelo. Libertar essa restrição reduziu o erro, e o F quantificou essa melhoria tendo em conta a complexidade adicional e o erro residual.

A partir daqui vamos mudar a pergunta de investigação, mas não vamos mudar de filosofia.

2 Um Modelo, Vários Nomes

Dois Grupos: O Primeiro Contraste

Pergunta de investigação: há uma diferença média de massa corporal entre os dois grupos de pinguins que estamos a comparar?

Agora queremos duas médias. Em vez de inventar uma matemática nova, acrescentamos uma coluna ao modelo que diga a que grupo pertence cada observação.

Vamos começar com a codificação que usaremos ao longo deste capítulo:

Grupo Contraste \(X\)
grupo A -1
grupo B +1

O modelo é

\[ \hat Y=\beta_0+\beta_1X. \]

Façamos as duas substituições.

Para o grupo A, \(X=-1\):

\[ \hat Y_A=\beta_0-\beta_1. \]

Para o grupo B, \(X=+1\):

\[ \hat Y_B=\beta_0+\beta_1. \]

Agora os coeficientes deixam de ser símbolos abstractos:

  • \(\beta_0\) fica exactamente a meio das duas médias;
  • \(\beta_1\) é metade da diferença entre B e A;
  • portanto, a diferença B - A é \(2\beta_1\).

Se \(\beta_1=0\), as duas previsões colapsam para a mesma média. É por isso que testar este coeficiente é testar a diferença entre os dois grupos.

O Factor Tem de Ser Traduzido em Números

No R, sex pode aparecer como palavras como female e male, mas o modelo acaba por precisar de uma coluna numérica. Essa tradução chama-se codificação por contrastes. Não é um detalhe de implementação escondido no computador: decide o que o intercepto e os coeficientes significam.

Para tornar a nossa comparação explícita:

penguins <- palmerpenguins::penguins
penguins <- subset(penguins, !is.na(body_mass_g) & !is.na(sex))
penguins$sex <- factor(penguins$sex, levels = c("female", "male"))

contrasts(penguins$sex) <- matrix(
  c(-1, 1),
  ncol = 1,
  dimnames = list(c("female", "male"), "male_vs_female")
)

contrasts(penguins$sex)
       male_vs_female
female             -1
male                1

E então ajustamos o modelo:

Código
m0 <- lm(body_mass_g ~ 1, data = penguins)
m1 <- lm(body_mass_g ~ sex, data = penguins)

anova(m0, m1)
Res.Df RSS Df Sum of Sq F Pr(>F)
332 215259666 NA NA NA NA
331 176380769 1 38878897 72.96099 0
Código

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

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 ***
sexmale_vs_female    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

Voltámos a acrescentar apenas um parâmetro. Por isso, neste caso de dois grupos, o F omnibus e o t desse coeficiente satisfazem \(F=t^2\).

A codificação de tratamento (treatment coding) usa normalmente 0 para um grupo de referência e 1 para o outro. Nesse caso, o intercepto é a média do grupo 0 e o declive é a diferença completa entre os grupos.

O modelo ajustado pode representar exactamente as mesmas duas médias. O que muda é a forma de descrever essas médias através dos coeficientes. Preferimos -1/+1 aqui porque prepara directamente a leitura de efeitos principais e interacções em delineamentos factoriais.

Com os erros-padrão habituais de lm(), y ~ grupo usa uma variância residual comum aos dois grupos. A equivalência é, portanto, com o teste-t de Student com variâncias agrupadas, não com o teste de Welch. O contraste continua a representar a diferença entre os grupos; muda a forma como quantificamos a sua incerteza.

Um Preditor Quantitativo: Regressão Linear

Pergunta de investigação: o número de mortes de peixes-boi aumenta à medida que aumenta o número de barcos a motor registados?

Troquemos agora a coluna -1/+1 por uma variável que pode assumir muitos valores. A equação não muda:

\[\hat Y=\beta_0+\beta_1X.\]

A diferença é a interpretação. \(\beta_1\) deixa de representar uma única diferença entre dois grupos e passa a representar a mudança prevista em \(Y\) por cada unidade adicional de \(X\).

Código
manatees <- read.csv("../data/manatees.csv")

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

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
Código

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

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

A comparação pergunta se permitir um declive reduz o erro face ao modelo que prevê apenas a média. Outra vez, acrescentámos um parâmetro, portanto o F da comparação é o quadrado do t do declive.

Correlação como Caso Simples

Numa regressão simples com intercepto, se padronizarmos \(X\) e \(Y\), o declive é a correlação de Pearson \(r\). Além disso, nesse caso simples, \(R^2=r^2\).

Isto é útil pedagogicamente porque mostra que «correlação» e «regressão» não vivem em universos matemáticos separados. Continuam, contudo, a responder a formas de apresentação diferentes e nenhuma delas, por si, estabelece causalidade.

Duas Tradições, o Mesmo Modelo Linear

É frequente a regressão ser ensinada com palavras como previsão, declive, intercepto, resíduo e variável explicativa. Essa linguagem é natural em muitas aplicações quantitativas, incluindo engenharia, calibração, previsão e modelação de processos. A tradição da ANOVA cresceu muito ligada ao delineamento experimental, onde falamos mais de factores, níveis, efeitos, interacções e decomposição da variância. O desenvolvimento de ANOVA e do delineamento experimental ficou particularmente associado ao trabalho de Fisher no início do século XX.

Estas tradições levaram a vocabulários diferentes, e esses vocabulários podem ser úteis porque destacam aspectos diferentes do problema. Matematicamente, porém, grande parte da ANOVA introdutória pode ser escrita como regressão com preditores categóricos. Não precisamos de escolher entre «ser pessoa de ANOVA» e «ser pessoa de regressão»: precisamos de perceber o modelo.

A história completa é muito mais longa do que «regressão veio da engenharia e ANOVA veio da Psicologia». Regressão e mínimos quadrados têm raízes em astronomia, geodesia e biometria; ANOVA foi desenvolvida no contexto da estatística matemática e do delineamento experimental, com enorme influência dos trabalhos de Fisher. O ponto pedagógico aqui é apenas que as tradições de uso deixaram linguagens diferentes para a mesma família de modelos.

Vários Preditores: Até Aqui, Tudo É Aditivo

Pergunta de investigação: que características da água ainda ajudam a prever o oxigénio quando as restantes já estão no modelo?

Se acrescentarmos preditores quantitativos,

\[\hat Y=\beta_0+\beta_1X_1+\beta_2X_2+\beta_3X_3+\ldots\]

cada declive descreve uma associação condicionada aos restantes preditores. É isto que normalmente queremos dizer por «controlar estatisticamente para» as outras variáveis. Não significa que tenhamos removido confundimento do mundo real nem que a associação se tenha tornado causal.

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

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

Call:
lm(formula = Oxygen ~ pH + Temp + NTU + Conductivity, data = water)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.55259 -0.17601  0.03673  0.26211  1.62329 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -19.635139   1.591156 -12.340  < 2e-16 ***
pH             2.626421   0.257724  10.191  < 2e-16 ***
Temp          -0.016087   0.024754  -0.650    0.516    
NTU           -0.417950   0.053411  -7.825 3.09e-14 ***
Conductivity   0.032833   0.002212  14.846  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.4602 on 495 degrees of freedom
Multiple R-squared:  0.6894,    Adjusted R-squared:  0.6869 
F-statistic: 274.7 on 4 and 495 DF,  p-value: < 2.2e-16

Até aqui todos os termos são aditivos. O declive de NTU é o mesmo em todos os valores de Temp; o declive de Temp é o mesmo em todos os valores de NTU. Mais à frente vamos libertar precisamente essa restrição.

O Teste Global e os Testes Parciais

Há pelo menos duas perguntas diferentes:

  1. O conjunto dos preditores melhora o modelo de uma só média?
  2. Cada preditor continua a acrescentar informação quando todos os outros já estão no modelo?

A primeira compara Oxygen ~ 1 com o modelo completo. A segunda exige uma comparação diferente para cada termo.

Código
m_intercept <- lm(Oxygen ~ 1, data = water)
anova(m_intercept, 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

Um Factor com Mais de Dois Níveis: ANOVA

Pergunta de investigação: a massa corporal média é a mesma nas três espécies de pinguins, ou precisamos de permitir médias diferentes?

Com dois grupos bastava uma coluna -1/+1. Com três grupos já não conseguimos representar todas as diferenças com uma única coluna. Um factor com \(k\) níveis precisa de \(k-1\) contrastes independentes.

Para três espécies, por exemplo, podemos usar contrastes que somam zero:

Espécie \(C_1\) \(C_2\)
Adelie 1 0
Chinstrap 0 1
Gentoo -1 -1

O modelo torna-se

\[ \hat Y=\beta_0+\beta_1C_1+\beta_2C_2. \]

Neste caso:

  • \(\beta_0\) é a média das três médias de espécie;
  • \(\beta_1\) diz quanto a média de Adelie se afasta dessa média geral;
  • \(\beta_2\) diz quanto a média de Chinstrap se afasta dela;
  • o desvio de Gentoo fica determinado pelos outros dois, porque os desvios têm de somar zero.

Isto é mais importante do que parece. O factor não entra no modelo como a palavra species. Entra através destas colunas numéricas. É por isso que os contrastes fazem parte da definição do modelo.

No R podemos tornar essa escolha explícita:

penguins <- palmerpenguins::penguins
penguins <- subset(penguins, !is.na(body_mass_g) & !is.na(species))
penguins$species <- factor(
  penguins$species,
  levels = c("Adelie", "Chinstrap", "Gentoo")
)

contrasts(penguins$species) <- contr.sum(3)
contrasts(penguins$species)
          [,1] [,2]
Adelie       1    0
Chinstrap    0    1
Gentoo      -1   -1

Depois ajustamos exactamente o modelo que já conhecemos:

Código
m0 <- lm(body_mass_g ~ 1, data = penguins)
m1 <- lm(body_mass_g ~ species, data = penguins)

anova(m0, m1)
Res.Df RSS Df Sum of Sq F Pr(>F)
341 219307697 NA NA NA NA
339 72443483 2 146864214 343.6263 0

O teste omnibus pergunta se as duas colunas de contraste, em conjunto, acrescentam informação. Se ambas puderem ser tratadas como zero, as três médias previstas colapsam para uma só.

Depois do omnibus vem a pergunta que realmente costuma interessar: que comparação entre espécies queríamos fazer? Podemos comparar pares, mas também podemos escrever contrastes planeados que correspondam melhor à hipótese. O factor dá-nos o espaço de comparações; não nos obriga a fazer todas as comparações possíveis.

Com codificação de tratamento, um nível é a referência e os coeficientes comparam os outros níveis com essa referência. Com contr.sum, o intercepto fica centrado nas médias dos níveis e os coeficientes obedecem a uma restrição de soma zero.

Quando as duas parametrizações descrevem o mesmo espaço de modelos, podem produzir as mesmas previsões ajustadas. O que muda é que comparação aparece directamente em cada coeficiente. Isso torna-se crucial quando há interacções e quando pedimos testes de efeitos de ordem inferior.

3 Interacções e ANCOVA

Dois Factores: A Interacção Já Faz Parte da Pergunta

Pergunta de investigação: a diferença associada a um factor mantém-se igual em todos os níveis do outro, ou muda consoante o contexto?

Suponhamos dois factores, A e B, ambos com dois níveis. Usando -1/+1, cada célula recebe também o produto dos dois contrastes:

A B A:B
-1 -1 +1
+1 -1 -1
-1 +1 -1
+1 +1 +1

O modelo completo é

\[ \hat Y=\beta_0+\beta_AA+\beta_BB+\beta_{AB}AB. \]

Esta tabela parece uma pequena convenção de codificação, mas contém a ideia da ANOVA factorial inteira.

  • \(2\beta_A\) é a diferença média entre os níveis de A, fazendo a média sobre B;
  • \(2\beta_B\) é a diferença média entre os níveis de B, fazendo a média sobre A;
  • \(4\beta_{AB}\) é a diferença de diferenças.

A escala dos coeficientes vem da escolha -1/+1. Se usássemos outra codificação, os números dos \(\beta\) podiam mudar; as médias previstas por célula e a comparação substantiva não precisavam de mudar.

O Modelo Sem Interacção É uma Restrição

O modelo

y ~ A + B

não significa apenas «ainda não testámos a interacção». Ele afirma que a diferença entre A=-1 e A=+1 é a mesma nos dois níveis de B. As duas partes somam-se de forma aditiva.

O modelo

y ~ A * B

expande para

y ~ A + B + A:B

e liberta essa restrição.

Por isso, quando a teoria permite ou prevê uma interacção, o modelo com interacção é o modelo que representa a teoria. Podemos comparar o modelo aditivo com ele para testar a interacção:

m_add <- lm(y ~ A + B, data = dados)
m_int <- lm(y ~ A * B, data = dados)

anova(m_add, m_int)

Mas não tratamos a interacção como uma decoração que só entra depois de termos garantido efeitos principais. Se a pergunta importante é uma diferença de diferenças, o delineamento e o tamanho da amostra também devem ser pensados para essa diferença de diferenças.

Interacção Entre Dois Factores

Uma interacção diz simplesmente:

a diferença associada a A depende do nível de B.

Num 2 × 2, podemos escrevê-la directamente:

\[ (\bar Y_{A_2B_2}-\bar Y_{A_1B_2})- (\bar Y_{A_2B_1}-\bar Y_{A_1B_1}). \]

Se este valor for zero, a diferença de A é a mesma nos dois níveis de B. Se não for, mudou.

É aqui que a tabela por células volta a ser mais útil do que a linguagem abstracta de «efeitos». Antes de interpretar um valor-p, escrevam ou desenhem as quatro médias. Perguntem: que diferença mudou, em que direcção e quanto?

ImportanteUm Efeito Principal É um Resumo

Com -1/+1, \(\beta_A\) resume A fazendo a média sobre os níveis de B. Isso pode ser uma pergunta substantiva perfeitamente legítima. Mas, quando a interacção é grande, esse resumo pode esconder duas diferenças muito distintas.

Por isso, uma tabela ANOVA não termina a interpretação de uma interacção. Precisamos normalmente das médias por célula, efeitos simples ou contrastes que digam exactamente onde está a diferença de diferenças.

ANCOVA: Factores e uma Covariável Quantitativa

Pergunta de investigação: depois de ter em conta uma variável quantitativa relevante, continuam a existir diferenças entre os grupos — e podemos realmente assumir que essa variável tem o mesmo declive em todos eles?

ANCOVA não exige uma nova matemática. É regressão linear com pelo menos um preditor categórico e um quantitativo.

A versão tradicional mais simples é:

y ~ grupo + x

Aqui x é a covariável. O modelo permite grupos com interceptos diferentes, mas obriga todos os grupos a partilhar o mesmo declive de x. É a famosa hipótese de declives paralelos.

ANCOVA Factorial Tradicional

Podemos ter dois factores categóricos e permitir a sua interacção, mantendo a covariável apenas como efeito aditivo:

y ~ A * B + x

Este modelo permite:

  • diferenças entre níveis de A;
  • diferenças entre níveis de B;
  • uma interacção A × B;
  • uma associação linear com x;
  • mas o mesmo declive de x em todas as células A × B.

Esta é uma forma muito comum de ANCOVA. É importante perceber que «não incluir interacções com a covariável» não é uma propriedade matemática obrigatória da ANCOVA. É uma restrição do modelo.

Quando Essa Convenção Pode Ser Má Ideia

Se houver razão científica para pensar que a relação entre x e a resposta muda entre grupos, impor um único declive pode esconder precisamente o padrão que nos interessa. Podemos então comparar o modelo anterior com um modelo que liberta essas interacções.

Por exemplo, para uma covariável e dois factores:

y ~ A * B * x

inclui A × x, B × x e A × B × x, além dos termos de ordem inferior. Não somos obrigados a incluir toda a estrutura se a pergunta justificar apenas uma parte, mas devemos respeitar a hierarquia dos termos e dizer claramente que modelo estamos a testar.

Aviso«Controlar pela Covariável» Não Resolve Causalidade

Acrescentar x altera a quantidade que estimamos. Pode reduzir variação residual e melhorar precisão; pode também introduzir ou agravar viés em certos cenários causais. Uma ANCOVA não transforma automaticamente dados observacionais numa experiência.

Interacção Entre um Factor e uma Variável Quantitativa

Consideremos:

y ~ grupo * x

Isto expande para

y ~ grupo + x + grupo:x

A interacção pergunta se o declive de x é igual entre grupos.

Com dois grupos, podemos imaginar:

\[ \hat Y_{grupo0}=\beta_0+\beta_1x \]

e

\[ \hat Y_{grupo1}=(\beta_0+\beta_2)+(\beta_1+\beta_3)x. \]

\(\beta_3\) é então a diferença entre declives. Se for zero, as rectas são paralelas. Se não for, a associação entre x e y depende do grupo.

Este é o mesmo conceito que vimos na interacção entre dois factores. Mudou a geometria da variável, não a ideia de «o efeito de um preditor depende do outro».

Interacção Entre Duas Variáveis Quantitativas

Finalmente:

\[ Y=\beta_0+\beta_1X+\beta_2Z+\beta_3XZ+e. \]

ou, em R:

y ~ x * z

O declive de x é

\[ \frac{\partial \hat Y}{\partial X}=\beta_1+\beta_3Z. \]

Portanto, quando existe uma interacção, não há um único declive de X. O declive depende do valor de Z. E, simetricamente, o declive de Z depende de X.

É por isso que interacções quantitativas costumam ser apresentadas através de declives simples em valores substantivamente úteis do moderador. Média, média ± 1 DP e quantis são convenções possíveis, não leis. Se um desses valores ficar fora da região onde há dados, a interpretação passa a depender de extrapolação.

Se substituirmos \(X\) por \(X-c\), o modelo continua a representar a mesma superfície quando os termos correspondentes são incluídos. O que muda é o ponto onde interpretamos o intercepto e o «efeito principal» de outro preditor. Centrar numa média pode tornar coeficientes mais interpretáveis, mas não é necessário para que uma interacção exista nem a torna estatisticamente válida.

4 Testes no Modelo Completo

A Nossa Convenção Para Tabelas ANOVA

As expressões Tipo I, Tipo II e Tipo III não são três botões que calculam a mesma pergunta de maneiras diferentes. Podem testar hipóteses diferentes.

Neste site vamos usar como regra geral Tipo III para modelos factoriais em que as interacções fazem parte do modelo científico, com contrastes centrados (soma zero) para os factores.

A razão é a mesma que já usamos em regressão múltipla: quando queremos testar um termo, os restantes termos do modelo podem funcionar como nuisance parameters para essa pergunta: não são o alvo do teste, mas continuam a ser estimados. Não precisamos de os declarar iguais a zero só porque não são o alvo daquele teste.

Se o modelo é

y ~ A * B

e estamos a perguntar por A, queremos saber se há uma diferença média associada a A enquanto o modelo continua a permitir diferenças de B e uma interacção A:B. Não queremos testar A num modelo provisório que primeiro obriga A:B a desaparecer.

É esta a ideia central da nossa preferência por Tipo III: cada termo é avaliado no contexto do modelo completo que decidimos que representa a hipótese científica.

Tipo III Como Teste no Modelo Completo

Num 2 × 2 com a codificação -1/+1 que acabámos de usar, a leitura é particularmente limpa:

  • a linha A testa \(\beta_A=0\): a diferença de A em média sobre B;
  • a linha B testa \(\beta_B=0\): a diferença de B em média sobre A;
  • a linha A:B testa \(\beta_{AB}=0\): a diferença de diferenças.

Ao testar A, por exemplo, \(\beta_B\) e \(\beta_{AB}\) continuam no modelo. Para essa pergunta funcionam como nuisance parameters: não são o alvo do teste, mas não os obrigamos a desaparecer.

A analogia com «controlar» outros termos numa regressão é útil, com uma cautela: a interacção não é um confundidor. É parte da estrutura que a nossa teoria permite. O ponto comum é simplesmente não fingir que os outros termos são zero para conseguir testar aquele que nos interessa.

No R, a escolha dos contrastes tem de ser explícita para que os testes de ordem inferior representem estas médias marginais:

contrasts(dados$A) <- matrix(c(-1, 1), ncol = 1)
contrasts(dados$B) <- matrix(c(-1, 1), ncol = 1)

m <- lm(y ~ A * B, data = dados)
car::Anova(m, type = "III")

Para factores com mais níveis, a mesma ideia aplica-se a um conjunto de contrastes. A linha do factor pergunta se esse conjunto pode ser zero, mantendo no modelo os restantes termos que especificámos.

Porque os Contrastes Encaixam Tão Bem Nesta Lógica

No 2 × 2 equilibrado, as colunas -1/+1 para A e B e a coluna produto A:B são ortogonais entre si e ao intercepto. As quatro médias das células ficam decompostas em quatro peças muito claras:

  1. média geral;
  2. contraste de A;
  3. contraste de B;
  4. contraste A × B, isto é, a diferença de diferenças.

É exactamente a decomposição que queremos interpretar.

Com factores de mais de dois níveis, a codificação soma-zero mantém o factor centrado relativamente ao intercepto. Se queremos decompor os vários graus de liberdade em comparações individuais, podemos escolher contrastes planeados ortogonais quando isso corresponde às hipóteses científicas; os produtos desses contrastes dão as decomposições correspondentes das interacções.

Há três ideias próximas que convém não fundir.

Células equilibradas tornam a geometria clássica da ANOVA especialmente limpa. Com o mesmo número de observações por célula e contrastes mutuamente ortogonais, as componentes do delineamento ficam ortogonais na matriz do modelo.

Contrastes soma-zero são o que precisamos aqui para centrar os factores e dar aos testes Tipo III de ordem inferior a interpretação marginal pretendida. contr.sum() cumpre esta condição, mas, com mais de dois níveis, as suas colunas não são necessariamente ortogonais entre si.

Tipo III não exige células equilibradas. Pelo contrário, uma das razões históricas para a sua utilização é produzir testes que não dependem da frequência relativa das células quando o delineamento é desequilibrado, desde que as funções relevantes continuem estimáveis. Células vazias tornam a situação mais delicada.

Portanto: equilíbrio + contrastes ortogonais dão a decomposição mais limpa; soma-zero dá a parametrização marginal que queremos; Tipo III define a hipótese a testar no modelo completo.

Esta convenção não apareceu aqui do nada. Tipo III é particularmente familiar em muitos fluxos de trabalho experimentais baseados em SPSS e SAS, incluindo muita da tradição de ANOVA ensinada em psicologia. O SPSS usa Tipo III por defeito no GLM e descreve-o como o tipo mais comum; o PROC GLM do SAS também apresenta testes Tipo III por defeito.

Isso explica parte da tradição. A razão para a mantermos aqui, contudo, é conceptual: encaixa bem com a nossa decisão de especificar primeiro o modelo científico completo e testar cada componente sem apagar provisoriamente os outros termos que a teoria permite.

Tipo I é sequencial. Os termos são testados pela ordem em que aparecem. Pode ser exactamente o que queremos quando há uma ordem científica clara, mas também pode tornar a resposta dependente de uma ordem arbitrária da fórmula.

Tipo II testa cada efeito principal depois dos outros efeitos principais, mas não condicionado a uma interacção que o contém. Isto é coerente quando o modelo substantivo é aditivo ou quando queremos explicitamente formular a pergunta supondo ausentes essas interacções de ordem superior.

A diferença para a nossa convenção fica então bastante concreta: Tipo II formula o teste do termo sem o condicionar às interacções que o contêm; Tipo III mantém essas interacções no modelo que define o teste.

A tabela omnibus continua a não saber qual é a pergunta mais interessante do estudo. Se tínhamos uma comparação planeada entre duas células, um efeito simples ou um contraste específico, estimem-no directamente e apresentem a estimativa e a incerteza. Tipo III é a nossa convenção para a tabela; os contrastes substantivos continuam a ser a forma mais directa de responder a hipóteses específicas.

5 Ligações e Extensões

O Que Une Tudo Isto?

A tabela abaixo é uma tradução entre nomes tradicionais e modelos lineares. Não é uma lista de testes que tenham de decorar.

Nome tradicional Modelo linear típico Pergunta central
teste-t de uma amostra y ~ 1 face a uma referência precisamos de estimar um intercepto diferente da referência?
teste-t de duas amostras y ~ grupo um contraste entre dois níveis reduz o erro?
regressão simples y ~ x precisamos de um declive de x?
correlação de Pearson regressão simples padronizada existe associação linear?
regressão múltipla y ~ x1 + x2 + ... que termos acrescentam informação condicionada aos outros?
ANOVA unifactorial y ~ factor o conjunto de contrastes do factor melhora o modelo?
ANOVA factorial aditiva y ~ A + B diferenças médias de A e B chegam?
ANOVA factorial com interacção y ~ A * B a diferença associada a A depende de B?
ANCOVA y ~ A + x ou y ~ A * B + x diferenças entre grupos condicionadas à covariável, assumindo declive comum?
ANCOVA com declives diferentes y ~ A * x a relação com a covariável depende do grupo?
moderação quantitativa y ~ x * z o declive de x depende de z?

A vantagem desta tradução é podermos crescer sem trocar de filosofia a cada capítulo. Acrescentamos parâmetros, libertamos restrições e perguntamos sempre que comparação responde à hipótese.

E Quando a Resposta Não É Bem Representada por uma Recta?

Até aqui trabalhámos com respostas quantitativas que podemos modelar directamente numa escala linear. Mas uma probabilidade não pode descer abaixo de 0 nem subir acima de 1, e uma contagem não pode prever -3 animais.

A lógica dos contrastes, interacções e comparação de modelos não desaparece. O que muda é a forma como representamos a resposta e a sua variabilidade.

Isso merece mais do que quatro secções técnicas no fim de um capítulo já longo. Continuem em Modelos Lineares Generalizados, onde começamos pelos problemas concretos que obrigam a mudar de família ou de ligação e só depois introduzimos a matemática.

E as Medidas Repetidas?

Até aqui tratámos cada linha como uma observação independente. Quando a mesma pessoa ou unidade aparece várias vezes, precisamos de representar essa estrutura de dependência. A lógica de comparação de modelos continua, mas o termo de erro deixa de ser tão simples.

No capítulo de medidas repetidas e modelos mistos vamos construir primeiro uma ANOVA de medidas repetidas a partir de modelos lineares que incluem o participante, perceber exactamente o que essa manobra faz e ver por que se torna incómoda em delineamentos maiores. Só depois passamos às funções que automatizam essa contabilidade e aos modelos mistos.

E os Testes «Não-Paramétricos»?

Até aqui evitámos de propósito organizar os métodos em «paramétricos» e «não-paramétricos». Isto não é falta de rigor. É uma escolha pedagógica: primeiro queremos perceber bem a estrutura comum — uma pergunta, um modelo, parâmetros, restrições e uma comparação — antes de acrescentar mais uma árvore de decisões para decorar.

Isso não quer dizer que os testes baseados em postos vivam noutro universo. Muitos encaixam surpreendentemente bem nesta linguagem quando pensamos em modelos ajustados aos postos (ranks). Mann–Whitney, Wilcoxon, Kruskal–Wallis e Spearman podem ser ligados desta forma a modelos que já conhecemos. A página de Jonas Lindeløv, Common statistical tests are linear models, mostra essas ligações de forma especialmente clara.

Não precisamos de abrir já essa caixa. A equivalência nem sempre é exacta em todos os detalhes do teste, e transformar os dados em postos também muda a quantidade que estamos a analisar.

Por agora, o objectivo é não perder a estrutura conceptual por estarmos a decorar famílias de testes. Depois voltaremos aos pressupostos com bastante mais detalhe.

E aí não vamos ensinar apenas a receita clássica:

«testar pressupostos → se algum teste der significativo → trocar para um teste não-paramétrico».

Essa sequência aparece muitas vezes em manuais e aulas introdutórias porque é simples. Mas é apenas uma simplificação. Um teste formal de normalidade, por exemplo, não decide sozinho se um modelo é adequado, e «não-paramétrico» não quer dizer «sem pressupostos».

Vamos antes perguntar que pressuposto importa para que conclusão, olhar para o delineamento e para diagnósticos relevantes, e estudar o que acontece quando o modelo simplifica demasiado. Isso inclui transformações quando fazem sentido, modelos para outras distribuições, erros-padrão robustos, bootstrap, estimadores robustos, métodos baseados em postos e análises de sensibilidade.

Ou seja: adiar esta discussão não a torna menos rigorosa. Permite-nos tratá-la depois com mais rigor, quando já temos linguagem suficiente para perceber o que cada alternativa muda e porquê.

Também estamos, portanto, a adiar deliberadamente uma pergunta inevitável:

quando é que cada modelo é uma representação suficientemente boa dos dados?

Falaremos disso em Pressupostos: Do Ritual ao Diagnóstico. Primeiro construímos o mapa. Depois aprendemos onde ele pode falhar — e o que fazer quando falha.

6 Síntese

Se retiverem apenas uma ideia deste capítulo, retenham esta: os nomes tradicionais mudam mais depressa do que a matemática subjacente. Uma média, uma diferença entre grupos, um declive, um factor com vários níveis e uma interacção podem ser representados como parâmetros do mesmo modelo linear.

Fazer inferência passa então por saber:

  1. que parâmetros o modelo contém;
  2. que restrição representa a hipótese que queremos avaliar;
  3. que modelo compacto corresponde a essa restrição;
  4. que variabilidade deve servir de termo de erro;
  5. que estimativa e incerteza respondem realmente à pergunta científica.

O tutorial prático de modelos lineares parte daqui, mas faz deliberadamente o contrário deste capítulo: começa pela pergunta de investigação e mostra apenas o código necessário para a responder.

7 Exercícios

  1. Escreva os modelos compacto e aumentado para testar se uma média difere de

    1. Explique por que y ~ 0 sozinho não testa directamente essa hipótese se a resposta não tiver sido recentrada.
  2. Num contraste -1/+1 para dois grupos, derive as duas médias previstas. Explique por que o intercepto fica a meio das médias e por que o coeficiente do contraste é metade da diferença entre elas.

  3. Num modelo y ~ x1 + x2 + x3, explique o que significa testar x2 mantendo os outros dois preditores no modelo. Compare essa pergunta com o teste global do conjunto dos três preditores.

  4. Um factor tem quatro níveis. Quantos coeficientes são necessários para o representar além do intercepto? Por que o teste omnibus não corresponde a um único teste-t?

  5. Compare y ~ A + B com y ~ A * B. Que restrição sobre as diferenças entre níveis está presente no primeiro e é libertada no segundo?

  6. Explique a diferença entre y ~ A * B + x e y ~ A * B * x. Qual deles obriga a covariável a ter o mesmo declive em todas as células?

  7. Num modelo y ~ x * z, qual é o declive de x quando z = 0? E quando z = 2? O que representa o coeficiente de x:z?

  8. Num 2 × 2 com uma interacção teoricamente importante, explique por que neste material preferimos uma tabela Tipo III com contrastes centrados. Depois dê um caso em que uma comparação sequencial Tipo I seria substantivamente mais natural.

8 Recursos