################################## ## Analyza dat v R ## ################################## ## Lekce 5 - Zakladni testy ## ################################## # v teto lekci si sami vyzkousime zakladní testy, ktere jsem si ukazali na prednasce # vycisteni databaze objektu pred nasi praci rm(list=ls()) # nastaveni pracovniho adresare # toto si musite nastavit sami podle sebe, podle toho, kam jste si soubory ulozili setwd("/home/marketa/prednasky/Motol-seminar/podklady/") # nacte tabulku, predpoklada umisteni souboru v pracovnim adresari # nacte tabulku, predpoklada umisteni souboru v pracovnim adresari Kojeni = read.csv2("Kojeni.csv",sep=",") attach(Kojeni) # zpristupneni novych udaju ####################### # Jednovyberove testy # ####################### # podíváme se, jak je to s hmotností dětí v 6. měsících # podle tabulek je prumer u chlapcu 7.8 a u divek 7.2kg hist(hmotnost[hoch==1]) hist(hmotnost[hoch==0]) summary(hmotnost[hoch==1]) # testujeme, zda je prumerna hmotnost chlapcu v souboru 7.8kg ######## # T-TEST # podoba testu # t.test(x, y = NULL, # alternative = c("two.sided", "less", "greater"), # mu = 0, paired = FALSE, var.equal = FALSE, # conf.level = 0.95, ...) # do polozky "mu" muzeme zadat tu testovanou hodnotu # defaultni je oboustranna alternativa a ocekavana stredni hodnota 0 tt = t.test(hmotnost[hoch==1],mu=7800) print(tt) str(tt) # nezamitame hypotezu, ze nasi chlapci pochazji ze stejne populace jako "tabulkova" # jak je to pro divky tt = t.test(hmotnost[hoch==0],mu=7200) print(tt) # hypotezu zamitame, nase divky jsou tezsi nez tabulkove # overme si jeste, ze rozdeleni lze povazovat na normalni shapiro.test(hmotnost[hoch==1]) shapiro.test(hmotnost[hoch==0]) # lze, u divek i chlapcu, uziti t-testu je tedy ok ################# # ZNAMENKOVY TEST # pomocna promenna "hmotnost chlapcu - testovana hodnota pom.hmot = hmotnost[hoch==1] - 7800 hist(pom.hmot) # pocet vsech nenulovych pozorovani N = sum(pom.hmot!=0,na.rm=T) # je jich 47, chlapcu je 49, u dvou se porodni hmotnost rovnala 7800g # pocet vsech zapornych k = sum(pom.hmot<0,na.rm=T) # test pomoci testu proporci (nebo ekvivalentne binomicky test) # prop.test(x, n, p = NULL, alternative = c("two.sided", "less", "greater"), conf.level = 0.95, correct = TRUE) prop.test(k,N,p = 0.5) binom.test(k,N,p = 0.5) # nezamitame # a pro divky pom.hmot = hmotnost[hoch==0] - 7200 N = sum(pom.hmot!=0,na.rm=T) # je jich 47, chlapcu je 49, u dvou se porodni hmotnost rovnala 7800g k = sum(pom.hmot<0,na.rm=T) prop.test(k,N,p = 0.5) # ani zde nezamitame ################# # WILCOXONUV TEST # wilcox.test(x, y = NULL, # alternative = c("two.sided", "less", "greater"), # mu = 0, paired = FALSE, exact = NULL, correct = TRUE, # conf.int = FALSE, conf.level = 0.95, ...) # vicemene stejna struktura jako u t-testu wilcox.test(hmotnost[hoch==1],mu=7800) wilcox.test(hmotnost[hoch==0],mu=7200) # rozhodovani stejne jako u t-testu # exaktni test nelze pouzit, protoze v datech jsou "ties", tedy hodnoty, ktere obsazuji stejne misto v poradi, # a take "nuly", tedy osoby s hmotnosti stejnou je je mu ################ # Parove testy # ################ # testujeme, zda je vyska otce o ocekavanych 10cm (populacni informace) vyssi nez vyska matek # H0: otcove jsou o 10 cm vyssi nez matky, H1: oboustranna # propojeni v paru skrze dite ######## # T-TEST # podminka paroveho t-testu je normalita rozdilu shapiro.test(vyska.o-vyska.m) # ok # prohledneme si to qqnorm(vyska.o-vyska.m);qqline(vyska.o-vyska.m) t.test(vyska.o, vyska.m, paired = T,mu = 10) # ve skutecnosti jsou muzi jeste vyssi nez o 10 cm ################# # ZNAMENKOVY TEST N = sum((vyska.o-vyska.m-10)!=0,na.rm=T) k = sum((vyska.o-vyska.m-10)<0,na.rm=T) # pouze cca tretina otcu ma vysku o mene nez 10 cm oproti matkam prop.test(k,N,p = 0.5) # zamitame H0 ################# # testujeme, zda jsou pary vekove stejne, tedy zda je stejny vek otce a vek matky ditete shapiro.test(vek.o-vek.m) # neni normalne rozdelene # prohledneme si to qqnorm(vek.o-vek.m);qqline(vek.o-vek.m) N = sum((vek.o-vek.m)!=0,na.rm=T) k = sum((vek.o-vek.m)<0,na.rm=T) # pouze cca tretina otcu ma vysku o mene nez 10 cm oproti matkam prop.test(k,N,p = 0.5) # zamitame - jen malo otcu je mladsich nez matka ################# # WILCOXONUV TEST # predpoklada symetrii!!!!!! hist(vyska.o-vyska.m) # na vysku lze pouzit hist(vek.o-vek.m) # na vek uz to fungovat nebude wilcox.test(vyska.o, vyska.m, paired = T,mu = 10) ###################### # Dvouvyberove testy # ###################### # testujeme, zda se lisi hmotnost sestimesicnich chlapcu a divek ######## # T-TEST # podminku normality uz mame overenou, pro kazdou podskupinu t.test(hmotnost[hoch==1],hmotnost[hoch==0]) # defaultne predpoklada, ze rozptyly nejsou stejne bartlett.test(hmotnost,hoch) # podle tohoto testu jsou t.test(hmotnost[hoch==1],hmotnost[hoch==0],var.equal=T) # alternativni zapis t.test(hmotnost~hoch,var.equal=T) # neprilis velky rozdil # testujeme, zda se lisi vek matky podle toho, zda planovala tehotenstvi # promenna "Plan" # vek matky zjevne neni normalne rozdeleny (nutno po kategoriich!) tapply(vek.m,Plan,shapiro.test) ################# # WILCOXONUV TEST wilcox.test(vek.m[Plan==0],vek.m[Plan==1]) wilcox.test(vek.m~Plan) #wilcox.test(x, y = NULL, # alternative = c("two.sided", "less", "greater"), # mu = 0, paired = FALSE, exact = NULL, correct = TRUE, # conf.int = FALSE, conf.level = 0.95, ...) # chceme videt rozdil v poloze wilcox.test(vek.m~Plan,conf.int=T) ggplot(Kojeni,aes(x=vek.m,fill=factor(Plan)))+geom_histogram(position = "dodge") ggplot(Kojeni,aes(x=vek.m,color=factor(Plan)))+ stat_ecdf(geom = "step") # neplanovana tehitenstvi viditelne v mladsim veku ########################### # KOLMOGOROV-SMIRNOVUV TEST ks.test(vek.m[Plan==0],vek.m[Plan==1]) ###################### # Analyza rozptylu # ###################### # vypnu stara data a natahnu nova detach(Kojeni) kenton = read.csv2("kentonfood.csv",sep="\t") kenton$design = factor(kenton$design) # potrebju, aby to byl faktor attach(kenton) # zpristupneni novych udaju # jde o data z ucebnice linearnich modelu Neter & al. # ctyri ruzne designy plechovek se prodavaly v peti prodejnach a merily se prodeje ggplot(kenton,aes(x=design,y=sales,fill=design))+geom_boxplot() # hypoteza: na designu plechovky nezalezi # overeni predpokladu (malo pozorovani, tak to asi projde) bartlett.test(sales,design) # rozptyly OK anova(lm(sales~factor(design))) fit = lm(sales~design) # normalitu je treba overit na reziduech!!!!! shapiro.test(resid(fit)) # normalita OK # rozhodne tam rozdily jsou # ale kde? Ktere kusy jsou stejne a ktere jine? # nasledujici testy umi jen vystup z aov a ne s lm fit = aov(sales~design) TukeyHSD(fit) plot(TukeyHSD(fit)) # rozdily jsou hlavne mezi 4 a zbytkem # 1,2,3 se mezi sebou nelisi # je to neco jako t-testy pro jednotlive kousky, ale se "spolecnou" p-hopdnotou ########################### # KOLMOGOROV-SMIRNOVUV TEST # vratime se ke Kojeni #detach(kenton) attach(kojeni) Kojeni$Vzdelani = relevel(Vzdelani,"zakladni") # abychom meli "zakladni" na zacatku ggplot(Kojeni,aes(x=Vzdelani,y=vek.m,fill=factor(Vzdelani)))+geom_boxplot() # neni normalne rozdeleno tapply(vek.m,Vzdelani,shapiro.test) # pouzijme rozsireni Wilcoxonova testu - poradovy test kruskal.test(vek.m,Vzdelani) # neumime ale rozhodovat automatizovane rozdily # jak by to vyslo na plechovky kruskal.test(sales,design)