tests/testthat/test-tte.R

In beastt: Bayesian Evaluation, Analysis, and Simulation Software Tools for Trials

################################################################################
# Code to test various functions of beastt - TTE endpoint
################################################################################
skip_on_cran()
### Source R script with Stan code and compile Stan models
source("code_for_tte_tests_Stan.R")
# compile Stan model - power prior
stan_mod_pp <- rstan::stan_model(model_code = BDB_stan_pp)
# compile Stan model - posterior
stan_mod_post <- rstan::stan_model(model_code = BDB_stan_post)

################################################################################
# calc_power_prior_weibull
################################################################################

##### Test for calc_power_prior_weibull() with:
#####    (1) prop_scr_obj (not unweighted external data)
#####    (2) Laplace approximation
test_that("calc_power_prior_weibull returns correct values for first case", {
  ## Define values to be used for both beastt code and comparison code
  init_norm_prior <- dist_normal(0, 10)   # initial prior for intercept parameter
  init_hn_scale <- 50            # scale of the half-normal initial prior for the shape parameter
  init_prior <- dist_normal(mu = 0.5, sigma = 10)   # proper initial prior for power prior
  ps_obj <- calc_prop_scr(internal_df = filter(int_tte_df, trt == 0),
                          external_df = ex_tte_df,
                          id_col = subjid,
                          model = ~ cov1 + cov2 + cov3 + cov4)

  ## beastt code
  pwr_prior <- calc_power_prior_weibull(external_data = ps_obj,
                                        response = y,
                                        event = event,
                                        intercept = init_norm_prior,
                                        shape = init_hn_scale,
                                        approximation = "Laplace")
  pp_mean_beastt <- parameters(pwr_prior)$mu[[1]]
  pp_cov_beastt <- parameters(pwr_prior)$sigma[[1]]

  ## Comparison code
  ipws <- ps_obj$external_df$`___weight___`
  calc_neg_log_dens <- function(x, y_vec, event_vec, ipw_vec, beta0_mean, beta0_sd, alpha_scale){

    # Extract elements of x and save as parameters log(alpha) and beta0
    log_alpha <- x[1]
    beta0 <- x[2]

    # Log density
    log_dens <-
      sum( ipw_vec * event_vec * dweibull(x = y_vec, shape = exp(log_alpha), scale = exp(-beta0), log = TRUE) -
             ipw_vec * (1 - event_vec) * (y_vec / exp(-beta0))^exp(log_alpha) ) +   # log IWP likelihood of external data
      dnorm(x = beta0, mean = beta0_mean, sd = beta0_sd, log = TRUE) +       # log density of normal prior on beta0
      .5 * log(2/pi) - log(alpha_scale) - exp(2 * log_alpha) / (2 * alpha_scale^2)  # log dens of half norm prior on alpha
    -log_dens

  }
  inits <- c(0, 0)   # initial parameters for log(alpha) and beta0, respectively
  optim_pp_ctrl <- optim(par = inits,
                         fn = calc_neg_log_dens,
                         method = "Nelder-Mead",
                         y_vec = ex_tte_df$y,
                         event_vec = ex_tte_df$event,
                         ipw_vec = ipws,
                         beta0_mean = parameters(init_norm_prior)$mu,
                         beta0_sd = parameters(init_norm_prior)$sigma,
                         alpha_scale = init_hn_scale,
                         hessian = TRUE)
  pp_mean_comp <- optim_pp_ctrl$par              # mean vector
  pp_cov_comp <- solve(optim_pp_ctrl$hessian)    # covariance matrix

  ## Check that means and SDs of the normal power priors are equal using both methods
  expect_equal(pp_mean_beastt, pp_mean_comp)
  expect_equal(pp_cov_beastt, pp_cov_comp)

  ## Check that a distribution is returned
  expect_s3_class(pwr_prior, "distribution")
})

