Dlaczego ANOVA jest równoważna regresji liniowej?

Czytałem, że ANOVA i regresja liniowa są tym samym. Jak to możliwe, biorąc pod uwagę, że wynikiem ANOVA jest pewna wartość i pewna wartość p, na podstawie których wnioskujesz, czy próbka oznacza, że ​​różne próbki są takie same lub różne.$F$$p$

Ale zakładając, że średnie nie są równe (odrzuć hipotezę zerową), ANOVA nie mówi ci nic o współczynnikach modelu liniowego. Więc w jaki sposób regresja liniowa jest taka sama jak ANOVA?

ANOVA i regresja liniowa są równoważne, gdy dwa modele testują te same hipotezy i stosują identyczne kodowanie. Modele różnią się podstawowym celem: ANOVA zajmuje się przede wszystkim prezentowaniem w danych różnic między średnimi kategorii, podczas gdy regresja liniowa dotyczy głównie oszacowania średniej odpowiedzi próbki i związanego z nią .$\sigma^2$

Nieco aforyzmicznie można opisać ANOVA jako regresję ze zmiennymi fikcyjnymi. Łatwo zauważyć, że dzieje się tak w prostej regresji ze zmiennymi kategorialnymi. Zmienna kategoryczna zostanie zakodowana jako matryca wskaźnikowa (matryca 0/1zależna od tego, czy pacjent jest częścią danej grupy, czy nie), a następnie zastosowana bezpośrednio do rozwiązania układu liniowego opisanego regresją liniową. Zobaczmy przykład z 5 grupami. Dla celów argumentu przyjmuję, że średnia group1równa 1, średnia group2równa 2, ... i średnia group5równa 5. (Używam MATLAB, ale dokładnie to samo jest równoważne w 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

Jak widać w tym scenariuszu, wyniki były dokładnie takie same. Drobna różnica liczbowa wynika z faktu, że projekt nie jest idealnie wyważony, a także z procedury szacowania podkładu; ANOVA kumuluje błędy numeryczne nieco bardziej agresywnie. W tym względzie pasujemy do punktu przecięcia LRbetas(1); moglibyśmy dopasować model bez przechwytywania, ale nie byłby to „standardowa” regresja liniowa. (W takim przypadku wyniki byłyby jeszcze bliższe ANOVA.)

$F$

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

Wynika to z tego, że procedury sprawdzają tę samą hipotezę, ale używając różnych sformułowań: ANOVA sprawdzi jakościowo, czy „ stosunek jest wystarczająco wysoki, aby zasugerować, że żadne grupowanie jest niewiarygodne ”, podczas gdy regresja liniowa sprawdzi jakościowo, czy „ stosunek jest wystarczająco wysoki, aby sugerować jedynie przechwycenie model jest prawdopodobnie nieodpowiedni ”.
(Jest to dość dowolna interpretacja „ możliwości zobaczenia wartości równej lub większej niż ta obserwowana w ramach hipotezy zerowej ” i nie jest to definicja podręcznika).

$\beta$, itd. Oczywiście, gdy zaczyna się dodawanie wielu zmiennych towarzyszących w jego modelu regresji, prosta jednokierunkowa ANOVA nie ma bezpośredniej równoważności. W takim przypadku powiększa się informacje wykorzystane do obliczenia średniej odpowiedzi regresji liniowej o informacje, które nie są bezpośrednio dostępne dla jednoczynnikowej ANOVA. Wierzę, że można ponownie wyrazić rzeczy w kategoriach ANOVA, ale jest to w większości ćwiczenie akademickie.

Interesującym artykułem na ten temat jest artykuł Gelmana z 2005 r. Zatytułowany: Analiza wariancji - dlaczego jest ważniejszy niż kiedykolwiek . Podniesiono kilka ważnych kwestii; Nie w pełni popieram ten artykuł (myślę, że osobiście dużo bardziej zgadzam się z poglądem McCullacha), ale może to być konstruktywna lektura.

Na koniec: fabuła gęstnieje, gdy masz mieszane modele efektów . Istnieją różne koncepcje dotyczące tego, co można uznać za uciążliwe lub aktualne informacje dotyczące grupowania danych. Te kwestie są poza zakresem tego pytania, ale myślę, że warto je skinąć.

Pozwól, że dodam trochę koloru do idei, że OLS z kategorycznymi ( sztucznie kodowanymi ) regresorami jest równoważny czynnikom z ANOVA. W obu przypadkach istnieją poziomy (lub grupy w przypadku ANOVA).

W regresji OLS najczęściej regresory mają również ciągłe zmienne. Logicznie modyfikują one związek w modelu dopasowania między zmiennymi kategorialnymi a zmienną zależną (DC). Ale nie do tego stopnia, że ​​równolegle nie można go rozpoznać.

Na podstawie mtcarszestawu danych możemy najpierw wizualizować model lm(mpg ~ wt + as.factor(cyl), data = mtcars)jako nachylenie określone przez zmienną ciągłą wt(wagę), a różne punkty przecięcia rzutują efekt zmiennej jakościowej cylinder(cztery, sześć lub osiem cylindrów). Jest to ostatnia część, która tworzy równoległość z jednokierunkową ANOVA.

Zobaczmy to graficznie na wykresie po prawej stronie (trzy wykresy po lewej stronie są uwzględnione w celu porównania bocznego z modelem ANOVA omówionym natychmiast potem):

Każdy silnik cylindrowy jest oznaczony kolorami, a odległość między dopasowanymi liniami z różnymi punktami przechwytywania a chmurą danych jest równoważna wariacji wewnątrzgrupowej w analizie ANOVA. Zauważ, że przecięcia w modelu OLS ze zmienną ciągłą ( weight) nie są matematycznie takie same, jak wartość różnych średnich wewnątrzgrupowych w ANOVA, ze względu na wpływ weighti różne macierze modelu (patrz poniżej): średnia mpgdla samochodów 4-cylinder, jest na przykład mean(mtcars$mpg[mtcars$cyl==4]) #[1] 26.66364, podczas gdy OLS „bazowy” osią (odzwierciedlające umownie cyl==4(od najniższego do najwyższego cyframi porządkuje R)) jest znacząco odmienne: summary(fit)$coef[1] #[1] 33.99079. Nachylenie linii jest współczynnikiem zmiennej ciągłej weight.

Jeśli spróbujesz stłumić efekt weightmentalnego prostowania tych linii i przywrócenia ich do linii poziomej, skończysz na wykresie ANOVA modelu aov(mtcars$mpg ~ as.factor(mtcars$cyl))na trzech wykresach podrzędnych po lewej stronie. weightREGRESSOR jest obecnie na zewnątrz, ale związek z punktami do różnych przechwytuje jest grubsza zachowane - jesteśmy po prostu obraca się w kierunku przeciwnym do ruchu wskazówek zegara i rozkładanie wcześniej nakładających się wykresy dla każdego innego poziomu (ponownie, tylko jako urządzenie wizualnej „zobaczyć” związek; nie jako matematyczna równość, ponieważ porównujemy dwa różne modele!).

cylinder$\small 20$$x$

I poprzez sumę tych pionowych segmentów możemy ręcznie obliczyć resztki:

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)

