### MATEMATICKA STATISTIKA 1 ### R cviceni c. 4 (cviceni c. 10) ### ### Dvouvyberovy problem s kvantitativni odezvou ##### Nastaveni pracovniho adresare ##### ================================================= setwd("H:/nmsa331") rm(list=ls()) ### alpha pro testy a intervaly spolehlivosti alpha <- 0.05 ############ Iq = read.table("Iq2.txt",header=T) summary(Iq) attach(Iq) # vizualizace kategorialnich velicin plot(FPohlavi) # jine moznosti pie(table(FPohlavi),col=c("lightblue","pink")) barplot(table(FPohlavi),col=c("lightblue","pink")) ### Zakladni popisne statistiky pro tento problem ### ----------------------------------------------- summary(IQ) sd(IQ) # 3. IQchlapci=IQ[FPohlavi=="chlapec"] IQdivky=IQ[FPohlavi=="divka"] # char. polohy dle pohlavi summary(IQchlapci) summary(IQdivky) # nebo rychleji pomoci 1 prikazu tapply(IQ,FPohlavi,summary) # podminene sd tapply(IQ,FPohlavi,sd) # co usuzujete z popisnych statistik? # 4. # graficke znazorneni problemu: boxplot(IQ~FPohlavi) # totez jako boxplot(IQchlapci,IQdivky,names=c("chlapci","divky")) # muzeme vybarvit a doplnit popisky: boxplot(IQ~FPohlavi,ylab="IQ",col=c("lightblue","pink")) # 5. # uvazime normalitu zkoumanych dvou rozdeleni: par(mfrow=c(1,2)) hist(IQchlapci,prob=T,xlim=c(60,150)) hist(IQdivky,prob=T,xlim=c(60,150)) qqnorm(IQchlapci,main="Chlapci") qqline(IQchlapci) qqnorm(IQdivky,main="Divky") qqline(IQdivky) par(mfrow=c(1,1)) ###### Kolmogorovuv-Smirnovuv dvouvyberovy test # 6. ### ks.test(IQdivky, IQchlapci, exact = FALSE) # 7. e1=ecdf(IQchlapci) e2=ecdf(IQdivky) max(abs(e1(IQ)-e2(IQ))) IQ[which.max(abs(e1(IQ)-e2(IQ)))] source("figKS.R") figKS(IQ, FPohlavi) ##### ========================================================================================== ##### ##### A. Dvouvyberovy t-test se shodou rozptylu ##### Testuje rovnost strednich hodnot. ##### ##### Oznacme mu1 = E(IQ | divka) ##### mu2 = E(IQ | chlapec) ##### ##### Testujme H0: mu1 = mu2 ##### H1: mu1 != mu2 ##### ========================================================================================== # 8. t.test(IQ ~ FPohlavi, var.equal = TRUE) ### "vlnkovy" zapis t.test(IQchlapci,IQdivky, var.equal = TRUE) ### totez zadanim dvou vektoru obsahujicich dva porovnavane vybery # vyzkousejte t.test(IQdivky,IQchlapci, var.equal = TRUE) # 9. # vypocet testove statistiky xmean=mean(IQchlapci) ymean=mean(IQdivky) nx=length(IQchlapci) ny=length(IQdivky) S2=1/(nx+ny-2)*((nx-1)*var(IQchlapci)+(ny-1)*var(IQdivky)) (T.stat=(xmean-ymean)/(sqrt(S2*(1/nx+1/ny)))) # kriticka hodnota qt(1-alpha/2,df=nx+ny-2) # p-hodnota 2 * pt(-abs(T.stat), df = nx+ny-2) # interval spolehlivosti xmean-ymean + c(-1,1)*(sqrt(S2*(1/nx+1/ny)))*qt(1-alpha/2,df=nx+ny-2) ##### 11. ##### B. Welchuv dvouvyberovy t-test (bez predpokladu shodnych rozptylu) ##### t.test(IQ ~ FPohlavi) # nebo t.test(IQchlapci,IQdivky) # manualne testova statistika (T.W=(xmean-ymean)/(sqrt(var(IQchlapci)/nx+var(IQdivky)/ny))) # kriticka hodnota # stupne volnosti viz https://en.wikipedia.org/wiki/Welch%27s_t-test # t.test(IQ~FPohlavi)$parameter qt(1-alpha/2,df=t.test(IQ~FPohlavi)$parameter) # p-hodnota 2*pt(abs(T.W),df=t.test(IQ~FPohlavi)$parameter, lower=FALSE) #### # 12. # rozmyslete si, proc zadavame prave takto t.test(IQdivky,IQchlapci,mu=5,alternative="greater") t.test(IQ ~ FPohlavi, mu=-5,alternative="greater") # testova statistika (T.stat = (ymean-xmean- 5)/(sqrt(var(IQchlapci)/nx+var(IQdivky)/ny))) # 13. # porovnani sily opak=1000 n1=100 n2=150 p.s=numeric(opak) p.w=numeric(opak) for(i in 1:opak){ x=rnorm(n1,0.2,sd=1) y=rnorm(n2,0,sd=1) p.s[i]=t.test(x,y,var.equal=TRUE)$p.val p.w[i]=t.test(x,y)$p.val } mean(p.s<=0.05) mean(p.w<=0.05) op = par(mfrow=c(1,2)) hist(p.s,prob=TRUE) abline(h=1,col=2,lwd=2) hist(p.w,prob=TRUE) abline(h=1,col=2,lwd=2) par(op) ##### 14. ##### Wilcoxonuv rank-sum test = Mannuv-Whitneyuv test ##### ##### Jak presne vypadaji testovane hypotezy? ##### Na zaklade jakeho pravdepodobnostniho modelu lze tento test pouzit? ##### ============================================================================================ ## asymptoticky test: wilcox.test(IQ~FPohlavi, correct = FALSE) # pokus o presny test wilcox.test(IQ~FPohlavi, exact=TRUE) # manualni vypocet wilcox.test(IQ~FPohlavi,correct=FALSE)$stat # 15. r=rank(IQ) (W1=sum(r[FPohlavi=="chlapec"])) (W2=sum(r[FPohlavi=="divka"])) # 16. ##### Spocitame Mann-Whitneyovu statistiku # vypocet pres for cyklus: U1=0;U2=0 for(i in 1:nx) for(j in 1:ny){ U1=U1+I(IQchlapci[i]IQdivky[j])+1/2*I(IQchlapci[i]==IQdivky[j]) } U1 U2 ## to same rychleji sum(outer(IQchlapci,IQdivky,function(x,y) (xy)*1 + (x==y)*.5)) # pres matice x1 <- matrix(rep(IQchlapci, ny), ncol = ny) # ve sloupcich IQchlapci x2 <- matrix(rep(IQdivky, nx), ncol = ny, byrow = TRUE) # v radcich IQdivky head(x1) head(x2) # (U1 <- sum(x1 < x2) + 0.5 * sum(x1 == x2)) ### Pocet konkordantnich paru plus polovina poctu paru se shodou (U2 <- sum(x1 > x2) + 0.5 * sum(x1 == x2)) ### Pocet diskordantnich paru plus polovina poctu paru se shodou # 17. # vztah mezi Wilcoxonovou statistikou a M-W statistikou: (nx * ny + 0.5 * nx * (nx + 1) - W1) ### U1 (nx * ny + 0.5 * ny * (ny + 1) - W2) ### U2 # 18. ## simulace (funkce simuluj z "figKS.R", nactena prikazem source("figKS.R")) simuluj simuluj(n=10, m=30, scenar="normal-posun") simuluj(n=100, m=300, scenar="normal-posun") simuluj(n=10, m=30, scenar="normal-hetero") simuluj(n=100, m=300, scenar="normal-hetero") simuluj(n=10, m=30, scenar="normal-exp") simuluj(n=100, m=300, scenar="normal-exp") simuluj(n=10, m=30, scenar="normal-exp2") simuluj(n=100, m=300, scenar="normal-exp2") ##### ##### ##### ### Samostatna prace ##### ##### # 1. grafy op = par(mfrow=c(1,2)) hist(ZN8[FPohlavi=="chlapec"],main="Znamky chlapci") hist(ZN8[FPohlavi=="divka"],main="Znamky divky") par(op) boxplot(ZN8~FPohlavi) # vyber testu a jeho provedeni -> samostatne #### 3. # p-hodnota Wilcoxonova testu # zkusime, zda sedi alespon p-hodnota: nx=length(IQchlapci) ny=length(IQdivky) EW=nx*(nx+ny+1)/2 # varW.bez=nx*ny*(nx+ny+1)/12 #vypocet korekce t=table(r) kor=sum(t^3-t)/((nx+ny)*(nx+ny-1)) varW.kor=nx*ny*(nx+ny+1-kor)/12 U=(W1-EW)/sqrt(varW.kor) 2*pnorm(-abs(U)) # porovname: wilcox.test(IQ~FPohlavi,correct=FALSE)$p.val