get.rejection.rates: Get rejection rates of the p value and fragility index based...

In brb225/FragilityTools: Toolbox For Patient Count Fragility Measures

Description

This function is used to calculate a monte carlo estimate of the size and power of the p value based and fragility index based statistical test.

Usage

get.rejection.rates(
  get.p.val,
  get.sample.null = NULL,
  get.sample.alt = NULL,
  phi,
  n,
  alpha = 0.05,
  algorithm = "greedy",
  cl = NULL,
  nsim = 1000
)

Arguments

get.p.val	a function which inputs a contingency table and outputs a p value
get.sample.null	a function which inputs a sample size and outputs a sample under the null hypothesis. By default NULL, which indicates to skip calculations under the null hypothesis
get.sample.alt	a function which inputs a sample size and outputs a sample under the alternative hypothesis. By default NULL, which indicates to skip calculations under the alternative hypothesis
phi	a numeric for the fragility index quantile
n	a numeric for the sample size
alpha	a numeric for the significance cutoff
algorithm	a string indicating the algorithm to use to calculate the fragility index. By default 'greedy'. See bin.fi for more options.
cl	a cluster from the parallel package, not currently supported
nsim	a numeric for the number of monte carlo simulations.

Value

a length 4 vector with the size and power of both tests.

Examples

get.p.val <- function(mat) {
  pearson.test.stata <- function(mat) {
    n1 <- sum(mat[1, ])
    n2 <- sum(mat[2, ])
    p1 <- mat[1, 1] / n1
    p2 <- mat[2, 1] / n2

    pbar <- (n1 * p1 + n2 * p2) / (n1 + n2)
    ts <- (p1 - p2) / sqrt(pbar * (1 - pbar) * (1 / n1 + 1 / n2))
    p_value <- 2 * pnorm(abs(ts), lower.tail = FALSE)
    return(ifelse(is.nan(p_value), 1, p_value))
  }
  suppressWarnings(p <- pearson.test.stata(mat))

  return(ifelse(is.na(p), 1, p))
}

get.sample.null <- function(ss) draw.binom(ss, theta1 = .14, theta2 = .14, row.prop = 1 / 2, matrix = TRUE)
get.sample.alt <- function(ss) draw.binom(ss, theta1 = .14, theta2 = .08, row.prop = 1 / 2, matrix = TRUE)

get.rejection.rates(get.p.val, get.sample.null, get.sample.alt, phi=2, n=100, algorithm = "walsh")

brb225/FragilityTools documentation built on Jan. 21, 2022, 1:26 a.m.