Wieloczynnikowa regresja wielokrotna w R.

Mam 2 zmienne zależne (DV), na które na każdy wynik może mieć wpływ zestaw 7 zmiennych niezależnych (IV). DV są ciągłe, podczas gdy zestaw IV składa się z kombinacji zmiennych ciągłych i binarnie kodowanych. (W kodzie poniżej zmienne ciągłe są pisane dużymi literami, a zmienne binarne małymi literami.)

Celem badania jest odkrycie, jak zmienne IV wpływają na te DV. Zaproponowałem następujący model wielowymiarowej regresji wielokrotnej (MMR):

my.model <- lm(cbind(A, B) ~ c + d + e + f + g + H + I)

Aby zinterpretować wyniki, nazywam dwie instrukcje:

1. summary(manova(my.model))
2. Manova(my.model)

Wyniki obu połączeń są wklejone poniżej i znacznie się różnią. Czy ktoś może wyjaśnić, które oświadczenie należy wybrać, aby właściwie podsumować wyniki MMR i dlaczego? Wszelkie sugestie będą mile widziane.

Dane wyjściowe za pomocą summary(manova(my.model))instrukcji:

> summary(manova(my.model))
           Df   Pillai approx F num Df den Df    Pr(>F)    
c           1 0.105295   5.8255      2     99  0.004057 ** 
d           1 0.085131   4.6061      2     99  0.012225 *  
e           1 0.007886   0.3935      2     99  0.675773    
f           1 0.036121   1.8550      2     99  0.161854    
g           1 0.002103   0.1043      2     99  0.901049    
H           1 0.228766  14.6828      2     99 2.605e-06 ***
I           1 0.011752   0.5887      2     99  0.556999    
Residuals 100                                              
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Dane wyjściowe za pomocą Manova(my.model)instrukcji:

> library(car)
> Manova(my.model)

Type II MANOVA Tests: Pillai test statistic
  Df test stat approx F num Df den Df    Pr(>F)    
c  1  0.030928   1.5798      2     99   0.21117    
d  1  0.079422   4.2706      2     99   0.01663 *  
e  1  0.003067   0.1523      2     99   0.85893    
f  1  0.029812   1.5210      2     99   0.22355    
g  1  0.004331   0.2153      2     99   0.80668    
H  1  0.229303  14.7276      2     99 2.516e-06 ***
I  1  0.011752   0.5887      2     99   0.55700    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 

Krótko mówiąc, to dlatego, że baza-R manova(lm())wykorzystuje kolejnych modeli porównań dla tak zwanego typu I suma kwadratów, natomiast car„s Manova()domyślnie wykorzystuje model do porównań typu II sumy kwadratów.

Zakładam, że znasz podejście porównywania modeli do analizy ANOVA lub analizy regresji. Podejście to definiuje te testy poprzez porównanie modelu ograniczonego (odpowiadającego hipotezie zerowej) z modelem nieograniczonym (odpowiadającym hipotezie alternatywnej). Jeśli nie znasz tego pomysłu, polecam doskonałe „Projektowanie eksperymentów i analizowanie danych” Maxwell & Delaney (2004).

Dla SS typu I, ograniczonym modelem w analizie regresji dla twojego pierwszego predyktora cjest model zerowy, który używa tylko terminu bezwzględnego: lm(Y ~ 1)gdzie Yw twoim przypadku byłaby wielowymiarowa DV zdefiniowana przez cbind(A, B). Nieograniczony model dodaje następnie predyktor c, tj lm(Y ~ c + 1).

W przypadku SS typu II nieograniczonym modelem w analizie regresji dla pierwszego predyktora cjest pełny model, który obejmuje wszystkie predyktory oprócz ich interakcji, tj lm(Y ~ c + d + e + f + g + H + I). Model ograniczony usuwa predyktor cz modelu nieograniczonego, tj lm(Y ~ d + e + f + g + H + I).

Ponieważ obie funkcje opierają się na różnych porównaniach modeli, prowadzą do różnych wyników. Trudno odpowiedzieć na pytanie, które z nich jest lepsze - tak naprawdę zależy od twoich hipotez.

Poniżej założono, że znasz się na tym, jak obliczane są statystyki testowe na wielu odmianach, takie jak Pillai-Bartlett Trace, na podstawie modelu zerowego, pełnego modelu i pary modeli z ograniczeniami nieograniczonymi. Dla zwięzłości rozważam tylko predyktory ci testuję Htylko c.

N <- 100                             # generate some data: number of subjects
c <- rbinom(N, 1, 0.2)               # dichotomous predictor c
H <- rnorm(N, -10, 2)                # metric predictor H
A <- -1.4*c + 0.6*H + rnorm(N, 0, 3) # DV A
B <-  1.4*c - 0.6*H + rnorm(N, 0, 3) # DV B
Y <- cbind(A, B)                     # DV matrix
my.model <- lm(Y ~ c + H)            # the multivariate model
summary(manova(my.model))            # from base-R: SS type I
#           Df  Pillai approx F num Df den Df  Pr(>F)    
# c          1 0.06835   3.5213      2     96 0.03344 *  
# H          1 0.32664  23.2842      2     96 5.7e-09 ***
# Residuals 97                                           

Dla porównania, wynik z car„s Manova()funkcję używając SS typu II.

library(car)                           # for Manova()
Manova(my.model, type="II")
# Type II MANOVA Tests: Pillai test statistic
#   Df test stat approx F num Df den Df  Pr(>F)    
# c  1   0.05904   3.0119      2     96 0.05387 .  
# H  1   0.32664  23.2842      2     96 5.7e-09 ***

