UP | HOME

R scripts for "Random effects structure for confirmatory hypothesis testing with mixed-effects models: Keep it maximal"

General Information

Below are R scripts that can be used to replicate the Monte Carlo simulations reported in the manuscript. Please note:

  • You will need to install the simgen package to be able to run the simulations.
  • The output from the fitstepwise function for model selection has four-dimensional data stored in two dimensional format. The four dimensions are:
    1. run (each line in the output file is a single Monte Carlo run)
    2. alpha level (the level at which inclusion of the random slope is tested)
    3. variable (the statistic output from the fitting procedure; e.g., \(t\), \(\chi^2\))
    4. the direction of model building (forward, backward)
  • To reassemble the data into a four dimensional array for ease of processing, using the following code.
x <- read.csv("somefilename.csv", header=TRUE)

ax <- array(as.matrix(x), dim=c(nrow(x), 10, 6, 2),
            dimnames=list(NULL, crit=c(.01,.05,seq(.1,.8,.1)),
              stat=c("fm","t","chi","pt","pchi","eff"),
              dir=c("forward","backward")))

Standard runs (no model selection)

library(simgen)

# number of Monte Carlo runs (nexp)
# set to small value for testing purposes
# NOTE: change this to 100000 for full replication
pmx <- createParamMx(nexp=50, h0=TRUE)

# ----------- Type I Error (H0 true) ------------

pmx[,"eff"] <- 0

