Package {mcstatsim}

Calculate Bias and Bias Monte Carlo Standard Error

Description

This function computes the bias and the Monte Carlo Standard Error (MCSE) of the bias for a set of estimates relative to a true parameter value. The bias is the difference between the mean of the estimates and the true parameter. The MCSE of the bias is calculated as the square root of the variance of the estimates divided by the number of estimates, providing a measure of the precision of the bias estimate.

Usage

calc_bias(estimates, true_param, digits = 4)

Arguments

Value

A list with two components: 'bias', the calculated bias of the estimates, and 'bias_mcse', the Monte Carlo Standard Error of the bias, indicating the uncertainty associated with the bias estimate.

Examples

estimates <- rnorm(100, mean = 50, sd = 10)
true_param <- 50
bias_info <- calc_bias(estimates, true_param)
print(bias_info)

Calculate Coverage Probability and its Monte Carlo Standard Error

Description

Computes the coverage probability of a confidence interval, defined as the proportion of times the true parameter value falls within the calculated lower and upper bounds across a set of simulations. Additionally, calculates the Monte Carlo Standard Error (MCSE) of the coverage probability to assess the uncertainty associated with this coverage estimate. This function is useful for evaluating the accuracy and reliability of confidence intervals generated by statistical models or estimation procedures.

Usage

calc_coverage(lower_bound, upper_bound, true_param, digits = 4)

Arguments

Value

A list with two components: 'coverage', the calculated coverage probability of the confidence intervals, and 'coverage_mcse', the Monte Carlo Standard Error of the coverage. This MCSE provides a measure of the precision of the coverage probability estimate.

Examples

set.seed(123) # For reproducibility
estimates <- rnorm(100, mean = 50, sd = 10)
ci_lower <- estimates - 1.96 * 10
ci_upper <- estimates + 1.96 * 10
true_param <- 50
coverage_info <- calc_coverage(ci_lower, ci_upper, true_param)
print(coverage_info)

Calculate Mean Squared Error and its Monte Carlo Standard Error

Description

Computes the Mean Squared Error (MSE) of a set of estimates relative to a true parameter value, along with the Monte Carlo Standard Error (MCSE) for the MSE. The MCSE takes into account the variance, skewness, and kurtosis of the estimates to provide a more accurate measure of uncertainty. This function is useful for assessing the accuracy of simulation or estimation methods by comparing the squared deviations of estimated values from a known parameter.

Usage

calc_mse(estimates, true_param, digits = 4)

Arguments

Value

A list with two components: 'mse', the calculated Mean Squared Error of the estimates, and 'mse_mcse', the Monte Carlo Standard Error of the MSE, offering insight into the reliability of the MSE calculation.

Examples

estimates <- rnorm(100, mean = 50, sd = 10)
true_param <- 50
mse_info <- calc_mse(estimates, true_param)
print(mse_info)

Calculate Rejection Rate and its Monte Carlo Standard Error

Description

Computes the rejection rate of hypotheses tests based on a vector of p-values and a specified significance level (alpha). The rejection rate is the proportion of p-values that are lower than alpha, indicating significant results. Additionally, the function calculates the Monte Carlo Standard Error (MCSE) for the rejection rate, which quantifies the uncertainty associated with the estimated rejection rate. This function is useful for assessing the overall type I error rate or the power of a statistical test across multiple simulations or experimental replications.

Usage

calc_rejection_rate(p_values, alpha = 0.05, digits = 4)

Arguments

Value

A list with two components: 'rejection_rate', the proportion of tests that resulted in rejection of the null hypothesis, and 'rejection_rate_mcse', the Monte Carlo Standard Error of the rejection rate, providing an estimate of its variability.

Examples

set.seed(123) # For reproducibility
p_values <- runif(100, min = 0, max = 1) # Simulated p-values
rejection_info <- calc_rejection_rate(p_values)
print(rejection_info)

Calculate Relative Bias and its Monte Carlo Standard Error

Description

Computes the relative bias of a set of estimates with respect to a true parameter value, along with the Monte Carlo Standard Error (MCSE) of the relative bias. Relative bias is the ratio of the mean of the estimates to the true parameter, providing a scale-independent measure of bias. This function is particularly useful for evaluating the accuracy of estimates in situations where the magnitude of the true parameter is crucial to the interpretation of bias. The function gracefully handles cases where the true parameter is zero by returning 'NA' for both relative bias and its MCSE, avoiding division by zero errors.

