Por que ANOVA é equivalente a regressão linear?

Eu li que ANOVA e regressão linear são a mesma coisa. Como pode ser isso, considerando que a saída da ANOVA é um valor e algum valor p com base no qual você conclui se a média da amostra nas diferentes amostras é igual ou diferente.$F$$p$

Mas, assumindo que as médias não são iguais (rejeite a hipótese nula), a ANOVA não diz nada sobre os coeficientes do modelo linear. Então, como a regressão linear é igual à ANOVA?

ANOVA e regressão linear são equivalentes quando os dois modelos testam contra as mesmas hipóteses e usam uma codificação idêntica. Os modelos diferem em seu objetivo básico: a ANOVA se preocupa principalmente em apresentar diferenças entre as médias das categorias nos dados, enquanto a regressão linear se preocupa principalmente em estimar uma resposta média da amostra e um associado .$\sigma^2$

Um tanto aforicamente, pode-se descrever a ANOVA como uma regressão com variáveis ​​fictícias. Podemos ver facilmente que esse é o caso da regressão simples com variáveis ​​categóricas. Uma variável categórica será codificada como uma matriz indicadora (uma matriz que 0/1depende se um sujeito faz parte de um determinado grupo ou não) e, em seguida, usada diretamente para a solução do sistema linear descrita por uma regressão linear. Vamos ver um exemplo com 5 grupos. Por uma questão de argumento, assumirei que a média de group1igual a 1, a média de group2igual a 2, ... e a média de group5igual a 5. (eu uso MATLAB, mas exatamente a mesma coisa é equivalente em R.)

rng(123);               % Fix the seed
X = randi(5,100,1);     % Generate 100 random integer U[1,5]
Y = X + randn(100,1);   % Generate my response sample
Xcat = categorical(X);  % Treat the integers are categories

% One-way ANOVA
[anovaPval,anovatab,stats] = anova1(Y,Xcat);
% Linear regression
fitObj = fitlm(Xcat,Y);

% Get the group means from the ANOVA
ANOVAgroupMeans = stats.means
% ANOVAgroupMeans =
% 1.0953    1.8421    2.7350    4.2321    5.0517

