From f0a4423fd4e1a2317b5cf60edfe47b5b885d64e9 Mon Sep 17 00:00:00 2001 From: Andree Valle Campos Date: Thu, 3 Oct 2024 17:55:44 +0100 Subject: [PATCH] fixed line breaks + added package namespace --- .../reproduction_number_cluster_size.qmd | 40 ++++++++------ .../reproduction_number_serological_data.qmd | 16 +++--- .../estimate_infections.qmd | 54 +++++++++++-------- 3 files changed, 62 insertions(+), 48 deletions(-) diff --git a/analyses/quantify_transmission/reproduction_number_cluster_size.qmd b/analyses/quantify_transmission/reproduction_number_cluster_size.qmd index a3a5b92..76d486b 100644 --- a/analyses/quantify_transmission/reproduction_number_cluster_size.qmd +++ b/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 diff --git a/analyses/quantify_transmission/reproduction_number_serological_data.qmd b/analyses/quantify_transmission/reproduction_number_serological_data.qmd index eaba079..8379eea 100644 --- a/analyses/quantify_transmission/reproduction_number_serological_data.qmd +++ b/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] diff --git a/analyses/reconstruct_transmission/estimate_infections.qmd b/analyses/reconstruct_transmission/estimate_infections.qmd index 82d6de7..9b6905d 100644 --- a/analyses/reconstruct_transmission/estimate_infections.qmd +++ b/analyses/reconstruct_transmission/estimate_infections.qmd @@ -50,17 +50,12 @@ 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", @@ -68,12 +63,12 @@ incidence_data <- incidence2::covidregionaldataUK |> 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 @@ -81,43 +76,56 @@ 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(