cps.count: Simulation-based power estimation for cluster-randomized...

In nickreich/clusterPower: Power Calculations for Cluster-Randomized and Cluster-Randomized Crossover Trials

Description

\loadmathjax

This function uses Monte Carlo methods (simulations) to estimate power for cluster-randomized trials with integer-valued outcomes. Users can modify a variety of parameters to suit the simulations to their desired experimental situation.

Users must specify the desired number of simulations, number of subjects per cluster, number of clusters per treatment arm, between-cluster variance, and two of the following three parameters: mean event rate per unit time in one group, the mean event rate per unit time in the second group, and/or the mean difference in event rates between groups. Default values are provided for significance level, analytic method, whether progress updates are displayed, and whether the simulated data sets are retained.

Note that if all units have the same observation time, you can use the mean count instead of the "mean event per unit time" in the preceding paragraph.

Usage

cps.count(
  nsim = NULL,
  nsubjects = NULL,
  nclusters = NULL,
  c1 = NULL,
  c2 = NULL,
  cDiff = NULL,
  sigma_b_sq = NULL,
  sigma_b_sq2 = NULL,
  family = "poisson",
  negBinomSize = 1,
  analysis = "poisson",
  method = "glmm",
  alpha = 0.05,
  quiet = FALSE,
  allSimData = FALSE,
  irgtt = FALSE,
  seed = NA,
  nofit = FALSE,
  poorFitOverride = FALSE,
  lowPowerOverride = FALSE,
  timelimitOverride = TRUE,
  optmethod = "Nelder_Mead"
)

Arguments

nsim	Number of datasets to simulate; accepts integer. Required.
nsubjects	Number of subjects per cluster; accepts either a scalar (implying equal cluster sizes for the two groups), a vector of length two (equal cluster sizes within arm), or a vector of length sum(nclusters) (unequal cluster sizes within arm). If a vector of > 2 is provided in nsubjects, sum(nclusters) must match the nsubjects vector length. Required.
nclusters	Number of clusters per treatment group; accepts a single integer (if there are the same number of clusters in each arm) or a vector of 2 integers (if there are not). Required. At least 2 of the following 3 arguments must be specified:
c1	The mean event rate per unit time in the first arm.
c2	The mean event rate per unit time in the second arm.
cDiff	Expected difference in mean event rates between groups, defined as cDiff = c1 - c2.
sigma_b_sq	Between-cluster variance; if sigma_b_sq2 is not specified, between-cluster variances are assumed to be equal in the two arms. Accepts numeric. Required.
sigma_b_sq2	Between-cluster variance for clusters in the second arm. Only required if between-cluster variances differ between treatment arms.
family	Distribution from which responses are simulated. Accepts Poisson ('poisson') or negative binomial ('neg.binom'); default = 'poisson'. Required.
negBinomSize	Only used when generating simulated data from the negative binomial (family = 'neg.binom'), this is the target for number of successful trials, or the dispersion parameter (the shape parameter of the gamma mixing distribution). Must be strictly positive but need not be integer. Defaults to 1.
analysis	Family used for data analysis; currently only applicable when method = 'glmm'. Accepts c('poisson', 'neg.binom'); default = 'poisson'. Required.
method	Data analysis method, either generalized linear mixed effects model (GLMM) or generalized estimating equations (GEE). Accepts c('glmm', 'gee'); default = 'glmm'. Required.
alpha	The level of significance of the test, the probability of a Type I error. Default = 0.05.
quiet	When set to FALSE, displays simulation progress and estimated completion time. Default = FALSE.
allSimData	Option to include a list of all simulated datasets in the output object. Default = FALSE.
irgtt	Logical. Default = FALSE. Is the experimental design an individually randomized group treatment trial? For details, see ?cps.irgtt.count.
seed	Option to set the seed. Default is NA.
nofit	Option to skip model fitting and analysis and instead return a dataframe with the simulated datasets. Default = FALSE.
poorFitOverride	Option to override stop() if more than 25% of fits fail to converge.
lowPowerOverride	Option to override stop() if the power is less than 0.5 after the first 50 simulations and every ten simulations thereafter. On function execution stop, the actual power is printed in the stop message. Default = FALSE. When TRUE, this check is ignored and the calculated power is returned regardless of value.
timelimitOverride	Logical. When FALSE, stops execution if the estimated completion time is more than 2 minutes. Defaults to TRUE.
optmethod	Option to fit with a different optimizer. Defaults to Nelder_Mead.

Details

If family = 'poisson', the data generating model is: \mjsdeqny_ij \sim Poisson(e^c_1 + b_i) for the first group or arm, where \mjseqnb_i \sim N(0,\sigma_b^2), while for the second group,

\mjsdeqn

y_ij \sim Poisson(e^c_2 + b_i) where \mjseqnb_i \sim N(0,\sigma_b_2^2); if \mjseqn\sigma_b_2^2 is not specified, then the second group uses \mjseqnb_i \sim N(0,\sigma_b^2).