% Get the beta coefficients from the linear regression
LRbetas = [fitObj.Coefficients.Estimate'] 
% LRbetas =
% 1.0953    0.7468    1.6398    3.1368    3.9565

% Rescale the betas according the intercept
scaledLRbetas = [LRbetas(1) LRbetas(1)+LRbetas(2:5)]
% scaledLRbetas =
% 1.0953    1.8421    2.7350    4.2321    5.0517

% Check if the two results are numerically equivalent
abs(max( scaledLRbetas - ANOVAgroupMeans)) 
% ans =
% 2.6645e-15

Como pode ser visto neste cenário, os resultados foram exatamente os mesmos. A diferença numérica diminuta se deve ao fato de o projeto não estar perfeitamente equilibrado, bem como ao procedimento de estimativa subjacente; a ANOVA acumula erros numéricos um pouco mais agressivamente. A esse respeito, ajustamos um intercepto LRbetas(1); poderíamos ajustar um modelo sem interceptação, mas isso não seria uma regressão linear "padrão". (Os resultados estariam ainda mais próximos da ANOVA nesse caso.)

$F$

abs( fitObj.anova.F(1) - anovatab{2,5} )
% ans =
% 2.9132e-13 

Isso ocorre porque os procedimentos testam a mesma hipótese, mas com palavras diferentes: a ANOVA verificará qualitativamente se " a proporção é alta o suficiente para sugerir que nenhum agrupamento é implausível " enquanto a regressão linear verificará qualitativamente se " a proporção é alta o suficiente para sugerir apenas uma interceptação modelo é possivelmente inadequado ".
(Essa é uma interpretação um tanto livre da " possibilidade de ver um valor igual ou superior ao observado sob a hipótese nula " e não pretende ser uma definição de livro-texto.)

$\beta$, etc. Claramente quando alguém começa a adicionar múltiplas covariáveis ​​em seu modelo de regressão, uma ANOVA unidirecional simples não tem uma equivalência direta. Nesse caso, aumentamos as informações usadas para calcular a resposta média da regressão linear com informações que não estão diretamente disponíveis para uma ANOVA de sentido único. Acredito que se pode re-expressar as coisas em termos de ANOVA mais uma vez, mas é principalmente um exercício acadêmico.

Um artigo interessante sobre o assunto é o artigo de Gelman, de 2005, intitulado: Analysis of Variance - Por que é mais importante do que nunca . Alguns pontos importantes levantados; Não sou totalmente favorável ao artigo (acho que pessoalmente me alino muito mais à visão de McCullach), mas pode ser uma leitura construtiva.

Como nota final: o enredo engrossa quando você tem modelos de efeitos mistos . Lá você tem conceitos diferentes sobre o que pode ser considerado um incômodo ou informações reais sobre o agrupamento de seus dados. Essas questões estão fora do escopo desta questão, mas acho que merecem um aceno.

Deixe-me colocar um pouco de cor na idéia de que o OLS com regressores categóricos ( codificados por modelos) é equivalente aos fatores da ANOVA. Nos dois casos, existem níveis (ou grupos no caso da ANOVA).

Na regressão OLS, é mais comum também ter variáveis ​​contínuas nos regressores. Eles modificam logicamente o relacionamento no modelo de ajuste entre as variáveis ​​categóricas e a variável dependente (CD). Mas não ao ponto de tornar o paralelo irreconhecível.

Com base no mtcarsconjunto de dados, podemos primeiro visualizar o modelo lm(mpg ~ wt + as.factor(cyl), data = mtcars)como a inclinação determinada pela variável contínua wt(peso) e as diferentes intercepções projetando o efeito da variável categórica cylinder(quatro, seis ou oito cilindros). É esta última parte que forma um paralelo com uma ANOVA unidirecional.

Vamos vê-lo graficamente na subparcela à direita (as três subparcelas à esquerda estão incluídas para comparação lado a lado com o modelo ANOVA discutido imediatamente depois):

Cada mecanismo de cilindro é codificado por cores, e a distância entre as linhas ajustadas com diferentes interceptações e a nuvem de dados é equivalente à variação dentro do grupo em uma ANOVA. Observe que as interceptações no modelo OLS com uma variável contínua ( weight) não são matematicamente iguais ao valor das diferentes médias dentro do grupo na ANOVA, devido ao efeito weighte às diferentes matrizes do modelo (veja abaixo): a média mpgpara carros de 4 cilindros, por exemplo, é mean(mtcars$mpg[mtcars$cyl==4]) #[1] 26.66364, ao passo que o OLS intercepção "linha de base" (reflectindo por convenção cyl==4(menor para o maior numerais ordenação em R)) é marcadamente diferente: summary(fit)$coef[1] #[1] 33.99079. A inclinação das linhas é o coeficiente para a variável contínua weight.

Se você tentar suprimir o efeito de weightendireitar mentalmente essas linhas e retorná-las à linha horizontal, você terminará com o gráfico ANOVA do modelo aov(mtcars$mpg ~ as.factor(mtcars$cyl))nas três sub plotagens à esquerda. O weightregressor agora está fora, mas a relação dos pontos para as diferentes interceptações é praticamente preservada - estamos simplesmente girando no sentido anti-horário e espalhando as plotagens anteriormente sobrepostas para cada nível diferente (novamente, apenas como um dispositivo visual para "ver" a conexão; não como uma igualdade matemática, pois estamos comparando dois modelos diferentes!).

cylinder$\small 20$$x$

E é através da soma desses segmentos verticais que podemos calcular manualmente os resíduos:

mu_mpg <- mean(mtcars$mpg)                      # Mean mpg in dataset
TSS <- sum((mtcars$mpg - mu_mpg)^2)             # Total sum of squares
SumSq=sum((mtcars[mtcars$cyl==4,"mpg"]-mean(mtcars[mtcars$cyl=="4","mpg"]))^2)+
sum((mtcars[mtcars$cyl==6,"mpg"] - mean(mtcars[mtcars$cyl=="6","mpg"]))^2)+
sum((mtcars[mtcars$cyl==8,"mpg"] - mean(mtcars[mtcars$cyl=="8","mpg"]))^2)

O resultado: SumSq = 301.2626e TSS - SumSq = 824.7846. Comparado a:

Call:
   aov(formula = mtcars$mpg ~ as.factor(mtcars$cyl))

Terms:
                as.factor(mtcars$cyl) Residuals
Sum of Squares               824.7846  301.2626
Deg. of Freedom                     2        29

Exatamente o mesmo resultado que testar com uma ANOVA o modelo linear com apenas o categórico cylindercomo regressor:

fit <- lm(mpg ~ as.factor(cyl), data = mtcars)
summary(fit)
anova(fit)

Analysis of Variance Table

Response: mpg
               Df Sum Sq Mean Sq F value    Pr(>F)    
as.factor(cyl)  2 824.78  412.39  39.697 4.979e-09 ***
Residuals      29 301.26   10.39 

O que vemos, então, é que os resíduos - a parte da variação total não explicada pelo modelo - e a variação são os mesmos, independentemente de você chamar um OLS do tipo lm(DV ~ factors) ou uma ANOVA ( aov(DV ~ factors)): quando removemos o valor modelo de variáveis ​​contínuas, acabamos com um sistema idêntico. Da mesma forma, quando avaliamos os modelos globalmente ou como uma ANOVA abrangente (não nível por nível), obtemos naturalmente o mesmo valor-p F-statistic: 39.7 on 2 and 29 DF, p-value: 4.979e-09.

Isso não significa que o teste de níveis individuais produzirá valores p idênticos. No caso do OLS, podemos invocarsummary(fit) e obter:

lm(formula = mpg ~ as.factor(cyl), data = mtcars)

                Estimate Std. Error t value                           Pr(>|t|)    
(Intercept)      26.6636     0.9718  27.437                           < 2e-16 ***
as.factor(cyl)6  -6.9208     1.5583  -4.441                           0.000119 ***
as.factor(cyl)8 -11.5636     1.2986  -8.905                           8.57e-10 ***

$p$ avaliações de valor , precisamos executar um teste de diferença significativa honesta de Tukey, que tentará reduzir a possibilidade de um erro do tipo I como resultado da realização de múltiplas comparações pareadas (daí, " p adjusted"), resultando em uma saída completamente diferente:

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = mtcars$mpg ~ as.factor(mtcars$cyl))

