tests/testthat/test-oc2S.R

In RBesT: R Bayesian Evidence Synthesis Tools

## test the analytical OC function via brute force simulation
set.seed(12354)

prior1 <- mixnorm(c(0.3, -0.2, 2), c(0.7, 0, 50), sigma = 1)
prior2 <- mixnorm(c(1.0, 0, 50), sigma = 1)

N1 <- 10
N2 <- 20

## type I error fairly large to 20% to make it easier to test (less
## simulations needed for accurate results)
pcrit <- 0.80
qcrit <- 0

## theta2 set such that we have about 75% power under this truth
theta1 <- 0
theta2 <- 0.5

Nsim <- 1e4

run_on_cran <- function() {
  if (identical(Sys.getenv("NOT_CRAN"), "true")) {
    return(FALSE)
  }
  return(TRUE)
}

oc2S_normal_MC <- function(
  prior1,
  prior2,
  N1,
  N2,
  theta1,
  theta2,
  pcrit = 0.975,
  qcrit = 0
) {
  mean_sd1 <- sigma(prior1) / sqrt(N1)
  mean_sd2 <- sigma(prior2) / sqrt(N2)

  mean_prior1 <- prior1
  sigma(mean_prior1) <- mean_sd1
  mean_prior2 <- prior2
  sigma(mean_prior2) <- mean_sd2

  mean_samp1 <- rnorm(Nsim, theta1, mean_sd1)
  mean_samp2 <- rnorm(Nsim, theta2, mean_sd2)

  dec <- rep(NA, Nsim)

  for (i in 1:Nsim) {
    post1 <- postmix(mean_prior1, m = mean_samp1[i], se = mean_sd1)
    post2 <- postmix(mean_prior2, m = mean_samp2[i], se = mean_sd2)
    dec[i] <- as.numeric(pmix(RBesT:::mixnormdiff(post1, post2), qcrit) > pcrit)
  }

  mean(dec)
}

Voc2S_normal_MC <- Vectorize(oc2S_normal_MC, c("theta1", "theta2"))

## first test that the analytic difference distribution for normal
## works as expected

test_that("Analytical convolution of normal mixture matches numerical integration result", {
  skip_on_cran()

  pdiff <- RBesT:::mixnormdiff(prior1, prior2)
  x <- seq(-20, 20, length = 21)
  d1 <- dmix(pdiff, x)
  d2 <- dmixdiff(prior1, prior2, x)
  dres <- abs(d1 - d2)
  expect_equal(sum(dres > 1e-5), 0)
  p1 <- pmix(pdiff, x)
  p2 <- pmixdiff(prior1, prior2, x)
  pres <- 100 * abs(p1 - p2)
  expect_equal(sum(pres > 2), 0)
})

## test that the type I error is matching, i.e. is not off by more than 2%
test_that("Type I error is matching between MC and analytical computations in the normal mixture case", {
  skip_on_cran()

  x <- c(-2, 0)
  alpha <- oc2S(
    prior1,
    prior2,
    N1,
    N2,
    decision2S(pcrit, qcrit),
    sigma1 = sigma(prior1),
    sigma2 = sigma(prior2)
  )(x, x)
  alphaMC <- Voc2S_normal_MC(prior1, prior2, N1, N2, x, x, pcrit, qcrit)
  res <- 100 * abs(alpha - alphaMC)
  expect_equal(sum(res > 2), 0)
})


## test that the power is matching, i.e. is not off by more than 2%
test_that("Power is matching between MC and analytical computations in the normal mixture case", {
  skip_on_cran()

  power <- oc2S(
    prior1,
    prior2,
    N1,
    N2,
    decision2S(pcrit, qcrit),
    sigma1 = sigma(prior1),
    sigma2 = sigma(prior2)
  )(theta1, theta2)
  powerMC <- oc2S_normal_MC(
    prior1,
    prior2,
    N1,
    N2,
    theta1,
    theta2,
    pcrit,
    qcrit
  )
  res <- 100 * abs(power - powerMC)
  expect_equal(sum(res > 2), 0)
})

## further test by cross-checking with Gsponer et. al, "A practical
## guide to Bayesian group sequential designs", Pharmaceut. Statist.
## (2014), 13 71-80, Table 1, Probability at interim

