Code used to determine power for GxSES in Bates et al (2016) doi http://dx.doi.org/10.1016/j.intell.2016.02.003

simulate GxE power.R

# Updated for umx 1.7+ 2017-06-12 04:28PM
# Notes: If you're on Mac or Unix, install the parallelOpenMx to get parallel (4x speedup or more)
# source('http://openmx.psyc.virginia.edu/getOpenMx.R')
library(umx)
umx_set_optimizer("NPSOL") # good optimizer for these data
umx_set_cores(detectCores()) # Max cores for speed

nSimulations = 1000                  # Number of simulations
nMZpairs     = nDZpairs = 500        # Number of twin pairs
pvalues      = rep(NA, nSimulations) # placeholder for the p-values from mxCompare-ing the 2 models you are testing

# ===============================================================
# = Set mean and lower and upper h2 expected across 4SDs of SES =
# = Alter these to explore the parameter space.                 =
# ===============================================================
avgA = .5; minA = .43; maxA = .57

for (i in 1:nSimulations) {
	umx_msg(i); # i = 1
	simData  = umx_make_TwinData(nMZpairs, nDZpairs, AA = c(avg = avgA, min = minA, max = maxA), CC = 0, EE = .1);
	# build base model
	m1 = umxGxE(selDVs = "var", selDefs = "M", suffix = "_T", dzData = simData[[2]], mzData = simData[[1]], autoRun = FALSE);
	# build AE model with A' (m2), and AE without A' (m3)
	m2 = umxModify(m1, update = c("c_r1c1", "cm_r1c1", "em_r1c1"), name = "AE plus am"); # drop C and  c- and e-moderation
	m3 = umxModify(m2, update = c("am_r1c1"), name = "AE model, no am"); # now drop a moderation
	# umxCompare(m2, m3)
	# store the p-value for dropping A' from an AE model with a moderation
	pvalues[i] = mxCompare(m2, m3)["p"][2,1]
	# can add other models, of course, e.g. pAm_vs_Cm etc.
}

notNA = na.omit(pvalues)
message("Power was ", (length(notNA[notNA < .05])/length(notNA))*100, "% given mean 'a' of ", .5, " (min and max " , minA, " and ", maxA, " respectively). Based on ", nSimulations, " simulations, of which ", nSimulations - length(notNA), " were NA.")

# Power was >99.9% given mean 'a' of 0.5 (min and max 0.3 and 0.7 respectively). Based on 1000 simulations, of which 0 failed

# Power was 98% given mean 'a' of 0.5 (min and max 0.4 and 0.6 respectively). Based on 50 simulations, of which 0 were NA.
# Power was 51.2% given mean 'a' of 0.5 (min and max 0.45 and 0.55 respectively). Based on 1000 simulations.

# Based on 1000 simulations, power was 84.9% given a mean path-coefficient 'a' sqrt(h^2) of 0.5, swinging from a minimum of 0.43 at -2 SDs of SES, and rising to 0.57 at +SD above mean SES.