sim_power {mpower}

Power analysis using Monte Carlo simulation

Description

Power analysis using Monte Carlo simulation

Usage

sim_power(
  xmod,
  ymod,
  imod,
  s = 100,
  n = 100,
  cores = 1,
  file = NULL,
  errorhandling = "stop",
  snr_iter = 10000,
  cluster_export = c()
)

Arguments

xmod

A MixtureModel object.

ymod

An OutcomeModel object.

imod

An InferenceModel object.

s

An integer for number of Monte Carlo simulations.

n

An integer for sample size in each simulation.

cores

An integer for number of processing cores. When cores > 1, parallelism is automatically applied.

file

A string, a file name with no extension to write samples to periodically. By default write to an RDS file.

errorhandling

A string "remove", "stop", or "pass". If an error occurs in any iteration, remove that iteration (remove), return the error message verbatim in the output (pass), or terminate the loop (stop). Default is "remove". See R package 'foreach' for more details.

snr_iter

An integer for number of Monte Carlo samples to estimate SNR

cluster_export

A vector of functions to pass to the parallel-processing clusters

Value

A PowerSim object. Attributes:

s

a number of simulations.

snr

a real number for SNR of the OutcomeModel.

n

a number of sample sizes.

xmod

the MixtureModel used.

ymod

the OutcomeModel used.

imod

the InferenceModel used.

sims

an output matrices.

Examples

data("nhanes1518")
chems <- c("URXCNP", "URXCOP", "URXECP", "URXHIBP", "URXMBP", "URXMC1",
"URXMCOH", "URXMEP","URXMHBP", "URXMHH", "URXMHNC", "URXMHP", "URXMIB",
"URXMNP", "URXMOH", "URXMZP")
chems_mod <- mpower::MixtureModel(nhanes1518[, chems], method = "resampling")
bmi_mod <- mpower::OutcomeModel(f = "0.2*URXCNP + 0.15*URXECP +
0.1*URXCOP*URXECP", family = "binomial")
logit_mod <- mpower::InferenceModel(model = "glm", family = "binomial")
logit_out <- mpower::sim_power(xmod=chems_mod, ymod=bmi_mod, imod=logit_mod,
s=100, n=2000, cores=2, snr_iter=2000)
logit_df <- summary(logit_out, crit="pval", thres=0.05, how="lesser")
plot_summary(logit_out, crit="pval", thres=0.05, how="lesser")