##### Test for calc_power_prior_weibull() with:
#####    (1) prop_scr_obj (not unweighted external data)
#####    (2) MCMC approximation
test_that("calc_power_prior_weibull returns correct values for second case", {
  ## Define values to be used for both beastt code and comparison code
  init_norm_prior <- dist_normal(0, 10)   # initial prior for intercept parameter
  init_hn_scale <- 50            # scale of the half-normal initial prior for the shape parameter
  init_prior <- dist_normal(mu = 0.5, sigma = 10)   # proper initial prior for power prior
  ps_obj <- calc_prop_scr(internal_df = filter(int_tte_df, trt == 0),
                          external_df = ex_tte_df,
                          id_col = subjid,
                          model = ~ cov1 + cov2 + cov3 + cov4)

  ## beastt code
  set.seed(123)
  pwr_prior <- calc_power_prior_weibull(external_data = ps_obj,
                                        response = y,
                                        event = event,
                                        intercept = init_norm_prior,
                                        shape = init_hn_scale,
                                        approximation = "MCMC")
  pp_mean_beastt <- parameters(pwr_prior)$mu[[1]]
  pp_cov_beastt <- parameters(pwr_prior)$sigma[[1]]

  ## Comparison code
  ipws <- ps_obj$external_df$`___weight___`
  data_input <- list(
    N = nrow(ex_tte_df),     # sample size of the external control arm
    y = ex_tte_df$y,         # N x 1 vector of observed times
    e = ex_tte_df$event,     # N x 1 vector of event indicators (1: event, 0: no event)
    wgt = ipws,              # N x 1 vector of IPWs
    beta_mean = parameters(init_norm_prior)$mu,    # mean for normal prior on the regression parameters
    beta_sd = parameters(init_norm_prior)$sigma,   # sd for normal prior on the regression parameters
    shape_mean = 0,                                # mean of half normal prior on weibull shape parameter
    shape_sd = init_hn_scale                       # sd of half normal prior on weibull shape parameter
  )
  set.seed(123)
  stan_fit <- sampling(stan_mod_pp, data = data_input, pars = c("beta", "shape"),
                       iter = 26000, warmup = 1000, chains = 4)
  pp_draws <- as.matrix(stan_fit)           # posterior samples of each parameter
  pp_mean_comp <- c(mean(log(pp_draws[,"shape"])), mean(pp_draws[,"beta"]))
  pp_cov_comp <- cov(cbind(log(pp_draws[,"shape"]), pp_draws[,"beta"]))
  names(pp_mean_comp) <- c("log_alpha", "beta0")
  dimnames(pp_cov_comp)[[1]] <- dimnames(pp_cov_comp)[[2]] <- c("log_alpha", "beta0")

  ## Check that means and SDs of the normal power priors are equal using both methods
  expect_equal(all(abs(pp_mean_beastt-pp_mean_comp) < 0.0025), TRUE)
  expect_equal(all(abs(pp_cov_beastt-pp_cov_comp) < 0.001), TRUE)

  ## Check that a distribution is returned
  expect_s3_class(pwr_prior, "distribution")
})

