\begin{Schunk}
\begin{Sinput}
 ##' Monte Carlo simulation for 2 group t test
 ##' @param conds = a conditions data frame
 ##' @param var.equal: should the t.test use the equal variance 
 ##' assumption or the Welch corrected calculation (if FALSE).
 ##' Note default TRUE is different from R base.
 ##' @return a matrix summarizing rejection rates
 tPowerSim <- function(conds, var.equal = TRUE){
      ## Creates data by parsing N, M and SD strings
      getTdata <- function(rep, N, M, SD) {
          Nvec <- as.numeric(unlist(strsplit(N, ":")))
          Mvec <- as.numeric(unlist(strsplit(M, ":")))
          SDvec <- as.numeric(unlist(strsplit(SD, ":")))
          
          dat <- data.frame(first = c(rep(0, times = Nvec[1]),
                                      rep(1, times = Nvec[2])))
          dat$IQ <- rnorm(sum(Nvec), m = Mvec[(dat$first + 1)], 
                          sd = SDvec[(dat$first + 1)])
          dat$IQ <- round(dat$IQ)
          attr(dat, "rep") <- rep
          attr(dat, "parms") <- c(N = N, M = M, SD = SD)
          dat
      }
      ## conducts T test, keeps only p value
      conductTtest <- function (dframe, y = "IQ", x = "first", var.equal){
          t.test(formula(paste(y, "~", x)),
                 data = dframe, var.equal = var.equal)$p.value
      }
      ## orchestrates the data pull, analysis, and summary
      runOneSim <- function(rep, N, M, SD, var.equal){
          dframe <- getTdata(rep, N = N, M = M, SD = SD)
          reslt <- conductTtest(dframe, var.equal = var.equal)
          parms <- attr(dframe, "parms")
          dframe2 <- data.frame(rep = attr(dframe, "rep"), 
                                pvalue = reslt,
                                reject.05 = if (reslt <= 0.05) 1 else 0,
                                reject.1 = if (reslt <= 0.10) 1 else 0,
                                N = parms["N"], M = parms["M"], SD = parms["SD"])
          dframe2
      }
  
      # Reads the condition matrix, runs one row from it
      runOneCondition <- function(i, conds, var.equal){
          x <- conds[i, ]
          result.list <- lapply(1:x$nReps, runOneSim, 
                                N = x$N, M = x$M, SD = x$SD, var.equal = var.equal)
          do.call("rbind", result.list)
      }
       # Run all of the rows in the condition matrix
      listofresults <- lapply(1:NROW(conds), runOneCondition, conds, var.equal = var.equal)
      stackedResults <- do.call(rbind, listofresults)
      output <- aggregate(stackedResults[, c("reject.05", "reject.1")],
                          by = list(N = stackedResults$N, SD = stackedResults$SD, M = stackedResults$M),
                          mean)
      names(output) <- c("N", "SD", "M", "reject.05.mean", "reject.1.mean")
      output
  }
\end{Sinput}
\end{Schunk}
