Statistick´ e metody v ekonomii: Teoretick´ a v´ ychodiska, Jednofaktorov´ a a dvoufaktorov´ a anal´ yza rozptlylu Ing. Michael Rost, Ph.D.
Co je vlastnˇ e c´ılem? C´ılem statistick´ eho zpracov´ an´ı dat je pod´ an´ı informace o vlastnostech, povaze a z´ akonitostech projevuj´ıc´ıch se na pozorovan´ ych datech. Statistika zahrnuje z´ısk´ av´ an´ı, anal´ yzu a objektivn´ı interpretaci z´ıskan´ ych dat. Tj. zaˇ c´ın´ a jiˇ z pˇ red samotn´ ym proveden´ım experimentu!
c Rost 2007 °
Nˇ ekolik doporuˇ cen´ı:
• Definujte sv˚ uj probl´ em, kter´ y hodl´ ate ˇ reˇ sit (a to co moˇ zn´ a nejjednoduˇ seji), a vytvoˇ rte si sv´ e pracovn´ı hypot´ ezy.
• Urˇ cete, co budete mˇ eˇ rit a jak to budete mˇ eˇ rit.
• Vˇ enujte dostateˇ cnou pˇ r´ıpravu vaˇ semu experimentu. Je d˚ uleˇ zit´ a.
• Randomizujte, tj. zn´ ahodˇ nujte. M˚ uˇ zete se tak vyhnout systematick´ ym chyb´ am (napˇ r´ıklad pˇ ri odeˇ c´ıt´ an´ı hodnot z pˇ r´ıstroje).
c Rost 2007 °
Nˇ ekolik doporuˇ cen´ı: • Pˇ ri anal´ yze dat je nutno rozliˇ sovat s jak´ ymi typy dat pracujeme. Mˇ ejte na pamˇ eti, ˇ ze ne vˇ sechny metody jsou vhodn´ e pro vˇ sechny typy dat.
• Pokud vyuˇ z´ıv´ ate deskriptivn´ı statistiky k popisu analyzovan´ eho souboru, uvˇ edomte si, ˇ ze aritmetick´ y pr˚ umˇ er nemus´ı b´ yt vˇ zdy ”tou pravou” charakteristikou polohy. Ud´ avejte i dalˇ s´ı statistiky.
• Vizualizujte! Pokud to jde, pouˇ z´ıvejte spolu s ˇ c´ısly i grafickou reprezentaci, m˚ uˇ ze pomoci (pouˇ z´ıvejte vˇ sak vhodn´ e grafick´ e vyj´ adˇ ren´ı).
c Rost 2007 °
Statistick´ e zhodnocen´ı ve 3 ´ urovn´ıch V pr˚ ubˇ ehu statistick´ e anal´ yzy dat postupujte v nˇ ekolika ´ urovn´ıch:
♣ Explorativn´ı ´ uroveˇ n (EDA). Popisn´ e statistiky a grafick´ e zhodnocen´ı.
♣ Form´ aln´ı statistick´ y pˇ r´ıstup. Testy a testov´ an´ı hypot´ ez ˇ ci r˚ uzn´ e statistick´ e metody.
♣ Diagnostika. Zhodnocen´ı, zda byly dodrˇ zeny pˇ redpoklady pro pouˇ zit´ı v´ ami vybran´ ych statistick´ ych metod. Pˇ r´ıpadnˇ e proved’te r˚ uzn´ a n´ apravn´ a opatˇ ren´ı – transformace atd... c Rost 2007 °
´ Uvod do probl´ emu V ekonomick´ e praxi se m˚ uˇ zeme setkat se situac´ı, ve kter´ e potˇ rebujete simult´ annˇ e otestovat shodu k stˇ redn´ıch hodnot, kde k je vˇ etˇ s´ı neˇ z 2. Jak to prov´ est, to bude n´ apln´ı t´ eto pˇ redn´ aˇ sky ; −)
c Rost 2007 °
Pˇ r´ıklad Klinickou studi´ı byla sledov´ ana z´ avislost mezi dietou a dobou, za kterou dojde ke koagulaci krve. Byly namˇ eˇ reny tyto hodnoty: Typ diety A
B
C
D
62 60 63 59 64 65 62 62
63 67 71 64 65 66 58 59
68 66 71 67 68 68 63 63
56 62 60 68 63 64 63 59
yij
Maj´ı r˚ uzn´ e diety vliv na dobu, za kterou dojde ke koagulaci krve? c Rost 2007 °
Pˇ r´ıklad
c Rost 2007 °
V´ ychoz´ı situace Situaci kter´ eˇ cel´ıme, lze obecnˇ e popsat n´ asleduj´ıc´ım zp˚ usobem (pˇ redpokl´ ad´ ame, ˇ ze m´ ame k soubor˚ u): ˇ c´ıslo v´ ybˇ eru
poˇ cet prvk˚ u
zjiˇ stˇ en´ e hodnoty sledovan´ eho znaku
pr˚ umˇ er
rozptyl
1 2 ... i ... k
n1 n2 ... ni ... nk
y11, y12, · · · , y1j , · · · , y1n1 y21, y22, · · · , y2j , · · · , y2n2 ... yi1, yi2, · · · , yij , · · · , yini ... yk1, yk2, · · · , ykj , · · · , yknk
y¯1 y¯2 ... y¯i ... y¯k
s2 1 s2 ...2 s2 ...i s2 k
c Rost 2007 °
Obecn´ y postup
1. Zformulovat hypot´ ezy: H0 vs. HA. 2. Ovˇ eˇ ren´ı pˇ redpoklad˚ u
Pozn´ amka : Uvˇ edomte si, ˇ ze pokud prov´ ad´ıte form´ aln´ı anal´ yzu ˇ ci statistick´ e testov´ an´ı, pˇ ri kter´ em vyuˇ z´ıv´ ate p-value, vych´ az´ıte z´ aroveˇ n z jist´ ych pˇ redpoklad˚ u. Ty vˇ sak nemus´ı b´ yt splnˇ eny. Stupeˇ n validity z´ıskan´ eho p-value z´ aleˇ z´ı na tom jakou shodu vykazuj´ı naˇ se data s teoretick´ ymi rozdˇ elen´ımi. Proto kaˇ zdop´ adnˇ e ovˇ eˇ rujte pˇ redpoklady vaˇ sich model˚ u! c Rost 2007 °
Ovˇ eˇ ren´ı pˇ redpoklad˚ u Pˇ red vlastn´ı anal´ yzou rozptylu je nutno odpovˇ edˇ et na nˇ ekolik ot´ azek:
• Poch´ az´ı jednotliv´ e v´ ybˇ ery z norm´ aln´ıho rozdˇ elen´ı?
• Jsou jednotliv´ e v´ ybˇ ery nez´ avisl´ e?
• Lze se domn´ıvat, ˇ ze v´ ybˇ ery maj´ı shodn´ e rozptyly?
Prvn´ı a tˇ ret´ı poˇ zadavek lze ovˇ eˇ rit prostˇ rednictv´ım r˚ uzn´ ych test˚ u.
c Rost 2007 °
Obecn´ y postup pokraˇ cov´ an´ı
2. Stanovit hodnotu α, nejˇ castˇ eji vol´ıme α = 0, 05 nebo α = 0, 01.
3. Zvolit adekv´ atn´ı testov´ e krit´ erium a stanovit hodnotu testov´ eho krit´ eria
4. Zjistit zda F ∈ K nebo zda p-value ≤ α
5. Z´ avˇ er
c Rost 2007 °
Specifikace nulov´ e a alternativn´ı hypot´ ezy V pˇ r´ıpadˇ eˇ ze potˇ rebujeme simult´ annˇ e otestovat shodu k stˇ redn´ıch hodnot, je nulov´ a a alternativn´ı hypot´ eza specifikov´ ana jako: H0 : µ1 = µ2 = . . . = µk−1 = µk HA : non H0 . Z´ aroveˇ n vˇ sak testujeme jeˇ stˇ e homoskedasiticitu: 2 = σk2 H0 : σ12 = σ22 = . . . = σk−1
HA : non H0 , a tu testujeme zpravidla jako prvn´ı!!!
c Rost 2007 °
Testov´ an´ı homoskedasticity
c Rost 2007 °
Testy homoskedasticity Pˇ redpoklad homoskedasticity (shody rozptyl˚ u) je moˇ zno otestovat napˇ r´ıklad prostˇ rednictv´ım tzv. Bartlettova testu. Bartlett˚ uv test je univerz´ aln´ım testem v tom smyslu, ˇ ze jej lze vyuˇ z´ıt k hodnocen´ı homoskedasticity u vyv´ aˇ zen´ ych i nevyv´ aˇ zen´ ych soubor˚ u. Testujeme hypot´ ezu: H0 : σ12 = σ22 = · · · = σk2 , HA : non H0 . Testov´ ym krit´ eriem Bartlettova testu je veliˇ cina B, kter´ a je definov´ ana jako B = [(n − k)ln s2 −
k X i=1
(ni − 1)ln s2 i ]/C . c Rost 2007 °
Testy homoskedasticity ribliˇ znˇ e plat´ı B ∼ χ2(k − 1). Plat´ı-li H0 a je-li ni ≥ 6, pak pˇ Testovanou hypot´ ezu zam´ıt´ ame pokud plat´ı B ≥ χ2 1−α (k − 1) . Jednotliv´ e symboly vyuˇ zit´ e pˇ ri v´ ypoˇ ctu testov´ e statistiky lze definovat takto: n
i X 1 2 si = (yij − y¯i)2 ni − 1 j=1
i = 1, · · · , k ,
celkov´ y rozptyl s jako n
k X i 1 X s= (yij − y¯i)2 , n − 1 i=1 j=1
a konstantu C
C =1+
k X
1 1 − /3(k − 1) . n − 1 n − k i=1 i c Rost 2007 °
Hartley˚ uv test Dalˇ s´ım testem je tzv. Hartley˚ uv test homoskedasticity. Testovac´ı statistika m´ a v pˇ r´ıpadˇ e Hartleyova testu tvar: max s2 i . Fmax = min s2 i Ke stanoven´ı kritick´ eho oboru je nutno vyuˇ z´ıt speci´ alnˇ e sestrojen´ ych tabulek, nebot’ testovan´ a dvojice rozptyl˚ u nen´ı n´ ahodnˇ e zvolena. Nulovou hypot´ ezu o shodˇ e rozptyl˚ u zam´ıt´ ame na hladinˇ e v´ yznamnosti α, pokud testovac´ı statistika Fmax pˇ rekroˇ c´ı jistou kritickou hodnotu.
c Rost 2007 °
Cochran˚ uv test Dalˇ s´ım testem pro ovˇ eˇ ren´ı homoskedasticity je tzv. Cochran˚ uv test. V pˇ r´ıpadˇ e jeho pouˇ zit´ı zam´ıt´ ame H0, hypot´ ezu pokud hodnota testov´ eho krit´ eria s2 max C= 2 2 + . . . + s s1 + s2 2 k pˇ rekroˇ c´ı kritickou hodnotu Cochranovy statistiky. Jin´ ymi slovy, pokud hodnota C bude n´ aleˇ zet do kritick´ eho oboru, kter´ y je definov´ an jako K = {C ≥ C1−α(k, n − 1)} , pak zam´ıt´ ame hypot´ ezu o shodˇ e rozptyl˚ u.
c Rost 2007 °
Leven˚ uv test homogenity rozptyl˚ u Leven˚ uv test v podstatˇ e prov´ ad´ı anal´ yzu rozptylu na rezidu´ıch. Vyuˇ z´ıv´ a pˇ ritom promˇ ennou zij = |yij − y¯i| pro i = 1, 2, · · · , k a j = 1, 2, · · · , ni. V´ ysledn´ a hodnota testov´ e statistiky F je porovn´ av´ ana s kritickou hodnotou F -rozdˇ elen´ı s k − 1 a n − k stupni volnosti. Pro jist´ e pˇ r´ıpady jsou navrˇ zeny i modifikace Levenova testu. V pˇ r´ıpadˇ eˇ sikmosti souboru lze vyuˇ z´ıt m´ısto y¯i. medi´ anu. V pˇ r´ıpadˇ e v´ yrazn´ e ˇ spiˇ catosti souboru je pak m´ısto y¯i. doporuˇ cov´ an 10 % oˇ rezan´ y pr˚ umˇ er.
c Rost 2007 °
Anal´ yza rozptylu - ANOVA jednofaktorov´ a
• Necht’ Yij je n´ ahodnou veliˇ cinou oznaˇ cuj´ıc´ı j-t´ e pozorov´ an´ı v r´ amci i-t´ e skupiny. Symbol yij pak bude pˇ redstavovat pozorovanou hodnotu veliˇ ciny Yij z´ıskanou proveden´ım experimentu.
• Symbolem ni oznaˇ c´ıme poˇ cet pozorov´ an´ı v i-t´ e skupinˇ e. Pr˚ umˇ ery v jednotliv´ ych skupin´ ach tj. y¯1, y¯2, · · · , y¯k z´ısk´ ame jako n
1 Xi y¯i = yij . ni j=1
c Rost 2007 °
Anal´ yza rozptylu - ANOVA jednofaktorov´ a • Rozptyly uvnitˇ r jednotliv´ ych skupin oznaˇ c´ıme jako s2 i , kde i = 1, 2, · · · , k. Je zˇ rejm´ e, ˇ ze: n
i X 1 (yij − y¯i)2 s2 i = ni − 1 j=1
• Vnitroskupinov´ a tzv. pr˚ umˇ ern´ a rezidu´ aln´ı suma ˇ ctverc˚ u: n
k X i 1 X M SSr = (yij − y¯i)2 n − k i=1 j=1
c Rost 2007 °
Anal´ yza rozptylu - ANOVA jednofaktorov´ a • Celkov´ y pr˚ umˇ er oznaˇ c´ıme jako y¯ kde n
k X i 1 X y¯ = yij n i=1 j=1
n=
k X
ni .
i=1
• Pr˚ umˇ ern´ a suma ˇ ctverc˚ u vlivem r˚ uzn´ ych ´ urovn´ı faktor˚ u (skripta: Rozptyl mezi tˇ r´ıdami): k 1 X M SSA = ni(¯ yi − y¯)2 k − 1 i=1
c Rost 2007 °
Testov´ e krit´ erium
• Testov´ e krit´ erium M SSA F = . M SSr Kde pro testov´ e krit´ erium F za platnosti nulov´ e hypot´ ezy plat´ı: F ∼ F (v1 = k − 1; v2 = n − k) • Pokud symbolem F oznaˇ c´ıme hodnotu testov´ eho krit´ eria F urˇ cenou na z´ akladˇ e proveden´ eho experimentu, pak lze p-value definovat takto: P (F(k − 1; n − k) > F )
c Rost 2007 °
Tabulka anal´ yzy rozptylu V podstatˇ e je tato testovac´ı statistika zaloˇ zena na pomˇ eru pr˚ umˇ ern´ ych meziskupinov´ ych a vnitroskupinov´ ych souˇ ct˚ uˇ ctverc˚ u. V´ ysledky anal´ yzy rozptylu se zapisuj´ı do tzv. tabulky anal´ yzy rozptylu. Ta mˇ ela v minulosti sv˚ uj v´ yznam z hlediska v´ ypoˇ ct˚ u. V nejjednoduˇ sˇ s´ım pˇ r´ıpadˇ e m´ a n´ asleduj´ıc´ı podobu: Zdroj variability
Souˇ cet ˇ ctverc˚ u
Poˇ cet stupˇ n˚ u volnosti
Pr˚ umˇ ern´ y ˇ ctverec
Faktor
SSA
k−1
M SSA =
SSA k−1
Rezidu´ aln´ı
SSr
n−k
M SSr =
SSr n−k
Celkov´ y
SST
n−1
F
F =
Dosaˇ zen´ a hladina p
M SSA M SSr
p
c Rost 2007 °
Jak na to v R Pokraˇ cujme v ˇ reˇ sen´ı motivaˇ cn´ıho pˇ r´ıkladu. Nejprve data naˇ cteme a zn´ azorn´ıme. data<-read.table("dieta.txt",header=TRUE) data names(data) plot(koagulace~dieta,data,ylab="Cas koagulace krve",col="red") with(data,stripchart(koagulace~dieta,vertical=TRUE,method="stack", xlab="Dieta",ylab="Cas koagulace krve"))
c Rost 2007 °
60
60
A B C
dieta D 65
Cas koagulace krve
65
Cas koagulace krve
70
70
Vytvoˇ ren´ e grafy
A B C D
Dieta
c Rost 2007 °
Pokraˇ cov´ an´ı pr´ ace v R - Tabulka anal´ yzy rozptylu V´ ysledn´ a tabulka anal´ yzy rozptylu model<-lm(koagulace~dieta,data) anova(model) Analysis of Variance Table Response: koagulace Df Sum Sq Mean Sq F value Pr(>F) dieta 3 122.344 40.781 3.8823 0.01939 * Residuals 28 294.125 10.504 --Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1
Lze ˇ r´ıci, ˇ ze se doba koagulace krve liˇ s´ı v z´ avislosti na pod´ avan´ e dietˇ e.
c Rost 2007 °
pokraˇ cov´ an´ı pr´ ace v R cn´ı hladina Model yij = µ + αi + ²ij , kde α1 = 0 - referenˇ model<-lm(koagulace~dieta,data) summary(model) Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 62.125 1.146 54.216 < 2e-16 *** dietaB 2.000 1.621 1.234 0.22740 dietaC 4.625 1.621 2.854 0.00803 ** dietaD -0.250 1.621 -0.154 0.87850 --Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1 Residual standard error: 3.241 on 28 degrees of freedom Multiple R-Squared: 0.2938, Adjusted R-squared: 0.2181 F-statistic: 3.882 on 3 and 28 DF, p-value: 0.01939
c Rost 2007 °
... R Pro spr´ avn´ e pochopen´ı k´ odov´ an´ı jednotliv´ ych ´ urovn´ı m˚ uˇ zeme pouˇ z´ıt pˇ r´ıkaz model.matrix() (Intercept) dietaB dietaC dietaD 1 1 0 0 0 2 1 0 0 0 3 1 0 0 0 4 1 0 0 0 5 1 0 0 0 6 1 0 0 0 7 1 0 0 0 8 1 0 0 0 9 1 1 0 0 10 1 1 0 0 11 1 1 0 0 12 1 1 0 0 13 1 1 0 0 14 1 1 0 0 15 1 1 0 0 16 1 1 0 0 . . . c Rost 2007 °
Ovˇ eˇ ren´ı pˇ redpoklad˚ u homoskedasticity K ovˇ eˇ ren´ı pˇ redpokladu homoskedasticity m˚ uˇ zeme vyuˇ z´ıt napˇ r´ıklad Bartlett˚ uv test. bartlett.test(koagulace~dieta,data) Bartlett test of homogeneity of variances data: koagulace by dieta Bartlett’s K-squared = 4.1045, df = 3, p-value = 0.2504
c Rost 2007 °
Dalˇ s´ı diagnostika Lze vyuˇ z´ıt rezidua k posouzen´ı normality qqnorm(residuals(model)) plot(jitter(fitted(model)), residuals(model), xlab="Teoretick´ e z´ ıskan´ e modelem", ylab="Residua")
c Rost 2007 °
Diagnostick´ e grafy
6 4 2 −6
−4
−2
0
Residua
2 0 −2 −4 −6
Sample Quantiles
4
6
Normal Q−Q Plot
−2
−1
0 Theoretical Quantiles
1
2
62
63
64
65
66
Teoretické získané modelem
c Rost 2007 °
Dvoufaktorov´ a anal´ yza rozptylu
c Rost 2007 °
Dvoufaktorov´ a anal´ yza rozptylu V nˇ ekter´ ych pˇ r´ıpadech je nutn´ e uvaˇ zovat vliv souˇ casn´ eho p˚ usoben´ı dvou faktor˚ u. Situace se tak m´ırnˇ e komplikuje. Hovoˇ r´ıme pak o tzv. anal´ yze rozptylu pˇ ri dovojn´ em tˇ r´ıdˇ en´ı, neboli o dvoufaktorov´ e anal´ yze rozptylu. Uvaˇ zujme dva faktory: Faktor A necht’ m´ a I ´ urovn´ı i = 1, 2, · · · , I). Faktor B necht’ m´ a J ´ urovn´ı j = 1, 2, · · · , J). D´ ale uvaˇ zujme pouze pˇ r´ıpady, kdy m´ ame tzv. vyv´ aˇ zen´ e tˇ r´ıdˇ en´ı, tj. kdy vˇ sechny ˇ cetnosti nij jsou stejn´ e a jsou rovny nˇ ejak´ emu ˇ c´ıslu.
c Rost 2007 °
ˇ ed´ S a teorie Matematick´ y model lze (obecnˇ e pro tento typ ´ uloh) formulovat takto: Yijp = µ + αi + βj + ²ijp
I X
J X
αi = 0
i
βj = 0
j
Potˇ rebn´ e v´ ypoˇ cty: I 1 X 2 − 1Y 2 SA = Yi•• JP i=1 n •••
ST =
J 1 X 2 − 1Y 2 SB = Y•j• IP j=1 n •••
I X J X P X i
j
p
2 − Yijp
1 2 Y••• n c Rost 2007 °
ˇ ed´ S a teorie Se = ST − SA − SB
c Rost 2007 °
Troˇ sku to uspoˇ r´ ad´ ame V´ ypoˇ cty lze opˇ et uspoˇ r´ adat do pˇ rehledn´ e tabulky anal´ yzy rozptylu: Zdroj variability
Souˇ cet ˇ ctverc˚ u
Stupnˇ e volnosti
Pr˚ umˇ ern´ y souˇ cet ˇ ctverc˚ u
Testov´ a statistika
p-value
Faktor A
SA
dfA = I − 1
M SA =
SA dfA
FA =
M SA M Se
1 − F (FA |dfA ; dfe )
Faktor B
SB
dfB = J − 1
M SB =
SB dfB
FB =
M SB M Se
1 − F (FB |dfB ; dfe )
Residuum
Se
dfe = IJP − I − J + 1
M Se =
Celkov´ y
ST
dfT = IP J − 1
Se dfe
c Rost 2007 °
Pˇ r´ıklad Bylo sledov´ ano, zda ˇ cas potˇ rebn´ y k v´ yrobˇ e urˇ cit´ e souˇ c´ astky z´ avis´ı na denn´ı dobˇ e a na hluˇ cnosti v okol´ı. Experimentem, jemuˇ z se podrobilo dvan´ act pracovn´ık˚ u firmy byly zjiˇ stˇ eny n´ asleduj´ıc´ı hodnoty [v min]:
R´ ano V poledne Veˇ cer
Ticho 6 8 7
Hudba 7 5 6
Rozhlas 8 10 12
Hluk 6 5 7
Je doba v´ yroby z´ avisl´ a na denn´ı dobˇ e a na hluˇ cnosti okol´ı? Testovan´ e hypot´ ezy: H0A : α1 = α2 = α3 = 0 H0B : β1 = β2 = β3 = β4 = 0 c Rost 2007 °
Jak to sp´ ach´ ame v R
cas<-c("rano","poledne","vecer") prostredi<-c("ticho","hudba","hra","hluk") odezva<-c(6,8,7,7,5,6,8,10,12,6,5,7) mat<-expand.grid(cas=cas,prostredi=prostredi) data<-data.frame(odezva,mat) data odezva cas prostredi 1 6 rano ticho 2 8 poledne ticho 3 7 vecer ticho 4 7 rano hudba 5 5 poledne hudba 6 6 vecer hudba 7 8 rano hra 8 10 poledne hra 9 12 vecer hra 10 6 rano hluk 11 5 poledne hluk 12 7 vecer hluk c Rost 2007 °
pokraˇ cov´ an´ı v R
model<-lm(odezva~cas+prostredi,data=data) anova(model) Analysis of Variance Table Response: odezva Df Sum Sq Mean Sq F value Pr(>F) cas 2 3.50 1.75 1.0000 0.42188 prostredi 3 32.25 10.75 6.1429 0.02926 * Residuals 6 10.50 1.75 --Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1
Z v´ ysledk˚ u je zˇ rejm´ e, ˇ ze vliv denn´ı doby na ˇ cas potˇ rebn´ y k vyroben´ı souˇ c´ astky nebyl na hladinˇ e α = 0,05 na z´ akladˇ e pozorovan´ ych dat prok´ az´ an. Naopak lze zam´ıtnout hypot´ ezu, ·ˇ ze hluˇ cnost v okol´ı neovlivˇ nuje dobu potˇ rebnou k v´ yrobˇ e. c Rost 2007 °
Odhad efekt˚ u
c Rost 2007 °
Odhad efekt˚ u summary(model) Residuals: Min 1Q Median -1.5000 -0.7500 -0.1250
3Q 0.6875
Max 1.5000
Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 6.5000 0.9354 6.949 0.000441 *** caspoledne 0.2500 0.9354 0.267 0.798217 casvecer 1.2500 0.9354 1.336 0.229893 prostredihudba -1.0000 1.0801 -0.926 0.390259 prostredihra 3.0000 1.0801 2.777 0.032104 * prostredihluk -1.0000 1.0801 -0.926 0.390259 --Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1 Residual standard error: 1.323 on 6 degrees of freedom Multiple R-Squared: 0.773, Adjusted R-squared: 0.5838 F-statistic: 4.086 on 5 and 6 DF, p-value: 0.0581
c Rost 2007 °
Motivaˇ cn´ı pˇ r´ıklad - dvoufaktorov´ a ANOVA s interakcemi Byl proveden pokus, pˇ ri nˇ emˇ z byly sledov´ any v´ ynosy jist´ e zemˇ edˇ elsk´ e plodiny v t.ha−1 v z´ avislosti na typu p˚ udy (kysel´ a p˚ uda, norm´ aln´ı p˚ uda) a typu hnojen´ı (kontrola bez hnojen´ı, chl´ evsk´ a mrva, Ca hnojivo). Kaˇ zd´ a kombinace byla realizov´ ana dvakr´ at. Data:
Typ Norm´ aln´ı p˚ uda Kysel´ a p˚ uda
Kontrola 2,9 3,1 2,7 3,0
Mrva 3,6 3,9 3,5 3,4
Hnojen´ı Ca 3,2 3,5 3,8 4,2
M´ a druh p˚ udy ˇ ci typ hnojen´ı vliv na v´ ynos?
c Rost 2007 °
N´ aˇ s model Matematick´ y model lze (obecnˇ e pro tento typ ´ uloh) formulovat takto: Yijp = µ + αi + βj + (αβ)ij + ²ijp Naˇ se pracovn´ı (testovan´ a - nulov´ a) hypot´ ezy pˇ redpokl´ adaj´ı, ˇ ze sledovan´ e faktory, pˇ r´ıpadnˇ e jejich interakce, nemaj´ı na tvorbu v´ ynosu vliv, tj.: H0A : α1 = α2 = · · · = αI ; H 0 B : β 1 = β2 = · · · = βJ ; H0AB : (αβ)11 = (αβ)12 = · · · = (αβ)IJ . Pozn.: Pˇ ri samotn´ em ˇ reˇ sen´ı vyuˇ zity reparametrizaˇ cn´ı podm´ınky: I X i
αi = 0;
J X j
βj = 0;
I X i
(αβ)ij = 0 ∀ i;
J X
(αβ)ij = 0 ∀ j
j
c Rost 2007 °
Jak to zad´ an´ı dostaneme do Erka? #Nejprve vytvoˇ rı ´me sledovan´ e ´ urovnˇ e obou faktor˚ u a uspoˇ ra ´d´ ame je. #Uspoˇ r´ ad´ an´ ı nen´ ı nezbytnˇ e nutn´ e, ale pro jistotu. puda<-rep(c("normalni","kysela"),6) puda<-factor(puda,levels=c("normalni","kysela")) osetreni<-rep(c("kontrola","mrva","Ca"),c(4,4,4)) osetreni<-factor(osetreni,levels=c("kontrola","mrva","Ca")) #Vytvoˇ r´ ıme vektor obsahuj´ ıci jednotliv´ e hodnoty v´ ynos˚ u #pro dan´ e ´ urovnˇ e y<-c(2.9,2.7,3.1,3.0,3.6,3.5,3.9,3.4,3.2,3.8,3.5,4.2) #Vˇ se to uloˇ zı ´me do objektu data data<-data.frame(puda,osetreni,y) #Koukneme se na to, co jsme vytvoˇ rili data
c Rost 2007 °
Naˇ se data data puda 1 normalni 2 kysela 3 normalni 4 kysela 5 normalni 6 kysela 7 normalni 8 kysela 9 normalni 10 kysela 11 normalni 12 kysela
osetreni kontrola kontrola kontrola kontrola mrva mrva mrva mrva Ca Ca Ca Ca
y 2.9 2.7 3.1 3.0 3.6 3.5 3.9 3.4 3.2 3.8 3.5 4.2
trellis.device(theme="col.whitebg") xyplot(y~osetreni|puda,data=data,ylab="vynos v t/ha",pch=20,cex=1.2) xyplot(y~puda|osetreni,data=data,ylab="vynos v t/ha",pch=20,cex=1.2)
c Rost 2007 °
Vizualizujme? Ano! kontrola
normalni
mrva
Ca
kysela
Ca ● ●
4.0 ●
4.0
3.5
●
● ●
●
3.5
●
● ●
vynos v t/ha
vynos v t/ha
3.0 ●
kontrola
mrva
4.0
●
● ●
●
3.0
●
3.5
● ●
●
●
3.0
● ●
●
●
kontrola
mrva
normalni
Ca
osetreni
kysela
puda
c Rost 2007 °
Ted’ trochu v´ aˇ znˇ e - malinko to zamlˇ z´ıme - ˇ sed´ a, ˇ sed´ a je teorie Y••• =
I X J X P X i
j
Yijp
Yi•• =
p
j
I 1 X 2 1 2 SA = Yi•• − Y••• JP i=1 n
Se =
J X P I X X i
j
p
2 Yijp
J X P X
I J 1 XX 2 − Y P i j ij•
Yijp
Yij• =
p
P X
Yijp
p
J 1 X 2 1 2 SB = Y•j• − Y••• IP j=1 n
ST =
I X J X P X i
j
p
2 Yijp −
1 2 Y n •••
SAB = ST − SA − SB − Se
c Rost 2007 °
Troˇ sku to uspoˇ r´ ad´ ame V´ ypoˇ cty lze uspoˇ r´ adat do pˇ rehledn´ e tabulky anal´ yzy rozptylu: Zdroj variability
Souˇ cet ˇ ctverc˚ u
Stupnˇ e volnosti
Pr˚ umˇ ern´ y souˇ cet ˇ ctverc˚ u
Faktor A
SA
dfA = I − 1
M SA =
SA dfA
FA =
M SA M Se
1 − F (FA |dfA ; dfe )
Faktor B
SB
dfB = J − 1
M SB =
SB dfB
FB =
M SB M Se
1 − F (FB |dfB ; dfe )
Interakce
SAB
dfAB = IJ − I − J + 1
M SAB =
Residuum
Se
dfe = IJP − IJ
M Se =
Celkov´ y
ST
dfT = IP J − 1
SA B dfAB
Testov´ a statistika
FAB =
M SAB M Se
p-value
1 − F (FAB |dfAB ; dfe )
Se dfe
c Rost 2007 °
Katarze pˇ rich´ az´ı aneb zlat´ a praxe ;-) V´ ypoˇ cet provedeme v Erku. Syntaxe je velmi jednoduch´ a: model<-aov(y~puda*osetreni,data=data) summary(model) Df Sum Sq Mean Sq F value Pr(>F) puda 1 0.01333 0.01333 0.3333 0.584700 osetreni 2 1.36500 0.68250 17.0625 0.003344 ** puda:osetreni 2 0.52167 0.26083 6.5208 0.031285 * Residuals 6 0.24000 0.04000 --Signif. codes: 0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1 V´ ysledky lze struˇ cnˇ e shrnout. Vzhledem k hodnot´ am p-value lze ˇ r´ıci, ˇ ze se na z´ akladˇ e analyzovan´ ych dat nepodaˇ ril prok´ azat vliv p˚ udn´ıho typu na v´ ynos hospod´ aˇ rsk´ e plodiny (p-value 0, 5847). Naopak vliv oˇ setˇ ren´ı je statisticky v´ yznamn´ y p-value 0, 003344. Tot´ eˇ z plat´ı i pro interakci mezi p˚ udn´ım typem a oˇ setˇ ren´ım p-value 0, 031285. Pod´ıvejme se jeˇ stˇ e na graf interakc´ı. M˚ uˇ ze pomoci pˇ ri interpretaci.
c Rost 2007 °
4.0
Opˇ et graf - interakce
osetreni
3.4 3.2 3.0
výnos v t/ha
3.6
3.8
Ca mrva kontrola
normalni
kysela Pùdni typ
c Rost 2007 °
Zjist´ıme jeˇ stˇ e statisticky v´ yznamn´ e rozd´ıly mezi ´ urovnˇ emi faktor˚ u Pro faktor oˇ setˇ ren´ı (kontrola, chl´ evsk´ a mrva, v´ apnˇ en´ı): TukeyHSD(model,which="osetreni") Tukey multiple comparisons of means 95% family-wise confidence level Fit: aov(formula = y ~ puda * osetreni, data = data) $osetreni diff lwr upr mrva-kontrola 0.675 0.2410805 1.1089195 Ca-kontrola 0.750 0.3160805 1.1839195 Ca-mrva 0.075 -0.3589195 0.5089195
c Rost 2007 °
Tukey HSD test - grafy Graficky to m˚ uˇ zeme zn´ azornit napˇ r´ıklad takto: 95% family−wise confidence level
95% family−wise confidence level
kysela:kontrola−normalni:kontrola normalni:mrva−normalni:kontrola kysela:mrva−normalni:kontrola kysela−normalni normalni:Ca−normalni:kontrola kysela:Ca−normalni:kontrola normalni:mrva−kysela:kontrola −0.2
−0.1
0.0
0.1
0.2
0.3
kysela:mrva−kysela:kontrola
Differences in mean levels of puda normalni:Ca−kysela:kontrola
95% family−wise confidence level mrva−kontrola
kysela:Ca−kysela:kontrola kysela:mrva−normalni:mrva normalni:Ca−normalni:mrva kysela:Ca−normalni:mrva
Ca−kontrola normalni:Ca−kysela:mrva kysela:Ca−kysela:mrva kysela:Ca−normalni:Ca Ca−mrva 0.0
0.5
1.0
Differences in mean levels of osetreni
−1.0
−0.5
0.0
0.5
1.0
1.5
2.0
Differences in mean levels of puda:osetreni
c Rost 2007 °
Metoda zn´ ahodnˇ en´ ych blok˚ u Jedn´ a se o dosti ˇ cast´ y pˇ r´ıpad v poln´ıch pokusech. Princip metody je n´ asleduj´ıc´ı: Pozemek rozdˇ el´ıme na tolik blok˚ u, kolik m´ ame opakov´ an´ı. Uvnitˇ r bloku provedeme n´ aslednˇ e dˇ elen´ı do tolika parcel, kolik m´ a sledovan´ y faktor ´ urovn´ı. Na takto vznikl´ e parcely pˇ riˇ rad´ıme v r´ amci blok˚ u n´ ahodnˇ e jednotliv´ e ´ urovnˇ e faktoru. V´ yhodou metody je to, ˇ ze umoˇ zˇ nuje eliminovat vliv rozd´ılnost´ı mezi jednotliv´ ymi bloky. Mˇ ejme 8 odr˚ ud ovsa a kaˇ zdou z odr˚ ud budeme cht´ıt vys´ıt na 5 pokusn´ ych pozemc´ıch stejn´ e velikosti. Jak provedeme zn´ ahodnˇ en´ı? Oznaˇ cme jednotliv´ e odr˚ udy jako: a1; b2; c3; d4; e5; f6; g7; h8.
c Rost 2007 °
Baˇ sta pro Erko odrudy<-c("a1","b2","c3","d4","e5","f6","g7","h8") poradi<-c(sample(odrudy,rep=F),sample(odrudy,rep=F),sample(odrudy,rep=F), sample(odrudy,rep=F),sample(odrudy,rep=F)) barvy<-factor(poradi) barvy<-as.numeric(barvy) barva<-matrix(barvy,byrow=TRUE,5,8) plan<-matrix(poradi,byrow=TRUE,5,8) x<-1:5 y<-1:8 sit<-expand.grid(x,y) sit plot(sit$Var1,sit$Var2,type="n",xlab="bloky - opakovani", ylab="parcely - odrudy") text(sit$Var1,sit$Var2,plan,col=barva,cex=1.52) plan
c Rost 2007 °
8
e5
c3
c3
e5
h8
7
g7
a1
a1
c3
g7
6
h8
b2
e5
f6
a1
5
b2
d4
b2
a1
f6
4
f6
e5
f6
d4
e5
3
a1
h8
g7
g7
b2
2
d4
g7
h8
h8
c3
1
parcely − odrudy
Jak to tedy vlastnˇ e vypad´ a?
c3
f6
d4
b2
d4
1
2
3
4
5
bloky − opakovani
M´ ame pl´ an pokus˚ u. Ted’ uˇ z zb´ yv´ a jen prov´ est samotn´ y pokus ... c Rost 2007 °
Data z´ıskan´ e proveden´ım pokusu - jiˇ z uspoˇ r´ adan´ e ´ daje jsou v gramech. U Bloky Odr˚ uda 1 2 3 4 a1 296 357 340 331 b2 202 390 431 340 c3 437 334 426 320 d4 303 319 310 260 e5 469 405 442 487 f6 345 342 358 300 g7 324 339 357 352 h8 488 374 401 338
5 348 320 296 242 394 308 220 320
Naˇ cteme data do Erka a provedeme ANOVU. V pˇ r´ıpadˇ e metody zn´ ahodnˇ en´ ych blok˚ u postupujeme jako v pˇ r´ıpadˇ e dvoufaktorov´ e anal´ yzy rozptylu bez interakc´ı.
c Rost 2007 °
ˇ ed´ S a teorie Matematick´ y model lze (obecnˇ e pro tento typ ´ uloh) formulovat takto: Yijp = µ + αi + βj + ²ijp
I X
J X
αi = 0
i
βj = 0
j
Potˇ rebn´ e v´ ypoˇ cty: I 1 X 2 1 2 SA = Yi•• − Y••• JP i=1 n
ST =
J 1 X 2 1 2 SB = Y•j• − Y••• IP j=1 n
I X J X P X i
j
p
2 Yijp −
1 2 Y n •••
Se = ST − SA − SB c Rost 2007 °
Troˇ sku to uspoˇ r´ ad´ ame V´ ypoˇ cty lze opˇ et uspoˇ r´ adat do pˇ rehledn´ e tabulky anal´ yzy rozptylu: Zdroj variability
Souˇ cet ˇ ctverc˚ u
Stupnˇ e volnosti
Pr˚ umˇ ern´ y souˇ cet ˇ ctverc˚ u
Testov´ a statistika
p-value
ˇ´ R adky
SA
dfA = I − 1
M SA =
SA dfA
FA =
M SA M Se
1 − F (FA |dfA ; dfe )
Sloupce
SB
dfB = J − 1
M SB =
SB dfB
FB =
M SB M Se
1 − F (FB |dfB ; dfe )
Residuum
Se
dfe = IJP − I − J + 1
M Se =
Celkov´ y
ST
dfT = IP J − 1
Se dfe
c Rost 2007 °
Kanonenfutter . . . bloky<-factor(1:5) odrudy<-c("a1","b2","c3","d4","e5","f6","g7","h8") odr<-odrudy dat<-expand.grid(blok=bloky,odruda=odr) y<-c(296 , 357 , 340 , 331 , 348, 202 , 390 , 431 , 340 , 320, 437 , 334 , 426 , 320 , 296, 303 , 319 , 310 , 260 , 242, 469 , 405 , 442 , 487 , 394, 345 , 342 , 358 , 300 , 308, 324 , 339 , 357 , 352 , 220, 488 , 374 , 401 , 338 , 320) data<-data.frame(blok=dat$blok,odruda=dat$odruda,y) model2<-aov(y~odruda+blok,data=dat) summary(model2)
c Rost 2007 °
V´ ysledek proveden´ e anal´ yzy Z´ısk´ ame v´ yslednou tabulku anal´ yzy rozptylu. Df Sum Sq Mean Sq F value Pr(>F) 7 75534 10791 4.5214 0.001773 ** 4 25845 6461 2.7074 0.050411 . 28 66823 2387
odruda blok Residuals --Signif. codes:
0 ’***’ 0.001 ’**’ 0.01 ’*’ 0.05 ’.’ 0.1 ’ ’ 1
Nulovou hypot´ ezu o shodˇ e stˇ redn´ıch hodnot pro jednotliv´ e odr˚ udy m˚ uˇ zeme tedy zam´ıtnout, nebot’ p-value = 0,001773. Ovˇ eˇ rme jeˇ stˇ e hypot´ ezu o shodˇ e rozptyl˚ u mezi jednotliv´ ymi odr˚ udami. Vyuˇ zijme pro tento ´ uˇ cel Bartlett˚ uv test. bartlett.test(y,dat$odruda) Bartlett test for homogeneity of variances data: y and dat$odruda Bartlett’s K-squared = 10.5244, df = 7, p-value = 0.1608 c Rost 2007 °
