#' Simulation-based power estimation for cluster-randomized trials: Parallel Designs, Count Outcome with multiple arms
#'
#'
#'
#' @description
#' \loadmathjax
#'
#'
#' This function uses Monte Carlo methods (simulations) to estimate
#' power for cluster-randomized trials for integer-valued outcomes with two or more
#' trial conditions. 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, progress updates, 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.
#'
#'
#'
#' Users must specify the desired number of simulations, number of subjects per
#' cluster, number of clusters per treatment arm, group probabilities, and the
#' between-cluster variance. Significance level, analytic method, whether progress
#' updates are displayed, poor/singular fit override, and whether or not to return the
#' simulated data may also be specified.
#'
#' This user-friendly function calls an internal function; the internal function
#' can be called
#' directly by the user to return the fitted models rather than the power
#' summaries (see \code{?cps.ma.count.internal} for details).
#'
#' Users can spread the simulated data
#' generation and model fitting tasks across multiple cores using the
#' \code{cores} argument. Users should expect that parallel computing may make
#' model fitting faster than using a single core for more complicated models.
#' For simpler models, users may prefer to use single thread computing
#' (\code{cores}=1), as the processes involved in allocating memory and
#' copying data across cores also may take some time. For time-savings,
#' this function stops execution early if estimated power < 0.5 or more
#' than 25\% of models produce a singular fit or non-convergence warning
#' message, unless \code{poorFitOverride = TRUE}.
#'
#'
#'
#'
#'
#' @param nsim Number of datasets to simulate; accepts integer. Required.
#'
#' @param nsubjects Number of subjects per cluster; accepts an
#' integer (implying equal cluster sizes in all arms) if \code{narms}
#' and \code{nclusters} are provided. Alternately, a list with one integer per arm (if the
#' cluster sizes are the same within the arm), or a list of vectors where each vector represents an arm
#' and each entry in the vector is the number of subjects per cluster (if the cluster sizes are not the same
#' within the arms). Required.
#'
#' @param narms Number of trial arms; accepts integer. Required.
#'
#' @param 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 integers
#' representing the number of clusters in each arm (if nsubjects differs between arms).
#' If a list of vectors of cluster sizes is provided in \code{nsubjects}, then
#' the vector of cluster counts must match the length of the \code{nsubjects} vectors.
#' Required.
#'
#' @param counts Mean event per unit time for each arm; accepts a scalar
#' (if all arms have the same event rate) or
#' a vector of length \code{narms}. Required.
#'
#' @param family Distribution from which responses are simulated. Accepts Poisson
#' (\code{'poisson'}) or negative binomial (\code{'neg.binom'}); default = 'poisson'. Required.
#'
#' @param analysis Family used for data analysis; currently only applicable when \code{method = 'glmm'}.
#' Accepts c('poisson', 'neg.binom'); default = 'poisson'. Required.
#'
#' @param 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 positive and defaults to 1. Required when
#' family = 'neg.binom'.
#'
#' @param sigma_b_sq Between-cluster variance for each arm; accepts a scalar
#' (if all arms have the same between-cluster variance) or a vector of length
#' \code{narms}. Required.
#'
#' @param alpha The level of significance of the test, the probability of a
#' Type I error. Default = 0.05.
#'
#' @param quiet When set to \code{FALSE}, displays simulation progress and estimated
#' completion time. Default = \code{FALSE}.
#'
#' @param method Data analysis method, either generalized linear mixed effects model
#' (GLMM) or generalized estimating equations (GEE). Accepts \code{c('glmm', 'gee')};
#' default = \code{'glmm'}. Required.
#'
#' @param multi_p_method A string indicating the method to use for adjusting
#' p-values for multiple comparisons. Choose one of "holm", "hochberg",
#' "hommel", "bonferroni", "BH", "BY", "fdr", "none". The default is
#' "bonferroni". See \code{?p.adjust} for additional details.
#'
#' @param allSimData Option to include a list of all simulated datasets in the output object.
#' Default = \code{FALSE}.
#'
#' @param seed Option to set the seed. Default is NULL.
#'
#' @param cores Number of cores to be used for parallel computing. Accepts a
#' string ("all"), NA (no parallel computing), or scalar value indicating
#' the number of CPUs to use. Default = NA.
#'
#' @param tdist Logical value indicating whether cluster-level random effects
#' should be drawn from a \mjseqn{t} distribution rather than a normal distribution.
#' Default = \code{FALSE}.
#'
#' @param poorFitOverride Option to override \code{stop()} if more than 25\% of fits fail to converge;
#' default = \code{FALSE}.
#'
#' @param lowPowerOverride Option to override \code{stop()} if the power
#' is less than 0.5 after the first 50 simulations and every ten simulations
#' thereafter. On \code{stop}, the power calculated from the completed simulations is printed in the
#' stop message. Default = \code{FALSE}. When \code{TRUE}, this check is ignored and the
#' calculated power is returned regardless of value.
#'
#' @param timelimitOverride Logical. When FALSE, stops execution if the estimated
#' completion time is more than 2 minutes. Defaults to TRUE.
#'
#' @param return.all.models Logical; Returns all of the fitted models, the simulated data,
#' the overall model comparisons, and the convergence report vector. This is equivalent
#' to the output of cps.ma.count.internal(). See ?cps.ma.count.internal() for details.
#'
#' @param nofit Option to skip model fitting and analysis and return the simulated data.
#' Defaults to \code{FALSE}.
#'
#' @param opt Optimizer for model fitting, from the package \code{optimx} or \code{nloptwrap}.
#' Default is 'NLOPT_LN_BOBYQA'.
#'
#'
#'
#'
#' @return A list with the following components:
#' \describe{
#' \item{power}{Data frame with columns "power" (Estimated statistical power),
#' "lower.95.ci" (Lower 95\% confidence interval bound),
#' "upper.95.ci" (Upper 95\% confidence interval bound).}
#' \item{model.estimates}{Data frame with columns corresponding
#' to each arm with descriptive suffixes as follows:
#' ".Estimate" (Estimate of treatment effect for a given
#' simulation),
#' "Std.Err" (Standard error for treatment effect estimate),
#' ".zval" (for GLMM) | ".wald" (for GEE), and
#' ".pval" (the p-value estimate).}
#' \item{overall.power}{Table of F-test (when method="glmm") or chi^{2}
#' (when method="gee") significance test results.}
#' \item{overall.power.summary}{Summary overall power of treatment model
#' compared to the null model.}
#' \item{sim.data}{Produced when allSimData==TRUE. List of \code{nsim}
#' data frames, each containing:
#' "y" (simulated response value),
#' "trt" (indicator for treatment group or arm), and
#' "clust" (indicator for cluster).}
#' \item{model.fit.warning.percent}{Character string containing the percent
#' of \code{nsim} in which the glmm fit was singular or failed to converge,
#' produced only when method = "glmm" & allSimData = FALSE.
#' }
#' \item{model.fit.warning.incidence}{Vector of length \code{nsim} denoting
#' whether or not a simulation glmm fit triggered a "singular fit" or
#' "non-convergence" error, produced only when method = "glmm" &
#' allSimData=TRUE.
#' }
#' }
#' If \code{nofit = T}, a data frame of the simulated data sets, containing:
#' \itemize{
#' \item "arm" (Indicator for treatment arm)
#' \item "cluster" (Indicator for cluster)
#' \item "y1" ... "yn" (Simulated response value for each of the \code{nsim} data sets).
#' }
#'
#'
#' @details
#'
#' If \code{family = 'poisson'}, the data generating model is:
#' \mjsdeqn{y_{ijk} \sim \code{Poisson}(e^{c_k + b_{jk}}) }
#' for observation \mjseqn{i}, in cluster \mjseqn{j}, in treatment arm \mjseqn{k}, where \mjseqn{b_{jk}\sim N(0,\sigma^2_{b_{k}})}.
#'
#' If \code{family = 'neg.bin'}, the data generating model, using the
#' alternative parameterization of the negative binomial distribution
#' detailed in \code{stats::rnbinom}, is:
#'
#' \mjsdeqn{y_{ijk} \sim \code{NB}(\mu = e^{c_k + b_{jk}}, \code{size} = 1) }
#'
#' for observation \mjseqn{i}, in cluster \mjseqn{j}, in treatment arm \mjseqn{k}, where \mjseqn{b_{jk}\sim N(0,\sigma^2_{b_{k}})}.
#'
#'
#' Non-convergent models are not included in the calculation of exact confidence
#' intervals.
#'
#' For complicated models, we recommend using parallel processing with the \code{cores="all"} argument.
#' For simpler models, users may prefer to use single thread computing
#' (\code{cores}=1), as the processes involved in allocating memory and
#' copying data across cores also may take some time.
#'
#' By default, this function stops execution early if estimated power < 0.5 or if more
#' than 25\% of models produce a singular fit or non-convergence warning. In some cases, users
#' may want to ignore singularity warnings (see \code{?isSingular}) by setting \code{poorFitOverride = TRUE}.
#'
#'
#'
#' @author Alexandria C. Sakrejda (\email{acbro0@@umass.edu})
#' @author Alexander R. Bogdan
#' @author Ken Kleinman (\email{ken.kleinman@@gmail.com})
#'
#' @examples
#'
#' # For a 3-arm trial with 4, 4, and 5 clusters in each arm, respectively,
#' # specify the number of subjects in each cluster with 3 vectors in a list,
#' # each vector representing a study arm. For each cluster, in no particular
#' # order, denote the number of subjects. In this example, the first arm
#' # contains 150, 200, 50, and 100 subjects in each of the 4 clusters. The second
#' # arm contains 50, 150, 210, and 100 subjects in each of 4 clusters, while
#' # the third arm contains 70, 200, 150, 50, and 100 subjects in each of 5
#' # clusters. The expected outcomes for each arm are 10, 55, and 65, and
#' # the sigma_b_sq values are 1, 1, and 2, respectively. Assuming
#' # seed = 123, the overall power for this trial should be 0.81.
#'
#' \dontrun{
#' nsubjects.example <- list(c(150, 200, 50, 100), c(50, 150, 210, 100),
#' c(70, 200, 150, 50, 100))
#' counts.example <- c(10, 55, 65)
#' sigma_b_sq.example <- c(1, 1, 2)
#'
#' count.ma.rct.unbal <- cps.ma.count(nsim = 100,
#' narms = 3,
#' nsubjects = nsubjects.example,
#' counts = counts.example,
#' sigma_b_sq = sigma_b_sq.example,
#' alpha = 0.05, seed = 123)
#'}
#'
#' # For a different trial with 4 arms, each arm has 4 clusters which
#' # each contain 100 subjects. Expected counts for each arm are 30
#' # for the first arm, 35 for the second, 70 for the third, and 40
#' # for the fourth. Similarly, sigma_b_sq for each arm are 1
#' # for the first arm, 1.2 for the second, 1 for the third, and 0.9
#' # for the fourth. Assuming seed = 123, the overall power for this
#' # trial should be 0.84
#'
#' \dontrun{
#' count.ma.rct.bal <- cps.ma.count(nsim = 10, nsubjects = 100, narms = 4,
#' nclusters = 25, counts = c(30, 35, 70, 40),
#' sigma_b_sq = c(1, 1.2, 1, 0.9), seed = 123)
#'}
#' @export
cps.ma.count <- function(nsim = 1000,
nsubjects = NULL,
narms = NULL,
nclusters = NULL,
counts = NULL,
family = "poisson",
analysis = "poisson",
negBinomSize = 1,
sigma_b_sq = NULL,
alpha = 0.05,
quiet = FALSE,
method = 'glmm',
multi_p_method = "bonferroni",
allSimData = FALSE,
seed = NA,
cores = NA,
tdist = FALSE,
poorFitOverride = FALSE,
lowPowerOverride = FALSE,
timelimitOverride = TRUE,
return.all.models = FALSE,
nofit = FALSE,
opt = "NLOPT_LN_BOBYQA") {
converge <- NULL
# use this later to determine total elapsed time
start.time <- Sys.time()
if (is.null(narms)) {
stop("ERROR: narms is required.")
}
if (narms < 3) {
stop("ERROR: LRT significance not calculable when narms < 3. Use cps.count() instead.")
}
# create nclusters if not provided directly by user
if (isTRUE(is.list(nsubjects))) {
if (is.null(nclusters)) {
nclusters <- vapply(nsubjects, length, 0)
}
}
if (length(nclusters) == 1 & !isTRUE(is.list(nsubjects))) {
nclusters <- rep(nclusters, narms)
}
if (length(nclusters) > 1 & length(nsubjects) == 1) {
narms <- length(nclusters)
}
# input validation steps
if (!is.wholenumber(nsim) || nsim < 1 || length(nsim) > 1) {
stop("nsim must be a positive integer of length 1.")
}
if (isTRUE(is.null(nsubjects))) {
stop("nsubjects must be specified. See ?cps.ma.count for help.")
}
if (length(nsubjects) == 1 & !isTRUE(is.numeric(nclusters))) {
stop("When nsubjects is scalar, user must supply nclusters (clusters per arm)")
}
if (length(nsubjects) == 1 & length(nclusters) == 1 &
!isTRUE(is.numeric(narms))) {
stop("User must provide narms when nsubjects and nclusters are both scalar.")
}
# nclusters must be whole numbers
if (sum(is.wholenumber(nclusters) == FALSE) != 0 ||
nclusters < 1) {
stop("nclusters must be postive integer values.")
}
# nsubjects must be whole numbers
if (sum(is.wholenumber(unlist(nsubjects)) == FALSE) != 0 ||
unlist(nsubjects) < 1) {
stop("nsubjects must be positive integer values.")
}
# Generate nclusters vector when a scalar is provided but nsubjects is a vector
if (length(nclusters) == 1 & length(nsubjects) > 1) {
nclusters <- rep(nclusters, length(nsubjects))
}
# Create nsubjects structure from narms and nclusters
if (mode(nsubjects) == "list") {
str.nsubjects <- nsubjects
} else {
if (length(nsubjects) == 1) {
str.nsubjects <- lapply(nclusters, function(x)
rep(nsubjects, x))
} else {
if (length(nsubjects) != narms) {
stop("nsubjects must be length 1 or length narms if not provided in a list.")
}
str.nsubjects <- list()
for (i in 1:length(nsubjects)) {
str.nsubjects[[i]] <- rep(nsubjects[i], times = nclusters[i])
}
}
}
# allow entries to be entered as text for Shiny app
if (!is.numeric(counts)) {
counts <- as.numeric(unlist(strsplit(counts, split = ", ")))
}
if (!is.numeric(sigma_b_sq)) {
sigma_b_sq <- as.numeric(unlist(strsplit(sigma_b_sq, split = ", ")))
}
# allows for counts, sigma_b_sq to be entered as scalar
if (length(sigma_b_sq) == 1) {
sigma_b_sq <- rep(sigma_b_sq, narms)
}
if (length(counts) == 1) {
counts <- rep(counts, narms)
}
if (length(counts) != narms) {
stop(
"Length of counts must equal narms, or be provided as a scalar if counts for all arms are equal."
)
}
if (length(sigma_b_sq) != narms) {
stop(
"Length of variance parameters sigma_b_sq must equal narms, or be provided as a scalar
if sigma_b_sq for all arms are equal."
)
}
# Set warnings to OFF
# Note: Warnings will still be stored in 'warning.list'
options(warn = -1)
# run the simulations
count.ma.rct <- cps.ma.count.internal(
nsim = nsim,
str.nsubjects = str.nsubjects,
counts = counts,
sigma_b_sq = sigma_b_sq,
alpha = alpha,
quiet = quiet,
method = method,
all.sim.data = allSimData,
seed = seed,
poor.fit.override = poorFitOverride,
low.power.override = lowPowerOverride,
timelimitOverride = timelimitOverride,
tdist = tdist,
cores = cores,
family = family,
analysis = analysis,
negBinomSize = negBinomSize,
nofit = nofit,
opt = opt
)
# Set warnings to ON
# Note: Warnings will still be stored in 'warning.list'
options(warn = 0)
#option to return simulated data only
if (nofit == TRUE || return.all.models == TRUE) {
return(count.ma.rct)
}
# Create object containing summary statement
summary.message = paste0(
"Monte Carlo Power Estimation based on ",
nsim,
" Simulations: Parallel Design, Poisson Outcome, ",
narms,
" Arms."
)
# Create method object
long.method = switch(method, glmm = 'Generalized Linear Mixed Model',
gee = 'Generalized Estimating Equation')
# Create object containing group-specific variance parameters
armnames <- vector(mode = "character", length = narms)
for (i in 1:narms) {
armnames[i] <- paste0("Arm.", i)
}
var.parms = data.frame('sigma_b_sq' = sigma_b_sq,
"counts" = counts)
rownames(var.parms) <- armnames
# Create object containing group-specific cluster sizes
names(str.nsubjects) <- armnames
cluster.sizes <- str.nsubjects
models <- count.ma.rct[[1]]
#Organize output for GLMM
if (method == "glmm") {
Estimates = matrix(NA, nrow = nsim, ncol = narms)
std.error = matrix(NA, nrow = nsim, ncol = narms)
z.val = matrix(NA, nrow = nsim, ncol = narms)
p.val = matrix(NA, nrow = nsim, ncol = narms)
for (i in 1:nsim) {
Estimates[i, ] <- models[[i]][[10]][, 1]
std.error[i, ] <- models[[i]][[10]][, 2]
z.val[i, ] <- models[[i]][[10]][, 3]
p.val[i, ] <-
p.adjust(models[[i]][[10]][, 4], method = multi_p_method)
}
# Organize the row/col names for the model estimates output
keep.names <- rownames(models[[1]][[10]])
names.Est <- rep(NA, narms)
names.st.err <- rep(NA, narms)
names.zval <- rep(NA, narms)
names.pval <- rep(NA, narms)
names.power <- rep(NA, narms)
for (i in 1:length(keep.names)) {
names.Est[i] <- paste(keep.names[i], ".Estimate", sep = "")
names.st.err[i] <- paste(keep.names[i], ".Std.Err", sep = "")
names.zval[i] <- paste(keep.names[i], ".zval", sep = "")
names.pval[i] <- paste(keep.names[i], ".pval", sep = "")
names.power[i] <- paste(keep.names[i], ".power", sep = "")
}
colnames(Estimates) <- names.Est
colnames(std.error) <- names.st.err
colnames(z.val) <- names.zval
colnames(p.val) <- names.pval
#convergence
converge <- as.vector(rep(NA, times = nsim))
for (i in 1:nsim) {
converge[i] <-
ifelse(is.null(count.ma.rct[[1]][[i]]$optinfo$conv$lme4$messages),
TRUE,
FALSE)
}
cps.model.temp <- data.frame(converge, p.val)
colnames(cps.model.temp)[1] <- "converge"
cps.model.temp2 <-
dplyr::filter(cps.model.temp, converge == TRUE)
if (isTRUE(nrow(cps.model.temp2) < (.25 * nsim))) {
warning(paste0(
nrow(cps.model.temp2),
" models converged. Check model parameters."
),
immediate. = TRUE)
}
# Organize the LRT output
LRT.holder <-
matrix(
unlist(count.ma.rct[[2]]),
ncol = 3,
nrow = nsim,
byrow = TRUE,
dimnames = list(seq(1:nsim),
colnames(count.ma.rct[[2]][[1]]))
)
# Proportion of times P(>F)
LRT.holder <- cbind(LRT.holder, cps.model.temp["converge"])
LRT.holder <- LRT.holder[LRT.holder[, "converge"] == TRUE, ]
sig.LRT <- ifelse(LRT.holder[, 3] < alpha, 1, 0)
# Calculate and store power estimate & confidence intervals
sig.val <- ifelse(p.val < alpha, 1, 0)
pval.power <- apply(sig.val, 2, sum)
# Calculate and store power estimate & confidence intervals
power.parms <- cbind(confintCalc(
alpha = alpha,
multi = TRUE,
p.val = as.vector(cps.model.temp2[, 3:length(cps.model.temp2)])
),
multi_p_method)
# Store simulation output in data frame
ma.model.est <-
data.frame(Estimates, std.error, z.val, p.val, count.ma.rct[[3]])
ma.model.est <-
ma.model.est[, -grep('.*ntercept.*', names(ma.model.est))]
# performance messages
total.est <-
as.numeric(difftime(Sys.time(), start.time, units = 'secs'))
hr.est <- total.est %/% 3600
min.est <- total.est %/% 60
sec.est <- round(total.est %% 60, 0)
message(
paste0(
"Simulations Complete! Time Completed: ",
Sys.time(),
"\nTotal Runtime: ",
hr.est,
'Hr:',
min.est,
'Min:',
sec.est,
'Sec'
)
)
## Output objects for GLMM
# Create list containing all output (class 'crtpwr.ma') and return
if (allSimData == TRUE && return.all.models == FALSE) {
complete.output = structure(
list(
"overview" = summary.message,
"nsim" = nsim,
"power" =
prop_H0_rejection(
alpha = alpha,
nsim = nsim,
sig.LRT = sig.LRT
),
"beta" = power.parms['Beta'],
"overall.power2" = power.parms,
"overall.power" = LRT.holder,
"method" = long.method,
"alpha" = alpha,
"cluster.sizes" = cluster.sizes,
"n.clusters" = nclusters,
"variance.parms" = var.parms,
"counts" = counts,
"model.estimates" = ma.model.est,
"convergence" = count.ma.rct[[3]],
"sim.data" = count.ma.rct[[4]]
),
class = 'crtpwr.ma'
)
}
if (return.all.models == TRUE) {
complete.output = structure(
list(
"overview" = summary.message,
"nsim" = nsim,
"power" =
prop_H0_rejection(
alpha = alpha,
nsim = nsim,
sig.LRT = sig.LRT
),
"beta" = power.parms['Beta'],
"overall.power2" = power.parms,
"overall.power" = LRT.holder,
"method" = long.method,
"alpha" = alpha,
"cluster.sizes" = cluster.sizes,
"n.clusters" = nclusters,
"variance.parms" = var.parms,
"counts" = counts,
"model.estimates" = ma.model.est,
"convergence" = count.ma.rct[[3]],
"sim.data" = count.ma.rct[[4]],
"all.models" <- count.ma.rct
),
class = 'crtpwr.ma'
)
}
if (return.all.models == FALSE && allSimData == FALSE) {
complete.output = structure(
list(
"overview" = summary.message,
"nsim" = nsim,
"power" =
prop_H0_rejection(
alpha = alpha,
nsim = nsim,
sig.LRT = sig.LRT
),
"beta" = power.parms['Beta'],
"overall.power2" = power.parms,
"overall.power" = LRT.holder,
"method" = long.method,
"alpha" = alpha,
"cluster.sizes" = cluster.sizes,
"n.clusters" = nclusters,
"variance.parms" = var.parms,
"counts" = counts,
"model.estimates" = ma.model.est,
"convergence" = count.ma.rct[[3]]
),
class = 'crtpwr.ma'
)
}
} # end of GLMM options
#Organize output for GEE method
if (method == "gee") {
# Organize the output
Estimates = matrix(NA, nrow = nsim, ncol = narms)
std.error = matrix(NA, nrow = nsim, ncol = narms)
Wald = matrix(NA, nrow = nsim, ncol = narms)
Pr = matrix(NA, nrow = nsim, ncol = narms)
for (i in 1:nsim) {
Estimates[i, ] <- models[[i]]$coefficients[, 1]
std.error[i, ] <- models[[i]]$coefficients[, 2]
Wald[i, ] <- models[[i]]$coefficients[, 3]
Pr[i, ] <- models[[i]]$coefficients[, 4]
}
# Organize the row/col names for the output
keep.names <- rownames(models[[1]]$coefficients)
names.Est <- rep(NA, length(narms))
names.st.err <- rep(NA, length(narms))
names.wald <- rep(NA, length(narms))
names.pval <- rep(NA, length(narms))
names.power <- rep(NA, length(narms))
for (i in 1:length(keep.names)) {
names.Est[i] <- paste(keep.names[i], ".Estimate", sep = "")
names.st.err[i] <- paste(keep.names[i], ".Std.Err", sep = "")
names.wald[i] <- paste(keep.names[i], ".wald", sep = "")
names.pval[i] <- paste(keep.names[i], ".pval", sep = "")
names.power[i] <- paste(keep.names[i], ".power", sep = "")
}
colnames(Estimates) <- names.Est
colnames(std.error) <- names.st.err
colnames(Wald) <- names.wald
colnames(Pr) <- names.pval
# Organize the LRT output
LRT.holder <-
matrix(
unlist(count.ma.rct[[2]]),
ncol = 3,
nrow = nsim,
byrow = TRUE,
dimnames = list(seq(1:nsim),
c("Df", "X2", "P(>|Chi|)"))
)
# Proportion of times P(>F)
LRT.holder <- cbind(LRT.holder, count.ma.rct[["converged"]])
LRT.holder <- LRT.holder[LRT.holder[, "converged"] == TRUE, ]
sig.LRT <- ifelse(LRT.holder[, 3] < alpha, 1, 0)
# Calculate and store power estimate & confidence intervals
power.parms <- cbind(confintCalc(
alpha = alpha,
multi = TRUE,
p.val = Pr[, 2:narms]
),
multi_p_method)
# Store GEE simulation output in data frame
ma.model.est <- data.frame(Estimates, std.error, Wald, Pr)
ma.model.est <-
ma.model.est[, -grep('.*ntercept.*', names(ma.model.est))]
# performance messages
total.est <-
as.numeric(difftime(Sys.time(), start.time, units = 'secs'))
hr.est <- total.est %/% 3600
min.est <- total.est %/% 60
sec.est <- round(total.est %% 60, 0)
message(
paste0(
"Simulations Complete! Time Completed: ",
Sys.time(),
"\nTotal Runtime: ",
hr.est,
'Hr:',
min.est,
'Min:',
sec.est,
'Sec'
)
)
# Create list containing all output (class 'crtpwr.ma') and return
if (allSimData == TRUE & return.all.models == FALSE) {
complete.output = structure(
list(
"overview" = summary.message,
"nsim" = nsim,
"power" =
prop_H0_rejection(
alpha = alpha,
nsim = nsim,
sig.LRT = sig.LRT
),
"beta" = power.parms['Beta'],
"overall.power2" = power.parms,
"overall.power" = LRT.holder,
"method" = long.method,
"alpha" = alpha,
"cluster.sizes" = cluster.sizes,
"n.clusters" = nclusters,
"variance.parms" = var.parms,
"counts" = counts,
"model.estimates" = ma.model.est,
"sim.data" = count.ma.rct[[3]]
),
class = 'crtpwr.ma'
)
}
if (return.all.models == TRUE) {
complete.output = structure(
list(
"overview" = summary.message,
"nsim" = nsim,
"power" =
prop_H0_rejection(
alpha = alpha,
nsim = nsim,
sig.LRT = sig.LRT
),
"beta" = power.parms['Beta'],
"overall.power2" = power.parms,
"overall.power" = LRT.holder,
"method" = long.method,
"alpha" = alpha,
"cluster.sizes" = cluster.sizes,
"n.clusters" = nclusters,
"variance.parms" = var.parms,
"counts" = counts,
"model.estimates" = ma.model.est,
"all.models" <- count.ma.rct
),
class = 'crtpwr.ma'
)
}
if (return.all.models == FALSE && allSimData == FALSE) {
complete.output = structure(
list(
"overview" = summary.message,
"nsim" = nsim,
"power" =
prop_H0_rejection(
alpha = alpha,
nsim = nsim,
sig.LRT = sig.LRT
),
"beta" = power.parms['Beta'],
"overall.power2" = power.parms,
"overall.power" = LRT.holder,
"method" = long.method,
"alpha" = alpha,
"cluster.sizes" = cluster.sizes,
"n.clusters" = nclusters,
"variance.parms" = var.parms,
"counts" = counts,
"model.estimates" = ma.model.est
),
class = 'crtpwr.ma'
)
}
}
return(complete.output)
}# end of fxn
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.