fixed line breaks + added package namespace

analyses/quantify_transmission/reproduction_number_cluster_size.qmd

@@ -36,27 +36,28 @@ knitr::opts_chunk$set(
 # Load required packages
 library(epichains)
 library(MCMCpack)
-```
 
-```{r}
 # Define data
 mers_clusters = c(rep(1,27),c(2,2),c(3,4),c(4,3),c(5,2),7,13,26)
 
 # Show summary table of frequencies
-freq_df <- as.data.frame(table(mers_clusters)); names(freq_df) <- c("Cluster size", "Frequency")
+freq_df <- as.data.frame(table(mers_clusters))
+names(freq_df) <- c("Cluster size", "Frequency")
+
+# Frequencies of MERS Clusters
+freq_df
 
-# Create a table for the HTML document
-knitr::kable(freq_df, caption = "Frequencies of MERS Clusters")
 # Define likelihood function
 lik_function <- function(param) {
-  if (any(param <= 0)) return(-Inf) # Ensure positive parameters
+  # Ensure positive parameters
+  if (any(param <= 0)) return(-Inf)
   
   # Extract values of R and k
   r_val <- as.numeric(param[1])
   k_val <- as.numeric(param[2])
 
   # Define likelihood
-  log_likelihood <- likelihood(
+  log_likelihood <- epichains::likelihood(
     chains = mers_clusters,
     statistic = "size",
     offspring_dist = rnbinom,
@@ -79,21 +80,26 @@ n_burn <- 1e3
 # Initial guess for c(R,k):
 init_param <- c(R=0.5, k=0.5)
 # Run MCMC to estimate parameters
-result_mcmcpack <- MCMCmetrop1R(lik_function, 
-                                theta.init = init_param, 
-                                burnin = n_burn, 
-                                mcmc = n_iter, 
-                                thin = 1)
+result_mcmcpack <- MCMCpack::MCMCmetrop1R(
+  lik_function, 
+  theta.init = init_param, 
+  burnin = n_burn, 
+  mcmc = n_iter, 
+  thin = 1
+)
 
 # Calculate effective sample size (i.e. measure of MCMC mixing)
-ess_mcmcpack <- effectiveSize(result_mcmcpack)
+ess_mcmcpack <- coda::effectiveSize(result_mcmcpack)
 
 # Plot posterior estimates
-plot(result_mcmcpack)
-# Define helper function to calculate median and 95% credible interval from data.frame of MCMC samples
+# plot(result_mcmcpack)
+
+# Define helper function to calculate median and 95% credible interval
+# from data.frame of MCMC samples
 get_param <- function(x){
-  apply(x,2,function(y){val = signif(quantile(y,c(0.5,0.025,0.0975)),3);
-                        val_text <- paste0(val[1]," (95%: CrI: ",val[2],"-",val[3],")")})
+  apply(x,2,function(y){
+    val <- signif(quantile(y,c(0.5,0.025,0.0975)),3)
+    val_text <- paste0(val[1], " (95%: CrI: ", val[2], "-", val[3], ")")})
 }
 
 # Get posterior median and 95% CrI

analyses/quantify_transmission/reproduction_number_serological_data.qmd

@@ -34,31 +34,31 @@ knitr::opts_chunk$set(
 #| warning: false
 
 # Load required packages
-library(readr)
 library(finalsize)
 library(socialmixr)
-library(dplyr)
 library(MCMCpack)
-```
+library(readr)
+library(dplyr)
+library(ggplot2)
 
-```{r}
 # Load serological data
-hk_serology = read_csv("https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4100506/bin/rspb20140709supp2.csv")
+hk_serology <- readr::read_csv("https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4100506/bin/rspb20140709supp2.csv")
 
 # Allocate ages to bands
-age_limits = c(3,18,25,40,65)
+age_limits <- c(3,18,25,40,65)
 age_group <- sapply(hk_serology$age,function(x){sum(x>=age_limits)})
 hk_serology$age_group <- age_group
 
 # Load HK social mixing data
-contact_data <- contact_matrix(
+contact_data <- socialmixr::contact_matrix(
   polymod,
   countries = "United Kingdom",
   age.limits = age_limits,
   symmetric = TRUE
 )
 
-demography_vector = contact_data$demography$population
+demography_vector <- contact_data$demography$population
+
 # Define model
 model_function <- function(param){
   r0 <- param[1]

analyses/reconstruct_transmission/estimate_infections.qmd

@@ -50,74 +50,82 @@ library(ggplot2) # to generate plots
 # Set number of cores
 withr::local_options(list(mc.cores = 4))
 
-# - - - 
-# Example 1: 
-# reconstruct SARS-CoV-2 infection dynamics in the UK from
-# daily data on deaths, 2020
-# - - -
-
-
 # Extract data on UK COVID deaths and format for EpiNow2
-incidence_data <- incidence2::covidregionaldataUK |> 
+incidence_data <- incidence2::covidregionaldataUK %>% 
+  # convert to tibble format
+  dplyr::as_tibble() %>% # for simpler dataframe output
   # preprocess missing values
-  tidyr::replace_na(list(deaths_new = 0)) |>
+  tidyr::replace_na(list(deaths_new = 0)) %>%
   # compute the daily incidence
   incidence2::incidence(
     date_index = "date",
     counts = "deaths_new",
     count_values_to = "confirm",
     date_names_to = "date",
     complete_dates = TRUE
-  ) |>
+  ) %>%
   dplyr::select(-count_variable)
 
 # Focus on early 2020 period and sort by ascending date
-incidence_data <- incidence_data |> 
-                  dplyr::filter(date<"2020-07-01" & date>="2020-03-01")
+incidence_data <- incidence_data %>% 
+  dplyr::filter(date<"2020-07-01" & date>="2020-03-01")
 
 # Preview data
 incidence_data
 
 # Define parameters
 # Extract infection-to-death distribution (from Aloon et al)
 incubation_period_in <-
-  epiparameter::epidist_db(disease = "covid",epi_dist = "incubation",single_epidist = T)
+  epiparameter::epidist_db(
+    disease = "covid",
+    epi_dist = "incubation",
+    single_epidist = T
+  )
 
 # Summarise distribution and type
 print(incubation_period_in)
 
 # Get parameters and format for EpiNow2 using LogNormal input
-incubation_params <- epiparameter::get_parameters(incubation_period_in) # Get parameters
+incubation_params <- epiparameter::get_parameters(incubation_period_in)
 
 # Find the upper 99.9% range by the interval
 incubation_max <- round(quantile(incubation_period_in,0.999))
 
-incubation_period <- LogNormal(meanlog = incubation_params[["meanlog"]], 
-                               sdlog = incubation_params[["sdlog"]], 
-                               max = incubation_max)
+incubation_period <- EpiNow2::LogNormal(
+  meanlog = incubation_params[["meanlog"]], 
+  sdlog = incubation_params[["sdlog"]], 
+  max = incubation_max
+)
 
 ## Set onset to death period (from Linton et al)
 onset_to_death_period_in <-
-  epiparameter::epidist_db(disease = "covid",epi_dist = "onset to death",single_epidist = T)
+  epiparameter::epidist_db(
+    disease = "covid",
+    epi_dist = "onset to death",
+    single_epidist = T
+  )
 
 # Summarise distribution and type
 print(onset_to_death_period_in)
 
 # Get parameters and format for EpiNow2 using LogNormal input
-onset_to_death_params <- epiparameter::get_parameters(onset_to_death_period_in) # Get parameters
+onset_to_death_params <- epiparameter::get_parameters(onset_to_death_period_in)
 
 # Find the upper 99.9% range by the interval
 onset_to_death_max <- round(quantile(onset_to_death_period_in,0.999))
 
-onset_to_death_period <- LogNormal(meanlog = onset_to_death_params[["meanlog"]], 
-                               sdlog = onset_to_death_params[["sdlog"]], 
-                               max = onset_to_death_max)
+onset_to_death_period <- LogNormal(
+  meanlog = onset_to_death_params[["meanlog"]], 
+  sdlog = onset_to_death_params[["sdlog"]], 
+  max = onset_to_death_max
+)
 
 ## Combine infection-to-onset and onset-to-death
 infection_to_death <- incubation_period + onset_to_death_period
 
 # Plot underlying delay distributions
-plot(infection_to_death)
+# plot(infection_to_death)
+
 # Extract serial interval distribution distribution (from Yang et al)
 serial_interval_in <-
   epiparameter::epidist_db(