######################################## ## Analyza dat v R ## ######################################## ## Lekce 7 - Kontingenční tabulky ## ######################################## # v teto lekci si ukazeme, jak pracujeme s kvalitativnimi znaky # 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/kody/") # nacte tabulku, predpoklada umisteni souboru v pracovnim adresari ichs = read.csv2("ichs.csv",sep=";",fileEncoding = "WINDOWS-1250") attach(ichs) names(ichs) # pozor, data jsou oproti prednasce "zmensena" Kojeni = read.csv2("Kojeni.csv",sep=",") # nacte uzitecne knihovny library(ggplot2) library(scales) # procentova skala do grafu ##################################### ## Prohlizeni kvalitativnich znaku ## ##################################### summary(Vzdel) table(Vzdel) # jednoduchy barplot barplot(table(Vzdel)) #vidime, ze serazeno abecedne, prehodime si zakladni vzdelani na zacatek ichs$Vzdelani = relevel(Vzdel,"zákl.") barplot(table(ichs$Vzdelani)) # jednoduchy barplot barplot(table(Kour)) #vidime, ze serazeno abecedne, prehodime přehodime silneho nakonec a nekuraka na zacatek ichs$Koureni = factor(ichs$Kour,levels = levels(ichs$Kour)[c(2,1,4,3)]) barplot(table(ichs$Koureni)) # prezentovatelnejsi graf pomoci ggplotu ggplot(ichs,aes(x=Vzdelani))+geom_bar() #vybarveno podle kategorie ggplot(ichs,aes(x=Vzdelani,fill=Vzdelani))+geom_bar() # popisky, nadpisy atd ggplot(ichs,aes(x=Vzdelani,fill=Vzdelani))+geom_bar()+ ggtitle("Vzdělání mužů v souboru")+ scale_x_discrete("Vzdělání")+ scale_y_continuous("Absolutní četnost",breaks = seq(0,80,by=10))+ scale_fill_discrete("Vzdělání",labels = c("základní","maturita","vysokoškolské")) # krizeni s kourenim, absolutni pocty ggplot(ichs,aes(x=Vzdelani,fill=Koureni))+geom_bar() #obracene ggplot(ichs,aes(fill=Vzdelani,x=Koureni))+geom_bar() # vedle sebe ggplot(ichs,aes(fill=Vzdelani,x=Koureni))+geom_bar(position = "dodge") # ve sloupci a relativní cetnosti ggplot(ichs,aes(fill=Vzdelani,x=Koureni))+geom_bar(position = "fill") # popisky, nadpisy, procenta atd ggplot(ichs,aes(fill=Vzdelani,x=Koureni))+geom_bar(position = "fill")+ ggtitle("Vzdělání mužů v souboru podle kouření")+ scale_x_discrete("Status kouření")+ scale_y_continuous("Relativní četnost",breaks = seq(0,1,by=.10),labels=percent)+ scale_fill_discrete("Vzdělání",labels = c("základní","maturita","vysokoškolské")) ###################### ## Test dobre shody ## ###################### # priklad srovnani rozlozeni krevnich skupin ve vyberu se znamou hodnoutou v populaci, # ze ktere by mel vyber pochazet # namerene cetnosti ve vyberu 200 osob oo = c(56,72,54,18) # cetnosti v originálni populaci pp = c(0.35,0.35,0.2,0.1) #ocekavane cetnosti ee = pp*sum(xx) # rucni napocitani # prispevky jednotlivych kategorii chi.pomocny = ((oo-ee)^2)/ee # testova statistika a kriticka hodnota chi = sum(chi.pomocny) krit.chi = qchisq(0.95,df=length(oo)-1) # podivame se na vysledek c(chi,krit.chi) # rozhodnuti testu chi>krit.chi # p-hodnota 1-pchisq(chi,df=length(oo)-1) # uziti funkce z R, jeden vektor a je treba specifikovat ocekavane pravdepodobnosti (p = ...) chisq.test(oo,p=pp) #podivame se, co nam nabizi za vystupy tt = chisq.test(oo,p=pp) str(tt) # absolutní rozdíly pro každou z kategorií oproti očekávání # pokud umocnime na druhou, dostaneme jednotlive prispevky k chi statistice tt$residuals tt$residuals^2-chi.pomocny # standardizovane rozdily # prevazene velikosti kategorie, ukazuje, ktera skupina se nejvice lisi bez ohledu na velikost tt$stdr ###################### ## Test homogenity ## ###################### # testujeme, zda je rozlozeni ve dvou populacich stejne cc = c(121,120,79,33) dd = c(118,95,121,30) # tohle je spatne!!!! # tohle by predpokladalo, ze prvni pozice v kazdem vektoru odpovida jedne osobe, druha druhe osobe atd chisq.test(cc,dd) # tento zapis testu ve skutecnosti generuje tabulku, na ktere pak dela test table(cc,dd) # podivejte se na ten pocet stupnu volnosti, to vas varuje # toto je spravny zapis, zadavame primo tabulku 2x4 chisq.test(rbind(cc,dd)) ####################################### ## Test Hardy-Weinbergovy rovnovahy ## ####################################### # namerena data genotypu, odpovidaji AA, Aa, aa aa = c(18,17,6) # slo by srovnat s nejakou populaci # napriklad, kdyby nam nekdo rekl, ze "spravne" rozlozeni genotypu je c(0.4,0.45,0.15) # pak by to byl chisq.test se dvema stupni volnosti (tri kategorie) chisq.test(aa,p = c(0.4,0.45,0.15)) # toto ale neni ten pripad # nemam odhad (rozlozeni cetnosti) zadano zvnejsku a netestuju, zda to populaci odpovida # ve skutecnosti testuju, zda je dochazi k nezavislemu krizeni a nedochazi k nejakemu posunu od genotypu k fenotypu # delam predpoklad nezavislosti a na jeho zaklade stanovim odhad cetnosti A za vyse uvedeneho predpokladu n = sum(aa) theta = (2*aa[1]+aa[2])/(2*n) theta # dale postupujeme jako v prvnim prikladu ee = c(theta^2,2*theta*(1-theta),(1-theta)^2)*n chi = sum(((aa-ee)^2)/ee) p = 1-pchisq(chi,df=1) # muzeme si to napsat jako funkci # jinak je to soucasti rady baliku hwe = function(aa){ n = sum(aa) theta = (2*aa[1]+aa[2])/(2*n) ee = c(theta^2,2*theta*(1-theta),(1-theta)^2)*n chi = sum(((aa-ee)^2)/ee) p = 1-pchisq(chi,df=1) res = c(chisq = chi, df = 1,p.val=p) print(res) } hwe(aa) ######################## ## Test nezavislosti ## ######################## # vracime se k datum ichs # nejprve je odpojime a pak znova pripojime, protoze jsem mezitim pridali nejake nove promenne detach(ichs) attach(ichs) # kontingencni tabulka tt = table(Koureni,Vzdelani) # kontingencni tabulka s marginalnimi cetnostmi addmargins(tt) # kontingencni tabulka s relativnimi cetnostmi # sectou se vsechny do 1 prop.table(tt) # radkova procenta prop.table(tt,margin = 1) apply(prop.table(tt,margin = 1),1,sum) # kontrola souctu po radkach # sloupcova procenta prop.table(tt,margin = 2) apply(prop.table(tt,margin = 2),2,sum) # kontrola souctu po sloupcich # test primo na tabulku chisq.test(tt) # test primo na data chisq.test(Koureni,Vzdelani) # upozorneni, ze nekde by mol byt problem s malou cetnosti chi = chisq.test(Koureni,Vzdelani) chi$expected Problem je u slabych kuraku, ocekavane jsou vsechny pod 3 # mozna naprava # 1. pripojeni slabych kuraku ke kurakum Koureni1 = Koureni Koureni1[which(Koureni=="slabý")] = "silný" table(Koureni1,Vzdelani) # tak to nam nepomohlo, musime zrusit cely ten faktor, ne ho jen vyprazdnit Koureni1 = as.character(Koureni) # musime zrusit faktory, jinak to bude zlobit Koureni1[which(Koureni1=="slabý")] = "silný" Koureni1[which(Koureni1=="silný")] = "kuřák" table(Koureni1,Vzdelani) chisq.test(Koureni1,Vzdelani) # 2. pouzitivjineho testu, treba Fisherova fisher.test(Koureni,Vzdelani) ff = fisher.test(Koureni,Vzdelani) str(ff) # lze specifikovat i jednostrannou alternativu ###################### ## McNemaruv test ## ###################### # "vzdelani" otce # vygeneruje nahodne se stejnou cetnosti kategorii jako u matek Kojeni$Vzdelani = relevel(Kojeni$Vzdelani,"zakladni") pom = rmultinom(nrow(Kojeni),1, prob = table(Kojeni$Vzdelani)/nrow(Kojeni)) Kojeni$VzdelaniO = factor(apply(pom*c(1,2,3),2,sum),labels = c("základní","maturita","VŠ")) addmargins(table(Kojeni$Vzdelani,Kojeni$VzdelaniO)) mcnemar.test(Kojeni$Vzdelani,Kojeni$VzdelaniO) # generovali jsme náhodně, takže četnost "nevyvážených" dvojic by mela byt na obe strany stejna