$`as.factor(mtcars$cyl)`
          diff        lwr        upr                                      p adj
6-4  -6.920779 -10.769350 -3.0722086                                    0.0003424
8-4 -11.563636 -14.770779 -8.3564942                                    0.0000000
8-6  -4.642857  -8.327583 -0.9581313                                    0.0112287

Por fim, nada é mais tranquilizador do que dar uma olhada no motor sob o capô, que não é outro senão as matrizes do modelo e as projeções no espaço da coluna. Estes são realmente bastante simples no caso de uma ANOVA:

$$\small\begin{bmatrix} y_1 \\ y_2 \\ y_3 \\ \vdots \\\vdots\\\vdots\\.\\y_n \end{bmatrix} = \begin{bmatrix} \color{magenta} 1 & 0 & 0 \\ \color{magenta}1 & 0 & 0 \\ \vdots & \vdots & \vdots \\ \color{magenta} 0 & 1 & 0 \\ \color{magenta}0 & 1 & 0 \\ \vdots & \vdots & \vdots \\ .&.&.\\\color{magenta} 0 & 0 & 1 \\ \color{magenta}0 & 0 & 1 \\ \end{bmatrix} \begin{bmatrix} \mu_1\\ \mu_2\\ \mu_3 \end{bmatrix} +\begin{bmatrix} \varepsilon_1 \\ \varepsilon_2\\ \varepsilon_3\\ \vdots\\ \vdots\\ \vdots\\ .\\ \varepsilon_n \end{bmatrix}\tag 1$$

cyl 4cyl 6cyl 8$\small y_{ij} = \mu_i + \epsilon_{ij}$$\mu_i$$j$$i$$y_{ij}$

Por outro lado, a matriz do modelo para uma regressão OLS é:

$$\small\begin{bmatrix}y_1 \\ y_2 \\ y_3 \\ y_4 \\ \vdots \\ y_n \end{bmatrix} = \begin{bmatrix} 1 & x_{12} & x_{13}\\ 1 & x_{22} & x_{23} \\ 1 & x_{32} & x_{33} \\ 1 & x_{42} & x_{43} \\ \vdots & \vdots & \vdots \\1 & x_{n2} & x_{n3} \end{bmatrix} \begin{bmatrix} \beta_0 \\ \beta_1 \\ \beta_2 \end{bmatrix} + \begin{bmatrix} \varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_3 \\ \varepsilon_4 \\ \vdots \\ \varepsilon_n \end{bmatrix}$$

$\small y_i = \beta_0 + \beta_1\, x_{i1} + \beta_2\, x_{i2} + \epsilon_i$$\beta_0$$\beta_1$$\beta_2$weightdisplacement

