5-covid-covariates.R

rm(list=ls())
source('functions/ib_func.R')
source('functions/brease_func.R')
source('functions/brease.R')
source('functions/contour_func.R')
source('functions/helpers.R')
library(dplyr)
library(latex2exp)
library(rjags)
library(rstan)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = TRUE)

# stratified Pfizer COVID vaccine data (Polack et al., 2023)

# age groups (75+, 65-74, 56-64, 16-55)
groups <- c("75+", "65-74", "56-64", "16-55")
N0 <- c(785,3880-785,7543-3880,9955)
y0 <- c(5,19-5,48-19,114)
N1 <- c(774,3848-774,7500-3848,9897)
y1 <- c(0,1-0,3-1,5)

## Data for other analysis: remove comments to run.

## race groups (white / black / others)
# groups <- c("white", "black", "others")
# N0 <- c(14670,1486,1355)
# y0 <- c(146,7,9)
# N1 <- c(14504,1502,1405)
# y1 <- c(7,0,1)

## countries (Argentina, Brazil, USA)
# groups <- c("Argentina", "Brazil", "USA")
# N0 <- c(2521,1121,13506)
# y0 <- c(35,8,119)
# N1 <- c(2545,1129,13359)
# y1 <- c(1,1,6)


# number of strata
n_strata <- length(N0)
N <- N1 + N0

# probabilities for posterior summaries
probs <- c(0.5,0.025,0.975) # for posterior median and 95% CI estimation


# No pooling --------------------------------------------------------------

## No pooling analysis:
## Treats each strata as its own, with no sharing of information.

## hyperparameters
mu0 <- 1/2
n0 <- 2
mus <- mue <- .3
ns <- ne <- 1

## computes results for each study
nopool <- list()
for(j in 1:n_strata){
  nopool[[j]] <-  gibbs.brease(y0  = y0[j],  y1 = y1[j], 
                               N0  = N0[j],  N1 = N1[j], 
                               mu0 = mu0, n0 = n0,
                               mue = mue, ne = ne, 
                               mus = mus, ns = ns,
                               n.iter = 1e4, burn.in = 1e3)
}
ef.nopool <- do.call("rbind", lapply(nopool, function(x) quantile(1-x$theta1/x$theta0, probs)))
ef.nopool <- as.data.frame(round(ef.nopool*100, 2))
ef.nopool <- cbind(groups, N, ef.nopool)
ef.nopool


# Partial Pooling ---------------------------------------------------------

# hierarchical BREASE analysis

## hierarchical hyperparameters
lambda <- c(1/2, 1/2, 1/2)
nu     <- c(10, 10, 10)
a_beta <- lambda*nu
b_beta <-  (1-lambda)*nu
a_gamma <- rep(10, 3)
b_gamma <- rep(.1, 3)

## data for hierarchical model
bhm.data <- list(
  n_strata = n_strata,
  N0 = N0,
  N1 = N1,
  y0 = y0,
  y1 = y1,
  a_beta = a_beta,
  b_beta = b_beta,
  a_gamma = a_gamma,
  b_gamma = b_gamma
)

## fit hierarchical model in Stan.
## there should be no divergent transitions
fit <- stan(
  file = 'functions/bhm-code.stan',
  data = bhm.data,
  chains = 4, iter = 5000, warmup = 2500,
  control = list(adapt_delta = 0.95, max_treedepth = 15)
)

## extract posterior samples
sims <- rstan::extract(fit)

## computes ef for each study
ef <- 1-sims$theta1/sims$theta0
ef <- (apply(ef, 2, quantile, probs))
ef <- apply(ef*100, 2, round, 2)
ef <- as.data.frame(t(ef))
ef <- cbind(groups, N,  ef)
ef