Introduction to simulation for experimental design using R

quantitative phenotype differs significantly between European-derived and Asian-derived populations

Common genetic variants account for differences in gene expression among ethnic groups

Interestingly, the arrays used to measure expression for the CEU [European-derived] individuals were primarily processed from 2003 to 2004,whereas the arrays used to measure expression for the ASN [Asian-derived] individuals were all processed in 2005–2006.

.

microarray batch effects appear to be completely confounded with population effects.”

On the design and analysis of gene expression studies in human populations

reply

Accounts for variability

• replication
• randomization
• blocking
  • stratification
  • clustering
    

• Feasible!
• Results can be analyzed!

Task view

Park grass :

• What questions can you ask with the park grass experiment as it was designed?

Redesign park grass experiment:

1. Given the experimental units provided (22 rows x 4 columns + 11 cells= 99) which of the park grass experiment questions (effect of N, P, Na, etc.) would you address and how?
2. Given the questions park grass experiment, come up with a minimal design that addresses the historic set of questions/hypothesis

• Maybe start by listing the treatments and or associated questions. This could be done by specifying the complete model in R style formula syntax.

** WHat other famous failures of experimental design do you know? **

Analysis despite flaws in experimental design:
Determinants of Species Richness in the Park Grass Experiment

Park grass inspired experiment using proper experimental design:
Effects of grassland management on plant C:N:P stoichiometry:implications for soil element cycling and storage (example publication, dozens of related publications)

Create a spreadsheet simulation of an experiment testing :

• testing the effect of nitrogen, phosphorus and [CO2] on chlorophyll in mesocosms

and/or

• testing difference between species in how (supra-optimal or full range) temperature affects their running spead

and/or

• the effect of gear choice and location on fisherman's catch

What were issues you identified?

• limited types of distributions
• tedious setup
  • design (repetitive copy pasting)
  • multiple trials (very long sheet)
  • etc.
• challenge of scaling
• requires switching to R for downstream steps
• …

Throwing a six face die

Create the die

die <-seq(from=1, to=6, by=1)

Throw the die

sample(x=die, size=1, replace=T)

[1] 4

Throw the 1000 times and keep the results

die_results <- sample(x=die,
                      size=1000,
                      replace=T)

1. Throw a 10 faced die
2. Throw a coin
3. Throw cheaters die that produces 6 50% of the time

die10 <-seq(from=1, to=10, by=1)

coin <- c("head", "tails")

p_dice<-c(1,1,1,1,1,5)
die_results <- sample(x=die,
                      size=1000,
                      replace=T,
                      prob= p_dice)

Normal distribution

rnorm(n = number of observations,
      mean = vector of means,
      sd = vector of means)

1. Create a sample representing the number of dolphins seen per 10 hours of observation, for 100 observations (hint: assuming this is a number of independent events occurring in a fixed time)
2. Create a sample representing the frequency of species in a quadrat (hint: assuming a few species will be very common and many species will be very rare)

Additional challenges:

• Plot these distributions.
• Repeat your sampling while the distribution while varying key parameters
• Plot summary statistics of your sampling with with various parameters

2.

species <- rlnorm(n = 1000,
                  mean = 25,
                  sd=1)
hist(species)

• Import one of your spreadsheet simulation into R
• Create a plot in R of the data simulated in Excel
• In R, create a new column containing simulated values using a more sensible distribution of error

design setup “manually”

Example 2 factors, 2 levels, 3 replicates, full factorial

one_replicate <- data.frame(factorA=factor(rep(c("FacA_Level1",
                                     "FacA_Level2"),
                                     3)),
                            factorB=factor(rep(c("FacB_Level1",
                                     "FacB_Level2", "FacB_Level3"),
                                     each=2)))
head(one_replicate)                              

      factorA     factorB
1 FacA_Level1 FacB_Level1
2 FacA_Level2 FacB_Level1
3 FacA_Level1 FacB_Level2
4 FacA_Level2 FacB_Level2
5 FacA_Level1 FacB_Level3
6 FacA_Level2 FacB_Level3

Alternatives:

• if conducting full factorial, use expand.grid()
• create design in spreadsheet and import

function to create scaled numeric predictors from all factors

num_scale_fact <- function(data_frame){
  ind <- sapply(data_frame, is.factor)
  data_frame[paste0("num_",names(data_frame[ind]))] <- 
    lapply(data_frame[ind],
                                                              function(x)scale(as.numeric(x)))
  return(data_frame)
}

experiment1 <- NULL
for(i in LETTERS[1:3]){
  one_replicate$replicate <- factor(i)
  one_replicate$response <- with(one_replicate,
                        num_factorA*
                          num_factorB+rnorm(n = nrow(one_replicate), sd=2))
  experiment1 <- rbind(experiment1,one_replicate)
}