test_that("Gsponer et al. results match (normal end-point)", {
  skip_on_cran()

  ocRef <- data.frame(
    delta = c(0, 40, 50, 60, 70),
    success = c(1.1, 32.2, 50.0, 67.6, 82.2),
    futile = c(63.3, 6.8, 2.5, 0.8, 0.2)
  )
  sigmaFixed <- 88

  priorT <- mixnorm(c(1, 0, 0.001), sigma = sigmaFixed, param = "mn")
  priorP <- mixnorm(c(1, -49, 20), sigma = sigmaFixed, param = "mn")

  ## the success criteria is for delta which are larger than some
  ## threshold value which is why we set lower.tail=FALSE
  successCrit <- decision2S(c(0.95, 0.5), c(0, 50), FALSE)
  ## the futility criterion acts in the opposite direction
  futilityCrit <- decision2S(c(0.90), c(40), TRUE)

  nT1 <- 20
  nP1 <- 10

  oc <- data.frame(delta = c(0, 40, 50, 60, 70))

  ## Note that due to the fact that only a single mixture component is
  ## used, the decision boundary is a linear function such that only few
  ## evaluations of the boundary are needed to estimate reliably the
  ## spline function

  ## Table 1, probability for interim for success
  oc$success <- oc2S(
    priorP,
    priorT,
    nP1,
    nT1,
    successCrit,
    Ngrid = 1,
    sigma1 = sigmaFixed,
    sigma2 = sigmaFixed
  )(-49, -49 - oc$delta)

  ## Table 1, probability for interim for futility
  oc$futile <- oc2S(
    priorP,
    priorT,
    nP1,
    nT1,
    futilityCrit,
    Ngrid = 1,
    sigma1 = sigmaFixed,
    sigma2 = sigmaFixed
  )(-49, -49 - oc$delta)

  ## Table 1, first three columns, page 74
  oc[-1] <- lapply(100 * oc[-1], round, 1)

  resFutility <- abs(ocRef$futile - oc$futile)
  resSuccess <- abs(ocRef$success - oc$success)

  expect_equal(sum(resFutility > 2), 0, info = "futility")
  expect_equal(sum(resSuccess > 2), 0, info = "success")
})


## failure when doing repeated evaluations which came up in consulting
test_that("Ensure that repeated oc2S evaluation works for normal case", {
  skip_on_cran()

  samp_sigma <- 3

  n_ia <- 38
  n_final <- 2 * n_ia
  n_ia_to_final <- n_final - n_ia
  sem_ia <- samp_sigma / sqrt(n_ia)

  theta_ctl <- 0
  delta <- 1.04

  obs_P <- 0.11
  obs_T <- 1.28

  prior <- mixnorm(c(1, 0, 0.001), sigma = samp_sigma, param = "mn")
  postP_interim <- postmix(prior, m = obs_P, se = sem_ia)
  postT_interim <- postmix(prior, m = obs_T, se = sem_ia)

  successCrit <- decision2S(c(0.9), c(0), FALSE)

  interim_CP <- oc2S(
    postT_interim,
    postP_interim,
    n_ia_to_final,
    n_ia_to_final,
    successCrit,
    sigma1 = samp_sigma,
    sigma2 = samp_sigma
  )

  cpd_ia <- interim_CP(obs_T, obs_P)
  cpd_ia2 <- interim_CP(theta_ctl + delta, theta_ctl)

  expect_number(cpd_ia, lower = 0, upper = 1, finite = TRUE)
  expect_number(cpd_ia2, lower = 0, upper = 1, finite = TRUE)

  ## check that when calculating directly that the results
  ## are close enough
  interim_CPalt <- oc2S(
    postT_interim,
    postP_interim,
    n_ia_to_final,
    n_ia_to_final,
    successCrit,
    sigma1 = samp_sigma,
    sigma2 = samp_sigma
  )
  cpd_ia2alt <- interim_CPalt(theta_ctl + delta, theta_ctl)
  expect_number(
    abs(cpd_ia2 - cpd_ia2alt),
    lower = 0,
    upper = 1E-3,
    finite = TRUE
  )
})

