####################################################################################### ########################## Jednovýběrové testy ######################################## ######################### Martina Litschmannová ####################################### ####################################################################################### ################## Přiklad 1. (Albumin) ############################################### ####################################################################################### # Test o střední hodnotě / mediánu ## Načtení dat z xlsx souboru (pomoci balíčku readxl) library(readxl) albumin = read_excel("C:/Users/Martina/Výuka/Pravděpodobnost a statistika/DATA/aktualni/testy_jednovyberove.xlsx", sheet = "bilk_serum") colnames(albumin)="hodnoty" ## Explorační analýza boxplot(albumin$hodnoty) # Data neobsahují odlehlá pozorování. length(albumin$hodnoty)# sd zaokrouhlujeme na 3 platné cifry sd(albumin$hodnoty) # sd a míry polohy zaokrouhlujeme na tisíciny summary(albumin$hodnoty) # Pokud je třeba, nastavíme požadovaný formát vypočtených číselných charakteristik options(scipen = 100, # 100 (desetinné číslo), -100 (vědecká notace) digits = 7) # max. počet platných cifer mean(albumin$hodnoty) # Ověření normality qqnorm(albumin$hodnoty) qqline(albumin$hodnoty) hist(albumin$hodnoty) library(moments) skewness(albumin$hodnoty) moments::kurtosis(albumin$hodnoty)-3 # Vizualizace dat, šikmost i špičatost odpovídají norm. rozdělení. # Pro konečné rozhodnutí o normalitě dat použijeme test normality. # Předpoklad normality ověříme Shapirovovým - Wilkovovým testem. # H0: Data jsou výběrem z normálního rozdělení. # Ha: Data nejsou výběrem z normálního rozdělení. shapiro.test(albumin$hodnoty) # Na hl. významnosti 0,05 nelze předpoklad normality zamítnout # (Shapirův-Wilkův test, p-hodnota=0,236). # H0: mu = 35 g/l # Ha: mu <> 35 g/l t.test(albumin$hodnoty,mu=35, alternative = "two.sided",conf.level=0.95) # Na hl. významnosti 0,05 zamítáme nulovou hypotézu ve prospěch hypotézy alternativní # (t-test, t = -19,2, df = 217, p-hodnota << 0,001). # Střední hodnota albuminu se statisticky významně liší od 35 g/l. # Doplnění plynoucí z bodového a intervalového odhadu mu. # V populaci očekáváme střední hodnotu albuminu cca 34,487 g/l. # 95% intervalový odhad této střední hodnoty je # 34,434 g/l až 34,540 g/l.(I z tohoto lze vyvodit, že se # střední hodnota albuminu statisticky významně liší od 35 g/l. # Pro zájemce: Vizualizace klasického testu. n = length(albumin$hodnoty) x = seq(-20,4,0.01) f = dt(x,n -1) df = data.frame(x,f) pom = t.test(albumin$hodnoty,mu=35, alternative = "two.sided",conf.level=0.95) names(pom) library(ggplot2) ggplot(df,aes(x = x,y = f))+ geom_line()+ theme_bw()+ ylab("hustota p-sti nulového rozdělení")+ xlab("pozorovaná hodnota")+ stat_function(fun = dt, geom = "area", fill = "steelblue", xlim = c(min(x),qt(0.025,n-1)), alpha = 0.2, args = list(df = n) )+ geom_segment(aes(x = min(x), xend = qt(0.025,n-1), y= 0, yend= 0),color = "red", size = 0.5,show.legend = FALSE)+ stat_function(fun = dt, geom = "area", fill = "steelblue", xlim = c(qt(0.975,n-1),max(x)), alpha = 0.2, args = list(df = n) )+ geom_segment(aes(x = qt(0.975,n-1), xend = max(x), y= 0, yend= 0),color = "red", size = 0.5,show.legend = FALSE)+ annotate("text",x = -10,y=0.03,label = bquote("W"^"*"),size = 3,colour = "red",hjust = 0)+ geom_point(aes(x = pom$statistic,y = 0), shape = 3,size = 3) + annotate("text", x=pom$statistic, y=0.03, label= bquote("x"[obs]),colour = "black",hjust = 0.5, size = 3) ################## Přiklad 2. (Přežití) ############################################### ####################################################################################### # Test o střední hodnotě / mediánu ## Načtení dat z xlsx souboru (pomoci balíčku readxl) library(readxl) preziti = read_excel("C:/Users/Martina/Výuka/Pravděpodobnost a statistika/DATA/aktualni/testy_jednovyberove.xlsx", sheet = "preziti") colnames(preziti)="hodnoty" ## Explorační analýza boxplot(preziti$hodnoty) # Data obsahují odlehlá pozorování. Pomoci f-ce boxplot je umíme vypsat. pom=boxplot(preziti$hodnoty) pom # rozhodli-li jsme se pro odstranění odlehlých hodnot, pak preziti$hodnoty.bez=preziti$hodnoty # doporučujeme nepřepisovat původní data preziti$hodnoty.bez[preziti$hodnoty %in% pom$out]=NA ## Explorační analýza pro data bez odlehlých pozorování boxplot(preziti$hodnoty.bez) length(na.omit(preziti$hodnoty.bez)) # sd zaokrouhlujeme na 3 platné cifry sd(preziti$hodnoty.bez,na.rm=TRUE) # sd a míry polohy zaokrouhlujeme na desetiny summary(preziti$hodnoty.bez,na.rm=TRUE) # Ověření normality qqnorm(preziti$hodnoty.bez) qqline(preziti$hodnoty.bez) hist(preziti$hodnoty.bez) library(moments) skewness(preziti$hodnoty.bez,na.rm=TRUE) moments::kurtosis(preziti$hodnoty.bez,na.rm=TRUE)-3 # QQ - graf i histogram ukazují, že výběr pravděpodobně není výběrem z normálního rozdělení. # Šikmost i špičatost odpovídá norm. rozdělení. Pro konečné rozhodnutí o normalitě dat použijeme # test normality. # Předpoklad normality ověříme Shapirovovým . Wilkovovým testem. shapiro.test(preziti$hodnoty.bez) # Na hl. významnosti 0.05 zamítáme předpoklad normality # (Shapirův-Wilkův test, p-hodnota<<0,001). # H0: median = 22,2 měsíců # Ha: median > 22,2 měsíců library(BSDA) SIGN.test(preziti$hodnoty.bez,md=22.2, alternative="greater",conf.level=0.95) # Na hl. významnosti 0,05 nelze zamítnout nulovou hypotézu # (Znaménkový test, s = 40, p-hodnota = 0,939). # Medián doby přežití není statisticky významně větší než 22,2 měsíců. # Doplnění plynoucí z bodového a intervalového odhadu mu. # V populaci očekáváme medián doby přežití cca 20,0 měsíců. # Dle 95% levostranného intervalového odhadu lze # medián doby přežití očekávat větší než 16,9 měsíce.(I z tohoto lze vyvodit, že se # medián doby přežití statisticky významně nepřevyšuje 22 měsíců. ################## Přiklad 3. (Kroužky) ############################################### ####################################################################################### # Test o směrodatné odchylce # Předpokládáme normalitu dat (dle zadání) n = 80 # rozsah souboru s = 0.04 # mm .... výběrová směrodatná odchylka (bodový odhad sm. odchylky) # H0: sigma = 0.05 mm # Ha: sigma < 0.05 mm x.obs = (n-1)*s^2/0.05^2 x.obs p.hodnota = pchisq(x.obs,n-1) p.hodnota # Na hladině významnosti 0,05 zamítáme nulovou hypotézu ve prospěch alternativní # (test o rozptylu, X2 = 50,56,df = 79, p-hodnota = 0,005). # Směrodatná odchylka průměru kroužku je statisticky významně menší než 0,05 mm. ################## Přiklad 4. (Kroužky - data) ######################################## ####################################################################################### ## Načtení dat z xlsx souboru (pomoci balíčku readxl) library(readxl) krouzky = read_excel("C:/Users/Martina/Výuka/Pravděpodobnost a statistika/DATA/aktualni/testy_jednovyberove.xlsx", sheet = "krouzky") colnames(krouzky)="hodnoty" ## Explorační analýza boxplot(krouzky$hodnoty) # Data obsahují odlehlá pozorování. Pomoci f-ce boxplot je umíme vypsat. pom=boxplot(krouzky$hodnoty) pom # rozhodli-li jsme se pro odstranění odlehlých hodnot, pak krouzky$hodnoty.bez = krouzky$hodnoty krouzky$hodnoty.bez[krouzky$hodnoty %in% pom$out]=NA ## Explorační analýza pro data bez odlehlých pozorování boxplot(krouzky$hodnoty.bez) length(na.omit(krouzky$hodnoty.bez))# sd zaokrouhlujeme na 3 platné cifry sd(krouzky$hodnoty.bez,na.rm=TRUE) # sd a míry polohy zaokrouhlujeme na tisíciny summary(krouzky$hodnoty.bez,na.rm=TRUE) # Pokud je třeba, nastavíme požadovaný formát vypočtených číselných charakteristik options(scipen = 100, # 100 (desetinné číslo), -100 (vědecká notace) digits = 7) # max. počet platných cifer mean(krouzky$hodnoty.bez,na.rm=TRUE) # Ověření normality qqnorm(krouzky$hodnoty.bez) qqline(krouzky$hodnoty.bez) hist(krouzky$hodnoty.bez) library(moments) skewness(krouzky$hodnoty.bez,na.rm=TRUE) moments::kurtosis(krouzky$hodnoty.bez,na.rm=TRUE)-3 # Šikmost i špičatost odpovídá norm. rozdělení. Pro konečné rozhodnutí o normalitě dat použijeme # test normality. # Předpoklad normality ověříme Shapirovovým . Wilkovovým testem. shapiro.test(krouzky$hodnoty.bez) # Na hl. významnosti 0,05 nelze předpoklad normality zamítnout # (Shapirův-Wilkův test, W = 0,970, p-hodnota=0,063). # H0: sigma = 0,05 mm # Ha: sigma < 0,05 mm library(EnvStats) varTest(krouzky$hodnoty.bez,sigma.squared = 0.05^2, alternative = "less") # Na hladině významnosti 0,05 zamítáme H0 ve prospěch Ha # (test o rozptylu, X = 19,01, df = 77, p-hodnota < 0,001). # Jak najít 95% intervalový odhad směrodatné odchylky? pom = varTest(krouzky$hodnoty.bez,sigma.squared = 0.05^2, alternative = "less") names(pom) sqrt(pom$estimate) sqrt(pom$conf.int) # Doplnění plynoucí z bodového a intervalového odhadu mu. # V populaci očekáváme směrodatnou odchylku průměru kroužků menší než cca 0,025 mm. # Dle 95% pravostranného intervalového odhadu lze # směrodatnou odchylku průměru kroužků očekávat menší než 0,029 mm.(I z tohoto lze vyvodit, že se # směrodatná odchylka průměru kroužků je statisticky významně menší než 0,04 mm. # Pro zájemce: Vizualizace klasického testu. n = length(krouzky$hodnoty) x = seq(0.01,150,0.1) f = dchisq(x,n-1) df = data.frame(x,f) pom = varTest(krouzky$hodnoty.bez,sigma.squared = 0.05^2, alternative = "less") library(ggplot2) ggplot(df,aes(x = x,y = f))+ geom_line()+ theme_bw()+ ylab("hustota p-sti")+ xlab("pozorovaná hodnota")+ stat_function(fun = dchisq, geom = "area", fill = "steelblue", xlim = c(min(x),qchisq(0.05,n-1)), alpha = 0.2, args = list(df = n) )+ geom_segment(aes(x = min(x), xend = qchisq(0.05,n-1), y= 0, yend= 0),color = "red", size = 0.5,show.legend = FALSE)+ annotate("text",x = 30,y=0.003,label = bquote("W"^"*"),size = 3,colour = "red",hjust = 0)+ geom_point(aes(x = pom$statistic,y = 0), shape = 3,size = 3) + annotate("text", x=pom$statistic, y=0.003, label= bquote("x"[obs]),colour = "black",hjust = 0.5, size = 3) ################## Přiklad 5. (Rezistory) ############################################# ####################################################################################### n = 1000 # rozsah výběru x = 15 # počet "úspěchů" p = x/n # relativní četnost (bodový odhad pravděpodobnosti) p # Ověření předpokladů n > 9/(p*(1-p)) # Dále předpokládáme n/N < 0.05, tj. že daná populace (rezistorů) má rozsah # alespoň 1000/0.05 = 1000*20 = 20 000 rezistorů, což je asi vcelku reálný předpoklad :-) ## Clopperův - Pearsonův (exaktní) test ## H0: pi = 0.01 ## Ha: pi > 0.01 binom.test(x,n,0.01,alternative="greater",conf.level=0.95) # Na hladině významnosti 0,05 nezamítáme H0 # (Clopperův - Pearsonův test, p-hodnota = 0,082). # Nelze očekávat, že podíl vadných rezistorů ve výrobě statisticky významně # převyšuje 1 %. # Doplnění plynoucí z bodového a intervalového odhadu pi. # V produkci firmy TT očekáváme cca 1,5 % vadných rezistorů. # Dle 95% levostranného Clopperova - Pearsonova intervalového odhadu lze # v produkci firmy TT očekávat více než 0,9 % vadných rezistorů.(I z tohoto lze vyvodit, že # nelze tvrdit, že firma TT má ve své produkci více než 1 % vadných rezistorů.