####################################################################################### ################ Statistická laboratoř - 11. workshop ################################# ############### Adéla Vrtková, Martina Litschmannová ################################## ####################################################################################### ### Skript je určen pro účastníky workshopů pořádaných Statistickou laboratoří. ### V případě, že jej, popř. jeho části, chcete v nezměněné či mírně upravené podobě ### používat v rámci vlastních publikací nebo ve výukových kurzech realizovaných ### na vašem pracovišti, budeme rády, když nás o tom budete informovat a v příslušných ### materiálech uvedete poděkování Statistické laboratoři (http://k470.vsb.cz/statlab/). ### Za akceptování této prosby předem děkujeme! ###################################################################################### #### Akumulátory - kontingenční tabulky ############################################## ## Už umíme - načtení dat, převod do standartního datového formátu ,vynechání chybějících hodnot, atd. data=read.csv("http://am-nas.vsb.cz/lit40/DATA/aku.csv",sep=";", dec=",", header=TRUE) dataS = reshape(data, direction="long", varying=list(c("A5","B5","C5","D5"),c("A100","B100","C100","D100")), v.names=c("kap5","kap100"), times=c("A","B","C","D"), timevar="vyrobce") dataS=dataS[,-4] dataS=na.omit(dataS) dataS$pokles=dataS$kap5-dataS$kap100 ## cílem je posouzení relativního poklesu - zda byl pokles větší než o 10 % # definování nové proměnné - relativní pokles dataS$pokles.rel=dataS$pokles/dataS$kap5 # definování dichotomické proměnné - několik možností jak se k ní dostaneme dataS$pokles.dich=rep("NE",length(dataS$pokles.rel)) # nová proměnná, všude NE - rozumíme NE nebyl větší než o 10 % dataS$pokles.dich[which(dataS$pokles.rel>0.1)]="ANO" # tam, kde je rel. pokles nad 10 % dáme ANO # nebo dataS$pokles.dich[dataS$pokles.rel>0.1]="ANO" # tam, kde je rel. pokles nad 10 % dáme ANO # jinak - pomocí ifelse dataS$pokles.dich= ifelse( # pokud dataS$pokles.rel>0.1, # relativní pokles větší než 10% "ANO", # označ Ano "NE") # jinak označ Ne ### Prezentace kategoriální proměnné pomocí známé funkce table ################################ ## funkci table už známe z explorační analýzy kvalitativní proměnné - používali jsme ji pro jednu proměnnou table(dataS$pokles.dich) table(dataS$vyrobce) ## pomocí funkce table získáme kontingenční tabulku tab=table(dataS$vyrobce,dataS$pokles.dich) # pozor na pořadí!!! Dodržujme "příčina"-"důsledek" tab prop.table(tab) # relativní četnosti prop.table(tab)*100 # relativní četnosti v % prop.table(tab,1) # řádkové relativní četnosti prop.table(tab,2) # sloupcové relativní četnosti # grafická prezentace pomocí shlukového sloupcového grafu (angl. grouped barplot) barplot(tab, beside = TRUE) # nastavení typu grafu na shlukový # nevhodný graf - "záměna příčiny a následku" # oprava: tab2 = table(dataS$pokles.dich,dataS$vyrobce) barplot(tab2, beside = TRUE, # nastavení typu grafu na shlukový col = topo.colors(2), # nastavení barev výplní sloupců xlab="Výrobce", ylab="počet akumulátorů") legend(9,85, # souřadnice pro umístění legendy legend=c("ANO","NE"), # text legendy fill = topo.colors(2), # barva výplně horiz = TRUE) # nastavení stylu legendy (vedle sebe vs. pod sebou) # grafická prezentace pomocí skládaného sloupcového grafu (angl. stacked barplot) barplot(tab2, beside = FALSE, # nastavení typu grafu na skládaný col = topo.colors(2), xlab="Výrobce", ylab="počet akumulátorů") legend(3.6,86, # pozice legendy legend=c("ANO","NE"), # text legendy fill=topo.colors(2), # barva výplně horiz=TRUE) # nastavení stylu legendy (vedle sebe vs. pod sebou) # grafická prezentace pomocí 100% skládaného sloupcového grafu (angl. 100% stacked barplot) ptab2 = prop.table(tab2,2) # sloupcové relativní četnosti pro tab2 ptab2 ptab2 = prop.table(tab2,2)*100 # sloupcové relativní četnosti (%) pro tab2 ptab2 barplot(ptab2, beside = FALSE, # nastavení typu grafu na skládaný col = topo.colors(2), xlab="Výrobce", ylab="rel. četnost akumulátorů (%)") legend(4.15,100, # pozice legendy legend=c("ANO","NE"), # text legendy fill=topo.colors(2), # barva výplně horiz=FALSE) # nastavení stylu legendy (vedle sebe vs. pod sebou) ## grafická prezentace pomocí mozaikového grafu mosaicplot(tab) mosaicplot(tab, color = topo.colors(2), main = "", cex=1.5) # nastavení velikosti písma ## Kvízová otázka - v čem je mozaikový graf lepší? par(mfrow=c(2,2)) #graf 1 mosaicplot(tab, color = topo.colors(2), main = "") # graf 2 barplot(ptab2, beside = FALSE, # nastavení typu grafu na skládaný col = topo.colors(2), xlab="Výrobce", ylab="rel. četnost akumulátorů (%)") legend(3.8,100, # pozice legendy legend=c("ANO","NE"), # text legendy fill=topo.colors(2), # barva výplně horiz=FALSE) # nastavení stylu legendy (vedle sebe vs. pod sebou) # graf 3 barplot(tab2, beside = FALSE, # nastavení typu grafu na skládaný col = topo.colors(2), xlab="Výrobce", ylab="počet akumulátorů") legend(3,86, # pozice legendy legend=c("ANO","NE"), # text legendy fill=topo.colors(2), # barva výplně horiz=TRUE) # nastavení stylu legendy (vedle sebe vs. pod sebou) #graf 4 barplot(tab2, beside = TRUE, # nastavení typu grafu na shlukový col = topo.colors(2), # nastavení barev výplní sloupců xlab="Výrobce", ylab="počet akumulátorů") legend(8,85, # souřadnice pro umístění legendy legend=c("ANO","NE"), # text legendy fill = topo.colors(2), # barva výplně horiz = TRUE) # nastavení stylu legendy (vedle sebe vs. pod sebou) # # # # ## Odpověď: Mozaikový graf zohledňuje rel. zastoupení výrobců (marginální rel. četnosti) i rel. zastoupení typu ## akumulátorů pro jednotlivé výrobce (řádkové relativní četnosti). ###################################### ## Výrobce B a C má velmi špatné výsledky - vyloučíme je ze sledování a porovnáme pouze výrobce A a D tab tab3=tab[c(1,4),] tab3 # jiný způsob zadání tab3 (ručně) tab3 = matrix(c(44,42,8,71),nrow=2,byrow=TRUE) tab3 rownames(tab3) = c("A","D") colnames(tab3) = c("ANO","NE") tab3 ### Test nezávislosti v kont. tabulce a míra kontingence # H0: Výrobce a přítomnost poklesu kapacity o více než 10 % jsou statisticky nezávislé veličiny. # HA: Výrobce a přítomnost poklesu kapacity o více než 10 % nejsou statisticky nezávislé veličiny. chisq.test(tab3) # pozor - funkce sama nezkontroluje předpoklady testu! pom=chisq.test(tab3) pom$expected # očekávané četnosti>=2, # alespoň 80% pak >5 # - splněno -> lze testovat korektně -> dívám se na p-hodnotu pom$p.value ## Na hladině významnosti 0.05 zamítám H0, tzn. výrobce a přítomnost poklesu kapacity o více než 10 % # nejsou statisticky nezávislé veličiny. Tj. záleží na tom, od jakého výrobce akumulátory pocházejí. install.packages("lsr") library(lsr) cramersV(tab3) # míra kontingence Cramerovo V naznačuje středně silnou závislost ##### Relativní riziko a poměr šancí # Z grafu to vypadá, že výrobce D je na tom lépe - lze kvantifikovat kolikrát a jestli vůbec stat. významně? tab3 install.packages("epiR") library(epiR) epi.2by2(tab3) ## Interpretace RR: Výrobce A má cca 5.1x větší riziko (pravděpodobnost), že jeho akumulátory ## budou mít pokles kapacity o více než 10 %, než výrobce D. ## Intervalový odhad RR napoví statistickou významnost. ## S 95% spolehlivostí lze očekávat, že výrobce A má 2.5 až 10.1 krát větší riziko, že jeho akumulátory ## budou mít pokles kapacity o více než 10 %, než výrobce D. ## Interpretace OR: Výrobce A má cca 9.3x větší šanci, že jeho akumulátory ## budou mít pokles kapacity o více než 10 %, než výrobce D. ## Intervalový odhad OR napoví statistickou významnost. ## S 95% spolehlivostí lze očekávat, že výrobce A má 4.0 až 21.6 krát větší šanci, že jeho akumulátory ## budou mít pokles kapacity o více než 10 %, než výrobce D.