Case Study 2: SARS-CoV-2 inference in The Gambia

R Code Description

This document describes the steps and components of the provided R code for creating and running serological models using simulated data. The code defines observational models, antibody kinetics models, and runs inference using MCMC methods.

Step 1: Load Required Libraries

library(devtools)
library(tidyr)
library(dplyr)
library(magrittr)
library(serojump)
#devtools::load_all()
library(patchwork)


# Using more than one core might fail on windows, 
if (.Platform$OS.type == "windows") {
  mc.cores <- 1
}else {
  mc.cores <- 2
} # use as many as available, preferably 4.}

Step 2: Read Simulated Data

This data has been generated using serosim.

wave_full <- readRDS(file = "wave2_export.RDS")
#wave_full <- readRDS(file = "./vignettes/wave2_export.RDS")

sero_data_w2 <- wave_full$sero_data
known_inf_w2 <- wave_full$known_inf
inf_prior_w2 <- wave_full$inf_prior

# Check the entries are sensible
check_sero_no_single_entries(sero_data_w2)

## No single entries in data_sero!,

check_sero_timings(sero_data_w2)

## No individuals with less than 15 days in the study!

wave_full <- readRDS(file = "wave3_export.RDS")
#wave_full <- readRDS(file = "./vignettes/wave3_export.RDS")

sero_data_w3 <- wave_full$sero_data
known_inf_w3 <- wave_full$known_inf
inf_prior_w3 <- wave_full$inf_prior

# Check the entries are sensible
check_sero_no_single_entries(sero_data_w3)

## No single entries in data_sero!,

check_sero_timings(sero_data_w3)

## No individuals with less than 15 days in the study!

#known_inf_w2 <- known_inf_w2_check %>% filter(key_exposure == "inferrable")

Step 3: Define Likelihood and Kinetics Functions

obsLogLikelihood <- function(titre_val, titre_est, pars) {
    ll <- dnorm(titre_val, titre_est, pars[1], log = TRUE)
}

noInfSerumKinetics <- function(titre_est, timeSince, pars) {
    titre_est_log <- titre_est - (pars[1] * titre_est) * (timeSince)
    titre_est_log <- max(0, titre_est_log)
    titre_est_log
}


infTuenisPower2016 <- function(titre_est, timeSince, pars) {

    y1 <- pars[1]
    t1 <- pars[2]
    r <- pars[3]
    alpha <- pars[4]

    v <- 0.001
    mu <- 1 / t1 * y1

    if (timeSince < t1) {
        titre_est_boost <- exp(mu * timeSince)
    } else {
        titre_est_boost <- exp(y1) * (1 + (r - 1) * exp(y1)^{r - 1} * v * (timeSince - t1)) ^ {-1 / (r - 1)}
    }

    titre_est_log <- titre_est + log(titre_est_boost) * max(0, 1 - titre_est * alpha)
    titre_est_log
}

Step 4: Define the Observational and Kinetics Models

# Define the biomarkers and exposure types in the model
biomarkers <- c("spike", "NCP")
exposureTypes <- c("none", "delta", "vax", "predelta")
exposureFitted <- "delta"

# Define the observational model
observationalModel <- list(
    names = c("spike", "NCP"),
    model = makeModel(
        addObservationalModel("spike", c("sigma"), obsLogLikelihood),
        addObservationalModel("NCP", c("sigma_a"), obsLogLikelihood)),
    prior = bind_rows(
        addPrior("sigma", 0.0001, 2, "exp", 1, NA),
        addPrior("sigma_a", 0.0001, 2, "exp", 1, NA)
        ) # observational model,
)