Wynik: SumSq = 301.2626i TSS - SumSq = 824.7846. Porównać do:

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

Dokładnie taki sam wynik jak testowanie za pomocą ANOVA modelu liniowego z jedynie kategorialnym cylinderregresorem:

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 

Widzimy zatem, że reszty - część całkowitej wariancji nie wyjaśniona przez model - a także wariancja są takie same, niezależnie od tego, czy wywołasz OLS typu lm(DV ~ factors), czy ANOVA ( aov(DV ~ factors)): kiedy usuwamy model zmiennych ciągłych kończy się na tym samym systemie. Podobnie, gdy oceniamy modele globalnie lub jako zbiorczą ANOVA (nie poziom po poziomie), naturalnie otrzymujemy tę samą wartość p F-statistic: 39.7 on 2 and 29 DF, p-value: 4.979e-09.

Nie oznacza to, że testowanie poszczególnych poziomów da identyczne wartości p. W przypadku OLS możemy wywołać summary(fit)i uzyskać:

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$p adjusted

  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

Ostatecznie nic nie jest bardziej uspokajające niż zerknięcie na silnik pod maską, który jest niczym innym jak matrycami modelu i rzutami w przestrzeni kolumny. Są to w rzeczywistości dość proste w przypadku 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}$

Z drugiej strony matryca modelu dla regresji OLS to:

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

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 i usεr11852 mieli bardzo dobrą odpowiedź. Odpowiem na twoje pytanie dotyczące kodowania perspektywy R.

ANOVA nie mówi nic o współczynnikach modelu liniowego. Więc w jaki sposób regresja liniowa jest taka sama jak ANOVA?

W rzeczywistości możemy aovdziałać w Rmoże być używany tak samo jak lm. Oto kilka przykładów.

> 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

Jak widać, nie tylko możemy uzyskać współczynnik z modelu ANOVA, ale także możemy go wykorzystać do przewidywania, podobnie jak model liniowy.

Jeśli sprawdzimy plik pomocy pod kątem aovfunkcji, wyświetli się komunikat

Ten dostarcza opakowanie do lm do montażu modeli liniowych na symetryczne lub niesymetryczne wzorów doświadczalnych. Główną różnicą w stosunku do lm jest sposób, w jaki druk, podsumowanie itd. Radzą sobie z dopasowaniem: wyraża się to w tradycyjnym języku analizy wariancji, a nie w modelach liniowych.

Jeśli weźmiemy wszystkie wpisy danych i ułożymy je w jedną kolumnę Y, a pozostałe kolumny będą zmiennymi wskaźnikowymi 1 {i dane są elementem j-tej kolumny w oryginalnym układzie anova}, wówczas przyjmując prostą regresję liniową Y na w dowolnej innej kolumnie (powiedzmy kolumnie B), powinieneś uzyskać taką samą statystykę testową DF, SS, MS i F, jak w twoim problemie ANOVA.

Zatem ANOVA można „traktować jak” regresję liniową, zapisując dane ze zmiennymi binarnymi. Zauważ też, że współczynnik regresji dla, powiedzmy, regresji Y na B powinien być taki sam jak średnia. kolumny B, obliczonej na podstawie oryginalnych danych.