Teraz ręcznie sprawdź oba wyniki. Najpierw zbuduj macierz projektową i porównaj ją z macierzą projektową R.$X$

X  <- cbind(1, c, H)
XR <- model.matrix(~ c + H)
all.equal(X, XR, check.attributes=FALSE)
# [1] TRUE

Teraz zdefiniuj rzut ortogonalny dla pełnego modelu ( , używając wszystkich predyktorów). To daje nam macierz .$P_{f} = X (X'X)^{-1} X'$$W = Y' (I-P_{f}) Y$

Pf  <- X %*% solve(t(X) %*% X) %*% t(X)
Id  <- diag(N)
WW  <- t(Y) %*% (Id - Pf) %*% Y

Ograniczona i nieograniczone modele SS typu I oraz ich występy i , prowadzące do matrycy .$P_{rI}$$P_{uI}$$B_{I} = Y' (P_{uI} - P_{PrI}) Y$

XrI <- X[ , 1]
PrI <- XrI %*% solve(t(XrI) %*% XrI) %*% t(XrI)
XuI <- X[ , c(1, 2)]
PuI <- XuI %*% solve(t(XuI) %*% XuI) %*% t(XuI)
Bi  <- t(Y) %*% (PuI - PrI) %*% Y

Ograniczona i nieograniczone modele SS typu II oraz ich występy i , co prowadzi do matrycy .$P_{rI}$$P_{uII}$$B_{II} = Y' (P_{uII} - P_{PrII}) Y$

XrII <- X[ , -2]
PrII <- XrII %*% solve(t(XrII) %*% XrII) %*% t(XrII)
PuII <- Pf
Bii  <- t(Y) %*% (PuII - PrII) %*% Y

Ślad Pillai-Bartlett dla obu typów ss śladu .$(B + W)^{-1} B$

(PBTi  <- sum(diag(solve(Bi  + WW) %*% Bi)))   # SS type I
# [1] 0.0683467

(PBTii <- sum(diag(solve(Bii + WW) %*% Bii)))  # SS type II
# [1] 0.05904288

Zauważ, że obliczenia dla rzutów ortogonalnych naśladują wzór matematyczny, ale liczbowo są złym pomysłem. crossprod()Zamiast tego należy naprawdę używać dekompozycji QR lub SVD .

Cóż, wciąż nie mam wystarczającej liczby punktów, aby skomentować poprzednią odpowiedź i dlatego piszę ją jako osobną odpowiedź, więc proszę wybacz mi. (Jeśli to możliwe, popchnij mnie ponad 50 punktów rep;)

Oto 2 centy: Testowanie błędów typu I, II i III to zasadniczo warianty z powodu niezrównoważenia danych. (Defn niezrównoważony: brak równej liczby obserwacji w każdej z warstw). Jeśli dane są zrównoważone, testy błędów typu I, II i III dają dokładnie takie same wyniki.

Co się dzieje, gdy dane są niezrównoważone?

Rozważ model, który obejmuje dwa czynniki A i B; istnieją zatem dwa główne efekty i interakcja, AB. SS (A, B, AB) oznacza pełny model SS (A, B) oznacza model bez interakcji. SS (B, AB) wskazuje model, który nie uwzględnia efektów czynnika A i tak dalej.

Ta notacja ma teraz sens. Pamiętaj o tym.

SS(AB | A, B) = SS(A, B, AB) - SS(A, B)

SS(A | B, AB) = SS(A, B, AB) - SS(B, AB)

SS(B | A, AB) = SS(A, B, AB) - SS(A, AB)

SS(A | B)     = SS(A, B) - SS(B)

SS(B | A)     = SS(A, B) - SS(A)

Typ I, zwany także „sekwencyjną” sumą kwadratów:

1) SS(A) for factor A.

2) SS(B | A) for factor B.

3) SS(AB | B, A) for interaction AB.

Tak więc najpierw oceniamy główny efekt A, efekt B, biorąc pod uwagę A, a następnie oceniamy interakcję AB, biorąc pod uwagę A i B (w tym przypadku, gdy dane są niezrównoważone, pojawiają się różnice. Gdy najpierw oceniamy główny efekt, a następnie główny inny i następnie interakcja w „sekwencji”)

Typ II:

1) SS(A | B) for factor A.

2) SS(B | A) for factor B.

Test typu II znaczenie głównego efektu A po B i B po A. Dlaczego nie ma SS (AB | B, A)? Zastrzeżenie polega na tym, że metody typu II można użyć tylko wtedy, gdy już przetestowaliśmy interakcję jako nieistotną. Biorąc pod uwagę brak interakcji (SS (AB | B, A) jest nieistotny), test typu II ma lepszą moc niż typ III

Typ III:

1) SS(A | B, AB) for factor A.

2) SS(B | A, AB) for factor B.

Testowaliśmy więc interakcję podczas typu II i interakcja była znacząca. Teraz musimy użyć typu III, ponieważ uwzględnia on termin interakcji.

Jak już powiedział @caracal, kiedy dane są zrównoważone, czynniki są ortogonalne, a typy I, II i III dają takie same wyniki. Mam nadzieję, że to pomoże !

Ujawnienie: Większość nie jest moją pracą. Znalazłem link do tej doskonałej strony i miałem ochotę ją jeszcze bardziej ułożyć.