If family = 'neg.bin', the data generating model, using the alternative parameterization of the negative binomial distribution detailed in stats::rnbinom, is:

\mjsdeqn

y_ij \sim NB(\mu = e^c_1 + b_i, size = 1)

for the first group or arm, where \mjseqnb_i \sim N(0,\sigma_b^2), while for the second group,

\mjsdeqn

y_ij \sim NB(\mu = e^c_2 + b_i, size = 1) where \mjseqnb_i \sim N(0,\sigma_b_2^2); if \mjseqn\sigma_b_2^2 is not specified, then the second group uses \mjseqnb_i \sim N(0,\sigma_b^2).

Non-convergent models are not included in the calculation of exact confidence intervals.

Value

If nofit = F, a list with the following components

• Character string indicating total number of simulations, distribution of simulated data, and regression family

• Number of simulations

• Data frame with columns "Power" (Estimated statistical power), "lower.95.ci" (Lower 95 "upper.95.ci" (Upper 95 Note that non-convergent models are returned for review, but not included in this calculation.

• Analytic method used for power estimation

• Data frame containing families for distribution and analysis of simulated data

• Significance level

• Vector containing user-defined cluster sizes

• Vector containing user-defined number of clusters

• Data frame reporting \mjseqn\sigma_b^2 for each group

• Vector containing expected events per unit time and risk ratios based on user inputs

• Data frame with columns: "Estimate" (Estimate of treatment effect for a given simulation), "Std.err" (Standard error for treatment effect estimate), "Test.statistic" (z-value (for GLMM) or Wald statistic (for GEE)), "p.value", "converge" (Did model converge for that set of simulated data?)

• If allSimData = TRUE, a list of data frames, each containing: "y" (Simulated response value), "trt" (Indicator for treatment arm), "clust" (Indicator for cluster)

• Logical vector reporting whether models converged.

If nofit = T, a data frame of the simulated data sets, containing:

• "arm" (Indicator for treatment arm)

• "cluster" (Indicator for cluster)

• "y1" ... "yn" (Simulated response value for each of the nsim data sets).

Testing details

This function has been verified, where possible, against reference values from PASS11, CRTsize::n4incidence, clusterPower::cps.count, and clusterPower::cpa.count.

Author(s)

Alexander R. Bogdan, Alexandria C. Sakrejda ([e-mail]), and Ken Kleinman ([e-mail])

Examples

# Estimate power for a trial with 10 clusters in each arm with 20 subjects each,
# with sigma_b_sq = 0.1 in both arms. We expect mean event rates per unit time of
# 20 and 30 in the first and second arms, respectively, and we use 100 simulated
# data sets analyzed by the GEE method.

## Not run: 
count.sim = cps.count(nsim = 100, nsubjects = 20, nclusters = 10,
                      c1 = 20, c2 = 30, sigma_b_sq = 0.1,
                      family = 'poisson', analysis = 'poisson',
                      method = 'gee', alpha = 0.05, quiet = FALSE,
                      allSimData = FALSE, seed = 123)

## End(Not run)
# The resulting estimated power (if you set seed = 123) should be about 0.8.



# Estimate power for a trial with 10 clusters and 10 subjects per cluster in the
# first arm, 20 clusters and 20 subjects per cluster in the second, and
# sigma_b_sq = 0.1 in both arms. We expect mean event rates per unit time of
# 20 and 30 in the first and second arms, respectively, and we use 100 simulated
# data sets analyzed by the GLMM method.

## Not run: 
count.sim = cps.count(nsim = 100, nsubjects = c(10,20), nclusters = c(10,10),
                      c1 = 20, c2 = 30, sigma_b_sq = 0.1,
                      family = 'poisson', analysis = 'poisson',
                      method = 'glmm', alpha = 0.05, quiet = FALSE,
                      allSimData = FALSE, seed = 123)

## End(Not run)
# The resulting estimated power (if you set seed = 123) should be about 0.85.



# Estimate power for a trial with 5 clusters in the first arm, those clusters having
# 4, 5, 6, 7, and 7 subjects each, and 10 clusters in the second arm, those
# clusters having 5 subjects each, with sigma_b_sq = 0.1 in the first arm and
# sigma_b_sq2 = .05 in the second arm. We expect mean event rates per unit time
# of 20 and 30 in the first and second arms, respectively, and we use 100 simulated
# data sets analyzed by the GLMM method.

## Not run: 
count.sim = cps.count(nsim = 100, nsubjects = c(4, 5, 6, 7, 7, rep(5, times = 10)),
                      nclusters = c(5,10),
                      c1 = 20, c2 = 30,
                      sigma_b_sq = 0.1, sigma_b_sq2 = 0.05,
                      family = 'poisson', analysis = 'poisson',
                      method = 'glmm', alpha = 0.05, quiet = FALSE,
                      allSimData = FALSE, seed = 123)

## End(Not run)
# The resulting estimated power (if you set seed = 123) should be about 0.75.

nickreich/clusterPower documentation built on Feb. 3, 2021, 6:54 p.m.