Mam zestaw danych z 482 obserwacji.
data=Populationfull
Zamierzam przeprowadzić analizę asocjacji genotypu dla 3 SNP. Próbuję zbudować model do mojej analizy i używam aov (y ~ x, data = ...). Dla jednej cechy mam kilka ustalonych efektów i zmiennych towarzyszących, które zawarłem w modelu, na przykład:
Starts <- aov(Starts~Sex+DMRT3+Birthyear+Country+Earnings+Voltsec+Autosec, data=Populationfull) summary(Starts) Df Sum Sq Mean Sq F value Pr(>F) Sex 3 17.90 5.97 42.844 < 2e-16 *** DMRT3 2 1.14 0.57 4.110 0.017 * Birthyear 9 5.59 0.62 4.461 1.26e-05 *** Country 1 11.28 11.28 81.005 < 2e-16 *** Earnings 1 109.01 109.01 782.838 < 2e-16 *** Voltsec 1 12.27 12.27 88.086 < 2e-16 *** Autosec 1 8.97 8.97 64.443 8.27e-15 *** Residuals 463 64.48 0.14 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Odkryłem, że jeśli zmieniłem kolejność zmiennych w modelu, otrzymałem różne wartości p, patrz poniżej.
Starts2 <- aov(Starts~Voltsec+Autosec+Sex+DMRT3+Birthyear+Country+Earnings, data=Populationfull) summary(Starts2) Df Sum Sq Mean Sq F value Pr(>F) Voltsec 1 2.18 2.18 15.627 8.92e-05 *** Autosec 1 100.60 100.60 722.443 < 2e-16 *** Sex 3 10.43 3.48 24.962 5.50e-15 *** DMRT3 2 0.82 0.41 2.957 0.05294 . Birthyear 9 3.25 0.36 2.591 0.00638 ** Country 1 2.25 2.25 16.183 6.72e-05 *** Earnings 1 46.64 46.64 334.903 < 2e-16 *** Residuals 463 64.48 0.14 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Dlaczego otrzymuję różne wartości p w zależności od tego, w jakiej kolejności kodowane są zmienne / czynniki / zmienne towarzyszące / ustalone efekty (?)? Czy istnieje sposób, aby to „poprawić”? Czy to możliwe, że używam niewłaściwego modelu? Nadal jestem całkiem nowy w R, więc jeśli możesz mi w tym pomóc, proszę, zastanów się naprawdę, abym mógł zrozumieć odpowiedź hehe ... Dziękuję, mam nadzieję, że ktoś może mi pomóc to zrozumieć!
Populationfull
aby problem mógł być powtarzalny . Nie dzieje się tak w przypadku przykładu ze stronyaov()
pomocy.summary(aov(yield ~ block + N + P + K, npk)); summary(aov(yield ~ K + P + block + N , npk))
Earnings 1 109.01 109.01 782.838 < 2e-16 ***
twój drugi biegEarnings 1 46.64 46.64 334.903 < 2e-16 ***
. Twoje wyniki nie są takie same. Zacznij od sprawdzenia, czy nie zrobiłeś nic więcej niż zmiana kolejności zmiennych.car
zajrzeć do pakietu - implementuje on ANOVA typu II i typu III, które nie zależą od kolejności zmiennych, natomiastaov
ANOVA typu I.Odpowiedzi:
Problem pochodzi ze sposobu, w jaki
aov()
wykonuje domyślne testowanie istotności. Wykorzystuje analizę ANOVA o nazwie „Typ I”, w której testowanie odbywa się w kolejności, w której zmienne są określone w modelu. Tak więc w pierwszym przykładzie określa, ile wariancji jest wyjaśnionesex
i testuje jej znaczenie, a następnie jaka część pozostałej wariancji jest wyjaśnionaDMRT3
i testuje jej znaczenie pod względem tej pozostałej wariancji i tak dalej. W drugim przykładzieDMRT3
jest oceniany tylko poVoltsec
,Autosec
isex
w tej kolejności, więc pozostało mniej wariancjiDMRT3
do wyjaśnienia.Jeśli dwie zmienne predykcyjne są skorelowane, pierwsza wprowadzona do modelu uzyska pełny „kredyt”, pozostawiając mniejszą wariancję, którą można „wyjaśnić” drugą, co może wydawać się mniej „statystycznie znaczące” niż pierwsza, nawet jeśli jest nie funkcjonalnie. To pytanie i jego odpowiedź wyjaśniają różne typy analiz ANOVA.
Jednym ze sposobów, aby ominąć to, aby wyodrębnić się od zwężeń klasycznej analizy wariancji i używa prostego modelu liniowego, ze
lm()
w R, zamiastaov()
. To skutecznie analizuje wszystkie predyktory równolegle, „korygując” wszystkie predyktory jednocześnie. W takim przypadku dwa skorelowane predyktory mogą mieć duże błędy standardowe szacowanych współczynników regresji, a ich współczynniki mogą różnić się między różnymi próbkami z populacji, ale przynajmniej kolejność wprowadzania zmiennych do specyfikacji modelu nie będzie miała znaczenia.Jeśli twoja zmienna odpowiedzi jest pewnego rodzaju zmienną zliczającą, jak
Starts
sugeruje jej nazwa , prawdopodobnie i tak nie powinieneś używać ANOVA, ponieważ jest mało prawdopodobne, aby reszty były normalnie rozłożone, jak wymaga interpretacja wartości p . Zmienne zliczania są lepiej obsługiwane przy uogólnionych modelach liniowych (np.glm()
W R), które można traktować jako uogólnienielm()
dla innych rodzajów struktur błędów resztkowych.źródło