################################################################## ######### Testovani hypotez: Realna dat ################################################################## # # A hura na realna data! :-)) # # Soubor vysky.txt obsahuje udaje o vysce (cm) a pohlavi (1=devce, 2=chlapec) 27 nahodne vybranych stejne starych deti # Ukoly: 1. Nactete data do nejake promenne # Hint: napriklad pomoci prikazu vys <- read.table(...), prislusne parametry (...) doplnte vyuzitim help(read.table) # 2. Prohlednete si data (vypis nekolika hodnot, summary, obrazky - histogramy, krabice...) # Hint 1: k promennym uvnitr 'vys' se pristupuje primo po zavolani attach(vys) nebo bez volani attach() pomoci znaku $ - napriklad vys$vyska, # pripadne pomoci indexovani [], vypis promennych dostaneme pomoci names(vys) # Hint 2: pro zapsani zavislosti lze pouzit specialni operator '~', napr. boxplot(vys$vyska ~ vys$pohlavi) vykresli krabice pro vysku v zavislosti na pohlavi ditete # 3. Otestujte hypotezu H_0, ze prumerna vyska divek je rovna 143cm, proti oboustranne alternative # 4. Otestujte hypotezu H_0, ze prumerna vyska divek je stejna jako prumerna vyska chlapcu, proti oboustranne alternative # Hint: Jaky test pouzijete? Jsou splneny jeho predpoklady? Jak je overite? # 5. Testujte stejnou nulovou hypotezu jako v 4. proti jednostranne alternative, ze vyska divek je vetsi nez vyska chlapcu. # # # Vsimnete si, ze vysledek testu je (jako cokoli jineho v Rku) objekt, a proto ho muzeme ulozit do promenne - vyzkousejte tst <- t.test(...) # V promenne tst mame nyni k dispozici vse, co bylo drive jen na obrazovce. Vyvolame pres znak $, napr. tst$statistic nebo tst$p.value # Staci i cast nazvu, pokud jednoznacne identifikuje pozadovany vysledek: tst$p. jiz urcuje p-value, zatimco tst$p neni zatim nic (mohl by to byt i tst$parameter). # Pojmenovani jednotlivych casti vystupu viz help(t.test), sekce Value. # # Vyhoda takto ulozenych vysledku je zrejma - pokud chceme napr. napsat skript, ktery uzivateli poskytne rovnou rozhodnuti a ne "jenom p-value", staci za test pridat radek if (tst$p. < 0.05) print(...). # Podobne se muze hodit pro dalsi vypocty spocitana hodnota testove statistiky atd. # Jeste vice teto moznosti vyvolani pres $ vyuzijeme pozdeji pri lin. regresnich modelech # # ------------------------------------------------------------------- # # Jak overovat normalitu dat? # # 1. graficky - zejmena histogram a tzv. Q-Q plot (vyberove kvantily vs. teoreticke kvantily N(0,1)) # lambda <- 0.5 x <- rnorm(100, mean=10, sd=2) y <- rexp(100, rate=lambda) par(mfrow=c(2,2)) # rozdeleni graf. okna na 4 casti hist(x, main="Histogram N(10,4)") hist(y, main="Histogram Exp(0.5)") qqnorm(x, main="Q-Q plot N(10,4)"); qqline(x) qqnorm(y, main="Q-Q plot Exp(0.5)"); qqline(y) # # 2. existuje rada statistickych testu normality # a) Shapiruv - Wilkuv test, specialni test primo pro normalni rozdeleni shapiro.test(x) shapiro.test(y) # b) D'Agostinuv test - test zalozeny na vyberove sikmosti a spicatosti # Normovane normalni rozdeleni ma (jak zajiste kazdy vi) teoretickou sikmost = 0 a teoretickou spicatost = 3. S temito hodnotami se pak porovnavaji vyberove charakteristiky spocitane z dat. # Pro D'Agostinuv test existuje makro Doc. Zvary - soubor dagost.R, ktere najdete u me na webu. source("dagost.R") DAgostino.test(x) DAgostino.test(y) # c) Kolmogorovuv - Smirnovuv test pro shodu distribucnich funkci, # obecny test zalozeny na empiricke distribucni funkci, na rozdil od # predeslych dvou funguje pro vsechna spojita rozdeleni # nevyhoda je, ze pro teoreticke rozdeleni musime specifikovat parametry: ks.test(x, "pnorm", 10, 2) # tj N(10,4) ks.test(y, "pnorm", 2, 2) # tj N(2,4) #K-S test je zalozen na maximu D absolutnich hodnot rozdilu empiricke a #skutecne distribucni funkce par(mfrow=c(2,1)) nadpis1 <- expression(paste("Empiricka distribucni funkce vektoru X~N(10,4) vs. distribucni funkce N(10,4)")) plot(ecdf(x),do.points=FALSE,verticals=TRUE,main=nadpis1) grid <- min(x)+0:1000*(max(x)-min(x))/1000 lines(grid,pnorm(grid,10,2),col="red") nadpis2 <- sprintf("Empiricka distribucni funkce vektoru Y~Exp(%.1f) vs. distribucni funkce N(2,4)",lambda) plot(ecdf(y),do.points=FALSE,verticals=TRUE,main=nadpis2) grid <- 0:1000*max(y)/1000 lines(grid,pnorm(grid,2,2),col="red") # Pro overovani normalniho rozdeleni bych preferoval spise # Shapiro-Wilk nebo D'Agostino, K-S se vice hodi pro obecnejsi situace, #kdy napr. chceme rozhodnout, zda 2 vybery pochazeji ze stejneho rozdeleni. # # ------------------------------------------------------------------- # ##### Dvouvyberovy testy # # za predpokladu normality t.test(vyska ~ pohlavi, data=vys, var.equal=TRUE) # t-test (stejne rozptyly) t.test(vyska ~ pohlavi, data=vys) # Welchuv t-test # bez predpokladu normality wilcox.test(vyska ~ pohlavi, data=vys) # Wilcoxon rank-sum (Mann-Whitney) test # #### Parovy vs. Dvouvyberovy t-test # # v souborech quicksort1.txt a quicksort2.txt jsou ulozena data s velicanami pro tridici algoritmus quicksort # oba soubory obsahuji dva vybery, jedna se vsak o dve odlisne situace (proc?), ktere vyzaduji rozdilny statisticky pristup # jak budeme testovat nulovou hypotezu o rovnosti strednich hodnot v prvnim a jak v druhe pripade? # # q() # Koncime ...