simulate.plcp: Perform a simulation study using a 'study_parameters'-object
In powerlmm: Power Analysis for Longitudinal Multilevel Models

Description Usage Arguments Details See Also Examples

Perform a simulation study using a study_parameters-object

## S3 method for class 'plcp'
simulate(object, nsim, seed = NULL, formula = NULL,
  satterthwaite = FALSE, CI = FALSE, cores = 1, progress = FALSE,
  batch_progress = TRUE, ...)

## S3 method for class 'plcp_multi'
simulate(object, nsim, seed = NULL,
  formula = NULL, satterthwaite = FALSE, CI = FALSE, cores = 1,
  progress = FALSE, batch_progress = TRUE, ...)

`object`	An object created by `study_parameters`.
`nsim`	The number of simulations to run.
`seed`	Currently ignored.
`formula`	Model formula(s) used to analyze the data, see Details. Should be created using `sim_formula`. It is also possible to compare multiple models, e.g. a correct and a misspecified model, by combining the formulas using `sim_formula_compare`. See Examples. If `NULL` the formula is made automatically, using `create_lmer_formula`, which does not support objects with multiple simulation setups.
`satterthwaite`	Logical; if `TRUE` Satterthwaite's degrees of freedom approximation will be used when computing p-values. This is implemented using the `lmerTest`-package. See Details.
`CI`	Logical; if `TRUE` coverage rates for 95 % confidence intervals will be calculated. See Details.
`cores`	Number of CPU cores to use. If called from a GUI environment (e.g. RStudio) or a computer running Microsoft Windows, PSOCK clusters will be used. If called from a non-interactive Unix environment forking is utilized.
`progress`	`logical`; will display progress if `TRUE`. Currently ignored on Windows. Package `pbmclapply` is used to display progress, which relies on forking. N.B using a progress bar will noticeably increase the simulation time, due to the added overhead.
`batch_progress`	`logical`; if `TRUE` progress will be shown for simulations with multiple setups.
`...`	Optional arguments, see Saving in Details section.

See also vignette("simulations", package = "powerlmm") for a tutorial.

Model formula

If no data transformation is used, the available model terms are:

y the outcome vector, with potential missing data.
y_c the complete version of y, before dropout was simulated.
time the time vector.
treatment treatment indicator (0 = "control", 1 = "treatment").
subject subject-level id variable, from 1 to total number of subjects.
cluster for three-level models; the cluster-level id variable, from 1 to the total number of clusters.

See Examples and the simulation-vignette for formula examples. For objects that contain a single study setup, then the lmer formula can be created automatically using create_lmer_formula.

Satterthwaite's approximation, and CI coverage

To decrease the simulation time the default is to only calculate Satterthwaite's dfs and the CIs' coverage rates for the test of 'time:treatment'-interaction. This can be changed using the argument test in sim_formula.

Confidence intervals are both calculated using profile likelihood and by the Wald approximation, using a 95 % confidence level.

Saving intermediate results for multi-sims

Objects with multi-sims can be save after each batch is finished. This is highly recommended when many designs are simulated. The following additional arguments control saving behavior:

'save', logical, if TRUE each batch is saved as a RDS-file. Results are saved in your working directory, in the directory specified by save_folder.
'save_folder' a character indicating the folder name. Default is 'save'.
'save_folder_create', logical, if TRUE then save_folder will be created if it does not exist in your working directory.

sim_formula, sim_formula_compare, summary.plcp_sim, simulate_data

## Not run: 
# Two-level ---------------------------------------------------------------
p <- study_parameters(n1 = 11,
                      n2 = 25,
                      sigma_subject_intercept = 1.44,
                      sigma_subject_slope = 0.2,
                      sigma_error = 1.44,
                      effect_size = cohend(0.5))

f <- sim_formula("y ~ treatment * time + (1 + time | subject)")


res <- simulate(object = p,
                nsim = 1000,
                formula = f,
                satterthwaite = TRUE,
                progress = FALSE,
                cores = 1,
                save = FALSE)

summary(res)


# Three-level (nested) ------------------------------------------------------
p <- study_parameters(n1 = 10,
                      n2 = 20,
                      n3 = 4,
                      sigma_subject_intercept = 1.44,
                      icc_pre_cluster = 0,
                      sigma_subject_slope = 0.2,
                      icc_slope = 0.05,
                      sigma_error = 1.44,
                      effect_size = 0)

## compare correct and miss-specified model
f0 <- "y ~ treatment * time + (1 + time | subject)"
f1 <- "y ~ treatment * time + (1 + time | subject) + (0 + time | cluster)"
f <- sim_formula_compare("correct" = f1,
                         "wrong" = f0)

res <- simulate(object = p,
                nsim = 1000,
                formula = f,
                satterthwaite = TRUE,
                progress = FALSE,
                cores = 1,
                save = FALSE)

summary(res)

## Compare random effects using LRT,
## summarise based on best model from each sim
summary(res,
        model_selection = "FW",
        LRT_alpha = 0.1,
        para = "treatment:time")

# Partially nested design ---------------------------------------------------
p <- study_parameters(n1 = 11,
                      n2 = 10,
                      n3 = 4,
                      sigma_subject_intercept = 1.44,
                      icc_pre_cluster = 0,
                      sigma_subject_slope = 0.2,
                      cor_subject = -0.5,
                      icc_slope = 0.05,
                      sigma_error = 1.44,
                      partially_nested = TRUE,
                      effect_size = cohend(-0.5))

f <- sim_formula("y ~ treatment * time + (1 + time | subject) +
                  (0 + treatment:time | cluster)")

res <- simulate(object = p,
                nsim = 1000,
                formula = f,
                satterthwaite = TRUE,
                progress = FALSE,
                cores = 4,
                save = FALSE)

summary(res)

# Run multiple designs  -----------------------------------------------------
p <- study_parameters(n1 = 10,
                      n2 = 20,
                      n3 = c(2, 4, 6),
                      sigma_subject_intercept = 1.44,
                      icc_pre_cluster = 0,
                      sigma_subject_slope = 0.2,
                      icc_slope = 0.05,
                      sigma_error = 1.44,
                      effect_size = cohend(0.5))

f0 <- "y ~ treatment * time + (1 + time | subject)"
f1 <- "y ~ treatment * time + (1 + time | subject) + (0 + time | cluster)"
f <- sim_formula_compare("correct" = f1,
                         "wrong" = f0)

res <- simulate(object = p,
                nsim = 1000,
                formula = f,
                satterthwaite = TRUE,
                progress = FALSE,
                cores = 1,
                save = FALSE)

# Summarize 'time:treatment' results for n3 = c(2, 4, 6) for 'correct' model
summary(res, para = "time:treatment", model = "correct")

# Summarize cluster-level random slope  for n3 = c(2, 4, 6) for 'correct' model
summary(res, para = "cluster_slope", model = "correct")

## End(Not run)