Usage

calc_relative_bias(estimates, true_param, digits = 4)

Arguments

Value

A list with two components: 'rel_bias', the calculated relative bias of the estimates, and 'rel_bias_mcse', the Monte Carlo Standard Error of the relative bias. If 'true_param' is zero, both 'rel_bias' and 'rel_bias_mcse' will be 'NA'.

Examples

estimates <- rnorm(100, mean = 50, sd = 10)
true_param <- 50 # Non-zero true parameter
relative_bias_info <- calc_relative_bias(estimates, true_param)
print(relative_bias_info)

Calculate Relative Mean Squared Error and its Monte Carlo Standard Error

Description

Computes the Relative Mean Squared Error (Relative MSE) of a set of estimates with respect to a true parameter value, along with the Monte Carlo Standard Error (MCSE) of the Relative MSE. The Relative MSE is a normalized measure of error that scales the Mean Squared Error (MSE) by the square of the true parameter value, making it particularly useful for comparing the accuracy of estimates across different scales. The function also calculates the MCSE for the Relative MSE, taking into account the variance, skewness, and kurtosis of the estimates to provide a measure of uncertainty. The function returns 'NA' for both Relative MSE and its MCSE if the true parameter is zero, to avoid division by zero.

Usage

calc_relative_mse(estimates, true_param, digits = 4)

Arguments

Value

A list with two components: 'rel_mse', the calculated Relative Mean Squared Error of the estimates, and 'rel_mse_mcse', the Monte Carlo Standard Error of the Relative MSE. If 'true_param' is zero, both 'rel_mse' and 'rel_mse_mcse' will be 'NA'.

Examples

estimates <- rnorm(100, mean = 50, sd = 10)
true_param <- 50 # Non-zero true parameter
relative_mse_info <- calc_relative_mse(estimates, true_param)
print(relative_mse_info)

Calculate Relative Root Mean Squared Error and its Monte Carlo Standard Error

Description

Computes the Relative Root Mean Squared Error (Relative RMSE) of a set of estimates with respect to a true parameter value. The Relative RMSE is the square root of the Relative Mean Squared Error (see [calc_relative_mse()]), a scale-independent measure of error that facilitates comparisons across different scales of the true parameter. This function is especially useful for evaluating the accuracy of estimates when the magnitude of the true parameter varies significantly across different scenarios. The function gracefully handles cases where the true parameter is zero by returning 'NA' for both Relative RMSE and its MCSE, to avoid division by zero.

Usage

calc_relative_rmse(estimates, true_param, digits = 4)

Arguments

Value

A list with two components: 'rel_rmse', the calculated Relative Root Mean Squared Error of the estimates, and 'rel_rmse_mcse', the Monte Carlo Standard Error of the Relative RMSE, obtained from the MCSE of the Relative MSE by the delta method ('rel_mse_mcse / (2 * rel_rmse)'). Both are 'NA' when 'true_param' is zero, and 'rel_rmse_mcse' is 'NA' when 'rel_rmse' is zero.

Examples

estimates <- rnorm(100, mean = 50, sd = 10)
true_param <- 50 # Non-zero true parameter
relative_rmse_info <- calc_relative_rmse(estimates, true_param)
print(relative_rmse_info)

Calculate Root Mean Squared Error and its Monte Carlo Standard Error

Description

Computes the Root Mean Squared Error (RMSE) of a set of estimates relative to a true parameter value, along with the Monte Carlo Standard Error (MCSE) for the RMSE. The RMSE is a measure of the accuracy of the estimates, representing the square root of the average squared differences between the estimated values and the true parameter. The MCSE for the RMSE is calculated using jackknife estimates, providing an assessment of the uncertainty associated with the RMSE value.

Usage

calc_rmse(estimates, true_param, digits = 4)

Arguments

Value

A list with two components: 'rmse', the calculated Root Mean Squared Error of the estimates, and 'rmse_mcse', the Monte Carlo Standard Error of the RMSE. This MCSE is derived from jackknife estimates, offering insight into the reliability of the RMSE calculation.