##### Test for calc_power_prior_weibull() with:
#####    (1) unweighted external data (read in external df, not prop_scr_obj)
#####    (2) Laplace approximation
test_that("calc_power_prior_weibull returns correct values for third case", {
  ## Define values to be used for both beastt code and comparison code
  init_norm_prior <- dist_normal(0, 10)   # initial prior for intercept parameter
  init_hn_scale <- 50            # scale of the half-normal initial prior for the shape parameter
  init_prior <- dist_normal(mu = 0.5, sigma = 10)   # proper initial prior for power prior

  ## beastt code
  pwr_prior <- calc_power_prior_weibull(external_data = ex_tte_df,
                                        response = y,
                                        event = event,
                                        intercept = init_norm_prior,
                                        shape = init_hn_scale,
                                        approximation = "Laplace")
  pp_mean_beastt <- parameters(pwr_prior)$mu[[1]]
  pp_cov_beastt <- parameters(pwr_prior)$sigma[[1]]

  ## Comparison code
  calc_neg_log_dens <- function(x, y_vec, event_vec, ipw_vec, beta0_mean, beta0_sd, alpha_scale){

    # Extract elements of x and save as parameters log(alpha) and beta0
    log_alpha <- x[1]
    beta0 <- x[2]

    # Log density
    log_dens <-
      sum( ipw_vec * event_vec * dweibull(x = y_vec, shape = exp(log_alpha), scale = exp(-beta0), log = TRUE) -
             ipw_vec * (1 - event_vec) * (y_vec / exp(-beta0))^exp(log_alpha) ) +   # log IWP likelihood of external data
      dnorm(x = beta0, mean = beta0_mean, sd = beta0_sd, log = TRUE) +       # log density of normal prior on beta0
      .5 * log(2/pi) - log(alpha_scale) - exp(2 * log_alpha) / (2 * alpha_scale^2)  # log dens of half norm prior on alpha
    -log_dens

  }
  inits <- c(0, 0)   # initial parameters for log(alpha) and beta0, respectively
  optim_pp_ctrl <- optim(par = inits,
                         fn = calc_neg_log_dens,
                         method = "Nelder-Mead",
                         y_vec = ex_tte_df$y,
                         event_vec = ex_tte_df$event,
                         ipw_vec = rep(1, nrow(ex_tte_df)),
                         beta0_mean = parameters(init_norm_prior)$mu,
                         beta0_sd = parameters(init_norm_prior)$sigma,
                         alpha_scale = init_hn_scale,
                         hessian = TRUE)
  pp_mean_comp <- optim_pp_ctrl$par              # mean vector
  pp_cov_comp <- solve(optim_pp_ctrl$hessian)    # covariance matrix

  ## Check that means and SDs of the normal power priors are equal using both methods
  expect_equal(pp_mean_beastt, pp_mean_comp)
  expect_equal(pp_cov_beastt, pp_cov_comp)

  ## Check that a distribution is returned
  expect_s3_class(pwr_prior, "distribution")
})

##### Test for calc_power_prior_weibull() with:
#####    (1) unweighted external data (read in external df, not prop_scr_obj)
#####    (2) MCMC approximation
test_that("calc_power_prior_weibull returns correct values for fourth case", {
  ## Define values to be used for both beastt code and comparison code
  init_norm_prior <- dist_normal(0, 10)   # initial prior for intercept parameter
  init_hn_scale <- 50            # scale of the half-normal initial prior for the shape parameter
  init_prior <- dist_normal(mu = 0.5, sigma = 10)   # proper initial prior for power prior

  ## beastt code
  set.seed(123)
  pwr_prior <- calc_power_prior_weibull(external_data = ex_tte_df,
                                        response = y,
                                        event = event,
                                        intercept = init_norm_prior,
                                        shape = init_hn_scale,
                                        approximation = "MCMC")
  pp_mean_beastt <- parameters(pwr_prior)$mu[[1]]
  pp_cov_beastt <- parameters(pwr_prior)$sigma[[1]]

  ## Comparison code
  data_input <- list(
    N = nrow(ex_tte_df),    # sample size of the external control arm
    y = ex_tte_df$y,         # N x 1 vector of observed times
    e = ex_tte_df$event,     # N x 1 vector of event indicators (1: event, 0: no event)
    wgt = rep(1, nrow(ex_tte_df)),                 # N x 1 vector of 1s
    beta_mean = parameters(init_norm_prior)$mu,    # mean for normal prior on the regression parameters
    beta_sd = parameters(init_norm_prior)$sigma,   # sd for normal prior on the regression parameters
    shape_mean = 0,                                # mean of half normal prior on weibull shape parameter
    shape_sd = init_hn_scale                       # sd of half normal prior on weibull shape parameter
  )
  set.seed(123)
  stan_fit <- sampling(stan_mod_pp, data = data_input, pars = c("beta", "shape"),
                       iter = 26000, warmup = 1000, chains = 4)
  pp_draws <- as.matrix(stan_fit)           # posterior samples of each parameter
  pp_mean_comp <- c(mean(log(pp_draws[,"shape"])), mean(pp_draws[,"beta"]))
  pp_cov_comp <- cov(cbind(log(pp_draws[,"shape"]), pp_draws[,"beta"]))
  names(pp_mean_comp) <- c("log_alpha", "beta0")
  dimnames(pp_cov_comp)[[1]] <- dimnames(pp_cov_comp)[[2]] <- c("log_alpha", "beta0")

  ## Check that means and SDs of the normal power priors are equal using both methods
  expect_equal(all(abs(pp_mean_beastt-pp_mean_comp) < 0.002), TRUE)
  expect_equal(all(abs(pp_cov_beastt-pp_cov_comp) < 0.001), TRUE)

  ## Check that a distribution is returned
  expect_s3_class(pwr_prior, "distribution")
})

