# R for simulation study # t-test example out <- t.test(extra ~ group, data = sleep, var.equal = TRUE) names(out) out[["p.value"]] out[["p.value"]] < 0.05 # Aggregate function aggregate(extra ~ group, data = sleep, FUN = mean) aggregate(extra ~ group, data = sleep, FUN = sd) # Data generation set.seed(1234) y1 <- rnorm(100, 0, 1) y2 <- rnorm(100, 0, 1) y <- c(y1, y2) g <- c(rep(1, 100), rep(2, 100)) dat <- data.frame(g = g, y = y) out <- t.test(y ~ g, data = dat, var.equal = TRUE) out[["p.value"]] # for loop set.seed(1234) nrep <- 1000 result <- rep(NA, nrep) for(i in 1:nrep) { y1 <- rnorm(100, 0, 1) y2 <- rnorm(100, 0, 1) y <- c(y1, y2) g <- c(rep(1, 100), rep(2, 100)) dat <- data.frame(g = g, y = y) out <- t.test(y ~ g, data = dat, var.equal = TRUE) result[i] <- out[["p.value"]] } mean(result < 0.05) # Specify design conditions set.seed(1234) nrep <- 1000 n <- c(100, 20) v <- c(1, 10) result <- rep(NA, nrep) for(i in 1:nrep) { y1 <- rnorm(n[1], 0, sqrt(v[1])) y2 <- rnorm(n[2], 0, sqrt(v[2])) y <- c(y1, y2) g <- c(rep(1, n[1]), rep(2, n[2])) dat <- data.frame(g = g, y = y) out <- t.test(y ~ g, data = dat, var.equal = TRUE) result[i] <- out[["p.value"]] } mean(result < 0.05) # Nested for loop ns <- matrix(NA, 5, 2) ns[,1] <- c(20, 40, 60, 80, 100) ns[,2] <- c(100, 80, 60, 40, 20) vs <- matrix(1, 3, 2) vs[,2] <- c(1, 3, 10) ovresult <- NULL set.seed(1234) nrep <- 1000 for(k in 1:nrow(ns)) { n <- ns[k,] for(j in 1:nrow(vs)) { v <- vs[j,] result <- rep(NA, nrep) for(i in 1:nrep) { y1 <- rnorm(n[1], 0, sqrt(v[1])) y2 <- rnorm(n[2], 0, sqrt(v[2])) g <- c(rep(1, n[1]), rep(2, n[2])) dat <- data.frame(g = g, y = c(y1, y2)) out <- t.test(y ~ g, data = dat, var.equal = TRUE) result[i] <- out[["p.value"]] } result <- cbind(k, j, result) ovresult <- rbind(ovresult, result) } } head(ovresult) # Summarize the results colnames(ovresult) <- c("n", "v", "p") ovresult <- data.frame(ovresult, sig = ovresult[,"p"] < 0.05) ovresult[,"n"] <- factor(ovresult[,"n"], labels = c("20,100", "40,80", "60,60", "80,40", "100,20")) ovresult[,"v"] <- factor(ovresult[,"v"], labels = c("1,1", "1,3", "1,10")) aggregate(sig ~ n + v, data = ovresult, FUN = mean) # Exercise code y1 <- rbinom(60, 1, 0.2) y2 <- rbinom(60, 1, 0.2) yes1 <- sum(y1) yes2 <- sum(y2) total1 <- length(y1) total2 <- length(y2) out <- prop.test(x = c(yes1, yes2), n = c(total1, total2)) out[["p.value"]]