lm(mpg ~ wt + as.factor(cyl), data = mtcars)weight$\beta_0$weight$\beta_1$$\color{brown}1$cyl 4cyl 4$\color{brown}1$$\color{magenta}1$$(1),$cyl 6cyl 8

$$\small\begin{bmatrix}y_1 \\ y_2 \\ y_3 \\ y_4\\ y_5 \\ \vdots \\ y_n\end{bmatrix} = \begin{bmatrix} \color{brown}1 & x_1 \\ \color{brown}1 & x_2 \\\color{brown} 1 & x_3 \\ \color{brown}1 & x_4 \\ \color{brown}1 & x_5 \\ \vdots & \vdots \\\color{brown}1 & x_n \end{bmatrix} \begin{bmatrix} \beta_0 \\ \beta_1 \end{bmatrix}+ \begin{bmatrix}\color{red}1&0\\\color{red}1&0\\\color{red}1&0\\0&\color{blue}1\\0&\color{blue}1\\ \vdots & \vdots\\0&\color{blue}1\end{bmatrix} \begin{bmatrix} \tilde\mu_2 \\ \tilde\mu_3 \end{bmatrix} + \begin{bmatrix} \varepsilon_1 \\ \varepsilon_2 \\ \varepsilon_3 \\ \varepsilon_4 \\ \varepsilon_5\\ \vdots \\ \varepsilon_n \end{bmatrix}$$

$\color{red}1$$\tilde\mu_2.$$\tilde\cdot$ indica que, como no caso da interceptação de "linha de base" no modelo OLS, não sendo idêntica à média do grupo dos carros de 4 cilindros, mas refletindo isso, as diferenças entre os níveis no modelo OLS não são matematicamente as diferenças de grupos em médias:

fit <- lm(mpg ~ wt + as.factor(cyl), data = mtcars)
summary(fit)$coef[3] #[1] -4.255582 (difference between intercepts cyl==4 and cyl==6 in OLS)
fit <- lm(mpg ~ as.factor(cyl), data = mtcars)
summary(fit)$coef[2] #[1] -6.920779 (difference between group mean cyl==4 and cyl==6)

$\color{blue}1$$\tilde\mu_3$$\small y_i = \beta_0 + \beta_1\, x_i + \tilde\mu_i + \epsilon_i$

Antoni Parellada e usεr11852 tiveram uma resposta muito boa. Vou abordar sua pergunta para codificar a perspectiva com R.

A ANOVA não diz nada sobre os coeficientes do modelo linear. Então, como a regressão linear é igual à ANOVA?

De fato, podemos aovfuncionar Rpode ser usado da mesma forma que lm. Aqui estão alguns exemplos.

> lm_fit=lm(mpg~as.factor(cyl),mtcars)

> aov_fit=aov(mpg~as.factor(cyl),mtcars)

> coef(lm_fit)
    (Intercept) as.factor(cyl)6 as.factor(cyl)8 
      26.663636       -6.920779      -11.563636 

> coef(aov_fit)
    (Intercept) as.factor(cyl)6 as.factor(cyl)8 
      26.663636       -6.920779      -11.563636 

> all(predict(lm_fit,mtcars)==predict(aov_fit,mtcars))
[1] TRUE

Como você pode ver, não apenas podemos obter coeficientes do modelo ANOVA, mas também podemos usá-lo para previsão, assim como o modelo linear.

Se verificarmos o arquivo de ajuda quanto à aovfunção, ele diz

Isso fornece um invólucro a lm para ajustar modelos lineares a projetos experimentais balanceados ou desbalanceados. A principal diferença de lm está na maneira como a impressão, o resumo e assim por diante lidam com o ajuste: isso é expresso na linguagem tradicional da análise de variância, e não nos modelos lineares.

Se pegarmos todas as entradas de dados e as organizarmos em uma única coluna Y, com o restante das colunas sendo variáveis ​​indicadoras 1 {com dados é elemento da j-ésima coluna no arranjo original da anova}, tomando uma regressão linear simples de Y Em qualquer uma das outras colunas (por exemplo, coluna B), você deve obter a mesma estatística de teste DF, SS, MS e F do seu problema ANOVA.

Assim, a ANOVA pode ser 'tratada como' Regressão Linear escrevendo os dados com variáveis ​​binárias. Observe também que o coeficiente de regressão para, digamos, uma regressão de Y em B deve ser o mesmo que a média. da coluna B, calculada com os dados originais.