num.trt <- 34 mean.trt <- c(22.5,rep(32.5,8),rep(37.5,7),42.5,rep(47.5,7),rep(52.5,8), 55.5,55.5) sd.trt <- rep(1.42,num.trt) num.rep <- 5 t.out <- matrix(rep(0,2*num.trt*(num.trt-1)/2),ncol=2) pair.num <- 0 for (i1 in 1:(num.trt-1)) { for (i2 in (i1+1):num.trt) { pair.num <- pair.num + 1 y1 <- rnorm(num.rep,mean.trt[i1],sd.trt[i1]) y2 <- rnorm(num.rep,mean.trt[i2],sd.trt[i2]) mean1 <- mean(y1); sd1 <- sd(y1) mean2 <- mean(y2); sd2 <- sd(y1) sd_p <- sqrt((sd1^2+sd2^2)/2) t.out[pair.num,1] <- (mean1-mean2)/(sd_p*sqrt(2/num.rep)) t.out[pair.num,2] <- mean.trt[i1]-mean.trt[i2] } } t.out