Jednoduchá analýza kovariance
Osnova:
1. Výchozí situace v analýze kovariance
2. Příklady situací řešených analýzou kovariance
3. Model analýzy kovariance
4. Součty čtverců v analýze kovariance
5. Předpoklady v analýze kovariance a jejich ověřování
6. Testy hypotéz v analýze kovariance
7. Příklad na analýzu kovariance
1. Výchozí situace v analýze kovariance
Předpokládáme, že soubor n objektů se rozpadá do 2≥r skupin podle variant (úrovní)
nějakého faktoru A, přičemž v i-té skupině je in objektů, ri ,,1 K= , nn
r
i
i =∑=1
. Na těchto
objektech sledujeme hodnoty závisle proměnné veličiny Y a nezávisle proměnné veličiny X.
Jejich hodnoty u j-tého objektu z i-té skupiny označíme yij a xij. Veličina X se nazývá
doprovodná proměnná neboli kovariáta. Předpokládáme, že mezi Y a X existuje ve všech
skupinách lineární závislost, přičemž regresní přímky, které tuto závislost modelují, jsou
rovnoběžné.
Chceme zjistit, zda střední hodnoty proměnné Y jsou stejné ve všech skupinách
odpovídajících úrovním faktoru A, pokud eliminujeme vliv kovariáty.
Úkoly tohoto typu řeší analýza kovariance (ANCOVA). Představuje spojení regresní
analýzy a analýzy rozptylu (ANOVY). Při ANOVĚ uvažujeme v jednotlivých skupinách
kolísání hodnot veličiny Y kolem průměru, zatímco při ANCOVĚ je to kolísání kolem regresní
přímky.
2. Příklady situací řešených analýzou kovariance
Obor Závisle proměnná veličina Y Faktor A Kovariáta X
Medicína Systolický krevní tlak Pohlaví pacienta Věk pacienta
Pedagogika Počet bodů v testu
z matematiky
Metoda výuky Skóre v testu čtenářské
gramotnosti
Biologie Počet zvukových pulsů, které vydá
samec cvrčka za 1 s
Druh cvrčka Teplota vnějšího
prostředí
Pedologie Obsah arzénu v půdě Lokalita, odkud
půda pochází
Obsah hliníku v půdě
3. Model analýzy kovariance
Vzájemný vztah závisle proměnné veličiny Y, kovariáty X a faktoru A popisuje model
( ) ijijiij xxy ε+−β+α+µ= , injri ,,1,,,1 KK == , kde
yij … j-tá hodnota závisle proměnné veličiny Y v i-té skupině
µ … společná část střední hodnoty veličiny Y
αi … efekt i-té úrovně faktoru A na veličinu Y
β … směrnice regresní přímky popisující závislost Y na X v každé skupině
xij … j-tá hodnota kovariáty X v i-té skupině
x … průměr hodnot kovariáty X
εij ... j-tá hodnota náhodné odchylky od modelu v i-té skupině
Předpokládáme, že ( )2
,0~ σε Nij a že platí reparametrizační rovnice ∑=
=α
r
i
iin
1
0 .
V modelu ( ) ijijiij xxy ε+−β+α+µ= člen αi vyjadřuje analýzu rozptylu a člen ( )xxij −β
vyjadřuje regresní část modelu.
Známe-li regresní parametr β, můžeme místo původních hodnot yij pracovat s upravenými
hodnotami yij
*
, kde ( )xxyy ijijij −β−=
*
. Porovnáním této a původní rovnice získáme výsledný
tvar modelu analýzy kovariance: ijiijy ε+α+µ=
*
, injri ,,1,,,1 KK == .
Odtud je vidět, že analýza kovariance je vlastně analýza rozptylu aplikovaná na upravené
hodnoty závisle proměnné veličiny Y. Upravená hodnota je vlastně původní hodnota veličiny Y
přepočítaná pomocí regresního vztahu mezi Y a X na průměrnou hodnotu kovariáty X.
Ilustrace:
Y
X
yi
yi = α+ βxi
y*
i
x
4. Součty čtverců v analýze kovariance
Stejně jako v ANOVĚ je i v ANCOVĚ základem metody rozklad celkového součtu čtverců na
složky. Jeho vyjádření je však složitější než v ANOVĚ, protože kromě různých součtů čtverců
pro Y musíme uvažovat také součty čtverců pro kovariátu X a součty čtverců pro obě proměnné
X a Y. Rozlišujeme tedy tři typy součtů čtverců:
a) pro celkovou variabilitu:
( ) ( )∑∑= =
−=
r
i
n
j
ijT
i
yyyS
1 1
2
… variabilita jednotlivých pozorování veličiny Y kolem celkového
průměru
( ) ( )∑∑= =
−=
r
i
n
j
ijT
i
xxxS
1 1
2
… variabilita jednotlivých pozorování kovariáty X kolem celkového
průměru
( ) ( )( )∑∑= =
−−=
r
i
n
j
ijijT
i
yyxxyxS
1 1
, … společná variabilita jednotlivých pozorování veličin Y a X
kolem jejich celkových průměrů
b) pro meziskupinovou variabilitu
( ) ( ) ( )∑∑∑ == =
−=−=
r
i
ii
r
i
n
j
iA yynyyyS
i
1
2
1 1
2
… variabilita skupinových průměrů veličiny Y kolem
celkového průměru
( ) ( ) ( )∑∑∑ == =
−=−=
r
i
ii
r
i
n
j
iA xxnxxxS
i
1
2
1 1
2
… variabilita skupinových průměrů kovariáty X kolem
celkového průměru
( ) ( )( ) ( )( )∑∑∑ == =
−−=−−=
r
i
iii
r
i
n
j
iiA yyxxnyyxxyxS
i
11 1
, … společná variabilita skupinových
průměrů veličin Y a X kolem jejich celkových průměrů
c) pro vnitroskupinovou (reziduální) variabilitu
( ) ( )∑∑= =
−=
r
i
n
j
iijE
i
yyyS
1 1
2
… variabilita jednotlivých pozorování veličiny Y kolem příslušných
skupinových průměrů
( ) ( )∑∑= =
−=
r
i
n
j
iijE
i
xxxS
1 1
2
… variabilita jednotlivých pozorování kovariáty X kolem příslušných
skupinových průměrů
( ) ( )( )∑∑= =
−−=
r
i
n
j
iijiijE
i
yyxxyxS
1 1
, … společná variabilita jednotlivých pozorování veličin Y a X
kolem příslušných skupinových průměrů
Lze ukázat, že platí () () ()... EAT SSS += .
Parametr β v modelu ( ) ijijiij xxy ε+−β+α+µ= odhadneme metodou nejmenších čtverců:
( )
( )xS
yxS
E
E ,ˆ =β .
Dále zavedeme následující označení, které využijeme při testování hypotéz:
( ) ( )( )∑=
−−=
in
j
iijiijiE yyxxyxS
1
, , ,
( ) ( )∑=
−=
in
j
iijiE xxxS
1
2
, .
Platí tedy, že
( ) ( )∑=
=
r
i
iEE yxSyxS
1
, ,, ,
( ) ( )∑=
=
r
i
iEE xSxS
1
, .
5. Předpoklady v analýze kovariance a jejich ověřování
Analýzu kovariance lze korektně použít, jsou-li splněny následující předpoklady.
a) Jednotlivé skupiny jsou navzájem nezávislé.
Tento předpoklad musí být zajištěn vhodnou organizací experimentu, kterým se získají data.
b) Hodnoty veličiny Y jsou ve všech skupinách normálně rozloženy.
Ověření provedeme pomocí diagnostických grafů a testů normality. Při větších rozsazích výběrů
není ANCOVA citlivá na mírné porušení normality.
c) Hodnoty veličiny Y mají ve všech skupinách stejný rozptyl (předpoklad homoskedasticity).
Ověření provedeme pomocí krabicových diagramů a testů homogenity rozptylu. Mírné porušení
homoskedasticity příliš nevadí.
Poznámka: První tři předpoklady jsou shodné s předpoklady pro ANOVU. Další předpoklad je
specifický pro ANCOVU.
d) Regresní přímky modelující závislost Y na X jsou ve všech r skupinách rovnoběžné.
Ilustrace:
Nulová hypotéza H0: β=β==β=β :21 rK , alternativní hypotéza H1: aspoň jedna dvojice
regresních koeficientů se liší.
Pro testování rovnoběžnosti regresních přímek použijeme statistiku
( )
( )
( )
( )
( ) ( )
( )∑
∑
=
=
−
−
⋅
−
−
= r
i
xS
yxS
E
iE
iE
r
i
xS
yxS
iE
iE
iE
iE
yS
xS
yxS
r
rn
T
1
,
,
2
,
1
,
0
,
2
,
,
2
,
,
1
2
, která se za platnosti H0 řídí rozložením F(r-1, n-2r). H0
tedy zamítáme na hladině významnosti α, když ( )rnrFT 2,110 −−≥ α− .
V případě, že hypotézu o rovnoběžnosti regresních přímek zamítneme, nelze použít standardní
ANCOVU.
6. Testy hypotéz v analýze kovariance
Pomocí ANCOVY testujeme dvě hypotézy. První se týká regrese Y na X, druhá vztahu Y a
faktoru A.
a) Test nulovosti regrese
Nulová hypotéza 0:0 =βH , alternativní hypotéza 0:1 ≠βH .
Testová statistika pro test nulovosti regrese má tvar
( )
( )
( )
( ) ( )
( )xS
yxS
yS
xS
yxS
rnT
E
E
E
E
E
2
2
0
,
,
1
−
⋅−−=
a za platnosti H0 se řídí rozložením F(1, n-r-1). H0 tedy
zamítáme na hladině významnosti α, když ( )1,110 −−≥ α− rnFT .
V případě, že hypotézu o nulovosti regrese nezamítneme, nemá smysl provádět ANCOVU, stačí
ANOVA.
b) Test hypotézy o nevýznamnosti faktoru A
Nulová hypotéza H0: 021 =α==α=α rK , alternativní hypotéza H1: aspoň jeden efekt faktoru
A je nenulový.
Pro testování této hypotézy slouží statistika
( ) ( )
( )
( ) ( )
( )
( ) ( )
( )xS
yxS
yS
xS
yxS
yS
xS
yxS
yS
r
rn
T
E
E
E
E
E
E
T
T
T
2
22
0
,
,,
1
1
−
+−−
⋅
−
−−
=
, která se za platnosti H0 řídí rozložením
F(r-1, n-2r). H0 tedy zamítáme na hladině významnosti α, když ( )1,110 −−−≥ α− rnrFT .
Pokud zamítneme hypotézu o nevýznamnosti faktoru A, provedeme mnohonásobné
porovnávání, abychom zjistili, které dvojice skupin se liší na dané hladině významnosti.
7. Příklad na analýzu kovariance
Popis situace: Chovatel ústřic chtěl zjistit, zda poloha ústřic v různé výšce vodního sloupce
uměle ohřáté vody ovlivňuje jejich růst. Za tímto účelem provedl pilotní studii, která byla
organizována takto:
Náhodně vybral 200 ústřic a náhodně je rozdělil do 20 sáčků po 10 ústřicích. V nádrži, kam
přitéká chladicí voda z elektrárny, si vybral 5 lokací a do každé umístil 4 sáčky s ústřicemi.
Lokace jsou následující:
1 – chladné dno,
2 – chladný povrch,
3 – teplé dno,
4 – teplý povrch,
5 – střední hloubka, střední teplota (lokace 5 je považována za kontrolní).
Experiment probíhal po dobu jednoho měsíce. Na jeho počátku a na konci byl každý sáček
s ústřicemi zvážen a byla zaznamenána jeho hmotnost v gramech.
Roli faktoru A tedy hraje lokace (má 5 variant, v datovém souboru je označena jako proměnná
ID), závisle proměnnou veličinou Y je koncová hmotnost sáčku s deseti ústřicemi a
doprovodnou proměnnou – kovariátou X – je počáteční hmotnost.
Datový soubor:
1
ID
2
Y
3
X
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
chladne dno 924,21 771,12
chladne dno 1037,61 907,2
chladne dno 1068,795 935,55
chladne dno 878,85 759,78
chladny povrch 958,23 810,81
chladny povrch 898,695 759,78
chladny povrch 870,345 751,275
chladny povrch 861,84 759,78
teple dno 997,92 810,81
teple dno 824,985 635,04
teple dno 819,315 657,72
teple dno 856,17 691,74
teply povrch 992,25 830,655
teply povrch 765,45 618,03
teply povrch 1031,94 859,005
teply povrch 864,675 688,905
kontrola 697,41 578,34
kontrola 663,39 555,66
kontrola 859,005 711,585
kontrola 618,03 513,135
Úkolem je posoudit, zda lokace má vliv na koncovou hmotnost při eliminaci vlivu počáteční
hmotnosti. Problém budeme řešit pomocí analýzy kovariance.
7.1. Výpočet číselných charakteristik datového souboru
Rozkladová tabulka popisných statistik (chov_ustric.sta)
N=20 (V seznamu záv. prom. nejsou ChD)
ID Y
průměr
Y
N
Y
Sm.odch.
X
průměr
X
N
X
Sm.odch.
chladne dno 977,37 4 90,41 843,41 4 90,88
chladny povrch 897,28 4 43,58 770,41 4 27,23
teple dno 874,60 4 83,80 698,83 4 78,21
teply povrch 913,58 4 121,84 749,15 4 114,79
kontrola 709,46 4 104,87 589,68 4 85,65
Vš.skup. 874,46 20 123,18 730,30 20 114,47
Graf průměrů s 95% intervaly spolehlivosti
Graf průměrů a int. spolehlivosti (95,00%)
Y
X
chladne dno
chladny povrch
teple dno
teply povrch
kontrola
ID
300
400
500
600
700
800
900
1000
1100
1200
Hodnoty
Nejvyšších hodnot nabývá hmotnost na chladném dně, naopak nejnižší hodnoty vykazuje
kontrolní lokace.
7.2. Ověření předpokladů ANCOVY
a) Nezávislost daných pěti náhodných výběrů – splněno, zajištěno organizací experimentu.
b) Normalita hodnot veličiny Y v daných pěti náhodných výběrech – posoudíme pomocí NP
plotu a S-W testu.
Normální p-graf z Y; kategorizovaný ID
chov_ustric.sta 3v*20c
Pozorovaný kvantil
Oček.normál.hodnoty
ID: chladne dno
ID: chladny povrch
ID: teple dno
ID: teply povrch
ID: kontrola500 600 700 800 900 1000 1100
-1,2
-1,0
-0,8
-0,6
-0,4
-0,2
0,0
0,2
0,4
0,6
0,8
1,0
1,2
ID: chladne dno Y: SW-W = 0,9104; p = 0,4844
ID: chladny povrch Y: SW-W = 0,8837; p = 0,3545
ID: teple dno Y: SW-W = 0,777; p = 0,0669
ID: teply povrch Y: SW-W = 0,9379; p = 0,6416
ID: kontrola Y: SW-W = 0,8892; p = 0,3792
Ani v jednom případě nezamítáme normalitu na hladině významnosti 0,05.
c) Homogenita rozptylu hodnot veličiny Y v daných pěti náhodných výběrech – posoudíme
pomocí krabicového grafu a pomocí Levenova testu.
Kategoriz. krabicový graf: Y
Průměr
Průměr±SmCh
Průměr±1,96*SmCh
chladne dno
chladny povrch
teple dno
teply povrch
kontrola
ID
500
600
700
800
900
1000
1100
Y
Leveneův test homogenity rozpylů (chov_ustric.sta)
Označ. efekty jsou význ. na hlad. p < ,05000
Proměnná
SČ
efekt
SV
efekt
PČ
efekt
SČ
chyba
SV
chyba
PČ
chyba
F p
Y 9733,08 4 2433,270 25517,69 15 1701,179 1,430343 0,272145
Na hladině významnosti 0,05 nezamítáme hypotézu o homogenitě rozptylu.
d) Rovnoběžnost regresních přímek ve všech pěti lokacích
Jednorozm. výsledky pro každou záv. proměnnou (chov_ustric.sta)
Sigma-omezená parametrizace
Dekompozice typu II
Efekt
Stupně
volnosti
Y
SČ
Y
PČ
Y
F
Y
p
Abs. člen
ID
X
ID*X
Chyba
Celkem
1 15293457 15293457 67142,62 0,000000
4 1363 341 1,50 0,275182
1 125413 125413 550,60 0,000000
4 1116 279 1,22 0,360175
10 2278 228
19 288271
Zajímá nás řádek ID*X. Příslušná p-hodnota je 0,3602, tedy na hladině významnosti 0,05
nezamítáme hypotézu o shodě směrnic daných pěti regresních přímek.
Test můžeme ještě doplnit grafem:
Bodový graf z Y proti X; kategorizovaný ID
chov_ustric.sta 7v*20c
ID: chladne dno Y = 148,5897+0,9826*x
ID: chladny povrch Y = -259,3833+1,5014*x
ID: teple dno Y = 136,5891+1,0561*x
ID: teply povrch Y = 121,7847+1,0569*x
ID: kontrola Y = -12,2424+1,2239*x
X
Y
ID: chladne dno
ID: chladny povrch
ID: teple dno
ID: teply povrch
ID: kontrola450 500 550 600 650 700 750 800 850 900 950 1000
500
600
700
800
900
1000
1100
Vidíme, že směrnice regresních přímek popisujících závislost konečné hmotnosti na počáteční
hmotnosti v jednotlivých lokacích se pohybují od 0,9826 na chladném dně po 1,5014 na
chladném povrchu. Rozdíly směrnic regresních přímek však nejsou prokazatelné na hladině
významnosti 0,05.
7.3. Provedení ANCOVY
Nyní budeme testovat dvě hypotézy. První se týká regrese Y na X a tvrdí, že ve všech pěti
skupinách je směrnice regresní přímky nulová (test nulovosti regrese). Druhá se týká vztahu Y a
faktoru ID a tvrdí, že faktor ID je nevýznamný.
Jednorozm. výsledky pro každou záv. proměnnou (chov_ustric.sta)
Sigma-omezená parametrizace
Dekompozice typu II
Efekt
Stupně
volnosti
Y
SČ
Y
PČ
Y
F
Y
p
Abs. člen
X
ID
Chyba
Celkem
1 15293457 15293457 63092,26 0,000000
1 125413 125413 517,38 0,000000
4 9716 2429 10,02 0,000482
14 3394 242
19 288271
Vidíme, že p-hodnota na řádku X je blízká 0, tedy hypotézu o nulovosti regrese zamítáme na
hladině významnosti 0,05.
Na řádku ID je p-hodnota 0,000482, tedy vliv lokace na koncovou hmotnost ústřic je prokázán
na hladině významnosti 0,05.
Upozornění: Při provádění ANCOVY je zapotřebí zvolit Součet čtverců Typ II (parciální),
protože ne všechny úrovně faktoru jsou zastoupeny ve všech hodnotách kovariáty.
Provedené testy ještě doplníme o odhad regresního koeficientu β a výpočet upravených průměrů.
Odhady parametrů (chov_ustric.sta)
Sigma-omezená parametrizace
Efekt
Úroveň
Efekt
Sloupec Y
Param.
Y
Sm.Ch.
Y
t
Y
p
-95,00%
LmtSpol.
+95,00%
LmtSpol.
Y
Beta (ß)
Y
Sm.Ch. ß
-95,00%
LmtSpol.
+95,00%
LmtSpol.
Abs. člen
X
ID
ID
ID
ID
1 83,4139 34,95089 2,38660 0,031671 8,4517 158,3761
2 1,0832 0,04762 22,74608 0,000000 0,9810 1,1853 1,006667 0,044257 0,911745 1,101588
chladne dno 3 -19,6150 8,80317 -2,22818 0,042776 -38,4959 -0,7341 -0,103332 0,046375 -0,202796 -0,003867
chladny povrch 4 -20,6303 7,22004 -2,85736 0,012665 -36,1157 -5,1448 -0,108680 0,038035 -0,190257 -0,027103
teple dno 5 34,2278 7,12217 4,80581 0,000279 18,9523 49,5033 0,180312 0,037519 0,099840 0,260783
teply povrch 6 18,7021 7,02037 2,66397 0,018517 3,6449 33,7593 0,098522 0,036983 0,019201 0,177843
Na řádku X, ve sloupci Y Param. najdeme odhad 1,0832.
ID; Průměry MNČ (chov_ustric.sta)
Současný efekt: F(4, 14)=10,021, p=,00048
(Vypočteno pro průměry spoj. nezáv. prom.)
Č. buňky
ID Y
Průměr
Y
Sm.Ch.
Y
-95,00%
Y
+95,00%
N
1
2
3
4
5
chladne dno 854,8407 9,46656 834,5370 875,1445 4
chladny povrch 853,8255 8,01554 836,6339 871,0171 4
teple dno 908,6835 7,92750 891,6808 925,6863 4
teply povrch 893,1578 7,83617 876,3509 909,9647 4
kontrola 861,7712 10,26834 839,7478 883,7946 4
ID; Průměry MNČ
Současný efekt: F(4, 14)=10,021, p=,00048
(Vypočteno pro průměry spoj. nezáv. prom.)
Vertikály označují 0,95 intervaly spolehlivosti
chladne dno
chladny povrch
teple dno
teply povrch
kontrola
ID
820
830
840
850
860
870
880
890
900
910
920
930
940
YPrům. spoj. nez. prom.:
X: 730,296
Nejnižší upravený průměr konečné hmotnosti pozorujeme na lokaci chladný povrch, naopak
nejvyšší na lokaci teplé dno.
Na závěr provedeme mnohonásobné porovnávání, abychom identifikovala dvojice lokací, které
se liší na hladině významnosti 0,05.
Tukeyův HSD test; proměnná Y (chov_ustric.sta)
Přibližné pravděpodobnosti pro post hoc testy
Chyba: Between MSE = 242,40, sv = 14,000
Č. buňky
ID {1}
977,37
{2}
897,28
{3}
874,60
{4}
913,58
{5}
709,46
1
2
3
4
5
chladne dno 0,000179 0,000151 0,000504 0,000151
chladny povrch 0,000179 0,289479 0,590163 0,000151
teple dno 0,000151 0,289479 0,022974 0,000151
teply povrch 0,000504 0,590163 0,022974 0,000151
kontrola 0,000151 0,000151 0,000151 0,000151
7.4. Srovnání výsledků ANCOVY a ANOVY
Pokud budeme testovat hypotézu o nevýznamnosti vlivu lokace na konečnou hmotnost sáčku
ústřic pomocí jednofaktorové ANOVY, bez eliminace vlivu počáteční hmotnosti, dostaneme
tabulku
Analýza rozptylu (chov_ustric.sta)
Označ. efekty jsou význ. na hlad. p < ,05000
Proměnná
SČ
efekt
SV
efekt
PČ
efekt
SČ
chyba
SV
chyba
PČ
chyba
F p
Y 159464,2 4 39866,04 128806,6 15 8587,105 4,642547 0,012239
Testová statistika nabývá hodnoty 4,6426, odpovídající p-hodnota je 0,0122, tedy na hladině
významnosti 0,05 považuje ANOVA vliv lokace za prokázaný. Při použití ANCOVY je však
testová statistika 10,02 a p-hodnota pouze 0,0005.
Tukeyova metoda mnohonásobného porovnávání poskytne tabulku p-hodnot:
Tukeyův HSD test; proměn.:Y (chov_ustric.sta)
Označ. rozdíly jsou významné na hlad. p < ,05000
ID
{1}
M=977,37
{2}
M=897,28
{3}
M=874,60
{4}
M=913,58
{5}
M=709,46
chladne dno {1}
chladny povrch {2}
teple dno {3}
teply povrch {4}
kontrola {5}
0,739248 0,537759 0,862937 0,007414
0,739248 0,996604 0,999099 0,075213
0,537759 0,996604 0,973812 0,138138
0,862937 0,999099 0,973812 0,047617
0,007414 0,075213 0,138138 0,047617
Vidíme, že se liší pouze dvojice lokací (chladné dno, kontrola) a (teplý povrch, kontrola),
zatímco ANCOVA ukázala, že se neliší pouze dvojice lokací (chladný povrch, teplé dno) a
(chladný povrch, teplý povrch).