Examples

estimates <- rnorm(100, mean = 50, sd = 10)
true_param <- 50
rmse_info <- calc_rmse(estimates, true_param)
print(rmse_info)

Calculate Variance and its Monte Carlo Standard Error

Description

Computes the sample variance of a set of estimates and the Monte Carlo Standard Error (MCSE) for the variance. The MCSE is adjusted by the sample kurtosis to account for the shape of the distribution of the estimates. This function is particularly useful in simulation studies where understanding the variability of an estimator and the precision of this variability estimate is crucial.

Usage

calc_variance(estimates, digits = 4)

Arguments

Value

A list containing two elements: 'variance', the calculated sample variance of the estimates, and 'variance_mcse', the Monte Carlo Standard Error of the variance. The MCSE provides a measure of the uncertainty associated with the variance estimate, adjusted for kurtosis.

Examples

estimates <- rnorm(100, mean = 50, sd = 10)
variance_info <- calc_variance(estimates)
print(variance_info)

Calculate Average Width of Confidence Intervals and its Monte Carlo Standard Error

Description

This function computes the average width of a set of confidence intervals, represented by their lower and upper bounds, along with the Monte Carlo Standard Error (MCSE) of this average width. The average width provides insight into the precision of the estimation process, with narrower intervals typically indicating more precise estimates. The MCSE of the width offers a measure of the uncertainty associated with the average width calculation, useful for assessing the variability of interval estimates across simulations or bootstrap samples.

Usage

calc_width(lower_bound, upper_bound, digits = 4)

Arguments

Value

A list with two components: 'width', the average width of the confidence intervals, and 'width_mcse', the Monte Carlo Standard Error of the average width. This MCSE provides a quantification of the uncertainty in the average width estimate.

Examples

set.seed(123) # For reproducibility
estimates <- rnorm(100, mean = 50, sd = 10)
ci_lower <- estimates - 1.96 * 10
ci_upper <- estimates + 1.96 * 10
width_info <- calc_width(ci_lower, ci_upper)
print(width_info)

Flag Non-Finite or Degenerate Estimates in Simulation Output

Description

Scans the numeric columns of a replicate-level [runsim()] result for values that will quietly distort performance measures: 'NA', 'NaN', and 'Inf'. Reports the count and proportion per condition, so a handful of blown-up replications do not pass unnoticed into a summary table.

Usage

check_estimates(data, cols = NULL, by = NULL, digits = 4)

Arguments

Value

A dataframe with one row per (condition, column) that has at least one non-finite value: the 'by' columns, 'column', 'n_bad', 'n_total', and 'prop_bad', ordered by decreasing 'prop_bad'. A zero-row dataframe (with those columns) means every checked value was finite.

Examples

sim_func <- function(mu) {
  x <- rnorm(10, mu)
  data.frame(est = mean(x), ratio = mu / (mu - 1))  # Inf at mu == 1
}
grid <- expand.grid(mu = c(0, 1, 2))
res <- runsim(20, grid, sim_func, show_progress = FALSE, num_cores = 1)
check_estimates(res, by = "mu")

Combine Dataframes in a Nested List

Description

This function combines dataframes that are nested within a list of lists into a single dataframe.

Usage

combine_df(nested_list)

Arguments

Details

The function first checks if the input is a non-empty list of lists containing dataframes. It then iterates through each sublist, combining the dataframes and adding an ID column to indicate their source. Finally, all combined dataframes are row-bound into a single dataframe.

Value

A combined dataframe where each original dataframe is augmented with an ID column indicating its source list.

Examples

df1 <- data.frame(a = 1:3, b = 4:6)
df2 <- data.frame(a = 7:9, b = 10:12)
nested_list <- list(list(df1, df2), list(df1, df2))
combine_df(nested_list)

Summarise Failed Jobs from a Simulation Run

Description

Aggregates the per-job failure log attached to a [runsim()] result (when it was called with 'on_error = "warn"' or 'on_error = "omit"') into one row per simulation condition, with the number of failed replications, the failure rate, and an example error message.

Usage

failure_summary(x, n = NULL, digits = 4)

Arguments

Value