# Test for invalid external data
test_that("calc_power_prior_weibull handles invalid external data", {
  expect_error(calc_power_prior_weibull(external_data = c(1, 2, 3),
                                        response = y,
                                        event = event,
                                        intercept = dist_normal(0, 10),
                                        shape = 50,
                                        approximation = "Laplace"))
  expect_error(calc_power_prior_weibull(external_data = "abc",
                                        response = y,
                                        event = event,
                                        intercept = dist_normal(0, 10),
                                        shape = 50,
                                        approximation = "Laplace"))
})

# Test for invalid intercept
test_that("calc_power_prior_weibull handles invalid intercept", {
  expect_error(calc_power_prior_weibull(external_data = ex_tte_df,
                                        response = y,
                                        event = event,
                                        intercept = "a",
                                        shape = 50,
                                        approximation = "Laplace"))
  expect_error(calc_power_prior_weibull(external_data = ex_tte_df,
                                        response = y,
                                        event = event,
                                        intercept = c(1, 2, 3),
                                        shape = 50,
                                        approximation = "Laplace"))
})

# Test for invalid shape
test_that("calc_power_prior_weibull handles invalid shape", {
  expect_error(calc_power_prior_weibull(external_data = ex_tte_df,
                                        response = y,
                                        event = event,
                                        intercept = dist_normal(0, 10),
                                        shape = "a",
                                        approximation = "Laplace"))
  expect_error(calc_power_prior_weibull(external_data = ex_tte_df,
                                        response = y,
                                        event = event,
                                        intercept = dist_normal(0, 10),
                                        shape = c(1, 2, 3),
                                        approximation = "Laplace"))
})

# Test for invalid approximation method
test_that("calc_power_prior_weibull handles invalid approximation method", {
  expect_error(calc_power_prior_weibull(external_data = ex_tte_df,
                                        response = y,
                                        event = event,
                                        intercept = dist_normal(0, 10),
                                        shape = 50,
                                        approximation = 12))
  expect_error(calc_power_prior_weibull(external_data = ex_tte_df,
                                        response = y,
                                        event = event,
                                        intercept = dist_normal(0, 10),
                                        shape = 50,
                                        approximation = c(5, 10)))
  expect_error(calc_power_prior_weibull(external_data = ex_tte_df,
                                        response = y,
                                        event = event,
                                        intercept = dist_normal(0, 10),
                                        shape = 50,
                                        approximation = "zero"))
})

# Test for internal and external data with different response variable names
test_that("calc_power_prior_weibull handles different response variable names", {
  int <- int_tte_df
  ex <- ex_tte_df |>
    dplyr::rename(y2 = y)
  init_prior <- dist_normal(mu = 0, sigma = 10)
  init_hn_scale <- 50


  ps_obj <- calc_prop_scr(internal_df = filter(int, trt == 0),
                          external_df = ex,
                          id_col = subjid,
                          model = ~ cov1 + cov2 + cov3 + cov4)

  expect_error(calc_power_prior_weibull(external_data = ps_obj,
                                        response = y,
                                        event = event,
                                        intercept = init_norm_prior,
                                        shape = init_hn_scale,
                                        approximation = "MCMC"))
})

################################################################################
# calc_post_weibull
################################################################################

