####################################################################################### ########################### Vícevýběrové testy ######################################## ######################### Martina Litschmannová, Adéla Vrtková ######################## ####################################################################################### ################## Přiklad 1. (Doplnění tabulky ANOVA) ################################ ####################################################################################### n=c(40,40,42) # rozsahy výběrů prum=c(300,290,310) # průměry v jednotlivých skupinách / třídách s=c(33,34,31) # směrodatné odchylky v jednotlivých skupinách / třídách n.total=sum(n) # celkový rozsah výběrů k=3 # počet tříd df.b=k-1 # počet stupňů volnosti - meziskupinový df.e=n.total-k # počet stupňů volnosti - reziduální # celkový průměr prum.total=weighted.mean(prum,n) prum.total # meziskupinový součet čtverců ss.b=sum(n*(prum-prum.total)^2) ss.b # reziduální součet čtverců ss.e=sum((n-1)*s^2) ss.e # celkový součet čtverců ss.b+ss.e # rozptyl mezi skupinami / třídami ms.b=ss.b/df.b ms.b # rozptyl uvnitř skupin / tříd ms.e=ss.e/df.e ms.e # F-poměr F=ms.b/ms.e F # p-hodnota p=1-pf(F,df.b,df.e) p # Na hladině významnosti 0,05 zamítáme hypotézu o shodě středních hodnot (p-hodnota=0,024, ANOVA), # tj. střední hodnoty alespoň jedné dvojice skupin se statisticky významně liší. # odhady skupinových efektů efekt=prum-prum.total efekt # Oproti celkovému průměru vykazuje nejvíce podprůměrné výsledky skupina 2 # (o cca 10 jednotek nižší než celkový průměr). Naopak průměr skupiny 3 je # o cca 10 jednotek vyšší než celkový průměr. Průměrné výsledky skupiny 1 # odpovídají celkovému průměru. ################## Přiklad 2. (Kyselina listová) ###################################### ####################################################################################### setwd("C:/Users/Martina/OneDrive - VŠB-TU Ostrava/Výuka/Pravděpodobnost a statistika/DATA/aktualni/") library(readxl) kysel = read_excel("testy_vicevyberove.xlsx", sheet=1) colnames(kysel)=c("sk1","sk2","sk3") # přejmenování sloupců # převod do standardního datového formátu kysel.s=stack(kysel) colnames(kysel.s)=c("hodnoty","skupina") kysel.s=na.omit(kysel.s) par(mfrow = c(1,1)) boxplot(hodnoty~skupina,data = kysel.s) # nebo boxplot(kysel) # Data neobsahují odlehlá pozorování. # Ověření normality skupiny = c("sk1","sk2","sk3") par(mfrow=c(1,3)) for (i in 1:3){ qqnorm(kysel.s$hodnoty[kysel.s$skupina == skupiny[i]], main = skupiny[i], xlab = "norm. teoretické kvantily", ylab = "výběrové kvantily") qqline(kysel.s$hodnoty[kysel.s$skupina == skupiny[i]]) } library(moments) tapply(kysel.s$hodnoty,kysel.s$skupina,skewness,na.rm=TRUE) tapply(kysel.s$hodnoty,kysel.s$skupina,moments::kurtosis,na.rm=TRUE)-3 tapply(kysel.s$hodnoty,kysel.s$skupina,shapiro.test) # Na hladině významnosti 0,05 nezamítáme předpoklad normality. # Informace potřebné pro nastavení zaokrouhlování tapply(kysel.s$hodnoty,kysel.s$skupina,length)# sd zaokrouhlujeme na 2 platné cifry tapply(kysel.s$hodnoty,kysel.s$skupina,sd)# sd a míry polohy zaokrouhlujeme na celá čísla # Ověření shody rozptylů s2=tapply(kysel.s$hodnoty,kysel.s$skupina,var,na.rm=TRUE) s2 # výběrové rozptyly max(s2)/min(s2) # Dle krabicového grafu a informace o poměru největšího a nejmenšího rozptylů (<2) nepředpokládáme, # že se rozptyly statisticky významně liší # Předpoklad normality nebyl zamítnut -> Bartlettův test bartlett.test(hodnoty~skupina,data = kysel.s) # nebo bartlett.test(kysel) # Na hladině významnosti 0,05 nelze zamítnout předpoklad o shodě rozptylů # (Bartlettův test, x_OBS = 0,878, df = 2, p-hodnota = 0,645). # Chceme srovnávat stř. hodnoty nezávislých výběrů z normálních rozdělení se stejnými rozptyly -> ANOVA # příkaz aov() vyžaduje data ve standardním datovém formátu vysledky=aov(hodnoty~skupina,data = kysel.s) # POZOR! Nestačí použít příkaz aov(). Výstup příkazu musíme uložit do pomocné proměnné # a na tu aplikovat příkaz summary(). summary(vysledky) # Na hladině významnosti 0,05 zamítáme hypotézu o shodě středních hodnot (ANOVA, p-hodnota<<0,001) -> mnohonásobné porovnávání TukeyHSD(vysledky) # průměry prum=tapply(kysel.s$hodnoty,kysel.s$skupina,mean,na.rm=TRUE) prum # skupinové efekty n=tapply(kysel.s$hodnoty,kysel.s$skupina,length) prum.total=weighted.mean(prum,n) efekt=prum-prum.total efekt # Považujeme-li vysoký obsah kyseliny listové za pozitivní, pak statisticky významně nejlepších # výsledků dosáhli pacienti ze skupiny 1 (průměrný obsah kys. listové o cca 27 jednotek vyšší než # průměrný obsah kys. listové v krvi všech testovaných pacientů) a statisticky významně nejhorších výsledků # dosáhli pacienti ze skupiny 2 (průměrný obsah kys. listové o cca 26 jednotek nižší než # průměrný obsah kys. listové v krvi všech testovaných pacientů). Obsah kys. listové v krvi pacientů ze skupiny 3 # odpovídá celkovému průměru. Všechny tři skupiny pacientů jsou navzájem dle obsahu kys. listové # v krvi statisticky významně odlišné. # Další část skriptu již neobsahuje tak podrobné komentáře =) ################## Přiklad 3. (Králíci) ############################################### ####################################################################################### library(readxl) kralici = read_excel("testy_vicevyberove.xlsx", sheet=2) colnames(kralici)=c("viden","cesky","kalif") # přejmenování sloupců # převod do standardního datového formátu kralici.s=stack(kralici) colnames(kralici.s)=c("hodnoty","skupina") kralici.s=na.omit(kralici.s) par(mfrow=c(1,1)) boxplot(kralici) # Odstranění odlehlého pozorování pom=boxplot(kralici) pom kralici.s$hodnoty.bez=kralici.s$hodnoty kralici.s$hodnoty.bez[kralici.s$hodnoty.bez %in% pom$out]=NA # Krabicový graf boxplot(hodnoty.bez~skupina,data = kralici.s) # Ověření normality skupiny = c("viden","cesky","kalif") nadpisy = c("Vídeňský bílý\nmodrooký","Český \nstrakáč","Kalifornský") par(mfrow=c(1,3)) for (i in 1:3){ qqnorm(kralici.s$hodnoty.bez[kralici.s$skupina == skupiny[i]], main = nadpisy[i], xlab = "norm. teoretické kvantily", ylab = "výběrové kvantily") qqline(kralici.s$hodnoty.bez[kralici.s$skupina == skupiny[i]]) } library(moments) tapply(kralici.s$hodnoty.bez,kralici.s$skupina,skewness,na.rm=TRUE) tapply(kralici.s$hodnoty.bez,kralici.s$skupina,moments::kurtosis,na.rm=TRUE)-3 tapply(kralici.s$hodnoty.bez,kralici.s$skupina,shapiro.test) # Na hladině významnosti 0,05 nezamítáme předpoklad normality. # Informace potřebné pro nastavení zaokrouhlování tapply(kralici.s$hodnoty,kralici.s$skupina,length)# sd zaokrouhlujeme na 2 platné cifry tapply(kralici.s$hodnoty,kralici.s$skupina,sd)# sd a míry polohy zaokrouhlujeme na setiny (sjednocení napříč druhy králíků) # Ověření shody rozptylů s2=tapply(kralici.s$hodnoty.bez,kralici.s$skupina,var,na.rm=TRUE) s2 max(s2)/min(s2) # Dle krabicového grafu a informace o poměru největšího a nejmenšího rozptylů (blízký 2, avšak rozsah výběrů < 30) # je těžší odhadnout, zda lze předpokládat shodu rozptylů. Rozhodnout nám pomůže test. # Předpoklad normality nebyl zamítnut -> Bartlettův test bartlett.test(hodnoty.bez~skupina,data = kralici.s) # Na hladině významnosti 0,05 nelze zamítnout předpoklad o shodě rozptylů # (Bartlettův test, x_OBS = 3,1, df = 2, p-hodnota = 0,217). # Chceme srovnávat stř. hodnoty nezávislých výběrů z normálních rozdělení se stejnými rozptyly -> ANOVA # příkaz aov() vyžaduje data ve standardním datovém formátu vysledky=aov(hodnoty.bez~skupina,data = kralici.s) summary(vysledky) # Na hladině významnosti 0,05 zamítáme hypotézu o shodě středních hodnot (p-hodnota<<0,001, ANOVA) -> mnohonásobné porovnávání TukeyHSD(vysledky) # průměry prum=tapply(kralici.s$hodnoty.bez,kralici.s$skupina,mean,na.rm=TRUE) prum # skupinové efekty # rozsahy výběrů bez NA hodnot (POZOR! Funkce length() počítá i NA hodnoty.) library(dplyr) n=tapply(kralici.s$hodnoty.bez,kralici.s$skupina,n_distinct,na.rm = T) prum.total=weighted.mean(prum,n) efekt=prum-prum.total efekt ################## Přiklad 4. (Jakost) ################################################ ####################################################################################### library(readxl) jakost.s = read_excel("testy_vicevyberove.xlsx", sheet=3) colnames(jakost.s)=c("poradi","skupina") # přejmenování sloupců par(mfrow=c(1,1)) boxplot(poradi~skupina, data = jakost.s) # Ověření normality nemá smysl provádět - z povahy jde o diskrétní data (pořadí) # Informace potřebné pro nastavení zaokrouhlování tapply(jakost.s$poradi,jakost.s$skupina,length)# sd zaokrouhlujeme na 2 platné cifry tapply(jakost.s$poradi,jakost.s$skupina,sd)# sd a míry polohy zaokrouhlujeme na celá čísla # Ověření shody rozptylů s2=tapply(jakost.s$poradi,jakost.s$skupina,var,na.rm=TRUE) s2 max(s2)/min(s2) # Dle krabicového grafu a informace o poměru největšího a nejmenšího rozptylů (<2) # ze předpokládat shodu rozptylů. (Kruskalův - Wallisův test má větší sílu testu, jsou-li data homoskedasticitní.) # Jde o "pořadová" data, nemá smysl uvažovat o předpokladu normality -> Leveneův test library(car) leveneTest(poradi~skupina, data = jakost.s) jakost.s$skupina=as.factor(jakost.s$skupina) # Starší verze příkazu leveneTest vyžaduje, aby "skupina" byla kódována jako faktor, # nově je výpočet proveden, ale objevuje se doporučení k překódování proměnné "skupina". leveneTest(poradi~skupina, data = jakost.s) # Na hladině významnosti 0,05 nelze zamítnout předpoklad o shodě rozptylů # (Leveneho test, x_OBS = 0,4, df_num = 3, df_denom = 62, p-hodnota = 0,750). # Ověření symetrie par(mfrow=c(1,1)) boxplot(poradi~skupina, data = jakost.s) tapply(jakost.s$poradi,jakost.s$skupina,moments::skewness) skupiny = c("A","B","C","D") par(mfrow=c(1,1)) for (i in 1:4){ hist(jakost.s$poradi[jakost.s$skupina == skupiny[i]], main = skupiny[i], xlim = c(1,70), ylim = c(0,6), xlab = "pořadí", ylab = "četnost") } # Chceme srovnávat mediány nezávislých výběrů (data nelze považovat za výběry z norm. rozdělení) -> Kruskalův-Wallisův test kruskal.test(poradi~skupina,data = jakost.s) # Na hladině významnosti 0,05 nelze zamítnout hypotézu o shodě mediánů (Kruskalův-Wallisův test, x_OBS = 3,7, df = 3, p-hodnota=0,295). # Tj. statisticky významné rozdíly mezi výrobci (z hlediska pořadí výrobků v soutěži) neexistují. ################## Přiklad 5. (Trombin) ############################################### ####################################################################################### library(readxl) trombin.s = read_excel("testy_vicevyberove.xlsx", sheet=4, skip = 1) colnames(trombin.s)=c("hodnoty","skupina") # přejmenování sloupců par(mfrow=c(1,1)) boxplot(hodnoty~skupina,data = trombin.s) # Ověření normality skupiny = c("A","B","C") par(mfrow=c(1,3)) for (i in 1:3){ qqnorm(trombin.s$hodnoty[trombin.s$skupina == skupiny[i]], main = skupiny[i], xlab = "norm. teoretické kvantily", ylab = "výběrové kvantily") qqline(trombin.s$hodnoty[trombin.s$skupina == skupiny[i]]) } library(moments) tapply(trombin.s$hodnoty,trombin.s$skupina,skewness,na.rm=TRUE) tapply(trombin.s$hodnoty,trombin.s$skupina,moments::kurtosis,na.rm=TRUE)-3 tapply(trombin.s$hodnoty,trombin.s$skupina,shapiro.test) # Na hladině významnosti 0,05 zamítáme předpoklad normality. # Informace potřebné pro nastavení zaokrouhlování tapply(trombin.s$hodnoty,trombin.s$skupina,length)# sd zaokrouhlujeme na 2 platné cifry tapply(trombin.s$hodnoty,trombin.s$skupina,sd)# sd a míry polohy zaokrouhlujeme na setiny (sjednocení napříč skupinami) # Ověření shody rozptylů s2=tapply(trombin.s$hodnoty,trombin.s$skupina,var,na.rm=TRUE) s2 max(s2)/min(s2) # Dle krabicového grafu a informace o poměru největšího a nejmenšího rozptylů (>>2) # nelze předpokládat shodu rozptylů. Rozhodnout nám pomůže test. # Předpoklad normality byl zamítnut -> Leveneho test library(car) leveneTest(hodnoty~skupina,data = trombin.s) trombin.s$skupina=as.factor(trombin.s$skupina) leveneTest(hodnoty~skupina,data = trombin.s) # Na hladině významnosti 0,05 zamítáme předpoklad o shodě rozptylů # (Leveneho test, p-hodnota << 0,001). # Ověření symetrie par(mfrow=c(1,1)) boxplot(hodnoty~skupina, data = trombin.s) tapply(trombin.s$hodnoty,trombin.s$skupina,moments::skewness) skupiny = c("A","B","C") par(mfrow=c(1,1)) for (i in 1:3){ hist(trombin.s$hodnoty[trombin.s$skupina == skupiny[i]], main = skupiny[i], xlim = c(0.9*min(trombin.s$hodnoty),1.1*max(trombin.s$hodnoty)), ylim = c(1,5), xlab = "trombinový čas (s)", ylab = "četnost") } # Chceme srovnávat mediány nezávislých výběrů, která nemají normální rozdělení a mají různé rozptyly -> Kruskalův - Wallisův test kruskal.test(trombin.s$hodnoty,trombin.s$skupina) # Na hladině významnosti 0,05 zamítáme hypotézu o shodě mediánů (p-hodnota<<0,001, Kruskalův-Wallisův test). Tj. trombinový čas je statisticky významně # ovlivněn preparátem. -> mnohonásobné porovnávání library(dunn.test) dunn.test(trombin.s$hodnoty,trombin.s$skupina, method = "bonferroni", alpha = 0.05, # nastavíme-li hodnotu alfa, jsou st. významné rozdíly označeny * # (defaultně: alpha = 0.05) altp = T) # altP = T nastavuje p-hodnotu tak,aby se při rozhodování # o statistické významnosti srovnávala s alfa (defaultně: altp = F, pak srovnáváme s alfa/2) # mediány med=tapply(trombin.s$hodnoty,trombin.s$skupina,quantile,prob=0.5,na.rm=TRUE) med # skupinové efekty med.total=quantile(trombin.s$hodnoty,probs=0.5) efekt=med-med.total efekt