### MATEMATICKA STATISTIKA 1 ### R cviceni c. 3 (cviceni c. 10) ### Jednovyberove testy (znamenkovy test, Wilcoxonuv test), parove testy ### ### -------------------------------------------------- # Uvodni nastaveni setwd("H:/nmsa331") rm(list=ls()) ## vycisteni ### alpha pro testy a intervaly spolehlivosti alpha <- 0.05 Hosi=read.table("Hosi.txt",header=TRUE) head(Hosi) summary(Hosi) dim(Hosi) ####### # stejne jako minule si vybereme 100 hodnot porodnich hmotnosti # ale ted trochu jine hodnoty (jine seed) set.seed(04122018); # nastavuje generator (pseudo)-nahodnych cisel hmot100 <- sample(Hosi$por.hmot, 100) hmot=Hosi$por.hmot n=length(hmot100) # minule: jednovyberovy t-test # zajimalo nas, zda je porovni hmotnost rovna 3300 g t.test(hmot100,mu=3300) # nyni totez pomoci znamenkoveho testu #priprava dat (nLess <- sum(hmot100 < 3300)) # pocet deti s vahou mensi nez 3300 (nMore <- sum(hmot100 > 3300)) # pocet deti s vahou vetsi nez 3300 (n3300 <- sum(hmot100 == 3300)) # pocet deti s vahou prave 3300 (N <- nLess + nMore) ### Provedeni testu (Z = (nMore - N/2) / sqrt(N*1/2*1/2)) # testova statistika qnorm(1-alpha/2) # p- hodnota: 2*(1-pnorm(abs(Z))) # coz je totez jako 2*pnorm(-abs(Z)) # pomoci funkce prop.test prop.test(x=nMore, n=N,p=1/2,correct=FALSE) prop.test(x=sum(hmot100>3300),n=sum(hmot100!=3300),p=1/2,correct=F) # vizualizace situace Y <- hmot100[hmot100 != 3300]; # vybereme hmotnosti ruzne od 3300 Y <- factor(Y > 3300) # udelame z nich 0-1 veliciny a chapeme je jako "faktor" levels(Y) <- c("< 3300", "> 3300"); summary(Y) plot(Y,col=2:3); # totez jako: barplot(table(Y),col=2:3) # nebo barplot(c(nLess, nMore), col=2:3,legend=c("< 3300", "> 3300")) pie(table(Y),col=2:3) # Proporcionalne barplot(prop.table(as.matrix(c(nLess, nMore))), col=2:3, legend=c("< 3300", "> 3300")) #### ### ------------------------------------------------- ### 3. Wilcoxonuv test (signed-rank). ### Jak presne vypadaji testovane hypotezy? ### Na zaklade jakeho pravdepodobnostniho modelu lze tento test pouzit? ### Jaky je zaver? ### ------------------------------------------------- wilcox.test(hmot100, mu = 3300,correct=F) # hodnota testove statistiky: x <- hmot100 - 3300 x <- x[x != 0] ### nebudeme pouzivat porodni hmotnosti 3300 (R <- rank(abs(x))) #(Sign <- -1*(x < 0) + 1*(x > 0)) (Splus <- sum(R[x > 0])) ### = V ve vystupu (Sminus <- sum(R[x < 0])) # za H_0 by melo byt podobne (N <- length(x)) # Vypocet p-hodnoty manualne - samostatne cviceni (v) ## stredni hodnota testove statistiky za nulove hypotezy (ESplus <- N * (N + 1) / 4) ## rozptyl testove statistiky za nulove hypotezy - bez korekce varSplus_NoTies <- N * (N + 1) * (2 * N + 1) / 24; #vypocet korekce: (tx <- table(R)) # Pocet jednotlivych hodnot (varCorrectWithTies <- sum(tx^3 - tx) / 48) # Korekce v pripade shod (varSplus <- varSplus_NoTies - varCorrectWithTies) #bez korekce (U <- (Splus - ESplus) / sqrt(varSplus)) ### varianta s korekcemi na spojitost (U.corr <- (Splus - ESplus - 0.5) / sqrt(varSplus)) ### P-hodnoty (P <- 2*pnorm(-abs(U))) (P.corr <- 2*pnorm(-abs(U.corr))) # jak by to bylo bez korekce rozptylu: (U.bez.corr <- (Splus - ESplus) / sqrt(varSplus_NoTies)) (P.bez.corr <- 2*pnorm(-abs(U.bez.corr))) # srovname s wilcox.test(hmot100, mu = 3300, correct = FALSE)$p.val ### varianta bez korekci na spojitost wilcox.test(hmot100, mu = 3300)$p.val ### varianta s korekcemi na spojitost ########### # porovnani sily testu opak=1000 n=100 p.t=numeric(opak) p.z=numeric(opak) p.w=numeric(opak) for(i in 1:opak){ x=rnorm(n,mean=0.2,sd=1) p.t[i]=t.test(x,mu=0)$p.val p.w[i]=wilcox.test(x,mu=0)$p.val p.z[i]=prop.test(sum(x>0),n,p=1/2)$p.val } mean(p.t<=0.05) mean(p.w<=0.05) mean(p.z<=0.05) ############# ########################################################################### ## PAROVY PROBLEM ====================================================##### ## ===================================================================##### ########################################################################### ## Jsou nase data v souladu s timto tvrzenim, ze chlapci behem ## prvniho roku zivota priberou v prumeru alespon 6.5 kg? ## U kazdeho ditete mame k dispozici hmotnost pri narozeni a hmotnost ## v jednom roce. Nelze pochopitelne predpokladat, ze tyto hmotnosti ## jsou nezavisle. hmot<- Hosi$por.hmot hmot1 <- Hosi$hmotnost rozdil=hmot1-hmot # A) t-test: ## H_0: E (rozdil) <= 6500, H_1: E(rozdil)>6500 # vyhodnoceni normality: qqnorm(rozdil) qqline(rozdil) hist(rozdil,prob=T) curve(dnorm(x,mean=mean(rozdil),sd=sd(rozdil)),add=T,col="red") # test t.test(rozdil, alternative="greater",mu=6500) # Alternativne t.test(hmot1, hmot, alternative="greater",mu=6500, paired=T) #pozor na poradi! ##### # B) znamenkovy test n.more=sum(rozdil>6500) n=sum(rozdil!=6500) prop.test(n.more, n, alternative = "greater") ##### # C) Wilcoxon??v test wilcox.test(rozdil, mu = 6500, alternative = "greater") # Alternativni zpusob provedeni paroveho Wilcoxonova testu wilcox.test(hmot1, hmot, paired = TRUE, mu = 6500, alternative = "greater") wilcox.test(hmot, hmot1, paired = TRUE, mu = -6500, alternative = "less") # opet pozor na poradi! ################ rozdil=Hosi$vek.otce-Hosi$vek.matky summary(rozdil) hist(rozdil,prob=T) curve(dnorm(x,mean=mean(rozdil),sd=sd(rozdil)),add=T,col="red") boxplot(rozdil) qqnorm(rozdil) qqline(rozdil) t.test(rozdil,mu=2,alternative="greater") wilcox.test(rozdil,mu=2,alternative="greater") prop.test(sum(rozdil>2),sum(rozdil!=2),p=1/2,alternative="greater")