## test against Schmidli et. al, "Robust Meta-Analytic-Predictive
## Priors", Table 2, unif and beta case
test_that("Schmidli et al. results (binary end-point)", {
  skip_on_cran()

  ocRef_inf <- expand.grid(pc = seq(0.1, 0.6, by = 0.1), delta = c(0, 0.3))
  ocRef_inf$ref <- c(
    0,
    1.6,
    6.1,
    13.7,
    26.0,
    44.4, ## beta/delta=0
    81.6,
    87.8,
    93.4,
    97.9,
    99.6,
    100.0 ## beta/delta=0.3
  ) /
    100

  ocRef_uni <- expand.grid(pc = seq(0.1, 0.6, by = 0.1), delta = c(0, 0.3))
  ocRef_uni$ref <- c(
    1.8,
    2.3,
    2.4,
    2.6,
    2.8,
    2.6, ## unif/delta=0
    89.7,
    82.1,
    79.5,
    79.5,
    81.9,
    89.8 ## unif/delta=0.3
  ) /
    100
  dec <- decision2S(0.975, 0, lower.tail = FALSE)

  N <- 40

  prior_inf <- mixbeta(c(1, 4, 16))
  prior_uni <- mixbeta(c(1, 1, 1))

  N_ctl_uni <- N - round(ess(prior_uni, method = "morita"))
  N_ctl_inf <- N - round(ess(prior_inf, method = "morita"))

  design_uni <- oc2S(prior_uni, prior_uni, N, N_ctl_uni, dec)
  design_inf <- oc2S(prior_uni, prior_inf, N, N_ctl_inf, dec)

  res_uni <- design_uni(ocRef_uni$pc + ocRef_uni$delta, ocRef_uni$pc)
  res_inf <- design_inf(ocRef_inf$pc + ocRef_inf$delta, ocRef_inf$pc)

  expect_true(all(abs(100 * (res_uni - ocRef_uni$ref)) < 2.5))
  expect_true(all(abs(100 * (res_inf - ocRef_inf$ref)) < 2.5))
})

## some additional, very simple type I error tests and tests for the
## discrete case of correct critical value behavior

test_scenario <- function(oc_res, ref) {
  resA <- oc_res - ref
  expect_true(all(abs(resA) < eps))
}

expect_equal_each <- function(test, expected) {
  for (elem in test) {
    expect_equal(elem, expected)
  }
}

## design object, decision function, posterior function must return
## posterior after updatding the prior with the given value; we assume
## that the priors are the same for sample 1 and 2
test_critical_discrete <- function(boundary_design, decision, posterior, y2) {
  lower.tail <- attr(decision, "lower.tail")
  crit <- boundary_design(y2)
  post2 <- posterior(y2)
  if (lower.tail) {
    expect_equal(decision(posterior(crit - 1), post2), 1)
    expect_equal(decision(posterior(crit), post2), 1)
    expect_equal(decision(posterior(crit + 1), post2), 0)
  } else {
    expect_equal(decision(posterior(crit - 1), post2), 0)
    expect_equal(decision(posterior(crit), post2), 0)
    expect_equal(decision(posterior(crit + 1), post2), 1)
  }
}

## expect results to be 1% exact
eps <- 1e-2
alpha <- 0.05

dec <- decision2S(1 - alpha, 0, lower.tail = TRUE)
decB <- decision2S(1 - alpha, 0, lower.tail = FALSE)

## test binary case

beta_prior <- mixbeta(c(1, 1, 1))
if (!run_on_cran()) {
  design_binary <- oc2S(beta_prior, beta_prior, 100, 100, dec)
  boundary_design_binary <- decision2S_boundary(
    beta_prior,
    beta_prior,
    100,
    100,
    dec
  )
  design_binaryB <- oc2S(beta_prior, beta_prior, 100, 100, decB)
  boundary_design_binaryB <- decision2S_boundary(
    beta_prior,
    beta_prior,
    100,
    100,
    decB
  )
} else {
  design_binary <- function(...) {
    return(0.1)
  }
  design_binaryB <- function(...) {
    return(0.1)
  }
  boundary_design_binary <- function(...) {
    return(0.1)
  }
  boundary_design_binaryB <- function(...) {
    return(0.1)
  }
}
posterior_binary <- function(r) postmix(beta_prior, r = r, n = 100)
p_test <- 1:9 / 10
test_that("Binary type I error rate", {
  skip_on_cran()
  test_scenario(design_binary(p_test, p_test), alpha)
})
test_that("Binary crticial value, lower.tail=TRUE", {
  skip_on_cran()
  test_critical_discrete(boundary_design_binary, dec, posterior_binary, 30)
})
test_that("Binary crticial value, lower.tail=FALSE", {
  skip_on_cran()
  test_critical_discrete(boundary_design_binaryB, decB, posterior_binary, 30)
})
test_that("Binary boundary case, lower.tail=TRUE", {
  skip_on_cran()
  expect_numeric(
    design_binary(1, 1),
    lower = 0,
    upper = 1,
    finite = TRUE,
    any.missing = FALSE
  )
})
test_that("Binary boundary case, lower.tail=FALSE", {
  skip_on_cran()
  expect_numeric(
    design_binaryB(0, 0),
    lower = 0,
    upper = 1,
    finite = TRUE,
    any.missing = FALSE
  )
})