##### Test for calc_post_weibull() with:
#####    (1) MVN prior (not a mixture of two MVN priors)
#####    (2) one analysis time
test_that("calc_post_weibull returns correct values for first case", {
  ## Define values to be used for both beastt code and comparison code
  MVN_prior <- dist_multivariate_normal(mu = list(c(0.5, 0.6)),
                                        sigma = list(matrix(c(10, -.5, -.5, 10), nrow = 2)))   # normal prior
  analysis_time <- 12

  ## beastt code
  set.seed(123)
  post_dist <- calc_post_weibull(internal_data = filter(int_tte_df, trt == 0),
                                 response = y,
                                 event = event,
                                 prior = MVN_prior,
                                 analysis_time = analysis_time)
  surv_prob_beastt <- as.data.frame(extract(post_dist, pars = c("survProb")))[,1]
  post_mean_beastt <- mean(surv_prob_beastt)
  post_var_beastt <- var(surv_prob_beastt)

  ## Comparison code
  analysis_time_comp <- c(analysis_time, analysis_time)  # making dimension >1 so code will run
  data_input <- list(
    N = nrow(filter(int_tte_df, trt == 0)),     # sample size of the internal control arm
    y = filter(int_tte_df, trt == 0)$y,         # N x 1 vector of observed times
    e = filter(int_tte_df, trt == 0)$event,     # N x 1 vector of event indicators (1: event, 0: no event)
    T = length(analysis_time_comp),             # length of analysis times
    times = c(analysis_time_comp),              # analysis time
    cov_inf = parameters(MVN_prior)$sigma[[1]],      # informative prior component covariance matrix
    cov_vague = parameters(MVN_prior)$sigma[[1]],    # vague prior component covariance matrix
    mu_inf = parameters(MVN_prior)$mu[[1]],          # informative prior component mean vector
    mu_vague = parameters(MVN_prior)$mu[[1]],        # vague prior component mean vector
    w0 = 1                                           # prior mixture weight
  )
  set.seed(123)
  stan_fit <- sampling(stan_mod_post, data = data_input, pars = "survProb",
                       iter = 26000, warmup = 1000, chains = 4)
  post_draws <- as.matrix(stan_fit)           # posterior samples of each parameter
  post_mean_comp <- mean(post_draws[,"survProb[1]"])
  post_var_comp <- var(post_draws[,"survProb[1]"])

  ## Check that means and SDs of the normal power priors are equal using both methods
  expect_equal(abs(post_mean_beastt-post_mean_comp) < 0.001, TRUE)
  expect_equal(abs(post_var_beastt-post_var_comp) < 0.0001, TRUE)

  ## Check that a stanfit is returned
  expect_s4_class(post_dist, "stanfit")
})

##### Test for calc_post_weibull() with:
#####    (1) MVN prior (not a mixture of two MVN priors)
#####    (2) two analysis times
test_that("calc_post_weibull returns correct values for second case", {
  ## Define values to be used for both beastt code and comparison code
  MVN_prior <- dist_multivariate_normal(mu = list(c(0.5, 0.6)),
                                        sigma = list(matrix(c(10, -.5, -.5, 10), nrow = 2)))   # normal prior
  analysis_times <- c(6, 8)

  ## beastt code
  set.seed(123)
  post_dist <- calc_post_weibull(internal_data = filter(int_tte_df, trt == 0),
                                 response = y,
                                 event = event,
                                 prior = MVN_prior,
                                 analysis_time = analysis_times)
  surv_prob_beastt <- as.data.frame(extract(post_dist, pars = c("survProb")))
  post_means_beastt <- colMeans(surv_prob_beastt)
  post_vars_beastt <- apply(surv_prob_beastt, 2, var)

  ## Comparison code
  data_input <- list(
    N = nrow(filter(int_tte_df, trt == 0)),     # sample size of the internal control arm
    y = filter(int_tte_df, trt == 0)$y,         # N x 1 vector of observed times
    e = filter(int_tte_df, trt == 0)$event,     # N x 1 vector of event indicators (1: event, 0: no event)
    T = length(analysis_times),                 # length of analysis times
    times = analysis_times,                     # analysis times
    cov_inf = parameters(MVN_prior)$sigma[[1]],      # informative prior component covariance matrix
    cov_vague = parameters(MVN_prior)$sigma[[1]],    # vague prior component covariance matrix
    mu_inf = parameters(MVN_prior)$mu[[1]],          # informative prior component mean vector
    mu_vague = parameters(MVN_prior)$mu[[1]],        # vague prior component mean vector
    w0 = 1                                           # prior mixture weight
  )
  set.seed(123)
  stan_fit <- sampling(stan_mod_post, data = data_input, pars = "survProb",
                       iter = 26000, warmup = 1000, chains = 4)
  post_draws <- as.matrix(stan_fit)           # posterior samples of each parameter
  post_means_comp <- colMeans(post_draws[,c("survProb[1]", "survProb[2]")])
  post_vars_comp <- apply(post_draws[,c("survProb[1]", "survProb[2]")], 2, var)

  ## Check that means and SDs of the normal power priors are equal using both methods
  expect_equal(all(abs(post_means_beastt-post_means_comp) < 0.001), TRUE)
  expect_equal(all(abs(post_vars_beastt-post_vars_comp) < 0.0001), TRUE)

  ## Check that a stanfit is returned
  expect_s4_class(post_dist, "stanfit")
})

