Details
This is a long-running computation, taking roughly 3 days on a recent
macbook air. It is written to save data as it runs, one file per
simulation, so it can be continued if interrupted. To continue, just
rerun the loop: only simulation for which these files are missing are
run. It can also be run on a slurm cluster, and will do so if a
correctly configured template file batchtools.slurm.tmpl is present. We
include the one we use in the paper-results-details directory.
The loop gives each replication of a simulation its own RNG seed.
Within each replication, every call to a simulator, estimator, or
variance estimator uses this replication-specific seed: we explicitly
restore the seed before each call. As a result, changing the set of
simulations, estimators, or variance estimators considered does not
change the results we get for any specific
simulator/estimator/variance-estimator within a given replication. Nor
does interrupting and continuing the loop.
By using RNGSeq to choose replication-specific seeds for the
L’Ecuyer-CMRG RNG, we get streams of random numbers that are more or
less independent from replication to replication. See
help(RNGseq).
sim.replications = 1000
coverage.replications = 400
coverage.estimators = c('sdid','sc', 'did', 'difp')
coverage.sims = c('baseline', 'gun.law', 'abortion', 'random', 'hours', 'urate',
'T1.is.one', 'N1.is.one', 'blocksize.is.one',
'resample.200', 'resample.400',
'penn.democracy', 'penn.education', 'penn.random')
se.methods = c('bootstrap', 'jackknife', 'placebo')
cluster = file.exists('batchtools.slurm.tmpl')
# set up simulation grid
grid = expand.grid(ss = 1:length(simulators),
ee = 1:length(estimators),
rr = 1:sim.replications)
grid$estimate.se = names(simulators[grid$ss]) %in% coverage.sims &
names(estimators[grid$ee]) %in% coverage.estimators &
grid$rr <= coverage.replications
# associate a grid row to an output filename
output.file = function(row) {
sprintf('simulations/simulation-%s-%s-%d.rds',
names(simulators)[row$ss], names(estimators)[row$ee], row$rr)
}
rows.finished = function() {
output.files = sapply(1:nrow(grid), function(ii) { output.file(grid[ii,]) })
file.exists(output.files)
}
# set up RNG seeds for each replication.
# We need to pass a L'Ecuyer-type seed: a vector of 7 integers starting with 10407 like initial.seed below.
# If we pass a simple integer seed, the first call in a vanilla R session (R --vanilla) differs from subsequent ones.
# To get initial_seed I run: 'library(rngtools); ignore=RNGseq(n=1, seed=12345); initial.seed=RNGseq(n=1,seed=12345)'.
initial.seed = c(10407, -2132566924, 1638542565, 108172386, -1884566405, -1838154368, -250773631)
seeds = RNGseq(n=sim.replications, seed=initial.seed)
run.simulation = function(row) {
setRNG(seeds[[row$rr]])
setup = simulators[[row$ss]]$run()
estimate = estimators[[row$ee]](setup$Y, setup$N0, setup$T0)
se.estimates = sapply(se.methods, function(method) {
setRNG(seeds[[row$rr]])
if(row$estimate.se) { sqrt(vcov(estimate, method = method)) } else { NA }
})
cbind(data.frame(simulation = names(simulators)[row$ss],
estimator = names(estimators)[row$ee],
replication = row$rr,
estimate = c(estimate)),
t(as.data.frame(se.estimates)))
}
# kill warning that the futures library can't determine that we're
# using and deterministically seeding the L'Ecuyer RNG
options('future.options.rng.onMisuse', 'ignore')
## $future.options.rng.onMisuse
## NULL
##
## $ignore
## NULL
registerDoFuture()
# set up worker processes
if(cluster) {
# use multiple nodes on a Slurm Cluster
plan(batchtools_slurm, workers=1000, resources=list(ncpus=1, memory=1024))
} else {
# use multiple processes on this this computer
plan(multisession, workers = 8)
}
