8.7 What you can now do

Given the above code and workflow, you can now figure out how many subjects you might need to achieve 80% power, assuming a certain effect size (or a range of effect sizes as above) and assuming some specific values for the standard deviations and variance covariance matrices.

You can also study Type I error properties of your model as a function of whether the model is a maximal model or a varying intercepts only model.

Example: Compute power as a function of effect size and sample size. Note that the number of subjects has to be even, otherwise the simulation code will fail! One could put in a test for this in the code: if the number of subjects is divisible by 2, the modulo function should return 0:

10 %% 2

## [1] 0

11 %% 2

## [1] 1

## define a function for computing power
## (as a function of effect size and
## subject sample size:
compute_power <- function(b = NULL, nsubjects = 28) {
  if (nsubjects %% 2 != 0) {
    stop("No. of subjects must be divisible by 2.")
  }
  nsim <- 100
  sotvals <- rep(NA, nsim)
  failed <- rep(0, nsim)
  for (i in 1:nsim) {
    beta[2] <- b
    dat_sim <- gen_sim_lnorm2(
      nitem = 24,
      nsubj = nsubjects,
      beta = beta,
      Sigma_u = Sigma_u,
      Sigma_w = Sigma_w,
      sigma_e = sigma_e
    )

    ## no correlations estimated to avoid convergence problems:
    ## analysis done after log-transforming:
    m <- lmer(log(rt) ~ so + (so || subj) +
      (so || item), data = dat_sim, control = lmerControl(calc.derivs = FALSE))
    ## ignore failed trials
    if (any(grepl("failed to converge", m@optinfo$conv$lme4$messages))) {
      failed[i] <- 1
    } else {
      sotvals[i] <- summary(m)$coefficients[2, 3]
    }
  }

  ## power:
  paste("Power:", round(mean(abs(sotvals) > 2,
    na.rm = TRUE
  ), 2), sep = " ")
}

## usage: this will halt function
## with an error message
# compute_power(b=0.03,nsubjects=29)

compute_power(b = 0.03, nsubjects = 28)

[1] “Power: 0.04”