# Define the antibody kinetics model
abkineticsModel <- list(
    model = makeModel(
            addAbkineticsModel("none_s", "spike", "none",  c("wane"), noInfSerumKinetics),
            addAbkineticsModel("delta_s", "spike", "delta", c("y1_d", "t1_d", "r_d", "s"), infTuenisPower2016),
            addAbkineticsModel("vax_s", "spike", "vax", c("y1_vax", "t1_vax", "r_vax", "s"), infTuenisPower2016),
            addAbkineticsModel("predelta_s", "spike", "predelta", c("y1_pd", "t1_pd", "r_pd", "s"), infTuenisPower2016),
            addAbkineticsModel("none_ncp", "NCP", "none",  c("wane_a"), noInfSerumKinetics),
            addAbkineticsModel("delta_ncp", "NCP", "delta", c("y1_d_a", "t1_d_a", "r_d_a", "s_a"), infTuenisPower2016),
            addAbkineticsModel("vax_ncp", "NCP", "vax", c("y1_vax_a", "t1_vax_a", "r_vax_a", "s_a"), infTuenisPower2016),
            addAbkineticsModel("predelta_ncp", "NCP", "predelta", c("y1_pd_a", "t1_pd_a", "r_pd_a", "s_a"), infTuenisPower2016)
        ),
    prior = bind_rows(
        addPrior("y1_d", 0, 6, "unif",  0, 6), # ab kinetics
        addPrior("t1_d", 7, 21,  "unif",  7, 21), # ab kinetics
        addPrior("r_d", 1, 5, "unif", 1, 5), # ab kinetics 
        addPrior("y1_vax", 0, 6, "unif",  0, 6), # ab kinetics
        addPrior("t1_vax", 7, 21,  "unif",  7, 21), # ab kinetics
        addPrior("r_vax", 1, 5, "unif", 1, 5), # ab kinetics 
        addPrior("y1_pd", 0, 6, "unif",  0, 6), # ab kinetics
        addPrior("t1_pd", 7, 21,  "unif",  7, 21), # ab kinetics
        addPrior("r_pd", 1, 5, "unif", 1, 5), # ab kinetics 
        addPrior("s",  0, 1, "unif", 0, 1), # ab kinetics 
        addPrior("wane", 0.0, 0.01, "unif", 0.0, 0.01), # observational model
        addPrior("y1_d_a", 0, 6, "unif",  0, 6), # ab kinetics
        addPrior("t1_d_a", 7, 21,  "unif",  7, 21), # ab kinetics
        addPrior("r_d_a", 1, 5, "unif", 1, 5), # ab kinetics 
        addPrior("y1_vax_a", 0, 6, "unif",  0, 6), # ab kinetics
        addPrior("t1_vax_a", 7, 21, "unif",  7, 21), # ab kinetics
        addPrior("r_vax_a", 1, 5, "unif", 1, 5), # ab kinetics 
        addPrior("y1_pd_a", 0, 6, "unif",  0, 6), # ab kinetics
        addPrior("t1_pd_a", 7, 21,  "unif", 7, 21), # ab kinetics
        addPrior("r_pd_a", 1, 5, "unif", 1, 5), # ab kinetics 
        addPrior("s_a", 0, 1, "unif", 0, 1), # ab kinetics 
        addPrior("wane_a", 0.0, 0.01, "unif", 0.0, 0.01) # observational model
    )
)



model_w2 <- createSeroJumpModel(
    data_sero = sero_data_w2, 
    data_known = known_inf_w2, 
    biomarkers = biomarkers,
    exposureTypes = exposureTypes,
    exposureFitted = exposureFitted,
    observationalModel = observationalModel,
    abkineticsModel = abkineticsModel,
    exposurePriorTime = inf_prior_w2,
    exposurePriorTimeType = "empirical"
)

## OUTLINE OF INPUTTED MODEL
## There are  2  measured biomarkers:  spike, NCP 
## There are  4  exposure types in the study period:  none, delta, vax, predelta 
## The fitted exposure type is  delta 
## PRIOR DISTRIBUTIONS 
## Prior parameters of observationalModel are:  sigma, sigma_a 
## Prior parameters of abkineticsModel are:  y1_d, t1_d, r_d, y1_vax, t1_vax, r_vax, y1_pd, t1_pd, r_pd, s, wane, y1_d_a, t1_d_a, r_d_a, y1_vax_a, t1_vax_a, r_vax_a, y1_pd_a, t1_pd_a, r_pd_a, s_a, wane_a 
## No single entries in data_sero!, 
## No individuals with less than 15 days in the study!

## Warning in addExposurePrior_checkempirical(exp_prior, T_max): Warning: The number of rows in the exposure prior dataframe does not match the number of time points in the study
## Warning in addExposurePrior_checkempirical(exp_prior, T_max): Warning: The number of rows in the exposure prior dataframe does not match the number of time points in the study

## There are individuals which have exposure before their first bleed or before their last bleed, ids:  19, 23, 45, 92, 95, 147, 162, 205, 218, 6, 18, 24, 30, 63, 98, 142, 181, 192, 199, 208, 210, 214, 215, 245, 246, 254, 265 . It is hard to perform inference on these people, so we recommend you remove the exposures from the known inf, but, if included, the model ignores these exposures.

Step 4B: Prerun sanity plots

Before running the whole model it is good to check the data and the priors. This can be done using a suit of functions plotPriors function.

p1 <- plotSero(model_w2)
p2 <- plotPriorPredictive(model_w2)
p3 <- plotPriorInfection(model_w2)
p1 / p2

p3

Step 5: Define the Complete Model

# Define the biomarkers and exposure types in the model
biomarkers <- c("spike", "NCP")
exposureTypes <- c("none", "omicron", "vax")
exposureFitted <- "omicron"