# ANOVA
mcRun("fitanova", mcr.outfile="fitanova.h0.wsbi.12.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=12, wsbi=TRUE)
mcRun("fitanova", mcr.outfile="fitanova.h0.wsbi.24.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=24, wsbi=TRUE)
mcRun("fitanova", mcr.outfile="fitanova.h0.wswi.12.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=12, wsbi=FALSE)
mcRun("fitanova", mcr.outfile="fitanova.h0.wswi.24.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=24, wsbi=FALSE)

# lmer models, random intercept only
mcRun("fitlmer", mcr.outfile="fitlmerRI.h0.wsbi.12.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=12, wsbi=TRUE, ri.only=TRUE)
mcRun("fitlmer", mcr.outfile="fitlmerRI.h0.wsbi.24.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=24, wsbi=TRUE, ri.only=TRUE)
mcRun("fitlmer", mcr.outfile="fitlmerRI.h0.wswi.12.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=12, wsbi=FALSE, ri.only=TRUE)
mcRun("fitlmer", mcr.outfile="fitlmerRI.h0.wswi.24.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=24, wsbi=FALSE, ri.only=TRUE)

# lmer models, random intercept only, p-values from MCMC
mcRun("fitlmer.mcmc", mcr.outfile="fitlmerRImcmc.h0.wsbi.12.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=12, wsbi=TRUE)
mcRun("fitlmer.mcmc", mcr.outfile="fitlmerRImcmc.h0.wsbi.24.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=24, wsbi=TRUE)
mcRun("fitlmer.mcmc", mcr.outfile="fitlmerRImcmc.h0.wswi.12.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=12, wsbi=FALSE)
mcRun("fitlmer.mcmc", mcr.outfile="fitlmerRImcmc.h0.wswi.24.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=24, wsbi=FALSE)

# lmer models, no random correlation
mcRun("fitnocorr", mcr.outfile="fitnocorr.h0.wsbi.12.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=12, wsbi=TRUE)
mcRun("fitnocorr", mcr.outfile="fitnocorr.h0.wsbi.24.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=24, wsbi=TRUE)
mcRun("fitnocorr", mcr.outfile="fitnocorr.h0.wswi.12.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=12, wsbi=FALSE)
mcRun("fitnocorr", mcr.outfile="fitnocorr.h0.wswi.24.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=24, wsbi=FALSE)

# lmer models, no random correlation, p-values from MCMC
mcRun("fitnocorr.mcmc", mcr.outfile="fitnocorrMCMC.h0.wsbi.12.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=12, wsbi=TRUE)
mcRun("fitnocorr.mcmc", mcr.outfile="fitnocorrMCMC.h0.wsbi.24.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=24, wsbi=TRUE)
mcRun("fitnocorr.mcmc", mcr.outfile="fitnocorrMCMC.h0.wswi.12.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=12, wsbi=FALSE)
mcRun("fitnocorr.mcmc", mcr.outfile="fitnocorrMCMC.h0.wswi.24.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=24, wsbi=FALSE)

# lmer models, no within-unit slopes
mcRun("fitrsonly", mcr.outfile="fitrsonly.h0.wsbi.12.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=12, wsbi=TRUE)
mcRun("fitrsonly", mcr.outfile="fitrsonly.h0.wsbi.24.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=24, wsbi=TRUE)
mcRun("fitrsonly", mcr.outfile="fitrsonly.h0.wswi.12.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=12, wsbi=FALSE)
mcRun("fitrsonly", mcr.outfile="fitrsonly.h0.wswi.24.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=24, wsbi=FALSE)

# lmer models, maximal
mcRun("fitlmer", mcr.outfile="fitlmerRS.h0.wsbi.12.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=12, wsbi=TRUE)
mcRun("fitlmer", mcr.outfile="fitlmerRS.h0.wsbi.24.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=24, wsbi=TRUE)
mcRun("fitlmer", mcr.outfile="fitlmerRS.h0.wswi.12.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=12, wsbi=FALSE)
mcRun("fitlmer", mcr.outfile="fitlmerRS.h0.wswi.24.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=24, wsbi=FALSE)


# ----------- Power (H0 false) ------------

pmx[,"eff"] <- .8

# ANOVA
mcRun("fitanova", mcr.outfile="fitanova.h1.wsbi.12.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=12, wsbi=TRUE)
mcRun("fitanova", mcr.outfile="fitanova.h1.wsbi.24.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=24, wsbi=TRUE)
mcRun("fitanova", mcr.outfile="fitanova.h1.wswi.12.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=12, wsbi=FALSE)
mcRun("fitanova", mcr.outfile="fitanova.h1.wswi.24.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=24, wsbi=FALSE)

# lmer models, random intercept only
mcRun("fitlmer", mcr.outfile="fitlmerRI.h1.wsbi.12.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=12, wsbi=TRUE, ri.only=TRUE)
mcRun("fitlmer", mcr.outfile="fitlmerRI.h1.wsbi.24.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=24, wsbi=TRUE, ri.only=TRUE)
mcRun("fitlmer", mcr.outfile="fitlmerRI.h1.wswi.12.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=12, wsbi=FALSE, ri.only=TRUE)
mcRun("fitlmer", mcr.outfile="fitlmerRI.h1.wswi.24.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=24, wsbi=FALSE, ri.only=TRUE)

# lmer models, random intercept only, p-values from MCMC
mcRun("fitlmer.mcmc", mcr.outfile="fitlmerRImcmc.h1.wsbi.12.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=12, wsbi=TRUE)
mcRun("fitlmer.mcmc", mcr.outfile="fitlmerRImcmc.h1.wsbi.24.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=24, wsbi=TRUE)
mcRun("fitlmer.mcmc", mcr.outfile="fitlmerRImcmc.h1.wswi.12.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=12, wsbi=FALSE)
mcRun("fitlmer.mcmc", mcr.outfile="fitlmerRImcmc.h1.wswi.24.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=24, wsbi=FALSE)

# lmer models, no random correlation
mcRun("fitnocorr", mcr.outfile="fitnocorr.h1.wsbi.12.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=12, wsbi=TRUE)
mcRun("fitnocorr", mcr.outfile="fitnocorr.h1.wsbi.24.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=24, wsbi=TRUE)
mcRun("fitnocorr", mcr.outfile="fitnocorr.h1.wswi.12.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=12, wsbi=FALSE)
mcRun("fitnocorr", mcr.outfile="fitnocorr.h1.wswi.24.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=24, wsbi=FALSE)

# lmer models, no random correlation, p-values from MCMC
mcRun("fitnocorr.mcmc", mcr.outfile="fitnocorrMCMC.h1.wsbi.12.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=12, wsbi=TRUE)
mcRun("fitnocorr.mcmc", mcr.outfile="fitnocorrMCMC.h1.wsbi.24.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=24, wsbi=TRUE)
mcRun("fitnocorr.mcmc", mcr.outfile="fitnocorrMCMC.h1.wswi.12.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=12, wsbi=FALSE)
mcRun("fitnocorr.mcmc", mcr.outfile="fitnocorrMCMC.h1.wswi.24.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=24, wsbi=FALSE)

# lmer models, maximal
mcRun("fitlmer", mcr.outfile="fitlmerRS.h1.wsbi.12.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=12, wsbi=TRUE)
mcRun("fitlmer", mcr.outfile="fitlmerRS.h1.wsbi.24.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=24, wsbi=TRUE)
mcRun("fitlmer", mcr.outfile="fitlmerRS.h1.wswi.12.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=12, wsbi=FALSE)
mcRun("fitlmer", mcr.outfile="fitlmerRS.h1.wswi.24.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=24, wsbi=FALSE)

Model-selection runs

library(simgen)

# nexp set to small value for testing; change to 100000 for full runs
pmx <- createParamMx(nexp=10, h0=TRUE) 

# ----------- Type I Error (H0 true) ------------

pmx[,"eff"] <- 0

# WSBI design
mcRun("fitstepwise", mcr.outfile="fitstepwise.h0.wsbi.12.MS.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=12, wsbi=TRUE,
      mf=modSpace(wsbi=TRUE))
mcRun("fitstepwise", mcr.outfile="fitstepwise.h0.wsbi.24.MS.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=24, wsbi=TRUE,
      mf=modSpace(wsbi=TRUE))

# WSWI design
# Note: fitstepwise does both forward and backward selection
# forward subject slope first & backward item slope first
mf <- modSpace(FALSE)   # define the models to be tested
mcRun("fitstepwise", mcr.outfile="fitstepwise.h0.wswi.12.FSBI.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=12, wsbi=FALSE,
      mf=c(mf[3], mf[[2]][1], mf[1]))
mcRun("fitstepwise", mcr.outfile="fitstepwise.h0.wswi.24.FSBI.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=24, wsbi=FALSE,
      mf=c(mf[3], mf[[2]][1], mf[1]))
# forward item slope first & backward subject slope first
mcRun("fitstepwise", mcr.outfile="fitstepwise.h0.wswi.12.FIBS.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=12, wsbi=FALSE,
      mf=c(mf[3], mf[[2]][2], mf[1]))
mcRun("fitstepwise", mcr.outfile="fitstepwise.h0.wswi.24.FIBS.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=24, wsbi=FALSE,
      mf=c(mf[3], mf[[2]][2], mf[1]))
# 'best path' forward
mcRun("fitstepwise.bestpath", mcr.outfile="fitstepwise.h0.wswi.12.FB.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=12, wsbi=FALSE,
      forward=TRUE)
mcRun("fitstepwise.bestpath", mcr.outfile="fitstepwise.h0.wswi.24.FB.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=24, wsbi=FALSE,
      forward=TRUE)
# 'best path' backward
mcRun("fitstepwise.bestpath", mcr.outfile="fitstepwise.h0.wswi.12.BB.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=12, wsbi=FALSE,
      forward=FALSE)
mcRun("fitstepwise.bestpath", mcr.outfile="fitstepwise.h0.wswi.24.BB.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=24, wsbi=FALSE,
      forward=FALSE)

# ----------- Power (H0 false) ------------

pmx[,"eff"] <- .8

# WSBI design
mcRun("fitstepwise", mcr.outfile="fitstepwise.h1.wsbi.12.MS.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=12, wsbi=TRUE,
      mf=modSpace(wsbi=TRUE))
mcRun("fitstepwise", mcr.outfile="fitstepwise.h1.wsbi.24.MS.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=24, wsbi=TRUE,
      mf=modSpace(wsbi=TRUE))

# WSWI design
# Note: fitstepwise does both forward and backward selection
# forward subject slope first & backward item slope first
mf <- modSpace(FALSE)   # define the models to be tested
mcRun("fitstepwise", mcr.outfile="fitstepwise.h1.wswi.12.FSBI.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=12, wsbi=FALSE,
      mf=c(mf[3], mf[[2]][1], mf[1]))
mcRun("fitstepwise", mcr.outfile="fitstepwise.h1.wswi.24.FSBI.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=24, wsbi=FALSE,
      mf=c(mf[3], mf[[2]][1], mf[1]))
# forward item slope first & backward subject slope first
mcRun("fitstepwise", mcr.outfile="fitstepwise.h1.wswi.12.FIBS.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=12, wsbi=FALSE,
      mf=c(mf[3], mf[[2]][2], mf[1]))
mcRun("fitstepwise", mcr.outfile="fitstepwise.h1.wswi.24.FIBS.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=24, wsbi=FALSE,
      mf=c(mf[3], mf[[2]][2], mf[1]))
# 'best path' forward
mcRun("fitstepwise.bestpath", mcr.outfile="fitstepwise.h1.wswi.12.FB.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=12, wsbi=FALSE,
      forward=TRUE)
mcRun("fitstepwise.bestpath", mcr.outfile="fitstepwise.h1.wswi.24.FB.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=24, wsbi=FALSE,
      forward=TRUE)
# 'best path' backward
mcRun("fitstepwise.bestpath", mcr.outfile="fitstepwise.h1.wswi.12.BB.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=12, wsbi=FALSE,
      forward=FALSE)
mcRun("fitstepwise.bestpath", mcr.outfile="fitstepwise.h1.wswi.24.BB.csv",
      mcr.xdatFnc="mkDf", mcr.varying=pmx, nitem=24, wsbi=FALSE,
      forward=FALSE)

Author: Dale J. Barr, Roger Levy, Christoph Scheepers, and Harry J. Tily

Last Update: March 27, 2012