Avoid for loops in this context using the apltly named replicate

1. Plot experiment 1
2. Analyze experiment 1
3. You know that factorA and factorB have an an interactive effect, can you detect this from this experiment?
4. If you can detect an effect, what could you do to save money/time etc in this experiment?
5. If you can not detect an effect, what would you need to do to detect the known effect?

1.

require(ggplot2)
p <- qplot(data=experiment1,
           x=factorA,
           y=response,
           colour=factorB,
           geom="boxplot")
print(p)

Experimental setup using tools

require(AlgDesign)
experiment2 <- gen.factorial(levels=c(3,2,2),
                             nVars=3,
                             varNames=c("A","B","C"),
                             factors="all")

Experimental design task view

• create an expanded simulated data set and analyse each trial

• create the experimental design

experiment_design <- gen.factorial(levels=c(3,2,2),
                            nVars=3,
                            varNames=c("A","B","C"),
                            factors="all")

experiment_design <- num_scale_fact(experiment_design)

• add replication and trial (replication of the experiment)

trial_rep <- expand.grid(trial=as.factor(1:1000),
                         rep=as.factor(1:5))

expand.grid.df <- function(...) Reduce(function(...) merge(..., by=NULL), list(...))

rep_exp1 <- expand.grid.df(trial_rep,experiment_design)

https://stackoverflow.com/questions/11693599/alternative-to-expand-grid-for-data-frames

rep_exp1$response <- with(rep_exp1, num_A+num_B+num_B*num_C+num_A*num_B*num_C+rnorm(nrow(rep_exp1), sd=2))

• fit model to single experiment and use simulate

  simulate.lm : normal distribution with same MSE as original model

creat a base experiment

exp_d_rep <- do.call("rbind", replicate(3, experiment_design, simplify = FALSE))
exp_resp <- exp_d_rep
exp_resp$response <- with(exp_resp,
                       num_A+num_B+num_B*num_C+num_A*num_B*num_C+rnorm(nrow(exp_resp), sd=2))

fit <- lm(response~A*B*C, data=exp_resp)
exp_sim <- cbind(exp_d_rep,simulate(fit,nsim = 10))

require(tidyr)
exp_sim <- gather(exp_sim,contains("sim_"), key = trial ,value = response)

Sensitivity analysis of design: simulate a series of experimental design scenarios and see how it affects inference from the study.

Lets get super meta: setup an experimental design for your sensitivity analysis experiment of your experiment

sens_setup <- expand.grid(trial=as.factor(1:100),
                         replication=c(3:10),
                         error=c(0.1,0.5,1,5,10))

Setup a function

f_rep_error <- function(rep,error){
  rep_exp <- expand.grid.df(data.frame(rep=rep, replicate=1:rep),experiment_design)
  rep_exp$response <- with(rep_exp, num_A+num_B+num_B*num_C+num_A*num_B*num_C+rnorm(nrow(rep_exp), sd=error))
  return(rep_exp)
}

Apply the function to individual trials

trial_rep_error <- sens_setup %>% group_by(trial, replication, error) %>%
  do(with(.,f_rep_error(replication, error)))

Optimize the design using tools

For certain types of design, the allocation of experimental units can be optimized using established algorithsm.

For example,in some cases, a full factorial design can be reduced to a carefully selected subset of cases with only a marginal effect on estimated parameters.

see R vignette("AlgDesign")

• optBlock {AlgDesign}
• optFederov {AlgDesign}

• instead of gen.factorial use one of the optimize incomplete factorial designs using the same setup rep_exp1
  • can you detect the all the same effects?
  • are the parameter estimates similar?
• does the use of non-normaly distributed error (eg. poisson) affect your interpreation using aov?
• any difference with numeric rather than factor predictors?
• what are the effect of a non-linear interaction with one of your predictors?
• what are the effect of heteroscedasticity (eg. correlation of variance with one of the predictors)?

Dojo: deliberate practice

• pair/team coding

Alternative: pair/team coding on each other's experiments (turns focusing on each experiment in the team).

Produce an optimal design for the treatments in the Park Grass experiment

Design an experiment to:

• Detect the effect of biodiversity on ecosystem process
  • Include a species that has a high contribution to ecosystem function compared to the average of all other species
  • Include a pair of species that seperatly do not contribute to the process, but together have a strong impact on diversity
  • How many levels of biodiversity should your experiment contain?
  • How many replicate communities are needed at each level?
  • Can blocking help distinguish between composition effects (eg. sampling effect) and diversity effect?

Detecting the Effects of Diversity on Measures of Ecosystem Function: Experimental Design, Null Models and Empirical Observations

Redesign the experiment from your favourite recent publication

• Use the same number of experimental units and try to improve on the design
• Were the authors lucky? Insert problems into the study (eg. other distribution, loss of a replicate). Do you consistently get the same results.