# Define the observational model
observationalModel <- list(
    names = c("spike", "NCP"),
    model = makeModel(
        addObservationalModel("spike", c("sigma"), obsLogLikelihood),
        addObservationalModel("NCP", c("sigma_a"), obsLogLikelihood)),

    prior = bind_rows(
        addPrior("sigma", 0.0001, 2, "exp", 1, NA),
        addPrior("sigma_a", 0.0001, 2, "exp", 1, NA)
        ) # observational model,
)

# Define the antibody kinetics model
abkineticsModel <- list(
    model = makeModel(
            addAbkineticsModel("none_s", "spike", "none",  c("wane"), noInfSerumKinetics),
            addAbkineticsModel("omicron_s", "spike", "omicron", c("y1_o", "t1_o", "r_o", "s"), infTuenisPower2016),
            addAbkineticsModel("vax_s", "spike", "vax", c("y1_vax", "t1_vax", "r_vax", "s"), infTuenisPower2016),
            addAbkineticsModel("none_ncp", "NCP", "none",  c("wane_a"), noInfSerumKinetics),
            addAbkineticsModel("omicron_ncp", "NCP", "omicron", c("y1_o_a", "t1_o_a", "r_o_a", "s_a"), infTuenisPower2016),
            addAbkineticsModel("vax_ncp", "NCP", "vax", c("y1_vax_a", "t1_vax_a", "r_vax_a", "s_a"), infTuenisPower2016)
        ),
    prior = bind_rows(
        addPrior("y1_o", 0, 6, "unif",  0, 6), # ab kinetics
        addPrior("t1_o", 7, 21,"unif",  7, 21), # ab kinetics
        addPrior("r_o", 1, 5, "unif", 1, 5), # ab kinetics 
        addPrior("y1_vax", 0, 6, "unif",  0, 6), # ab kinetics
        addPrior("t1_vax", 7, 21, "unif",  7, 21), # ab kinetics
        addPrior("r_vax", 1, 5, "unif", 1, 5), # ab kinetics 
        addPrior("s",  0, 1, "unif", 0, 1), # ab kinetics 
        addPrior("wane", 0.0, 0.01, "unif", 0.0, 0.01), # observational model
        addPrior("y1_o_a", 0, 6, "unif",  0, 6), # ab kinetics
        addPrior("t1_o_a", 7, 21, "unif",  7, 21), # ab kinetics
        addPrior("r_o_a", 1, 5, "unif", 1, 5), # ab kinetics 
        addPrior("y1_vax_a", 0, 6, "unif",  0, 6), # ab kinetics
        addPrior("t1_vax_a", 7, 21, "unif",  7, 21), # ab kinetics
        addPrior("r_vax_a", 1, 5, "unif", 1, 5), # ab kinetics 
        addPrior("s_a", 0, 1, "unif", 0, 1), # ab kinetics 
        addPrior("wane_a", 0.0, 0.01, "unif", 0.0, 0.01) # observational model
    )
)


model_w3 <- createSeroJumpModel(
    data_sero = sero_data_w3, 
    data_known = known_inf_w3, 
    biomarkers = biomarkers,
    exposureTypes = exposureTypes,
    exposureFitted = exposureFitted,
    observationalModel = observationalModel,
    abkineticsModel = abkineticsModel,
    exposurePriorTime = inf_prior_w3,
    exposurePriorTimeType = "empirical"
)

## OUTLINE OF INPUTTED MODEL
## There are  2  measured biomarkers:  spike, NCP 
## There are  3  exposure types in the study period:  none, omicron, vax 
## The fitted exposure type is  omicron 
## PRIOR DISTRIBUTIONS 
## Prior parameters of observationalModel are:  sigma, sigma_a 
## Prior parameters of abkineticsModel are:  y1_o, t1_o, r_o, y1_vax, t1_vax, r_vax, s, wane, y1_o_a, t1_o_a, r_o_a, y1_vax_a, t1_vax_a, r_vax_a, s_a, wane_a 
## No single entries in data_sero!, 
## No individuals with less than 15 days in the study!

## Warning in addExposurePrior_checkempirical(exp_prior, T_max): Warning: The number of rows in the exposure prior dataframe does not match the number of time points in the study
## Warning in addExposurePrior_checkempirical(exp_prior, T_max): Warning: The number of rows in the exposure prior dataframe does not match the number of time points in the study

## There are individuals which have exposure before their first bleed or before their last bleed, ids:  111, 57, 227, 228, 248 . It is hard to perform inference on these people, so we recommend you remove the exposures from the known inf, but, if included, the model ignores these exposures.

Step 5B: Prerun sanity plots