# t-test for H0: muD=mux-muy=0 mux <- 0 muy <- 0.2 sigmax <- sigmay <- 1 n <- m <- 50 x <- rnorm(n,mean=mux,sd=sigmax) y <- rnorm(m,mean=muy,sd=sigmay) test <- t.test(x,y,paired=FALSE,alternative="two.sided",var.equal=TRUE) test #confidence interval for muD manually alpha <- 0.05 sp <- ((n-1)*var(x)+(m-1)*var(y))/(n+m-2) Delta <- qt(p=1-alpha/2,df=n+m-2)*sqrt(sp *(1/n+1/m)) conf.int <- c(mean(x)-mean(y) - Delta, mean(x)-mean(y) + Delta) test$conf.int[1:2] conf.int #check if the interval does what it should (contain the true parameter in 100*(1-alpha) % #systematic: R <- 10000 error <- rep(0,R) CI <- data.frame(lower=rep(0,R),upper=rep(0,R)) for(i in 1:R){ mux <- 0 muy <- 0.2 sigmax <- sigmay <- 1 n <- m <- 50 x <- rnorm(n,mean=mux,sd=sigmax) y <- rnorm(m,mean=muy,sd=sigmay) test <- t.test(x,y,paired=FALSE,alternative="two.sided",var.equal=TRUE) CI[i,] <- test$conf.int[1:2] } CI$contained <- ifelse(CI$lower<= mux- muy & CI$upper>= mux- muy,1,0) coverage <- mean(CI$contained) coverage #check the interrelation between CI and p-value of the two-sided test #systematic: R <- 10000 error <- rep(0,R) CI <- data.frame(lower=rep(0,R),upper=rep(0,R),p=rep(0,R)) for(i in 1:R){ mux <- 0 muy <- 0.2 sigmax <- sigmay <- 1 n <- m <- 50 x <- rnorm(n,mean=mux,sd=sigmax) y <- rnorm(m,mean=muy,sd=sigmay) test <- t.test(x,y,paired=FALSE,alternative="two.sided",var.equal=TRUE) CI[i,1:2] <- test$conf.int[1:2] CI$p[i] <- test$p.value } CI$contained <- ifelse(CI$lower<= 0 & CI$upper>= 0,1,0) CI$ok <- ifelse((CI$contained==1 & CI$p>=0.05)|(CI$contained==0 & CI$p<0.05),1,0) summary(CI) #bootstrap CI mux <- 0 muy <- 0.2 sigmax <- sigmay <- 1 n <- m <- 50 x <- rnorm(n,mean=mux,sd=sigmax) y <- rnorm(m,mean=muy,sd=sigmay) test <- t.test(x,y,paired=FALSE,alternative="two.sided",var.equal=TRUE) #just test$conf.int[1:2] boot.diff <- rep(0,R) for(i in 1:R){ x.boot <- sample(x,size = n,replace = TRUE) y.boot <- sample(y,size = m,replace = TRUE) boot.diff[i] <- mean(x.boot)-mean(y.boot) } ci.boot <- as.numeric(quantile(boot.diff,probs = c(alpha/2,1-alpha/2))) #compare the two intervals test$conf.int[1:2] ci.boot #systematical comparison of the two CIs outer.R <- 100 Results <- data.frame(lower.t=rep(0,outer.R),lower.boot=rep(0,outer.R),upper.t=rep(0,outer.R),upper.boot=rep(0,outer.R)) for(k in 1:outer.R){ mux <- 0; muy <- 0.2 sigmax <- sigmay <- 1 n <- m <- 50 x <- rnorm(n,mean=mux,sd=sigmax) y <- rnorm(m,mean=muy,sd=sigmay) test <- t.test(x,y,paired=FALSE,alternative="two.sided",var.equal=TRUE) #just Results[k,c(1,3)] <- test$conf.int[1:2] R <- 1000; boot.diff <- rep(0,R) for(i in 1:R){ x.boot <- sample(x,size = n,replace = TRUE) y.boot <- sample(y,size = m,replace = TRUE) boot.diff[i] <- mean(x.boot)-mean(y.boot) } Results[k,c(2,4)] <- as.numeric(quantile(boot.diff,probs = c(alpha/2,1-alpha/2))) } View(Results) #graphical overview library(ggplot2) library(reshape2) Results$i <- 1:nrow(Results) A1 <- subset(Results,select = c(1,3,5)) A1$type <- "ci.t" names(A1)[1:2] <- c("lower","upper") A2 <- subset(Results,select = c(2,4,5)) A2$type <- "ci.boot" names(A2)[1:2] <- c("lower","upper") A <- rbind(A1,A2) A <- A[order(A$i),] p <- ggplot(data=A, aes(x=i, y=upper, colour=type, group=type)) p <- p + geom_hline(yintercept = -0.2) p <- p + geom_errorbar(aes(ymin=lower, ymax=upper),width=1) p <- p + scale_color_manual(values=c("gray","magenta")) p <- p + scale_x_continuous(name="run") p <- p + scale_y_continuous(name="") p <- p + theme_bw() p