Cvičení z NSTP097 11. 1. 2010 Dvouvýběrové testy, analýza rozptylu Ze stránky www.karlin.mff.cuni.cz/~omelka/Vyuka_stp097.php si stáhněte soubor s daty cviceni4.RData. Otevřete si program R a pomocí File/Load_Workspace... si načtěte do R tento soubor. Příkazem ls() si vyzkoušejte, zda se načtení dat povedlo. Mezi vypsanými objekty byste měli vidět předdefinovanou funkci skewness a datový soubor plat. Jako obvykle doporučuji si psát použité příkazy do zvláštního textového souboru např. pomocí editoru zabudovaného ve windowsovském rozhraní R. Datový soubor Datový soubor plat obsahuje údaje z výběrového šetření 534 zaměstnaných osob žijících na různých místech v USA v roce 1985 (Berndt, ER. The Practice of Econometrics. 1991. NY: Addison-Wesley). Každý řádek obsahuje údaje pro jednoho respondenta šetření. Budeme vyšetřovat závislost výše mzdy na některých vybraných faktorech. Veličiny a jejich kódování Vzdělání (počet let ve škole) Geografická poloha (Jih/jinde) Pohlaví (Žena/Muž) Délka praxe (roky) Členství v odborech (Člen/Nečlen) Hodinová mzda v US $ Věk (roky) Rasa (Běloch/Hispanik/Jiná) Pracovní zařazení (Vedoucí prac./Obchod/Administrativa/Služby/Odborný prac./Ostatní) sektor Sektor (Stavebnictví/Průmysl/Ostatní) stav Rodinný stav (Ženatý–vdaná/Ostatní) vzdel jih pohl leta odbory mzda vek rasa zaraz
Zadání 1. Prostudujte si datový soubor plat dle následujících instrukcí. (a) Vypište si jména veličin příkazem names(plat) . (b) Prohlédněte si prvních dvacet pozorování plat[1:20,] . Tímto říkáme, že chceme vypsat řádky 1 až 20 a všechny sloupce. Obecně, chceme-li vypsat řádek i a sloupec j dat, napíšeme plat[i,j].
(c) Vypište si základní charakteristiky jednotlivých veličin souboru příkazem summary(plat) . Jak se liší výstupy pro veličiny mzda a rasa? Proč pro ně nedostaneme totéž? Zajistěte si přímý přístup k jednotlivým sloupcům (veličinám) attach(plat) . (d) Zobrazte histogram mezd hist(mzda) . Co byste řekli o rozdělení mezd? 2. Nejprve se budeme zabývat otázkou, zda jsou na Jihu v průměru stejné platy jako jinde ve Spojených státech. Ve značení používaném na přednášce máme dva nezávislé výběry X1 , . . . , Xn z nějaké distribuční funkce FX (platy na Jihu) a Y1 , . . . , Ym z obecně různé distribuční funkce FY (platy jinde). Zajímá nás rozdíl mezi R ∞distribučními funkcemi FX a FY , zejména pak rozdíl jejich středních hodnot EXi = −∞ x dFX (x) a R∞ EYj = −∞ y dFY (y). Jelikož hodinové mzdy nabývají mnoha kladných hodnot, můžeme rozdělení FX a FY považovat za spojitá s nosičem (0, ∞). (a) Zjistěte, čemu jsou rovny rozsahy výběru n a m, například příkazem table(jih) . (b) Zjistěte výběrové momenty rozdělení FX a FY : odhady středních hodnot pomocí tapply(mzda,jih,mean) , odhady rozptylů pak dostaneme stejným způsobem, nahradíme-li v předchozím příkazu mean za var. Pro odhady šikmosti napíšeme skewness namísto mean. Příkaz tapply(mzda,jih,mean) říká, že na vektor mzda chceme aplikovat funkci mean zvlášť pro jednotlivé kategorie faktorové veličiny jih. Tj. spočítáme např. průměr pro kategorii Jih a jinde zvlášť. Pomocí vypočtených statistik získáme hrubou představu o tom, zda se FX a FY liší střední hodnotou či rozptylem a zda jsou obě rozdělení přibližně symetrická. Připomeňte si definici šikmosti. Čemu je rovna šikmost symetrického rozdělení? (c) Porovnejte graficky rozdělení obou výběrů pomocí krabicových diagramů, tzv. box plotů) boxplot(mzda~jih) . Krabicový diagram mimo jiné znázorňuje vybrané výběrové kvantily: obvykle medián (tlustá čára uvnitř krabice) a kvartily (okraje krabice). Stejný obrázek bychom dostali, jestliže zadáme standardní plot(mzda~jih). Ověřte! (d) Porovnejte histogramy obou výběrů: par(mfrow=c(2,1)) hist(mzda[jih=="Jih"],breaks=5*(0:10),xlim=range(mzda)) hist(mzda[jih=="jinde"],breaks=5*(0:10),xlim=range(mzda)) par(mfrow=c(1,1))
Volba par(mfrow=c(2,1)) zajistí, že se histogramy nakreslí do jednoho panelu pod sebe. Chceme-li je mít vedle sebe, zadáme par(mfrow=c(1,2)). Volby xlim a breaks zajišťují porovnatelnost obou histogramů co do rozmezí osy x a šířky a rozmístění sloupců. Nakonec vrátíme grafiku do původního stavu par(mfrow=c(1,1)). (e) Nyní spočtěte empirické distribuční funkce FˆX,n a FˆY,m a nakreslete jejich grafy: a <- ecdf(mzda[jih=="Jih"])
2
b <- ecdf(mzda[jih=="jinde"]) oddo <- range(mzda)*c(0.9,1.1) par(col="blue") plot(a,xlim=oddo,main="Empiricka distribucni funkce mzdy",cex=0.3) par(col="black") lines(b,cex=0.3,lty=2) legend(30,0.4,lty=c(1,1),col=c("blue","black"),legend=c("Jih","Jinde"))
Do a, b si uložíme zmíněné empirické distribuční funkce. Do vektoru oddo si nastavíme rozsah na ose x pro kreslení obrázků. Příkaz par(col="barva") nastaví barvu, kterou se kreslí body a čáry, na barva. Děláme to kvůli lepšímu rozlišení funkcí a, b. Volba cex=0.3 zmenšuje zakreslované body. Příkaz legend přidá do obrázku legendu do místa s x-ovou souřadnicí 30 a y-ovou souřadnicí 0.4. Volba legend=c("Jih","Jinde") specifikuje popisky, zbytek potom formát čar v legendě, tak aby to odpovídalo obrázku. Co prozrazují empirické distribuční funkce o platech? Rozmyslete si následující věc: Jestliže pro náhodné veličiny U ∼ FU a V ∼ FV platí FU (x) ≤ FV (x) pro všechna x, která z veličin U, V nabývá „větších hodnotÿ? 3. Proveďte formální test hypotézy H0 : FX = FY proti alternativě H1 : ∃t : FX (t) 6= FY (t), dvouvýběrovým Kolmogorovovým-Smirnovovým testem: ks.test(mzda[jih=="Jih"],mzda[jih=="jinde"]) . Zamítáme na základě výsledku tohoto testu hypotézu o rovnosti distribucí? Poznámka: Všimněte si, že R píše Warning message: In ks.test(mzda[jih == "Jih"], mzda[jih == "jinde"]) : cannot compute correct p-values with ties
Funkce ks.test si stěžuje, že nelze spočítat p-hodnotu kvůli shodným pozorováním (porušení předpokladu o spojitém rozdělení). Příkazem sum(table(mzda)>1) zjistíme, kolik je ve vektoru mzda shodných pozorování (tj. pozorování, jejichž hodnota se vyskytuje v datech více než jednou). Kolik to je? Nyní tento problém ignorujme. 4. Nyní proveďte klasický dvouvýběrový t-test (za předpokladu shodnosti rozptylů). Připomeňte si předpoklady dvouvýběrového t-testu. Jak obecně interpretujeme výsledky ttestů při porušení normality? Která věta to zdůvodňuje? t.test(x=mzda[jih=="Jih"],y=mzda[jih=="jinde"],var.equal=T) ## nebo ekvivalentně, ale o něco jednodušeji t.test(mzda~jih, var.equal=T)
Volba var.equal=T nastavuje verzi t-testu, která předpokládá shodné rozptyly u obou skupin. Jakou hypotézu tento test testuje proti jaké alternativě? Zamítáte na základě výsledku testu nulovou hypotézu? Co je testová statistika? Jaké je za nulové hypotézy její rozdělení? 5. Vrátíme-li se k výsledkům bodu 2 (b). Můžeme si být jisti, že předpoklad rovnosti rozptylů není porušen? Co kdybychom provedli asymptotický z-test, který tento předpoklad nemá? 3
Spusťme t.test(mzda~jih, var.equal=F)
Jestliže napíšeme var.equal=F nebo nespecifikujeme nic, pak R provádí test, který nepředpokládá shodu rozptylů. Tento test se nazývá Welchův dvouvýběrový t-test, ale nejedná se o nic jiného než o testovou statistiku z-testu (viz přednáška) s kritickými hodnotami spočtenými z t-rozdělení s aproximovaným počtem stupňů volnosti. (Všimněte si hodnoty stupňů volnosti df ve výstupu.) Tato statistika nemá t-rozdělení, jedná se vskutku pouze o aproximaci. Pro velké m a n je to totéž jako použít limitní normální rozdělení. Liší se testová statistika a p-hodnota od klasického t-testu? Liší se rozhodnutí o zamítnutí H0 ? 6. Proveďme ještě dvouvýběrový Wilcoxonův test: wilcox.test(mzda[jih=="Jih"],mzda[jih=="jinde"]). (a) Na čem je založen Wilcoxonův test? Podívejte se ve výstupu na p-hodnotu. Zamítá se nulová hypotéza o tom, že obě srovnávaná rozdělení jsou totožná? Poznámka: To, že zkoušíme řadu různých testů stejné nebo podobné hypotézy na jedněch datech, má pouze pedagogické důvody. Chceme si to vyzkoušet. V praxi bychom si měli předem vybrat jeden určitý test a ten provést! Pokud však hledáme tak dlouho nějaký test, až konečně nulovou hypotézu zamítneme, pak se naše analýza podobá útrpnému právu. Data se nakonec téměř vždycky k něčemu přiznaji (zamítneme nulovou hypotézu), i kdyby byla nevinná (i kdyby nulová hypotéza platila). (b) Připomeňme, že Wilcoxonova testová statistika byla na přednášce definována jako Wn,m =
n X
Ri ,
i=1
kde Ri jsou pořadí pozorování náhodného výběru X1 , . . . , Xn mezi všemi X1 , . . . , Xn , Y1 , . . . , Ym . Testová statistika ve výstupu z R není přímo Wn,m , ale Wn,m − n(n + 1)/2. Ověřte to! Návod: Statistiku Wn,m spočítáte jako sum(rank(mzda)[jih=="Jih"]) Příkaz rank(mzda) nám dá pořadí pozorování ve vektoru mzda. Vezme z nich pouze ta, co odpovídají X’ům, tj. datům z jihu, a ty sečteme funkcí sum. 7. Problém různých rozptylů u dvouvýběrového t-testu lze řešit i transformací zkoumané veličiny nějakou funkcí. Opakujte body 2 (b), 4 a 5, kde místo mzdy budete analyzovat logaritmus mzdy o základu 10, tj. všude místo mzda pište log10(mzda). Zlepšila se rovnost rozptylů? Liší se po transformaci výsledky t-testu a z-testu více nebo méně než před transformací? Rozmyslete si: Jsou hypotézy testované t-testem a z-testem po transformaci ekvivalentní hypotézám testovaným na netransformovaných datech? Pokud ne obecně, za jakých předpokladů by byly? 8. Má smysl opakovat body 3 a 6 pro transformovaná data? Co by vyšlo? Proč? 9. Nyní přistoupíme k porovnání středních hodnot několika výběrů pomocí analýzy rozptylu. Připomeňte si, jaký problém řeší analýza rozptylu a jaké jsou předpoklady této metody. 4
Podívejme se třeba na vztah mezi rasou a mzdou. Rasa je veličina, která rozděluje pozorování na tři skupiny: Běloch/Hispanik/Jiná. Zajímá nás, zda mají tyto tři skupiny stejnou střední mzdu. (a) Jak velké (kolik pozorování) jsou zmíněné tři skupiny v našich datech? table(rasa) Znázorněte graficky pozorovaná rozdělení mezd v těchto třech skupinách: nejdříve pomocí box plotů boxplot(mzda~rasa) a pak pomocí histogramů par(mfrow=c(3,1)) hist(mzda[rasa=="Beloch"],breaks=5*(0:10),xlim=range(mzda)) hist(mzda[rasa=="Hispanik"],breaks=5*(0:10),xlim=range(mzda)) hist(mzda[rasa=="Jina"],breaks=5*(0:10),xlim=range(mzda)) par(mfrow=c(1,1))
(b) Porovnejte aritmetické průměry a odhady rozptylu mzdy ve třech rasových skupinách tapply(mzda,rasa,mean) , tapply(mzda,rasa,var) . Vypadají průměrné mzdy podobně? Proč je důležité porovnávat odhadnuté rozptyly? (c) Proveďme analýzu rozptylu. K tomu v R slouží funkce aov (Analysis Of Variance), jejíž výsledek však musíme ještě zpracovat funkcí summary. Takže zadejte summary(aov(mzda~rasa)) Na levé straně operátoru ~ je veličina obsahující měření Yij (viz značení na přednášce), na pravé straně je veličina definující skupiny. Rozhodněte, zda zamítáte či nezamítáte hypotézu o rovnosti středních hodnot. Čemu je rovna testová statistika a jaké je její rozdělení za nulové hypotézy? (d) Porovnejte výstup z funkce aov s tabulkou analýzy rozptylu z přednášky: Zdroj měnlivosti
Součet Stupně Podíl čtverců volnosti
Skupina
SSA
k−1
SSA k−1
Residuální
SSe
n−k
SSe n−k
Celkový
SSC
n−1
F SSA SSe / k−1 n−k
Poznáte kde co je, popř. co chybí a co přebývá? Dopočítejte celkový součet čtverců SSC =
ni k X X
(Yij − Y ·· )2
i=1 j=1
jednak pomocí výstupu funkce aov a jednak přímo z uvedené definice. Návod: sum((mzda-mean(mzda))^2) . (e) Hodnotu distribuční funkce F (x) v bodě x rozdělení Fn,m (F -rozdělení s n a m stupni volnosti) můžeme spočítat pomocí příkazu pf(x,n,m) a-kvantil téhož rozdělení pak dostaneme pomocí 5
qf(a,n,m) . Pomocí těchto dvou funkcí spočítejte kritickou hodnotu pro zamítání H0 u právě provedeného testu a ověřte p-hodnotu. Srovnejte s výstupem z R. 10. Opakujte body (b) a (c) po zlogaritmování mezd transformací log10(mzda). Změnil se výsledek testu? Která analýza je „spolehlivějšíÿ, ta s původními mzdami nebo po zlogaritmování? 11. Prozkoumejte, zda mzda souvisí s pracovním zařazením a se sektorem, v němž je respondent zaměstnán. Opakujte pouze body (a), (b), (c) a 10. 12. Na závěr je dobrou zvyklostí „propustitÿ náš datový soubor pomocí detach(plat) .
6