## check case where decision never changes due to prior being too
## strong

beta_prior1 <- mixbeta(c(1, 0.9, 1000), param = "mn")
beta_prior2 <- mixbeta(c(1, 0.1, 1000), param = "mn")
design_lower <- oc2S(beta_prior1, beta_prior2, 20, 20, dec) ## always 0
design_upper <- oc2S(beta_prior1, beta_prior2, 20, 20, decB) ## always 1
boundary_design_lower <- decision2S_boundary(
  beta_prior1,
  beta_prior2,
  20,
  20,
  dec
) ## always 0
boundary_design_upper <- decision2S_boundary(
  beta_prior1,
  beta_prior2,
  20,
  20,
  decB
) ## always 1

test_that("Binary case, no decision change, lower.tail=TRUE, critical value", {
  skip_on_cran()
  expect_equal_each(boundary_design_lower(0:20), -1)
})
test_that("Binary case, no decision change, lower.tail=FALSE, critical value", {
  skip_on_cran()
  expect_equal_each(boundary_design_upper(0:20), 21)
})
test_that("Binary case, no decision change, lower.tail=TRUE, frequency=0", {
  skip_on_cran()
  expect_equal_each(design_lower(p_test, p_test), 0.0)
})
test_that("Binary case, no decision change, lower.tail=FALSE, frequency=1", {
  skip_on_cran()
  expect_equal_each(design_upper(p_test, p_test), 1.0)
})


if (!run_on_cran()) {
  design_lower_rev <- oc2S(beta_prior2, beta_prior1, 20, 20, dec) ## always 1
  design_upper_rev <- oc2S(beta_prior2, beta_prior1, 20, 20, decB) ## always 0
  boundary_design_lower_rev <- decision2S_boundary(
    beta_prior2,
    beta_prior1,
    20,
    20,
    dec
  ) ## always 1
  boundary_design_upper_rev <- decision2S_boundary(
    beta_prior2,
    beta_prior1,
    20,
    20,
    decB
  ) ## always 0
} else {
  design_lower_rev <- function(...) {
    return(1)
  }
  design_upper_rev <- function(...) {
    return(0)
  }
  boundary_design_lower_rev <- function(...) {
    return(1)
  }
  boundary_design_upper_rev <- function(...) {
    return(0)
  }
}

