power_diff: Power simulation function for a two-group comparison

Description

There are two modes of operation:

power=NULL: simulate power based on given sample size n
n=NULL: search iteratively for a suitable sample size n for a given power

Usage

power_diff(
  distribution = c("exponential", "weibull"),
  twoPhase = FALSE,
  eff = stop("Provide parameters for both groups that reflect the effect!"),
  param = "delay1",
  test = c("bootstrap", "logrank", "logrank_pp", "LRT"),
  method = c("MPSE", "MLEw", "MLEc", "MLEn"),
  n = NULL,
  r = 1,
  sig.level = 0.05,
  power = NULL,
  nPowerSim = 1600,
  R = 200,
  nRange = c(5, 250),
  verbose = 0
)

Value

List of results of power simulation. Or NULL in case of errors.

Arguments

distribution: character. Which assumed distribution is used for the power calculation. Default is "exponential".
twoPhase: logical(1). Do we model two phases per group? Default is FALSE, i.e. a single delay phase per group.
eff: list of length 2. The two list elements must be numeric vectors that contain the model parameters (as understood by the delay-distribution functions provided by this package) for the two groups.
param: character. Parameter name(s) which are to be tested for difference and for which to simulate the power. Default value is 'delay1'. You can specify multiple parameters, by giving a vector or by concatenating them with a + in a single string.
test: character. Which test to use for this power estimation? Defaults to "bootstrap". Non-parametric logrank test is also possible (either "logrank" or "logrank_pp"). See also test_diff().
method: character. Which fitting method to use in case of a parametric test.
n: integer. Number of observations per group for the power simulation or NULL when n is to be estimated for a given power.
r: numeric. Ratio of both groups sizes, ny / nx. Default value is 1, i.e. balanced group sizes. Must be positive.
sig.level: numeric. Significance level. Default is 0.05.
power: numeric. NULL when power is to be estimated for a given sample size or a desired power is specified (and n is estimated).
nPowerSim: integer. Number of simulation rounds. Default value 1600 yields a standard error of 0.01 for power if the true power is 80%.
R: integer. Number of bootstrap samples for test of difference within each power simulation. It affects the resolution of the P-value for each simulation round. A value of around R=200 gives a resolution of 0.5% which might be enough for power analysis.
nRange: integer. Admissible range for sample size when power is specified and sample size is requested. The routine might not find the optimal sample size when this range is set too wide.
verbose: numeric. How many details are requested? Higher value means more details. 0=off, no details.

Details

In both cases, the distribution, the parameters that are tested for, the type of test and the effect size (eff=) need to be specified. Specify the effect as a list with two elements, each element holds the parameter vector describing the distribution of the outcome per group. The fitting method (method=) and the kind of significance test (test=) are handed down to test_diff(). The more power simulation rounds (parameter nPowerSim=) the more densely the space of data according to the specified model is sampled.

Note that estimating sample size n is computationally intensive. The iterative search uses heuristics to find a sample size n within the provided range (nRange=), and the estimated sample size might yield a slightly different power level. Hence, check the reported power in the output. The search algorithm comes to better results when the admissible range for sample size (nRange=) is chosen sensibly and not too wide. In case the estimated sample size and the achieved power is too high it might pay off to rerun the function with an adapted admissible range for the sample size by giving a narrower range in nRange=.

Examples

Run this code

# Simulate power for a given sample size:
# test for any difference in a delay-exponential model using logrank test
# the assumed effect is given in terms of model parameters for both groups
power_diff(
  eff = list(grA = c(delay1 = 5, rate1 = .09),
             grB = c(delay1 = 7, rate1 = .12)),
  test = "logrank",
  n = 16, power = NULL, nPowerSim = 300)

if (FALSE) {
# test for difference in delay in a delay-exponential model via a bootstrap test:
# the power is estimated based on nPowersim = 420 simulated datasets
# for real applications, use a higher nPowerSim (e.g. 1600) for more
# precise power estimation and a higher R (e.g. 400) for more precise
# P-value estimation in each simulation round
set.seed(123) # for reproducibility
power_diff(
  eff = list(grA = c(delay1 = 5.2, rate1 = .1),
             grB = c(delay1 = 7, rate1 = .12)),
  param = "delay1",
  test = "bootstrap", method = "MPSE",
  n = 16, power = NULL,
  nPowerSim = 420, R = 160)
}

if (FALSE) {
# test for difference in rate in a delay-exponential model via a bootstrap test:
# required sample size is estimated for a given power.
# provide a range for the sample size search, e.g. nRange = c(10, 30)
# this takes more time than the previous example, as the function
# iteratively searches for a sample size that yields the requested power
set.seed(1234) # for reproducibility
power_diff(
  eff = list(grA = c(delay1 = 5, rate1 = .07),
             grB = c(delay1 = 7, rate1 = .18)),
  param = "rate1",
  test = "bootstrap", method = "MPSE",
  n = NULL, power = 0.8,
  nPowerSim = 480, R = 150,
  nRange = c(18, 48))
}