Randomization demo

This is an R notebook that contains the code I used in the Randomization and sample size post in my series on Causal inference in ecology.

First, we set the default parameters for the simulation: the allele frequency at the locus influencing yield (p), the mean yield for high and low genotypes (hi, lo: high is dominant to low), the standard deviation of yield around the mean (sd), the number of simulations at each of the sample sizes (n_sim, and the P-value threshold for statistical significance (thresh).

## allele frequency
##
p <- sqrt(0.5)
## yield
##
hi <- 1.0
lo <- 0.5
## standard deviation
##
sd <- 0.1
## number of simulations at each sample size
##
n_sim <- 10000
## P-value threshold for significance
##
thresh <- 0.05

The we define the functions simulate() and report() that perform the simulation for one sample size and report the results, respectively.

simulate <- function(n_sim, n, p, hi, lo, sd, thresh) {
  hi_freq <- p^2 + 2*p*(1.0-p)
  lo_ct <- 0
  hi_ct <- 0
  ## n_hi is the total number of high-yield plants to be sampled across
  ## the two treatments
  ##
  ## 2*n because the sample size in each treatment is n
  ##
  n_hi <- rbinom(n_sim, 2*n, hi_freq)
  for (i in 1:n_sim) {
    ## select high-yield plants for high-N treatment
    ##
    n_hi_in_hi_N <- rhyper(1, n_hi[i], 2*n-n_hi[i], n)
    ## high-yield plants in low-N treatment is simply the number left
    ## out of n_hi
    ##
    n_hi_in_lo_N <- n_hi[i] - n_hi_in_hi_N
    ## low-yield plants in high-N
    ##
    n_lo_in_hi_N <- n - n_hi_in_hi_N
    ## low-yield plants in low-N
    ##
    n_lo_in_lo_N <- n - n_hi_in_lo_N
    hi_n <- c(rnorm(n_hi_in_hi_N, hi, sd),
              rnorm(n_lo_in_hi_N, lo, sd))
    lo_n <- c(rnorm(n_hi_in_lo_N, hi, sd),
              rnorm(n_lo_in_lo_N, lo, sd))
    test <- t.test(hi_n, lo_n)
    if ((test$statistic < 0.0) && (test$p.value < thresh)) {
      lo_ct <- lo_ct + 1
    } else if ((test$statistic > 0.0) && (test$p.value < thresh)) {
      hi_ct <- hi_ct + 1
    }
  }
  return(list(n=n,
              lo_ct=lo_ct,
              hi_ct=hi_ct))
}
report <- function(result, n_sim) {
  cat("Sample size: ", result$n, "\n",
      "         lo: ", result$lo_ct, "\n",
      "         hi: ", result$hi_ct, "\n",
      "    neither: ", n_sim - result$lo_ct - result$hi_ct, "\n")
}

With those in place, we’re ready to run the simulations and report the results. Notice that I’ve hardcoded small, medium, and large as 5, 10, and 20 (the second argument in the call to simulate).

small <- simulate(n_sim, 5, p, hi, lo, sd, thresh)
medium <- simulate(n_sim, 10, p, hi, lo, sd, thresh)
large <- simulate(n_sim, 20, p, hi, lo, sd, thresh)
report(small, n_sim)

Sample size:  5 
          lo:  141 
          hi:  135 
     neither:  9724 

report(medium, n_sim)

Sample size:  10 
          lo:  177 
          hi:  195 
     neither:  9628 

report(large, n_sim)

Sample size:  20 
          lo:  211 
          hi:  220 
     neither:  9569 

I produced the webpage you see here simply by selecting “Run->Restart R and Run All Chunks”.