A dataframe with the grid columns that describe each condition, plus 'n_failed', 'failure_rate', and 'example_error'. If 'x' has no '"errors"' attribute (no job failed, or 'on_error = "stop"' was used) a zero-row dataframe with those columns is returned.

Examples

sim_func <- function(mu) {
  if (mu < 0) stop("mu must be non-negative")
  data.frame(est = mean(rnorm(10, mu)))
}
params <- expand.grid(mu = c(-1, 0, 1))
res <- runsim(5, params, sim_func, show_progress = FALSE, num_cores = 1,
              on_error = "omit")
failure_summary(res)

Parallel Map Function using functionals::fmapn

Description

This function applies a given function over a list of parameters in parallel using multiple cores.

Usage

mcpmap(
  lists,
  func,
  num_cores = parallel::detectCores() - 1,
  show_progress = TRUE
)

Arguments

Details

The function ensures that all elements in the list have the same length and uses 'functionals::fmapn' for parallel processing. It sets the number of cores based on the operating system and then applies the function in parallel.

Value

A list of results from applying the function over the parameters.

Examples

params <- list(a = 1:3, b = 4:6)
mcpmap(params, function(a, b) a + b, num_cores = 2)

Replications Needed to Reach a Target Monte Carlo Standard Error

Description

Monte Carlo standard errors shrink like '1 / sqrt(n_sim)'. Given the MCSE actually achieved with 'n_current' replications, this returns the number of replications a future run would need to bring the MCSE down to 'target'.

Usage

mcse_target(mcse_current, n_current, target)

Arguments

Value

An integer vector ('ceiling') of the replications required to reach 'target'. Entries are 'NA' where 'mcse_current' is 'NA' or not positive.

Examples

# A coverage estimate with MCSE 0.011 from 1000 reps; how many reps for MCSE 0.005?
mcse_target(0.011, 1000, 0.005)

# Vectorised over a summarise_sim() column
# mcse_target(tab$coverage_mcse, tab$n_sim, target = 0.005)

Convergence Trace of a Monte Carlo Estimate

Description

Plots the running mean of a quantity against the number of replications, with a Monte Carlo standard error band, one panel per simulation condition. It shows whether a run used enough replications for the quantity of interest to settle.

Usage

plot_convergence(data, var, by = NULL, se_mult = 1)

Arguments

Value

A 'ggplot' object.

Examples

## Not run: 
sim_func <- function(mu) data.frame(est = mean(rnorm(30, mu)))
res <- runsim(500, expand.grid(mu = c(0, 2)), sim_func,
              show_progress = FALSE, num_cores = 1, seed = 1)
plot_convergence(res, "est", by = "mu")

## End(Not run)

Point-Range Plot of a Performance Measure Across Conditions

Description

Takes a summary table from [summarise_sim()] and plots one performance measure per condition as a point with a Monte Carlo error bar ('estimate +/- se_mult * MCSE'), the standard "does the method work across the design" picture.

Usage

plot_performance(
  summary,
  measure,
  x,
  colour = NULL,
  facet = NULL,
  se_mult = 2,
  ref = NULL
)

Arguments

Value

A 'ggplot' object.

Examples

## Not run: 
tab <- summarise_sim(res, by = c("mu", "n_obs"),
                     measures = list(coverage = ~ calc_coverage(lo, hi, true_val[1])))
plot_performance(tab, "coverage", x = "n_obs", colour = "mu", ref = 0.95)

## End(Not run)

Zip Plot of Confidence-Interval Coverage

Description

Draws the zip plot of Morris, White and Crowther (2019): within each condition the confidence intervals are sorted by how far their implied z-statistic sits from the target value, stacked vertically by that fractional centile, and coloured by whether they cover. A well-calibrated procedure fills the panel to the nominal line with the covering colour and shows the non-covering intervals only in the top tail.

Usage

plot_zip(data, lower, upper, true, by = NULL, conf_level = 0.95)

Arguments

Value

A 'ggplot' object.

References

Morris TP, White IR, Crowther MJ (2019). Using simulation studies to evaluate statistical methods. *Statistics in Medicine*, 38(11), 2074-2102.

Examples

