diff --git a/main/r/rstride/rmarkdown/social_contacts_usa2026.Rmd b/main/r/rstride/rmarkdown/social_contacts_usa2026.Rmd deleted file mode 100644 index 31d392c4..00000000 --- a/main/r/rstride/rmarkdown/social_contacts_usa2026.Rmd +++ /dev/null @@ -1,339 +0,0 @@ ---- -title: "USA Social Contact Matrices for STRIDE" -author: "Willem L" -date: "`r format(Sys.Date())`" -output: - html_document: - toc: true - toc_depth: 3 - number_sections: true ---- - -```{r setup, include=FALSE} -knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE, - fig.width = 6, fig.height = 4) -knitr::opts_knit$set(root.dir = '/Users/lwillem/opt/stride-751') -``` - -## Preamble - -This report prepares USA-specific social contact matrices for STRIDE. Age-specific contact rates by location (household, school, workplace, community) are derived from the Prem et al. (2020) synthetic contact matrices for the US, and conditioned on the school enrolment, employment and household structure observed in a synthetic Wisconsin (Dane County) population. The resulting age-specific contact rates are exported as an XML file for use by the STRIDE simulator. - -## Configuration parameters - -```{r params} -# load 'contactdata' package -suppressPackageStartupMessages(library('contactdata')) - -# load simid.rtools package -# require(devtools) -# devtools::install_github("lwillem/simid_rtools", force = F, quiet = T) -library('simid.rtools') - -# select country with ISO2 code -sel_country <- 'US' - -# select population file -pop_file <- 'data/pop_usa_wisconsin_dane474k_c1000.csv' - -# set output tag -cdata_tag <- paste0('contact_matrix_usa_conditional') - -# set output directory -output_dir <- smd_file_path('sim_output') -file_name <- file.path(output_dir, cdata_tag) -``` - -## Helper function: contact rates by location - -```{r cnt-data-function} -# define help function to load social data from Prem et al. by location -get_cnt_data <- function(location, country) { - cnt_matrix <- contact_matrix(country = country, - location = location, - data_source = c("2020")) - cnt_count <- rowSums(cnt_matrix) - all_ages <- 0:94 - approx(x = seq(0, 75, 5), - y = cnt_count, - xout = all_ages, - method = 'constant', - rule = 2:2)$y -} -``` - -## Load and process data - -### Contact data by location - -```{r load-contact-data} -# get data by location ---- -cnt_all <- get_cnt_data("all", country = sel_country) -cnt_home <- get_cnt_data("home", country = sel_country) -cnt_workplace <- get_cnt_data("work", country = sel_country) -cnt_school <- get_cnt_data("school", country = sel_country) -cnt_other <- get_cnt_data("other", country = sel_country) -``` - -```{r check-contact-data} -# check -cnt_all_list <- contact_df_countries(countries = sel_country, - location = "all", - geographic_setting = c("all"), - data_source = c("2020") - ) -sum(cnt_all_list$contact[cnt_all_list$age_from == "75_80"]) -cnt_all[76] # note: index 1 is age age 0 -``` - -### Additional contacts - -We accumulate "additional contacts" for cases where the contact data assumes contacts at a location (e.g. school) that the population data does not support (e.g. no enrolment), or where the number of contacts at home outnumbers the number of household members. - -```{r additional-contacts-init} -# define vector to accumulate "additional contact" -cnt_additional <- cnt_all * 0 -``` - -### Demography data - -```{r load-demography} -# load file -pop_usa <- read.csv(pop_file, header = T, sep = ',') - -# rename work_id to workplace_id -names(pop_usa) <- gsub('work', 'workplace', names(pop_usa)) - -# set 'missing' to NA -pop_usa$school_id[pop_usa$school_id == 0] <- NA -pop_usa$workplace_id[pop_usa$workplace_id == 0] <- NA - -# define age breaks for the age distribution -breaks_ages <- (0:95) - 0.5 - -# age distribution -age_counts <- hist(pop_usa$age, breaks = breaks_ages, plot = FALSE)$counts -``` - -### School enrolment - -```{r school-enrolment} -age_counts_school <- hist(pop_usa$age[!is.na(pop_usa$school_id)], breaks = breaks_ages, plot = FALSE)$counts -age_distr_school <- age_counts_school / age_counts -age_distr_school[is.na(age_distr_school)] <- 0 - -# calculate number of contacts conditional on being at school -cnt_school_conditional <- cnt_school * (7 / 5) # school is open 5 out of 7 -cnt_school_conditional[3:5] <- cnt_school_conditional[6] # extrapolate behaviour of age 5 to age 2 to 4 - -# set contact rates for ages not enrolled in school to zero -cnt_school_conditional[age_distr_school == 0] <- 0 - -# define school contacts for ages not at school as "additional" -cnt_additional <- cnt_school * (age_distr_school == 0) -cnt_additional[20] <- cnt_additional[21] # adjust artefact for age 19 (not enrolled in US data, but high number of contacts observed) -``` - -### Employment - -```{r employment} -age_counts_workplace <- hist(pop_usa$age[!is.na(pop_usa$workplace_id)], breaks = breaks_ages, plot = FALSE)$counts -age_distr_workplace <- age_counts_workplace / age_counts -age_distr_workplace[is.na(age_distr_workplace)] <- 0 - -# define conditional number of contacts for all ages as the average of a selection of the (most) active population -workplace_ages <- 18:69 -workplace_ages_select <- 30:49 -cnt_workplace_conditional <- cnt_workplace * 0 # start with zero's -cnt_workplace_conditional[workplace_ages + 1] <- mean(cnt_workplace[workplace_ages_select + 1]) # index = age + 1 - -# calculate number of contacts conditional on being at work -cnt_workplace_conditional <- cnt_workplace_conditional * (7 / 5) # account for working 5 days out of 7 -cnt_workplace_conditional <- cnt_workplace_conditional / mean(age_distr_workplace[workplace_ages_select + 1]) # account for employment rate -``` - -```{r workplace-size-check} -# optional: increase rates to account for small workplaces (n < mean number of contacts)? -workplace_size_count <- table(table(pop_usa$workplace_id)) -# number of people in workplaces with <= 7 people -num_people_workplace_leq7 <- sum(workplace_size_count[1:7] * 1:7) -# proportional to number of workers -num_people_workplace_leq7 / sum(!is.na(pop_usa$workplace_id)) -``` - -### Household contacts - -```{r household-contacts} -# define household sizes -hh_size <- data.frame(table(pop_usa$household_id)) -names(hh_size) <- c('household_id', 'household_size') - -# add household info to population matrix -pop_usa_edit <- merge(pop_usa, hh_size) - -# count number of households by member age and size -num_hh_age_size <- table(pop_usa_edit$age, pop_usa_edit$household_size) -num_age <- rowSums(num_hh_age_size) - -# define matrix to represent number of hh contacts by hh size -mat_cnt_size <- matrix(rep(1:ncol(num_hh_age_size) - 1, length(num_age)), - ncol = ncol(num_hh_age_size), - byrow = T) - -# assume fully connected households, get total number of contacts by age -hh_cnt_size <- num_hh_age_size * mat_cnt_size - -# get mean number of contacts at home by age if household is fully connected -cnt_home_fully_connected <- rowSums(hh_cnt_size) / num_age - -# impute/approximate missing ages -cnt_home_fully_connected <- approx(names(cnt_home_fully_connected), - cnt_home_fully_connected, - 0:(length(cnt_home) - 1))$y - -# set contacts with non-household-members as "additional" -cnt_additional <- cnt_additional + (cnt_home - cnt_home_fully_connected) -``` - -### Community - -```{r community-contacts} -# start from cnt_other and add "additional" -cnt_other_adj <- cnt_other + cnt_additional -``` - -## Explore: conditional vs. unconditional contact rates - -```{r explore-function} -# define function to explore (un)conditional contact rates -plot_conditional_contacts <- function(cnt_orig, cnt_conditional, pop_fraction, plot_main, xlim = c(0, 95)) { - - # define y_limit - ylim <- c(0, max(c(cnt_orig, cnt_conditional)) * (4 / 3)) - - # set margin - par(mar = c(5, 5, 2, 5)) - - # plot (un)conditional contact rates - plot(cnt_conditional, xlim = xlim, - ylim = ylim, - main = plot_main, - xlab = "age", - ylab = "mean number of contacts") - points(cnt_orig, col = 2) - legend('topleft', c('unconditional and 7d average', 'conditional'), fill = 2:1) - - # set population fraction scaling factor - if (!any(is.na(pop_fraction))) { - frac_schale <- max(cnt_conditional) * 3 / 4 - lines(pop_fraction * frac_schale, col = 4) - abline(h = frac_schale, lty = 3, col = 4) - axis(4, seq(0, frac_schale, length.out = 11), labels = seq(0, 10, 1) / 10, las = 2, col = 4) - mtext(paste('fraction enrolled in', plot_main), side = 4, padj = 4, col = 4, adj = 0.1) - } -} -``` - -```{r explore-school, fig.cap='Conditional vs. unconditional contact rates, school'} -plot_conditional_contacts(cnt_school, cnt_school_conditional, age_distr_school, 'school', xlim = c(0, 22)) -``` - -```{r explore-workplace, fig.cap='Conditional vs. unconditional contact rates, workplace'} -plot_conditional_contacts(cnt_workplace, cnt_workplace_conditional, age_distr_workplace, 'workplace') -``` - -```{r explore-household, fig.cap='Conditional vs. unconditional contact rates, household'} -plot_conditional_contacts(cnt_home, cnt_home_fully_connected, NA, 'household') -``` - -```{r explore-community, fig.cap='Conditional vs. unconditional contact rates, other/community'} -plot_conditional_contacts(cnt_other, cnt_other_adj, NA, 'other') -``` - -```{r explore-total, fig.cap='Contact rates, total'} -plot_conditional_contacts(cnt_all, cnt_all, NA, 'total') -``` - -## Save results - -### Store as RData (for R) - -```{r store-rdata} -# start with info on data and methods -cnt_data_meta <- list(data_source = 'USA social contact based on Prem et al (2017)', - method = "Reported average number of contacts by age group, conditional on presence", - author = Sys.info()['user'], - date = format(Sys.time()) - ) - -# add social contact data -social_cnt_data <- cnt_data_meta -social_cnt_data$regular_weekday <- cnt_all -social_cnt_data$regular_weekend <- cnt_all -social_cnt_data$household <- cnt_home_fully_connected -social_cnt_data$school <- cnt_school_conditional -social_cnt_data$workplace <- cnt_workplace_conditional -social_cnt_data$community_weekday <- cnt_other_adj -social_cnt_data$community_weekend <- cnt_other_adj - -save(social_cnt_data, file = paste0(file_name, '.RData')) -``` - -### Store as XML (for STRIDE) - -```{r store-xml} -library('XML') - -# create xml prefix -xml_prefix <- paste0(' This file is part of the Stride software [', format(Sys.time()), ']') - -cnt_matrices_lib <- list(regular_weekday = cnt_all, - regular_weekend = cnt_all, - household = cnt_home_fully_connected, - school = cnt_school_conditional, - workplace = cnt_workplace_conditional, - community_weekday = cnt_other_adj, - community_weekend = cnt_other_adj - ) -# setup XML doc (to add prefix) -xml_doc <- newXMLDoc() - -cnt_matrix_xml <- newXMLNode("matrices", doc = xml_doc) -cnt_matrix_meta <- newXMLNode("metadata", parent = cnt_matrix_xml) -smd_listToXML(cnt_matrix_meta, cnt_data_meta) - -for (i_context in 1:length(cnt_matrices_lib)) { - - # extract data - cnt_context_name <- names(cnt_matrices_lib)[i_context] - cnt_context_values <- cnt_matrices_lib[[i_context]] - - # add data to XML - cnt_context <- newXMLNode(cnt_context_name, parent = cnt_matrix_xml) - for (i in 1:length(cnt_context_values)) { - - participant <- newXMLNode("participant", parent = cnt_context) - - part_age <- newXMLNode("age", parent = participant) - xmlValue(part_age) <- paste(i) - - contacts <- newXMLNode("contacts", parent = participant) - - contact <- newXMLNode("contact", parent = contacts) - age <- newXMLNode("age", parent = contact) - xmlValue(age) <- 'all' - rate <- newXMLNode("rate", parent = contact) - xmlValue(rate) <- paste(cnt_context_values[i]) - } -} - -# create filename for xml output -out_filename <- paste0(file_name, '.xml') - -# save as XML, -# note: if we use an XMLdoc to include prefix, the line break dissapears... -# fix: http://r.789695.n4.nabble.com/saveXML-prefix-argument-td4678407.html -cat(saveXML(xml_doc, indent = TRUE, prefix = newXMLCommentNode(xml_prefix)), file = out_filename) -``` - -The final contact matrices for STRIDE are written to `r file.path(getwd(),out_filename)`.