#---------------------------------------------------------------------- # t-test for H0: muD=mux-muy<=0 mux <- 0 muy <- 0.2 sigmax <- sigmay <- 1 n <- 50 x <- rnorm(n,mean=mux,sd=sigmax) y <- rnorm(n,mean=muy,sd=sigmay) test <- t.test(x,y,paired=FALSE,alternative="greater") test #t.test zu fuss sp <- ((n-1)*var(x)+(n-1)*var(y))/(2*n-2) test.stat <- (mean(x)-mean(y))/(sqrt(sp *(1/n+1/n))) test.stat 1-pt(test.stat,df=2*n-2) #systematic: R <- 10000 error <- rep(0,R) for(i in 1:R){ mux <- 0 muy <- 0.2 sigmax <- sigmay <- 1 n <- 50 x <- rnorm(n,mean=mux,sd=sigmax) y <- rnorm(n,mean=muy,sd=sigmay) test <- t.test(x,y,paired=FALSE,alternative="greater") test if(test$p.value<0.05){error[i] <- 1} } mean(error) #permutation test for mu_D<=0 mux <- 0; muy <- 0.2 n <- 50; sigma <- 1 x.orig <- rnorm(n,mean=mux,sd=sigma) y.orig <- rnorm(n,mean=muy,sd=sigma) diff.orig <- mean(x.orig)-mean(y.orig) all.orig <- c(x.orig,y.orig) test <- t.test(x.orig,y.orig,paired=FALSE,alternative="greater") R <- 1000 diff.perm <- rep(0,R) for(i in 1:R){ all.perm <- sample(all.orig,size=2*n,replace = FALSE) x.perm <- all.perm[1:n] y.perm <- all.perm[(n+1):(2*n)] diff.perm[i] <- mean(x.perm)-mean(y.perm) } greater <- ifelse(diff.perm>=diff.orig,1,0) p.value <- mean(greater) p.value test$p.value #systematic outer.R <- 1000 outer.tp <- rep(0,outer.R) outer.p <- rep(0,outer.R) for(j in 1:outer.R){ mux <- 0; muy <- 0.2 n <- 50; sigma <- 1 x.orig <- rnorm(n,mean=mux,sd=sigma) y.orig <- rnorm(n,mean=muy,sd=sigma) diff.orig <- mean(x.orig)-mean(y.orig) all.orig <- c(x.orig,y.orig) outer.tp[j] <- t.test(x.orig,y.orig,paired=FALSE,alternative="greater")$p.value R <- 1000 dist.perm <- rep(0,R) for(i in 1:R){ all.perm <- sample(all.orig,size=2*n,replace = FALSE) x.perm <- all.perm[1:n] y.perm <- all.perm[(n+1):(2*n)] diff.perm[i] <- mean(x.perm)-mean(y.perm) } greater <- ifelse(diff.perm>=diff.orig,1,0) p.value <- mean(greater) outer.p[j] <- p.value } reject.rate.t <- mean(ifelse(outer.tp<0.05,1,0)) reject.rate.t reject.rate <- mean(ifelse(outer.p<0.05,1,0)) reject.rate A <- data.frame(p.value.t=outer.tp,p.value.perm=outer.p) library(ggplot2) p <- ggplot(data=A,aes(x=p.value.t,y=p.value.perm)) p <- p + geom_hline(yintercept = 0.05,color="magenta") p <- p + geom_vline(xintercept = 0.05,color="magenta") p <- p + geom_point() p <- p + theme_bw() p