Final size vs serology

main_script_UK.R

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
# Code to compare finalsize model outputs with ONS data in UK
# Adam Kucharski (2022)
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

library(devtools)
devtools::install_github("epiverse-trace/finalsize",force=T)
devtools::install_github("adamkucharski/quickfit")

library(finalsize)
library(rio)
library(tidyverse)
remotes::install_cran("socialmixr"); library(socialmixr)


# Define dates for BA4/5 wave in UK ------------------------------------------------
xrange <- as.Date(c("2022-06-01","2022-09-01"))

# Import ONS antibody data ------------------------------------------------
# Use higher threshold for cutoff (i.e. 800)

url_abs <- "https://www.ons.gov.uk/visualisations/dvc2323/fig1/datadownload.xlsx"
import_ons_abs_eng <- rio::import(file = url_abs,which = 6) # =2 is lower threshold

# Format data
data_ons_abs_eng <- import_ons_abs_eng[-(1:4),] # Remove header
names(data_ons_abs_eng) <- import_ons_abs_eng[4,] # Add col names
data_ons_abs_eng <- head(data_ons_abs_eng,-1) # Remove footer

age_columns <- seq(2,26,3)
age_labels <- c(16,25,35,50,60,65,70,75,80)

start_date <- as.Date("2020-12-07") # Define first date in dataset
date_list <- start_date + seq(0,7*(nrow(data_ons_abs_eng)-1),7) # Relabel with actual dates
data_ons_abs_eng <- data_ons_abs_eng |> add_column("date" = date_list)

# Extract start and end antibody levels
abs_ba4_start <- data_ons_abs_eng[data_ons_abs_eng$date>xrange[1],][1,]
abs_ba4_end <- data_ons_abs_eng[data_ons_abs_eng$date>xrange[2],][1,]

abs_ba4_start <- cbind(age_labels,as.numeric(abs_ba4_start[,age_columns])/100) # Antibodies pre-BA4 wave
abs_ba4_end <- cbind(age_labels,as.numeric(abs_ba4_end[,age_columns])/100) # Antibodies post-BA4 wave


# Import ONS incidence estimates ------------------------------------------------------------------
# from https://github.com/epiforecasts/inc2prev

#url_source <- "https://raw.githubusercontent.com/epiforecasts/inc2prev/master/outputs/cumulative_age.csv"
#import_cis_incidence <- rio::import(file = url_cis)

# Import locally owing to remote download error
data_c_incidence <- read_csv("~/Documents/COVID_data/cumulative_age.csv")

data_cis_inc <- data_c_incidence |> filter(name=="cumulative_exposure")

# Recode ages
data_cis_inc <- data_cis_inc |> mutate(variable=recode(variable, 
                                                       "2-10"=2,"11-15"=11,"16-24"=16,"25-34"=25,"35-49"=35,"50-69"=50,"70+"=70))

# Extract relevant values for start and end of wave
total_inf_ba4_start <- data_cis_inc |> filter(date==xrange[1]) |> select(variable,q50) |> arrange(variable)
total_inf_ba4_end <- data_cis_inc |> filter(date==xrange[2]) |> select(variable,q50) |> arrange(variable)

# Match age groups
#age_names <- c("2-10","11-15","16-24","25-34","35-49","50-69","70+")
age_labels_inc <- c(2,11,16,25,35,50,70)

# Plot comparison of total infections (circles) vs antibodies (triangles)
plot(total_inf_ba4_start$variable,total_inf_ba4_start$q50,pch=19,
     xlim=c(0,80),ylim=c(0,1),
     xlab="age",ylab="% 'immune' ")
points(total_inf_ba4_end$variable,total_inf_ba4_end$q50)
points(abs_ba4_start,pch=17,ylim=c(0,1)); points(abs_ba4_end,pch=2)

# Estimate attack rates based on each dataset
attack_abs <- (abs_ba4_end[,2]-abs_ba4_start[,2])/(abs_ba4_start[,2])
attack_prev <- (total_inf_ba4_end$q50-total_inf_ba4_start$q50)/(1-total_inf_ba4_start$q50)

# Plot attack rates based on infections (circles) vs antibodies (triangles)
plot(total_inf_ba4_start$variable,attack_prev,pch=19,
     xlim=c(0,80),ylim=c(0,0.5),
     xlab="age",ylab="attack rate")
points(abs_ba4_start[,1],attack_abs,pch=17,ylim=c(0,1))



# Load UK social mixing data ----------------------------------------------

contact_data <- contact_matrix(
  polymod,
  countries = "United Kingdom",
  age.limits = age_labels_inc, # Use age bands
  symmetric = TRUE
)


# Define model parameters -------------------------------------------------

demography_vector = contact_data$demography$population
#demography_vector <- demography_vector/sum(demography_vector)

contact_matrix = t(contact_data$matrix)
contact_matrix = contact_matrix / max(eigen(contact_matrix)$values)

names(demography_vector) = contact_data$demography$age.group
n_demo_grps = length(demography_vector)

contact_matrix = contact_matrix / demography_vector

alpha = 1 # Assume prior infection fully immunising

# Define susceptibility based on proportion infected
n_risk_grps = 1L
p_susceptibility = matrix(1, n_demo_grps, n_risk_grps)


# Run fitting -------------------------------

# Define likelihood function
log_l <- function(x,a,b){
  
  # Relative susceptibility
  susceptibility = matrix(
    data = (1-a*total_inf_ba4_start$q50), #c(1, alpha*rep(1,n_demo_grps-1)),
    n_demo_grps
  )
  
  # Outbreak size
  output_model <- finalsize::final_size(
    r0 = b,
    contact_matrix = contact_matrix,
    demography_vector = demography_vector,
    susceptibility = susceptibility,
    p_susceptibility = p_susceptibility,
    solver = "iterative"
  )
  
  cis_samples <- 1e4 # Assume 10,000 participants in each age group in ONS study
  
  n_positive <- round(x*cis_samples)
  n_total <- rep(cis_samples,length(x))
  
  likelihood_out <- dbinom(n_positive,
                           size = n_total,
                           prob = output_model$p_infected,
                           log = T
  )
  
  likelihood_out
  
}

# Estimate MLE for R0 (b) and cross-reactivity (a)
mle_estimate <- quickfit::estimate_MLE(log_l,data_in=attack_prev,n_param=2,a_initial=0.5,b_initial=3)
  
# Simulate MLE estimates
susceptibility = matrix(
  data = (1-mle_estimate$estimate[["a"]]*total_inf_ba4_start$q50), #c(1, alpha*rep(1,n_demo_grps-1)),
  n_demo_grps
)

output_model_mle <- finalsize::final_size(
  r0 = mle_estimate$estimate[["b"]],
  contact_matrix = contact_matrix,
  demography_vector = demography_vector,
  susceptibility = susceptibility,
  p_susceptibility = p_susceptibility,
  solver = "iterative"
)


# Plot outputs

plot(total_inf_ba4_start$variable,attack_prev,pch=1,
     xlim=c(0,80),ylim=c(0,0.5),
     xlab="age",ylab="attack rate")   
points(total_inf_ba4_start$variable,output_model_mle$p_infected,col="blue",pch=3,lwd=2)


# Additional reference parameters if needed -------------------------------

# Protection estimates from: https://www.sciencedirect.com/science/article/pii/S1473309922007290
# Previous infection with wildtype: 32·7% (27·7–37·4)
# Previous infection with alpha: 36·6% (32·9–40·1)
# Previous infection with delta: 52·4% (50·9–53·8)
# Previous infection with omicron: 59·3% (46·7-69·0)