##### Test for calc_post_weibull() with:
#####    (1) mixture prior (two MVN prior components)
#####    (2) one analysis time
test_that("calc_post_weibull returns correct values for third case", {
  ## Define values to be used for both beastt code and comparison code
  MVN_prior <- dist_multivariate_normal(mu = list(c(0.5, 0.6)),
                                        sigma = list(matrix(c(10, -.5, -.5, 10), nrow = 2)))   # normal prior
  mix_prior <- robustify_mvnorm(MVN_prior, n = sum(ex_tte_df$event), weights = c(.5, .5))
  analysis_time <- 12

  ## beastt code
  set.seed(123)
  post_dist <- calc_post_weibull(internal_data = filter(int_tte_df, trt == 0),
                                 response = y,
                                 event = event,
                                 prior = mix_prior,
                                 analysis_time = analysis_time)
  surv_prob_beastt <- as.data.frame(extract(post_dist, pars = c("survProb")))[,1]
  post_mean_beastt <- mean(surv_prob_beastt)
  post_var_beastt <- var(surv_prob_beastt)

  ## Comparison code
  analysis_time_comp <- c(analysis_time, analysis_time)  # making dimension >1 so code will run
  data_input <- list(
    N = nrow(filter(int_tte_df, trt == 0)),     # sample size of the internal control arm
    y = filter(int_tte_df, trt == 0)$y,         # N x 1 vector of observed times
    e = filter(int_tte_df, trt == 0)$event,     # N x 1 vector of event indicators (1: event, 0: no event)
    T = length(analysis_time_comp),             # length of analysis times
    times = c(analysis_time_comp),              # analysis time
    cov_inf = mix_sigmas(mix_prior)$informative,    # informative prior component covariance matrix
    cov_vague = mix_sigmas(mix_prior)$vague,        # vague prior component covariance matrix
    mu_inf = mix_means(mix_prior)$informative,      # informative prior component mean vector
    mu_vague = mix_means(mix_prior)$vague,          # vague prior component mean vector
    w0 = .5                                         # prior mixture weight
  )
  set.seed(123)
  stan_fit <- sampling(stan_mod_post, data = data_input, pars = "survProb",
                       iter = 26000, warmup = 1000, chains = 4)
  post_draws <- as.matrix(stan_fit)           # posterior samples of each parameter
  post_mean_comp <- mean(post_draws[,"survProb[1]"])
  post_var_comp <- var(post_draws[,"survProb[1]"])

  ## Check that means and SDs of the normal power priors are equal using both methods
  expect_equal(abs(post_mean_beastt-post_mean_comp) < 0.001, TRUE)
  expect_equal(abs(post_var_beastt-post_var_comp) < 0.0001, TRUE)

  ## Check that a stanfit is returned
  expect_s4_class(post_dist, "stanfit")
})

