tests/testthat/test-pos2S.R

In RBesT: R Bayesian Evidence Synthesis Tools

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

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

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

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

N1 <- 30
N2 <- 40

## we test here that the PoS indeed averages over the predictive of an
## informative prior weighting the conditional power respectively.

s <- 1

pcrit <- 0.80
qcrit <- 0

N_samp <- 1E4

dec <- decision2S(pcrit, qcrit)
decU <- decision2S(pcrit, qcrit, lower.tail = FALSE)

test_pos2S <- function(prior1, prior2, ia_dist1, ia_dist2, n1, n2, dec, decU) {
  skip_on_cran()

  ## the PoS is the expected value of the condition power integrated
  ## over the interim density which is what we check here
  cpo_analytic <- oc2S(prior1, prior2, n1, n2, dec)
  pos_analytic <- pos2S(prior1, prior2, n1, n2, dec)
  samp1 <- rmix(ia_dist1, N_samp)
  samp2 <- rmix(ia_dist2, N_samp)
  pos_mc <- mean(cpo_analytic(samp1, samp2))
  ## print(pos_mc)
  ## print(pos_analytic(ia_dist1, ia_dist2))
  expect_true(all(abs(pos_mc - pos_analytic(ia_dist1, ia_dist2)) < eps))
  lower.tail <- attr(dec, "lower.tail")
  if (lower.tail) {
    ## cat("Also testing lower.tail=FALSE\n")
    test_pos2S(prior1, prior2, ia_dist1, ia_dist2, n1, n2, decU)
  }
}

ia1 <- postmix(prior1, m = 0.2, se = s / sqrt(15))
ia2 <- postmix(prior2, m = 0, se = s / sqrt(15))

test_that("Normal PoS 2 sample function matches MC integration of CPO", {
  test_pos2S(
    prior1,
    prior2,
    ia1,
    ia2,
    N1,
    N2,
    dec,
    decU
  )
})

## also run a MC comparison
pos2S_normal_MC <- function(
  prior1,
  prior2,
  N1,
  N2,
  dtheta1,
  dtheta2,
  pcrit = 0.975,
  qcrit = 0
) {
  skip_on_cran()

  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

  pred_dtheta1 <- preddist(dtheta1, n = N1) ## , sigma=mean_sd1)
  pred_dtheta2 <- preddist(dtheta2, n = N2) ## , sigma=mean_sd1)

  ## mean_samp1 <- rnorm(Nsim, theta1, mean_sd1)
  ## mean_samp2 <- rnorm(Nsim, theta2, mean_sd2)
  mean_samp1 <- rmix(pred_dtheta1, N_samp)
  mean_samp2 <- rmix(pred_dtheta2, N_samp)

  dec <- rep(NA, N_samp)

  for (i in 1:N_samp) {
    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)
}

test_that("Normal PoS 2 sample function matches MC integration", {
  pos_mc <- pos2S_normal_MC(
    prior1,
    prior2,
    N1,
    N2,
    ia1,
    ia2,
    pcrit = 0.8,
    qcrit = 0
  )
  pos_analytic <- pos2S(prior1, prior2, N1, N2, dec)
  expect_true(all(abs(pos_mc - pos_analytic(ia1, ia2)) < eps))
})

beta_prior <- mixbeta(c(1, 1, 1))
beta_ia1 <- postmix(beta_prior, r = 20, n = 50)
beta_ia2 <- postmix(beta_prior, r = 30, n = 50)
test_that("Binomial PoS 2 sample function matches MC integration of CPO", {
  test_pos2S(
    beta_prior,
    beta_prior,
    beta_ia1,
    beta_ia2,
    N1,
    N2,
    dec,
    decU
  )
})


gamma_prior <- mixgamma(c(1, 1, 1), param = "mn")
alpha <- 0.05
dec_count <- decision2S(1 - alpha, 0, lower.tail = TRUE)
dec_countU <- decision2S(1 - alpha, 0, lower.tail = FALSE)
gamma_ia1 <- postmix(gamma_prior, m = 0.7, n = 60)
gamma_ia2 <- postmix(gamma_prior, m = 1.2, n = 60)
test_that("Poisson PoS 2 sample function matches MC integration of CPO", {
  test_pos2S(
    gamma_prior,
    gamma_prior,
    gamma_ia1,
    gamma_ia2,
    N1,
    N2,
    dec_count,
    dec_countU
  )
})

test_that("Binomial PoS 2 with IA returns results", {
  ## reported by user
  successCrit <- decision2S(c(0.9), c(0), lower.tail = FALSE)
  n0 <- 50
  n <- 100
  n_alt <- 140
  neutr_prior <- mixbeta(c(1, 1 / 3, 1 / 3))
  post_placeboIA <- postmix(neutr_prior, r = 13, n = n0)
  post_treatIA <- postmix(neutr_prior, r = 3, n = n0)
  # Criterion for PPoS at IA
  pos_final <- pos2S(post_treatIA, post_placeboIA, n - n0, n - n0, successCrit)
  pos_final_alt <- pos2S(
    post_treatIA,
    post_placeboIA,
    n_alt - n0,
    n_alt - n0,
    successCrit
  )
  # Predictive proba of success at the end
  expect_number(
    pos_final_alt(post_treatIA, post_placeboIA),
    na.ok = FALSE,
    lower = 0,
    upper = 1,
    finite = TRUE,
    null.ok = FALSE
  )
  expect_number(
    pos_final(post_treatIA, post_placeboIA),
    na.ok = FALSE,
    lower = 0,
    upper = 1,
    finite = TRUE,
    null.ok = FALSE
  )
})