## Not run: 
sim_func <- function(mu) {
  x <- rnorm(25, mu); est <- mean(x); se <- sd(x) / 5
  data.frame(lo = est - 1.96 * se, hi = est + 1.96 * se, mu = mu)
}
res <- runsim(300, expand.grid(mu = c(0, 1)), sim_func,
              show_progress = FALSE, num_cores = 1, seed = 1)
plot_zip(res, "lo", "hi", "mu", by = "mu")

## End(Not run)

Run Monte Carlo Simulations in Parallel

Description

This function executes a series of Monte Carlo simulations in parallel, providing detailed progress updates.

Usage

runsim(
  n,
  grid_params,
  sim_func,
  show_progress = TRUE,
  num_cores = parallel::detectCores() - 1,
  seed = NULL,
  on_error = c("stop", "warn", "omit"),
  save_path = NULL,
  checkpoint_every = 1L
)

Arguments

Details

The function first validates the input parameters. It then builds a single flat list of 'n * nrow(grid_params)' jobs (replication-major order: all conditions for replication 1, then all conditions for replication 2, and so on) and dispatches it in one load-balanced parallel pass via [mcpmap()]. Results are row-bound into a single dataframe, the condition's 'grid_params' columns are prepended (unless 'sim_func' already returned a column of that name), and the replication index is tagged in the 'ID' column.

If 'seed' is supplied, 'n * nrow(grid_params)' independent RNG streams are generated up front with [parallel::nextRNGStream()] and handed to the jobs by position. Because job positions are fixed (replication-major), a given '(replication, condition)' cell always draws from the same stream regardless of core count, so runs are bit-for-bit reproducible. Do not call 'set.seed()' inside 'sim_func' when using 'seed', as that would defeat the per-stream design.

Each job is wrapped so that an error in 'sim_func' is caught rather than lost in a worker process, and any warnings it emits are recorded. See 'on_error' for how failures are handled.

When 'save_path' is set the replications are run in chunks of 'checkpoint_every' and the checkpoint file is rewritten after each chunk. Per-job RNG streams are keyed by position, so a resumed run produces exactly the same results as an uninterrupted one.

Value

A combined dataframe of all simulation results. Every 'grid_params' column that 'sim_func' did not itself return is prepended, so each row identifies its own condition, and an 'ID' column gives the replication index (1 to 'n'). When 'seed' is supplied it is stored in 'attr(x, "seed")'. Any warnings emitted by 'sim_func' are collected in 'attr(x, "warnings")', and, when 'on_error' is not '"stop"', failed jobs are described in 'attr(x, "errors")' (see [failure_summary()]).

Examples

## Not run: 
library(mcstatsim)

# Define a simple simulation function
sim_function <- function(a, b) {
  Sys.sleep(0.1)  # Simulate a time-consuming process
  return(data.frame(result = a + b))
}

# Generate a grid of parameters
params <- expand.grid(a = 1:3, b = 4:6)

# Run simulations
results <- runsim(n = 1, grid_params = params, sim_func = sim_function)
print(results)

## End(Not run)

Summarise Replicate-Level Simulation Output into a Performance Table

Description

Collapses the raw, replicate-level dataframe returned by [runsim()] into one row per simulation condition, evaluating one or more performance measures (bias, coverage, rejection rate, ...) per condition. This is the bridge between 'runsim()' output and the 'calc_*' functions in this package.

Usage

summarise_sim(data, by = NULL, measures, digits = 4)

Arguments

Value

A dataframe with one row per group: the 'by' columns, then 'n_sim' (the number of distinct 'ID' values in the group, or the row count if there is no 'ID' column), then one column per performance-measure output.

Examples

sim_func <- function(mu, n_obs) {
  x <- rnorm(n_obs, mean = mu)
  est <- mean(x)
  se <- sd(x) / sqrt(n_obs)
  data.frame(est = est, lo = est - 1.96 * se, hi = est + 1.96 * se, true_val = mu)
}
grid <- expand.grid(mu = c(0, 1), n_obs = c(20, 50))
res <- runsim(100, grid, sim_func, show_progress = FALSE, num_cores = 1, seed = 1)

summarise_sim(
  res,
  by = c("mu", "n_obs"),
  measures = list(
    bias     = ~ calc_bias(est, true_val[1]),
    coverage = ~ calc_coverage(lo, hi, true_val[1])
  )
)