test_that("Binary case, no decision change (reversed), lower.tail=TRUE, critical value", {
  skip_on_cran()
  expect_equal_each(boundary_design_lower_rev(0:20), 20)
})
test_that("Binary case, no decision change (reversed), lower.tail=FALSE, critical value", {
  skip_on_cran()
  expect_equal_each(boundary_design_upper_rev(0:20), -1)
})
test_that("Binary case, no decision change (reversed), lower.tail=TRUE, frequency=0", {
  skip_on_cran()
  expect_equal_each(design_lower_rev(p_test, p_test), 1.0)
})
test_that("Binary case, no decision change (reversed), lower.tail=FALSE, frequency=1", {
  skip_on_cran()
  expect_equal_each(design_upper_rev(p_test, p_test), 0.0)
})
test_that("Binary case, log-link", {
  skip_on_cran()
  success <- decision2S(
    pc = c(0.90, 0.50),
    qc = c(log(1), log(0.50)),
    lower.tail = TRUE,
    link = "log"
  )
  prior_pbo <- mixbeta(
    inf1 = c(0.60, 19, 29),
    inf2 = c(0.30, 4, 5),
    rob = c(0.10, 1, 1)
  )
  prior_trt <- mixbeta(c(1, 1 / 3, 1 / 3))
  n_trt <- 50
  n_pbo <- 20
  design_suc <- oc2S(prior_trt, prior_pbo, n_trt, n_pbo, success)
  theta <- seq(0, 1, by = 0.1)
  expect_numeric(
    design_suc(theta, theta),
    lower = 0,
    upper = 1,
    finite = TRUE,
    any.missing = FALSE
  )
})
test_that("Binary case, logit-link", {
  skip_on_cran()
  success <- decision2S(
    pc = c(0.90, 0.50),
    qc = c(log(1), log(0.50)),
    lower.tail = TRUE,
    link = "logit"
  )
  prior_pbo <- mixbeta(
    inf1 = c(0.60, 19, 29),
    inf2 = c(0.30, 4, 5),
    rob = c(0.10, 1, 1)
  )
  prior_trt <- mixbeta(c(1, 1 / 3, 1 / 3))
  n_trt <- 50
  n_pbo <- 20
  design_suc <- oc2S(prior_trt, prior_pbo, n_trt, n_pbo, success)
  theta <- seq(0, 1, by = 0.1)
  expect_numeric(
    design_suc(theta, theta),
    lower = 0,
    upper = 1,
    finite = TRUE,
    any.missing = FALSE
  )
})

## check approximate method

beta_prior <- mixbeta(c(1, 1, 1))
design_binary_eps <- oc2S(beta_prior, beta_prior, 100, 100, dec, eps = 1E-3)
p_test <- seq(0.1, 0.9, by = 0.1)
test_that("Binary type I error rate", {
  skip_on_cran()
  test_scenario(design_binary_eps(p_test, p_test), alpha)
})

## 22 Nov 2017: disabled test as we trigger always calculation of the
## boundaries as of now.
## test_that("Binary results cache expands", {
##               design_binary_eps <- oc2S(beta_prior, beta_prior, 100, 100, dec, eps=1E-3)
##               design_binary_eps(theta1=0.99, theta2=0.8)
##               ## in this case the cache boundaries do not cover the
##               ## critical value
##               expect_true(is.na(design_binary_eps(theta1=0.99, y2=80)))
##               ## while now they do as theta1 is set to 0.1 and 0.9
##               ## internally which triggers recalculation of the
##               ## internal boundaries
##               expect_true(!is.na(design_binary_eps(y2=80)))
##           })

## test poisson case

gamma_prior <- mixgamma(c(1, 2, 2))

design_poisson <- oc2S(gamma_prior, gamma_prior, 100, 100, dec)
design_poissonB <- oc2S(gamma_prior, gamma_prior, 100, 100, decB)
boundary_design_poisson <- decision2S_boundary(
  gamma_prior,
  gamma_prior,
  100,
  100,
  dec
)
boundary_design_poissonB <- decision2S_boundary(
  gamma_prior,
  gamma_prior,
  100,
  100,
  decB
)
posterior_poisson <- function(m) postmix(gamma_prior, m = m / 100, n = 100)
lambda_test <- seq(0.5, 1.3, by = 0.1)
test_that("Poisson type I error rate", {
  skip_on_cran()
  test_scenario(design_poisson(lambda_test, lambda_test), alpha)
})
test_that("Poisson crticial value, lower.tail=TRUE", {
  skip_on_cran()
  test_critical_discrete(boundary_design_poisson, dec, posterior_poisson, 90)
})
test_that("Poisson crticial value, lower.tail=FALSE", {
  skip_on_cran()
  test_critical_discrete(boundary_design_poissonB, decB, posterior_poisson, 90)
})
## 22 Nov 2017: disabled test as we trigger always calculation of the
## boundaries as of now.
## test_that("Poisson results cache expands", {
##              design_poisson  <- oc2S(gamma_prior, gamma_prior, 100, 100, dec)
##              design_poisson(theta1=1, theta2=c(0.7,1))
##              expect_true(sum(is.na(design_poisson(y2=70:90)) ) == 4)
##              expect_true(sum(is.na(design_poisson(theta1=c(0.01, 1), y2=70:90)) ) == 0)
##          })