##### Test for calc_post_weibull() with:
#####    (1) mixture prior (two MVN prior components)
#####    (2) two analysis times
test_that("calc_post_weibull returns correct values for fourth case", {
  ## Define values to be used for both beastt code and comparison code
  MVN_prior <- dist_multivariate_normal(mu = list(c(0.5, 0.6)),
                                        sigma = list(matrix(c(10, -.5, -.5, 10), nrow = 2)))   # normal prior
  mix_prior <- robustify_mvnorm(MVN_prior, n = sum(ex_tte_df$event), weights = c(.5, .5))
  analysis_times <- c(6, 8)

  ## beastt code
  set.seed(123)
  post_dist <- calc_post_weibull(internal_data = filter(int_tte_df, trt == 0),
                                 response = y,
                                 event = event,
                                 prior = MVN_prior,
                                 analysis_time = analysis_times)
  surv_prob_beastt <- as.data.frame(extract(post_dist, pars = c("survProb")))
  post_means_beastt <- colMeans(surv_prob_beastt)
  post_vars_beastt <- apply(surv_prob_beastt, 2, var)

  ## Comparison code
  data_input <- list(
    N = nrow(filter(int_tte_df, trt == 0)),     # sample size of the internal control arm
    y = filter(int_tte_df, trt == 0)$y,         # N x 1 vector of observed times
    e = filter(int_tte_df, trt == 0)$event,     # N x 1 vector of event indicators (1: event, 0: no event)
    T = length(analysis_times),                 # length of analysis times
    times = analysis_times,                     # analysis times
    cov_inf = mix_sigmas(mix_prior)$informative,    # informative prior component covariance matrix
    cov_vague = mix_sigmas(mix_prior)$vague,        # vague prior component covariance matrix
    mu_inf = mix_means(mix_prior)$informative,      # informative prior component mean vector
    mu_vague = mix_means(mix_prior)$vague,          # vague prior component mean vector
    w0 = .5                                         # prior mixture weight
  )
  set.seed(123)
  stan_fit <- sampling(stan_mod_post, data = data_input, pars = "survProb",
                       iter = 26000, warmup = 1000, chains = 4)
  post_draws <- as.matrix(stan_fit)           # posterior samples of each parameter
  post_means_comp <- colMeans(post_draws[,c("survProb[1]", "survProb[2]")])
  post_vars_comp <- apply(post_draws[,c("survProb[1]", "survProb[2]")], 2, var)

  ## Check that means and SDs of the normal power priors are equal using both methods
  expect_equal(all(abs(post_means_beastt-post_means_comp) < 0.001), TRUE)
  expect_equal(all(abs(post_vars_beastt-post_vars_comp) < 0.0001), TRUE)

  ## Check that a stanfit is returned
  expect_s4_class(post_dist, "stanfit")
})

# Test for invalid internal data
test_that("calc_post_weibull handles invalid internal data", {
  expect_error(calc_post_weibull(internal_data = c(1, 2, 3),
                                 response = y,
                                 event = event,
                                 prior = dist_multivariate_normal(mu = list(c(0.5, 0.6)),
                                                                  sigma = list(matrix(c(10, -.5, -.5, 10), nrow = 2))),
                                 analysis_time = c(12, 24)))
  expect_error(calc_post_weibull(internal_data = "a",
                                 response = y,
                                 event = event,
                                 prior = dist_multivariate_normal(mu = list(c(0.5, 0.6)),
                                                                  sigma = list(matrix(c(10, -.5, -.5, 10), nrow = 2))),
                                 analysis_time = c(12, 24)))
})

# Test for invalid intercept
test_that("calc_post_weibull handles invalid prior", {
  expect_error(calc_post_weibull(internal_data = filter(int_tte_df, trt == 0),
                                 response = y,
                                 event = event,
                                 prior = c(1, 2, 3),
                                 analysis_time = c(12, 24)))
  expect_error(calc_post_weibull(internal_data = filter(int_tte_df, trt == 0),
                                 response = y,
                                 event = event,
                                 prior = "a",
                                 analysis_time = c(12, 24)))
})

# Test for invalid analysis time
test_that("calc_post_weibull handles invalid analysis_time", {
  expect_error(calc_post_weibull(internal_data = filter(int_tte_df, trt == 0),
                                 response = y,
                                 event = event,
                                 prior = dist_multivariate_normal(mu = list(c(0.5, 0.6)),
                                                                  sigma = list(matrix(c(10, -.5, -.5, 10), nrow = 2))),
                                 analysis_time = "a"))
  expect_error(calc_post_weibull(internal_data = filter(int_tte_df, trt == 0),
                                 response = y,
                                 event = event,
                                 prior = dist_multivariate_normal(mu = list(c(0.5, 0.6)),
                                                                  sigma = list(matrix(c(10, -.5, -.5, 10), nrow = 2))),
                                 analysis_time = dist_normal(0, 1)))
})