1
Analiza wariancji klasyfikacja prosta
Dane o przeżywalności chrząszczy hodowanych hodowlanych na czterech różnych pożywkach. Każda pożywka stanowi grupę po 5 pomiarów. Interesuje nas odpowiedz na pytanie, czy skład pożywki ma wpływ przeżywalności chrząszczy. Jest model to I analizy wariancji, ponieważ czynnik różnicujący grupy (skład pożywki) jest czynnikiem powtarzalnym znajdującym się pod kontrolą eksperymentatora, a nie zmienną losową. Sposób przeprowadzenia obliczeń jest taki sam dla modelu I i II. Polega on na oszacowaniu wariancji między grupami i wewnątrz grup.
X
MP0 58 60 51 66 62 59,4
MP5 65 70 64 75 68 68,4
MP2 69 62 70 63 65 65,8
MPR 63 68 68 60 66 65,0
Dla tych danych będziemy testować hipotezę zerową zakładającą, że zróżnicowanie przeżywalności między grupami nie jest większe niż wewnątrz grup, czyli miedzy różnymi pożywkami nie ma różnicy w przeżywalności chrząszczy.
Zgodnie z konwencją wskaźnikiem ij oznaczymy i-ty pomiar w j-tej grupie. W ten sposób drugi pomiar w trzeciej grupie (MP2) oznaczamy symbolem x2,3=62.
W naszych danych są a=4 grupy, w każdej grupie j mamy Nj=5 pomiarów, zatem we wszystkich grupach jest łącznie N=20 pomiarów.
Aby otrzymać ogólną (całkowitą) sumę kwadratów odchyleń posłużymy się wzorem:
a
j N
i
a
j N
i ij ij
j
j
N X X
2
2 (1)
Pierwszy składnik wzoru (1) otrzymujemy podnosząc do kwadratu każdy z pomiarów, a następnie sumując wszystkie wyniki
a
j N
i ij
j
X2 (58)2 (60)2 ... (66)2 84107
I II III IV V I2 II2 III2 IV2 V5
MP0 58 60 51 66 62 3364 3600 2601 4356 3844
MP5 65 70 64 75 68 4225 4900 4096 5625 4624
MP2 69 62 70 63 65 4761 3844 4900 3969 4225
MPR 63 68 68 60 66 3969 4624 4624 3600 4356
1293 84107
Drugi składnik wzoru (1), czyli wyraz poprawkowy obliczamy sumując wszystkie pomiary, podnosząc je do kwadratu, a następnie dzieląc przez liczbę wszystkich pomiarów N.
45 , 83592 20
/ 1671849 20
/ ) 1293 ( 20 / ) 66 ...
60 58 (
/ 2 2
2
N X
a
j N
i ij
j
Zatem całkowita (ogólna) suma kwadratów odchyleń wynosi 84107-83592,45=514,55
2
Międzygrupowa suma kwadratów odchyleń jest liczona według wzoru:
N X N
X
a
j N
i a ij
j j
N
i ij
j j
2 2
(2)
Drugi składnik tego wzoru jest identyczny jak w wzorze (1) do obliczenia całkowitej sumy kwadratów (jest to wyraz poprawkowy), został wcześniej wyliczony i wynosi 83592,45.
Pierwszy wyraz wzoru (2) otrzymujemy sumując wszystkie pomiary dla każdej grupy osobno. Następnie każdą z tych sum podnosimy do kwadratu i dzielimy przez liczbę pomiarów, na podstawie których została obliczona. Jeżeli liczba pomiarów w każdej grupie jest różna, to zgodnie ze wzorem (2) dla każdej grupy dzielimy przez liczbę pomiarów Nj w tej grupie
i i i i
i i
i
i X X X
X1 297, 2 342, 3 329, 4 325
8 , 83807 5
/ ) 325 ( 5 / ) 329 ( 5 / ) 342 ( 5 / ) 297 (
/ 2 2 2 2
a 3
j
j N
i
ij N
X
j
I II III IV V 2 2/Nj
MP0 58 60 51 66 62 297 88209 17641,8
MP5 65 70 64 75 68 342 116964 23392,8
MP2 69 62 70 63 65 329 108241 21648,2
MPR 63 68 68 60 66 325 105625 21125,0
1293 419039 83807,8
Zgodnie ze wzorem (2) międzygrupowa suma kwadratów odchyleń wynosi:
83807,8-83592,45=215,35
Suma kwadratów odchyleń (SK) równa się:
Ogólna SK = międzygrupowa SK + wewnątrzgrupowa SK
Wewnątrzgrupową suma kwadratów odchyleń (składnik błędu) oblicz się następująco:
Wewnątrzgrupowa SK = Ogólna SK - międzygrupowa SK Czyli dla przykładu: 514,55 – 215,35=299,20
Liczba stopni swobody dla całkowitej SK wynosi: df=N-1 = 20-1=19, dla międzygrupowej SK: df=a-1 = 4-1=3 a dla wewnątrzgrupowej SK:
df =
a
i
Nj 1) (5 1) (5 1) (5 1) (5 1) 16 (
Wzór ten pozwala obliczyć wewnątrzgrupowa liczbę stopni swobody nawet, gdy liczba pomiarów w poszczególnych grupach jest różna.
Całkowita df = międzygrupowa df + wewnątrzgrupowa df
Następnie uzyskane sumy kwadratów (SK) i stopnie swobody (df) zbieramy w tabelce
3
Oszacowanie wariancji między grupami i wewnątrz grup (niewyjaśnionej) otrzymujemy dzieląc sumy kwadratów odchyleń przez odpowiadające im stopnie swobody, odpowiednio 215,53/3=71,783 i 299,20/16=18,700
Tabela 1. Sumy kwadratów (SK), stopnie swobody (df), oszacowanie wariancji i stosunek F oszacowania wariancji międzygrupowej (miedzy pożywkami) do oszacowania wariancji wewnątrzgrupowej (błędu) dla przykładu.
Źródła zmienności SK df Oszacowanie
wariancji F
Całkowita (ogólna) Między pożywkami Błąd
514,55 215,35 299,20
19 3 16
71,783 18,700
3,839
Oszacowanie wariancji ogólnej w analizie wariancji można pominąć, ponieważ interesuje nas istnienie zmienności między grupami (pożywkami).
Aby ustalić, czy zróżnicowanie między grupami jest statystycznie istotne musimy obliczyć stosunek:
F = wariancja między grupami/wariancja w grupach Co w przykładzie daje F=71,783/18,700=3,839
W przypadku, gdy oszacowanie wariancji między grupami jest mniejsze niż w grupach, czyli gdy F 1, to możemy uznać, ze zebrane dane nie pozwalają na stwierdzenie zróżnicowania między grupami. W przeciwnym razie należy dokonać porównania z tabelą G w której podano wartości krytyczne rozkładu F. Tabela jest tak skonstruowana, że w główce tablicy podana jest liczba stopni swobody dla większego oszacowania wariancji (czyli między grupami), a w pierwszej kolumnie – dla mniejszego (czyli wewnątrz grup). Sprawdzamy najpierw wartości krytyczne dla poziomu istotności 0,05. Dla df=3 (większa wariancja) i df=16 (mniejsza wariancja) otrzymujemy krytyczny stosunek F0,05; 3; 16=3,24. Ponieważ otrzymany z obliczeń stosunek F=3,839 jest większy od krytycznego, to odrzucamy hipotezę zerową zakładającą, że poszczególne grupy nie różnią się między sobą. Skład pożywki ma zatem wpływ na przeżywalność chrząszczy. Odrzucając te hipotezę, akceptujemy prawdopodobieństwo popełnienia błędu I rodzaju P<0,05. Gdyby odczytany z tablicy G krytyczny stosunek F0,05; 3; 16=3,24 był większy od stosunku obliczonego, wówczas należałoby przyjąć hipotezę zerowa.
Po odrzuceniu hipotezy na poziomie istotności 0,05 należałoby sprawdzić, czy nie da się jej odrzucić z mniejszym błędem I rodzaju. Dlatego też sprawdzamy wartość krytyczną stosunku F dla poziomu istotności 0,025. Wynosi ona F0,025; 3; 16=4,08, czyli jest wyższa niż otrzymany stosunek F. Hipotezy zerowej przy tym poziomie istotności nie można odrzucić.
Wynika z tego, że hipotezę zerową należy odrzucić na poziomie istotności 0,05 (?0,025), zaś prawdopodobieństwa popełnienia błędu I rodzaju przy jej odrzuceniu zawiera się w przedziale 0,025<P<0,05.
Za pomocą analizy wariancji klasyfikacja prosta można testować hipotezę zerowa o braku różnic między dwoma tylko grupami. Taki test jest formalnie identyczny z testem t Studenta różnic między średnimi, gdy nie mamy do czynienia z parami związanymi i przy założeniu
Z oznaczeniem tym spotkamy się jeszcze przy obliczaniu NIR-u
4
jednorodności wariancji. Gdybyśmy te same dane sprawdzali oboma testami, to miedzy otrzymanym stosunkiem F a otrzymaną statystyka związek t2=F. Można zatem uważać test t Studenta za specjalny przypadek stosowania analizy wariancji do porównywania dwóch tylko grup.
Test a posteriori Tukeya (metoda T)
Poszukujemy tzw. najmniejszej istotnej różnicy (NIR) definiowanej jako:
NIR=(wartość krytyczna) (błąd standardowy)
Wartość krytyczna to rozstęp studentyzowany, który podany jest w tabeli H1. Rozstępy te (Q0,05; a; df) – podawany jest dla liczby zabiegów (grup) a i liczby stopni swobody df (liczby stopni swobody df =N-a, gdzie N – liczba wszystkich pomiarów) przy założonym poziomie istotności 0,05. NIR obliczamy na podstawie wariancji wewnątrz grup (zabiegów), czyli wariancji niewyjaśnionej zwanej też składnikiem błędu. Posługujemy się zatem stopniami swobody i oszacowaniem wariancji niewyjaśnionej (wewnątrzgrupowej).
Test Tukeya przeprowadzimy opierając się na danych o przeżywalności chrząszczy na różnych pożywkach. Dla poziomu istotności 0,05, liczby zabiegów a=4 i df=16 otrzymujemy Q0,05; 4; 16=4,05.
Błąd standardowy s obliczamy według wzoru ,
2
n sx s (3)
gdzie s2 oznacza oszacowanie wariancji niewyjaśnionej czyli składnika błędu, zaś n to liczba powtórzeń w jednym zabiegu.
W naszym przykładzie:
Składnik s2 został już obliczony i zamieszczony w tabeli 1.
NIR liczymy według wzoru
df x
a s
Q
NIR 0,05; ; ,(4) W rozpatrywanym przykładzie
934 , 5 1
7 ,
2 18 n sx s
NIR=4,05 1,934=7,8327
Dla każdego zabiegu (grupy) obliczamy średnią arytmetyczną X z pomiarów oraz dolny d i górny g zakres każdej średniej dla porównania ze średnimi pozostałych zabiegów. Zakresy te obliczamy ze wzorów
2 /
2 / NIR X g
NIR X
d (5)
1 Wartośći krytyczne są dostępne też w internecie. Szczególnie polecam ten pierwszy link.
http://academic.udayton.edu/gregelvers/psy216/tables/qtab.htm http://cse.niaes.affrc.go.jp/miwa/probcalc/s-range/
5 W naszym przykładzie:
MP0 MP5 MP2 MPR
X 59,4 68,4 65,8 65,0
d 55,48 64,48 61,88 61,08
g 63,32 72,32 69,72 68,92
Średnie i ich zakresy można przedstawić graficznie
Jeżeli zakresy średnich z dwóch zabiegów zachodzą na siebie, to znaczy, że między tymi zabiegami nie ma statystycznie istotnych różnic, jeżeli nie zachodzą, to znaczy, że są statystycznie istotne różnice i trzeba odrzucić hipotezę zerową o braku różnic w porównywanej parze zabiegów. Innymi słowy, hipotezę zerową dla dwóch średnich odrzucamy wówczas, gdy NIR jest mniejsza aniżeli różnica miedzy tymi średnimi.
Oprócz graficznego przedstawienia tych wyników można przedstawić wyniki porównania każdego zabiegu z każdym. Zabiegi należy ułożyć według wzrastającej średniej.
By dokonać interpretacji należy odjąć od siebie poszczególne średnie dla grup. Jeżeli:
1. NIR < różnica średnich – odrzucamy H0, są różnice statystycznie istotne 2. NIR różnica średnich – nie ma różnic
Dla naszego przykładu
MP0 MPR MP2 MP5
X 59,4 65,0 65,8 68,4
MP0 - - +
MPR - -
MP2 -
Z tego porównania wyraźnie widać, że różnice w przeżywalności są statystycznie istotne przy porównywaniu wyników zabiegu MP0 z zabiegiem MP5. Porównania pozostałych 5 par nie pozwalają na odrzucenie hipotezy zerowej, ponieważ, ich zakresy zachodzą na siebie.
Książki godne uwagi:
Łomnicki A. 2003. Wprowadzenie do statystyki dla przyrodników. PWN Warszawa.
Zieliński R. 1972. Tablice statystyczne. PWN Warszawa.
MP0
MP5
MP2 MPR
50 55 60 65 70 75
6
Przedstawione powyżej obliczenia można wykonać korzystając narzędzi zainstalowanych w programie Microsoft Excell. Ponieważ to narzędzie nie jest standardowo zainstalowane należy to zrobić samemu. W tym celu po uruchomieniu programu należy wejść w opcje Narzędzia, a następnie wybrać Dodatki. W okienku, które się pojawi należy zaznaczyć pierwsze trzy pozycje: Aktualizowanie łączy dodatków, Analysis ToolPak, oraz Analysis ToolPak-VBA. Wybór należy potwierdzić poprzez naciśnięcie przycisku OK.
Teraz można przystąpić do wprowadzenia danych. Dane mogą być wprowadzane w postaci wierszy lub kolumn. W naszym przykładzie dane są podane w wierszach.
MP0 58 60 51 66 62
MP5 65 70 64 75 68
MP2 69 62 70 63 65
MPR 63 68 68 60 66
Po wprowadzeniu danych ponownie rozwija się menu Narzędzia, a z niego wybiera się opcje Analiza Danych. W okienku, które się pojawia wybiera się Analiza wariancji:
jednoczynnikowa. Następnie pojawia się kolejne okno dialogowe. Jako Zakres wejściowy podaje się całość naszych danych (włącznie z nazwami), następnie należy wybrać sposób w jaki dane są podawane: wiersze albo kolumny (w naszym przykładzie wiersze). Ponieważ zaznaczyliśmy w zakresie wejściowym kolumnę z nazwami to w oknie dialogowym też to należy to zaznaczyć. Poziom istotności wybieramy, w zależności od potrzeb (zwykle 0,05 lub 0,01). Następnie potwierdzamy wybór przez przyciśnięcie przycisku OK. Ponieważ nic nie zmienialiśmy w opcjach wyjścia to wyniki pojawia się na nowym arkuszu w formie tabeli. Wygląda to następująco:
Analiza wariancji: jednoczynnikowa
PODSUMOWANIE
Grupy Licznik Suma Średnia Wariancja
MP0 5 297 59,4 30,8
MP5 5 342 68,4 19,3
MP2 5 329 65,8 12,7
MPR 5 325 65 12
ANALIZA WARIANCJI
Źródło wariancji SS df MS F Wartość-p Test F
Pomiędzy grupami 215,35 3 71,78333 3,838681 0,030278 3,238867
W obrębie grup 299,2 16 18,7
Razem 514,55 19
Proszę zwrócić uwagę, że w tabeli powyżej w kolumnie Test F podano odpowiednie krytyczne wartości F (strona 3) przy założonym poziomie istotności (w tym przypadku 0,05), co uniezależnia nas od tabeli G.
Dla porównania załączam poniżej tabelkę z wynikami, którą uprzednio sami sporządziliśmy.
Wytłuszczone dane są niezbędne do liczenia NIR-u:
Źródła zmienności SK df Oszacowanie
wariancji F
Całkowita (ogólna) Między pożywkami Błąd
514,55 215,35 299,20
19 3 16
71,783 18,700
3,839
7
Wielkość Wartość-p można obliczyć korzystając z funkcji Rozkład F wpisując: X to wartość, dla której ta funkcja ma być obliczona czyli F, Stopnie_swobody1 to licznik stopni swobody (df pomiędzy grupami), Stopnie_swobody2 to mianownik stopni swobody (df w obrębie grup). W naszym przypadku to odpowiednio 3,838681, 3 i 16.
Z kolei wielkość Test F można obliczyć też używając funkcji Rozkład F odwrócony wpisując: Prawdopodobieństwo to prawdopodobieństwo związane ze skumulowanym rozkładem F-Snedecora czyli poziom istotności, Stopnie_swobody1 to licznik stopni swobody (df pomiędzy grupami), Stopnie_swobody2 to mianownik stopni swobody (df w obrębie grup).
W naszym przypadku będzie to odpowiednio 0,05 oraz df (czyli 3 i 16).
W przypadku korzystania z pakietu Open Office sprawa się trochę komplikuje, gdyż nie jest dostępna funkcja, który by równie łatwo dokonywała wszystkich obliczeń. Należy wykorzystać funkcję: =ODCH.KWADRATOWE(xx:yy) w celu obliczenia wewnątrzgrupowej SK (strona 2). Jako zakres funkcji, czyli xx:yy podajemy dane dla poszczególnych grup. Następnie sumujemy wyniki uzyskane dla poszczególnych grup i uzyskujemy wewnątrzgrupową SK. Teraz należy tylko, tak jak na stronie 3, otrzymany wynik podzielić przez liczbę stopni swobody df (liczby stopni swobody df =N-a, gdzie N – liczba wszystkich pomiarów, a - liczba zabiegów (grup)).
MP0 MP5 MP2 MPR
58 65 69 63
60 70 62 68
51 64 70 68
66 75 63 60
62 68 65 66
Odchylenie kwadratowe 123,2 77,2 50,8 48
Suma odchyleń 299,2
N=20 a=4 df=16
Oszacowanie wariancji między grupami 18,7
Dalsze postępowanie zostało juz omówione wcześniej (strona 4). Przy opracowaniu części dotyczącej korzystania z Open Office korzystałem z pomocy dr Jacka Rożnowskiego, za co jestem mu serdecznie wdzięczny.
8
W przypadku użycia programu STATISTICA dane należy uporządkować w następujący sposób:
Rodzaj pożywki Ilość chrząszczy
1 MP0 58
2 MP0 60
3 MP0 51
4 MP0 66
5 MP0 62
6 MP5 65
7 MP5 70
8 MP5 64
9 MP5 75
10 MP5 68
11 MP2 69
12 MP2 62
13 MP2 70
14 MP2 63
15 MP2 65
16 MPR 63
17 MPR 68
18 MPR 68
19 MPR 60
20 MPR 66
By przeprowadzić analizę wariancji wybieramy:
Statystyka → ANOVA → jednoczynnikowa ANOVA. W pojawiającym się oknie dialogowym jako predykatory jakościowe wybieramy kolumnę 1 (Rodzaj pożywki), a jako listę zmiennych zależnych kolumnę 2 (Ilość chrząszczy). Zatwierdzamy wybór i w kolejnym oknie mamy szereg możliwości:
Wszystkie efekty
SS Stopnie -
swobody MS F p
Wyraz wolny 83592,45 1 83592,45 4470,184 0,000000
Rodzaj pożywki 215,35 3 71,78 3,839 0,030278
Błąd 299,20 16 18,70
Wygląda znajomo.
z kolei tzw. NIR można obliczyć wybierając kolejno zakładki Więcej czynników a następnie Post hoc. Następnie należy wybrać jeden z dostępnych testów (sugeruje zaznaczyć opcję jednorodne grupy i wybrany poziom istotności).