test_that("Normal OC 2-sample case works for n2=0, crohn-1", {
  crohn_sigma <- 88

  map <- mixnorm(c(0.6, -50, 19), c(0.4, -50, 42), sigma = crohn_sigma)

  ## add a 20% non-informative mixture component
  map_robust <- robustify(map, weight = 0.2, mean = -50, sigma = 88)

  poc <- decision2S(pc = c(0.95, 0.5), qc = c(0, -50), lower.tail = TRUE)

  weak_prior <- mixnorm(c(1, -50, 1), sigma = crohn_sigma, param = "mn")
  n_act <- 40
  ## n_pbo <- 20

  design_noprior_b <- oc2S(
    weak_prior,
    map,
    n_act,
    0,
    poc,
    sigma1 = crohn_sigma,
    sigma2 = crohn_sigma
  )

  expect_numeric(
    design_noprior_b(-20, -30),
    lower = 0,
    upper = 1,
    any.missing = FALSE
  )
})

test_that("Normal OC 2-sample case works for n2=0, crohn-2", {
  crohn_sigma <- 88

  map <- mixnorm(c(1.0, -50, 19), sigma = crohn_sigma)

  ## add a 20% non-informative mixture component
  map_robust <- robustify(map, weight = 0.2, mean = -50, sigma = 88)

  poc <- decision2S(pc = c(0.95, 0.5), qc = c(0, -50), lower.tail = TRUE)

  weak_prior <- mixnorm(c(1, -50, 1), sigma = crohn_sigma, param = "mn")
  n_act <- 40
  ## n_pbo <- 20

  design_noprior_b <- oc2S(
    weak_prior,
    map,
    n_act,
    0,
    poc,
    sigma1 = crohn_sigma,
    sigma2 = crohn_sigma
  )

  expect_numeric(
    design_noprior_b(-20, -30),
    lower = 0,
    upper = 1,
    any.missing = FALSE
  )
})

test_that("Normal OC 2-sample avoids undefined behavior, example 1", {
  skip_on_cran()

  sigma_ref <- 3.2
  ## map_ref <- mixnorm(c(0.51, -2.1, 0.39), c(0.42, -2.1, 0.995), c(0.06, -1.99, 2.32), sigma=sigma_ref)
  ## chagned so that weights sum to 1
  map_ref <- mixnorm(
    c(0.52, -2.1, 0.39),
    c(0.42, -2.1, 0.995),
    c(0.06, -1.99, 2.32),
    sigma = sigma_ref
  )
  prior_flat <- mixnorm(c(1, 0, 100), sigma = sigma_ref)
  alpha <- 0.05
  dec <- decision2S(1 - alpha, 0, lower.tail = FALSE)
  n <- 58
  k <- 2
  design_map <- oc2S(
    prior_flat,
    map_ref,
    n,
    n / k,
    dec,
    sigma1 = sigma_ref,
    sigma2 = sigma_ref
  )
  design_map_2 <- oc2S(
    prior_flat,
    map_ref,
    n,
    n / k,
    dec,
    sigma1 = sigma_ref,
    sigma2 = sigma_ref
  )

  x <- seq(-2.6, -1.6, by = 0.1)
  expect_numeric(design_map(x, x), lower = 0, upper = 1, any.missing = FALSE)
  expect_silent(design_map(-3, -4))
  expect_numeric(design_map(-3, -4), lower = 0, upper = 1, any.missing = FALSE)
  expect_numeric(design_map(-3, 4), lower = 0, upper = 1, any.missing = FALSE)
  expect_numeric(
    design_map(-1.6, -1.6),
    lower = 0,
    upper = 1,
    any.missing = FALSE
  )

  expect_numeric(
    design_map_2(-3, -4),
    lower = 0,
    upper = 1,
    any.missing = FALSE
  )
  expect_numeric(design_map_2(-3, 4), lower = 0, upper = 1, any.missing = FALSE)
  expect_numeric(
    design_map_2(-1.6, -1.6),
    lower = 0,
    upper = 1,
    any.missing = FALSE
  )
  expect_numeric(design_map_2(x, x), lower = 0, upper = 1, any.missing = FALSE)
})