This report provides R code to process the data and perform the analysis of antimicrobial-resistant E. coli in river water reported by:
Holcomb DA, Cantrell ME, Feistel DJ, Meyer M, Sterling J, Chen J, Krishnan K, Im S, Folster JP, Kim J, Burroughs M, Nobles S, Stahl A, Mattioli MC. (2026). Environmental monitoring for antimicrobial resistance in a metropolitan watershed. Journal of Water and Health. https://doi.org/10.2166/wh.2026.219
This analysis code and the study data are available in the manuscript repository. The whole-genome sequencing data are available from NCBI under BioProject PRJNA1423195 with accession numbers SRR55333857–SRR55334224.
Setup
Load packages and set options.
Code
## Optionsoptions(tidyverse.quiet =TRUE)options(dplyr.summarise.inform =FALSE)options(knitr.kable.NA ='')## Packages# Utilitieslibrary(readxl) # I/Olibrary(writexl) # I/Olibrary(lubridate) # date formattinglibrary(prettyunits) # number formattinglibrary(zoo) # time series datalibrary(units) # handle units in dataframes# Spatiallibrary(sf) # spatial data & operationslibrary(nngeo) # spatial operationslibrary(rmapshaper) # spatial operationslibrary(terra) # raster data & operationslibrary(exactextractr) # summarize raster data within polygonslibrary(FedData) # download US federal data sourceslibrary(tigris) # download US census shapefileslibrary(SSNbler) # prepare spatial stream networkslibrary(SSN2) # fit spatial stream network models# Modelslibrary(broom) # extract model estimateslibrary(fitdistrplus) # fit distributions to datalibrary(truncnorm) # draw from truncated distributionlibrary(DescTools) # gives correlation CIs# Figureslibrary(patchwork) # combine plotslibrary(ragg) # more consistent graphics outputlibrary(ggtext) # plot label formattinglibrary(scales) # axis formatting# Tidyverselibrary(knitr) # table renderinglibrary(kableExtra) # table headerslibrary(tidyverse) # load last to avoid conflicts
Data processing
Water samples
The variables measured in surface water grab samples, including microbial concentrations, turbidity, and optical brightener fluorescence, along with sampling site information, are provided in the file sample_data.xlsx. This file is structured with the data in the first sheet and a codebook in the second sheet. Each row in the data sheet corresponds to the measured E. coli concentration in a particular water sample from a single site and date, selected with a given antibiotic at a particular dose.
Read in the sample data, fix the dates, and save as RDS for later use:
# read sample datadf_samp <-read_xlsx("../data/sample_data.xlsx", sheet ="data") |>mutate(dt_collect =as_date(dt_collect))# save in RDS format for stabilitysaveRDS(df_samp, "../out/sample_data.rds")
Sampling site locations
Information about each of the sampling sites is contained in the waster sample data set. Extract the site information from the sample data, condense, and save for later:
Extract the coordinate reference system (CRS) from the site data and use it to convert the lat/lon coordinates for each of the sites into sf format to facilitate future spatial operations. Project the site locations into the NAD83(2011) / Georgia West (m) projection (EPSG:6446) and save for later use:
Polygons representing the perimeter of the drainage basin contributing to each sampling location were obtained using the USGS StreamStats Batch Processor v1.5.1. Sample site locations have already been manually adjusted in to coincide with the nearest feasible USGS StreamGrid raster cell prior to basin delineation. In addition to the sampling sites, a drainage basin was delineated for all areas upstream of the USGS monitoring station at Buford Dam (USGS02334430) to encompass Lake Lanier and all the Upper Chattahoochee watershed. This basin upstream of Buford Dam is used to truncate the drainage basins for all the sampling sites to restrict them to downstream of Lake Lanier, as all the river mainstem sites would otherwise share fully overlapping drainage basins upstream of the dam.
The raw drainage basin polygons output by USGS StreamStats are provided as the shapefile basin_raw.shp in the compressed folder basin.zip. After downloading and extracting the compressed file, read in the raw basins shapefile and perform spatial quality control operations, including filling holes in the polygons and validating the resulting geometries:
# raw basin polygons from streamstats shapefilebasin_raw <-read_sf("../gis/shape/basin/basin_raw.shp")# join to site info, fill holes, validatebasin_full <-readRDS("../out/site_data.rds") |># start with site info# add a row corresponding to the Buford Dam basin, upstream of all the other sites.add_case(site_num =1, site_id ="BUD", site_lab ="Buford Dam", site_type ="lake") |>left_join(basin_raw, by =c("site_id"="Name")) |># join the basin polygonsselect(site_num, site_id, site_lab, site_type, geometry) |>st_as_sf(crs =st_crs(basin_raw)) |># convert back to sf format using original crsst_transform(crs =6446) |># project as NAD83(2011) / Georgia Westst_remove_holes() |># remove holes from the polygons st_make_valid() |># validate geometry arrange(site_num) |># order by site number (roughly upstream -> downstream)mutate(site_id =fct_inorder(site_id),epsg ="EPSG:6446",proj ="NAD83(2011) / Georgia West")saveRDS(basin_full, "../gis/sf/basin_full.rds")
Separate the Buford Dam drainage basin (upstream-most) from the rest:
# upper basin as geometry so the attributes don't conflictbasin_upper <- basin_full |>filter(site_id =="BUD")# everything elsebasin_lower <- basin_full |>filter(site_id !="BUD")
Use the upper basin to clip the rest of the basins so they are all truncated at Buford Dam, then re-validate the geometries:
Land cover data were obtained from the USGS National Land Cover Database (NLCD) 2019 vintage via the FedData package. Both land cover classes and percent impervious surface data were obtained as separate raster files. For mapping purposes, the full study area scope (including the shared watershed upstream of Buford Dam) was used to select the raster coverage area.
Merge all the basins, including the basin above the dam, to obtain the extended scope for which to download land cover rasters:
Add a small (2 km) buffer to the scope so the retrieved rasters fully cover the study area extent and transform to WGS84 to work with the FedData package:
basin_scope_buff <- basin_scope_full |>st_buffer(dist =2000) |>st_transform(crs =4326) # transform to WGS84 to work with FedData
Download NLCD LULC 2019 and NLCD Impervious 2019 as separate rasters and save as GeoTIFF format:
Because re-projecting rasters requires a non-reversible transformation of the data, extract the existing CRS information from the downloaded rasters so it can be used to re-project the vector data for analyses:
# reload from saved TIFFnlcd_lulc <-rast("../gis/lulc/nlcd_2019_lulc.tif")# extract CRS as sfnlcd_crs <-st_crs(nlcd_lulc)saveRDS(nlcd_crs, "../gis/lulc/nlcd_crs_sf.rds")
Meteorology
Meteorological data are obtained from the NOAA GHCN-Daily records using the FedData package. The following code demonstrates how the data were originally obtained, but because the data available from NOAA may change over time, re-executing this code could return a slightly different data set than was used for the canonical analysis reported in the manuscript.
Set a search zone for weather stations within 20 km of the study area (defined as the truncated drainage basins), then transform to WGS84 to work with the FedData package:
template_fd <-readRDS("../gis/sf/basin_scope.rds") |>st_buffer(dist =20000) |># search within 20 km buffer of study areast_transform(crs =4326) # transform to WGS84 to work with FedData
Different stations report different weather variables, so download the GHCN data separately for precipitation and temperature to ensure broadest swath of stations are available for each for both full years of the study period (2018 & 2019):
The data cover septic system locations for 15 counties in the metropolitan Atlanta, GA area. These 15 counties fully encompass the study area, with the exception of a small portion of the Fairburn Rd drainage basin (downstream-most site), which extends into Carroll County, not covered in the septic system data. In order to associate septic systems with drainage basins, the Carroll County portion of the Fairburn Rd basin must be masked and the truncated basin area used in density calculations.
Download county boundary polygons from the US Census using the tigris package, project, and reformat:
Read and project septic system locations, then perform a spatial join with the basin polygons and save:
# read and project septic locationssf_septic <-readRDS("../gis/sf/septic_loc.rds") |>st_transform(crs =6446) |># project as NAD83(2011) / Georgia Westmutate(epsg ="EPSG:6446",proj ="NAD83(2011) / Georgia West")# spatial join to basinsseptic_loc <-st_intersection(sf_septic, basin_clip) |>select(septic_id,starts_with("site"),starts_with("wshed"),starts_with("dt"), source,septic_lat = lat, septic_lon = lon, epsg, proj, geometry)saveRDS(septic_loc, "../gis/sf/septic_pt.rds")
Variable generation
The variables used to represent E. coli outcomes and environmental factors in the analyses were created from the data sets obtained above through additional processing.
E. coli concentration: imputation of censored observations
Raw concentration data for both generic and resistant E. coli include values that are problematic for statistical analysis of microbial concentrations. First, the concentration of the target in samples where it was not detected was not necessarily zero, just less than the concentration that corresponds to the lower limit of detection (LOD). Because zeros are also problematic for logarithmic transformation, which is typically appropriate for analyzing microbial concentration data that scale geometrically rather than arithmetically, it is useful to consider non-detects as censored below the LOD. Similarly, samples in which the target is so abundant that all the wells in the most probable number (MPN) test are positive (i.e., TNTC: “too numerous to count”) are only known to have a concentration that is greater than that corresponding to the upper limit of quantification (ULOQ) of the test.
Censoring below the LOD and above the ULOQ correspond to left and right censoring, respectively, and can be addressed using imputation techniques to substitute informed values for the left- and right-censored observations to facilitate log-transformation and downstream analysis. Imputation is performed separately for each target (i.e., each for each selective antibiotic and dose) by fitting a probability distribution to the censored data and drawing a random value for each censored observation from the fitted distribution truncated below the LOD or above the ULOQ for left- and right-censored observations, respectively.
Prepare a data set of microbial concentrations for censored imputation processing. Because few meropenem-nonsusceptible E. coli (MEM-Ec) were detected and a relatively small number of samples were selected by cefotaxime at 4 μg/mL, restrict the data set to generic E. coli, cefotaxime-selected E. coli at 1 μg/mL (ESBL-Ec), and colistin-selected E. coli (CST-Ec). All concentrations will be log10-transformed for this analysis; for this reason, non-detect and TNTC observations have already been replaced with their respective LOD or ULOQ values, respectively, along with an indicator variables denoting whether they are left-, right-, or non-censored. These substituted values are simply for convenience and will be treated as missing in the imputation procedure.
For compatibility with the function from the fitdistrplus package for fitting censored distributions to these data, the concentration variables (mpn, mpn_lod, and mpn_uloq) along with this censored indicator variable (cens) need to be converted into a pair of variables named left and right that both contain the observed concentration when not censored and contain a missing value in corresponding column and the censoring value in the opposite column for censored observations. I.e., a left-censored observation should be missing in left and contain the LOD value in right, while a right-censored observation will contain the ULOQ value in left and be missing in right:
df_micro_cens <-readRDS("../out/sample_data.rds") |>filter(drug %in%c("cefotaxime", "colistin", "none"), drug_dose !=4, site_num >2,!is.na(mpn)) |>group_by(drug, drug_dose) |>mutate(nobs =n(),npos =sum(cens =="no"),ncens =sum(cens =="left"| cens =="right")) |>ungroup() |>mutate(# log10-transform the concentration and limitsacross(c(mpn, mpn_lod, mpn_uloq), ~log10(.x)),# create left-censored values variableleft =case_when(cens =="left"~NA_real_, cens =="no"~ mpn, cens =="right"~ mpn_uloq),# create right-censored value variableright =case_when(cens =="right"~NA_real_, cens =="no"~ mpn, cens =="left"~ mpn_lod))
Extract the number of censored observations by outcome:
Fit the censored normal distribution and extracte the estimated mean and standard deviation with their standard errors:
# fit distributionfit_cens_ctx <-fitdistcens(censdata = df_cens_ctx,distr ="norm")# extract parameter estimatesest_cens_ctx <-tibble(drug ="cefotaxime",nobs = fit_cens_ctx$n, # number of observations used to fit the distributionmu = fit_cens_ctx$estimate[1], # estimated mean of the normal distributionmu_se = fit_cens_ctx$sd[1], # standard error of the mean estimatesigma = fit_cens_ctx$estimate[2], # estimated standard deviation of the normal distributionsigma_se = fit_cens_ctx$sd[2]) # standard error of the SD estimate
The MLE means and standard deviations from the fitted censored normal distributions are used to parameter truncated normal distributions, from which suitable values are drawn to represent the censored observations. A random draw for each censored observation is used for analysis, but the expected value of the distribution truncated at the sample-specific LOD/ULOQ is also calculated for each censored observation for additional context.
Set a random seed for reproducibility, join the estimates to the censored observation data set, prepare the estimates for drawing from a truncated normal distribution. Define the lower bound a and the upper bound b for each sample to correspond with the censoring status of that observation. Imputes values for every observation (including the non-censored ones, which are then drawn from a normal distribution without truncation), retaining the imputed values only for censored observations:
# seed for reproducibilityset.seed(30333)# join estimates to data and imputeimp_cens <- df_micro_cens |>select(starts_with("site"),starts_with("drug"), wk_collect, dt_collect, samp_event, samp_id, test_id, lab_vol, uid, uid_gen, uid_wq, wq_turb, wq_optic, target, detect, mpn, cens, mpn_lod, mpn_uloq, left, right) |># merge distribution estimates left_join(est_cens, by ="drug") |># add bounds variablesmutate(# Define lower bound variablea =case_when(cens =="left"~-Inf, cens =="no"~-Inf, cens =="right"~ mpn_uloq),# Define upper bound variableb =case_when(cens =="left"~ mpn_lod, cens =="no"~Inf, cens =="right"~Inf)) |># perform imputation and only keep censoredmutate(imp_e =etruncnorm(a, b, mu, sigma), # expected value of truncated normal distimp_r =rtruncnorm(1, a, b, mu, sigma), # random draw from truncated normal distmpn_imp_e =ifelse(cens =="no", mpn, imp_e), # only keep expected value for censored obsmpn_imp_r =ifelse(cens =="no", mpn, imp_r)) # only keep random draw for censored obs
Clean up the data set and join the imputed values for generic E. coli to their corresponding resistance-selected (and non-selected [generic], for consistency) E. coli observations:
Basin land cover: distance-weighted, areal proportion
Variables representing the influence of land cover classes on each sampling location are generated by weighting each raster cell of a given class contained within a site’s drainage basin by the inverse of the squared Euclidian distance between that cell’s centroid and the sampling location. Summing these cells by class provides inverse distance squared weighted (IDSW) cell counts, which are divided by the total IDSW cell count across all classes to obtain the IDSW proportion of the basin represented by that land cover class. Similarly, the IDSW mean imperviousness is obtained as a weighted average of the percent imperviousness of each cell in the drainage basin.
Calculate distance and weights
Read in the site locations and basin polygons and reproject to match the NLCD projection extracted previously:
Write a function to calculate the distance between a sampling site and the centroid of each raster cell. The function should take as inputs a raster with the data to be extracted and an sf dataframe with the point locations to calculate distances from. The output will be a raster stack with a layer for each sampling location where the raster value in each cell records its distance from the sampling point. Include functionality to encode the resulting rasters with a weight based on the inverse of the distance (raised to an arbitrary power) rather than the raw distance itself.
site_dist <-function(x_rast, # raster to calculate distances for y_sf, # sf tibble with point locations to measure distances fromidw =FALSE, # return inverse distance weights instead of the distance?power =1# distance raised to power for calculating weights ){ n_feature <-nrow(y_sf) # number of sites to calculate distances to# loop through each site dist_list <-list() # Prepare empty list object for(i in1:n_feature){ y <- y_sf[i, ] # select site x <-distance(x_rast, y) # calculate distances# apply inverse distance weightingif(idw){ x_idw <-1/ (x^power) # raise distances to specified power dist_list <-c(dist_list, x_idw) # insert weights into list object } else { dist_list <-c(dist_list, x) # insert distances into list object } dist_stack <-rast(dist_list) # convert list to raster stack }return(dist_stack)}
Apply function to calculate inverse distance squared weights by specifying idw = TRUE and power = 2. Because the LULC and impervious surface rasters contain identical cell locations and extents, and differ only in the values recorded for each cell, the distance function need only be applied once to obtain the raster stack of weights that can be used for both:
# read in LULC raster (impervious would also work)nlcd_lulc <-rast("../gis/lulc/nlcd_2019_lulc.tif")# calculate inverse distance weightsnlcd_idsw <-site_dist(nlcd_lulc, site_nlcd,idw =TRUE, power =2)# name the raster layers by corresponding site idnames(nlcd_idsw) <- site_nlcd$site_id# save to drive as GeoTIFFwriteRaster(nlcd_idsw, "../gis/lulc/nlcd_site_idsw.tif",overwrite =TRUE)
Extract weighted land cover class proportions
Write a function for use with the exactextractr package to calculate the weighted counts of each class for a given basin, as well as basin-wide statistics. The exactextractr::exact_extract() function calculates fractional coverage for cells partially covered by the basin polygon and returns a dataframe with a row-per-cell with computed values given pre-specified names (e.g., coverage_fraction) that are accessed directly in the function before summarizing land cover class:
extract_idw <-function(df){ df |># calculate cell-level characteristicsmutate(cover_tot =sum(coverage_fraction), # total cells in basin, accounting for partial coveragefrac = coverage_fraction / cover_tot, # proportion of total in this cellcover_wt = coverage_fraction * weight, # cell's weighted coverage fractioncover_wt_tot =sum(coverage_fraction * weight), # total weighted cells in basinfrac_wt = cover_wt / cover_wt_tot, # proportion of weighted total in this cellwt_tot =sum(weight)) |># sum of all the weights across the basin# specify variables to retain after summarizing# only one site at a time so site and tot variables are static for all cells# `value` contains the class and is the only variable that separate summaries are calculated bygroup_by(site_num, site_id, value, cover_tot, cover_wt_tot, wt_tot) |># tally up cellssummarize(freq =sum(frac), # raw class proportion freq_wt =sum(frac_wt), # weighted class proportionwt =sum(weight), # sum of weights within classcount =sum(coverage_fraction), # raw class cell count in basincount_wt =sum(cover_wt), # weighted class cell count in basinarea =sum(area), # total area covered by classarea_wt =sum(area * weight)) # weighted area covered by class}
Write a second function to loop through each basin and its corresponding IDSW raster in the stack to apply the extraction function to each basin in turn. Include an option :
loop_idw <-function(x_rast, # raster with the values to extract y_poly, # sf with basin polygons (one per row) wt_rast, # raster stack with weights (must perfectly align with x_rast) include_cols =NULL# character vector of column names to retain in output ){ n_feature <-nrow(y_poly) # number of basins extract_df <-tibble() # prepare empty dataframe# loop through each basinfor(i in1:n_feature){ y <- y_poly[i, ] # select basin wt <- wt_rast[[i]] # select weights raster x <-exact_extract(x_rast, y,fun = extract_idw, # custom weighting functionweights = wt,summarize_df =TRUE, # produce summarized output dataframeinclude_cols = include_cols, # column names to retain (as character vector)include_area =TRUE, # calculate area within basinprogress =TRUE) # show progress indicator extract_df <-bind_rows(extract_df, x) # add basin summary to dataframe (inefficient) }return(extract_df)}
Apply the functions to the NLCD LULC raster to calculate IDSW land cover class proportions:
Each cell in the land cover raster contains a single, discrete value corresponding to the land cover class assigned to that cell (cell value is a categorical variable). By contrast, the imperviousness raster contains an integer value between 0 and 100 that represents the proportion of that cell that is covered by impervious surfaces. Rather than a weighted sum of the cells assigned to each class, a weighted mean imperviousness across all cells in the basin in calculated.
Weighted mean operations are built into exactextractr::exact_extract(), so revise the looping function to apply the built-in statistical functions (which require slightly different syntax in the arguments) to each basin in turn:
loop_idw_impv <-function(x_rast, # raster with the values to extract y_poly, # sf with basin polygons (one per row) wt_rast, # raster stack with weights (must perfectly align with x_rast) include_cols =NULL# character vector of column names to retain in output ){ n_feature <-nrow(y_poly) # number of basins extract_df <-tibble() # prepare empty dataframe# loop through each basinfor(i in1:n_feature){ y <- y_poly[i, ] # select basin wt <- wt_rast[[i]] # select weights raster x <-exact_extract(x_rast, y,fun =c("count","mean", "weighted_mean","stdev", "weighted_stdev"),weights = wt,append_cols = include_cols,progress =TRUE) extract_df <-bind_rows(extract_df, x) # add basin summary to dataframe }return(extract_df)}
Apply the functions to the NLCD LULC raster to calculate IDSW land cover class proportions:
The NLCD land cover classes are narrowly defined and are collapsed into broader categories for analysis:
Developed: 22, 23, 24
Forest: 41, 42, 43, 90
Open: a combined category encompassing
Agriculture: 81, 82
Barren: 31
Grassland: 52, 71, 72, 95
Open developed: 21
For convenience, a key is provided at nlcd_lulc_key.xlsx to map the numerical class values contained in the NLCD LULC raster to the condensed categories used in the analysis. Read in the key, join to the IDSW proportion estimates, and aggregate the class proportions up to the broader categories:
# land use classesnlcd_class <-read_xlsx("../data/nlcd_lulc_key.xlsx")# join and aggregatebasin_lulc_long <-readRDS("../out/basin_lulc_idsw.rds") |>left_join(nlcd_class, by =c("value")) |>mutate(method ="idsw") |>group_by(site_num, site_id, cat, method) |>summarise(freq =sum(freq),freq_wt =sum(freq_wt),cells =sum(count)) |>group_by(site_num, site_id, method) |>mutate(freq_tot =sum(freq),freq_wt_tot =sum(freq_wt),cells_tot =sum(cells)) |>ungroup()
Pivot the aggregated LULC category values wider to obtain separate variables for each class:
Daily weather: distance-weighted, rolling aggregate
Weather variables (precipitation in mm and maximum temperature in °C) are computed for each sampling location on each date of the study period (September 2018 - September 2019). Data from weather stations are weighted by the inverse of the squared distance of the station from the sample location to interpolate IDSW weather values for each site and date. Water samples are assigned the IDSW value for the day prior to that sample’s collection to ensure the weather values used reflect antecedent meteorological conditions. For precipitation, rolling sums over the preceding 2 days (“2-day rainfall”) and preceding 7 days, excluding the 2 most recent days (“3-7 day rainfall”), are computed to represent more- and less-recent cumulative antecedent precipitation, respectively.
Inverse distance squared interpolation
Load weather station and sampling site location data and re-project to match:
# weather station locationsloc_ghcn <-readRDS("../gis/sf/meteo_pt.rds") |>st_transform(crs =6446) # project as NAD83(2011) / Georgia West# sample site locationsloc_site <-readRDS("../gis/sf/site_pt.rds") |>st_transform(crs =6446) # project as NAD83(2011) / Georgia West
Compute the Euclidean distance matrix between weather stations and sampling sites:
Read in weather data, join distances, drop missing values, and add an indicator ranking values for a given weather type on each day on the basis of distance from the sampling site to assist with identifying “neighboring” stations:
df_meteo_join <-readRDS("../out/ghcn_data.rds") |>filter(!is.na(meteo_val)) |># drop missing meteo valuesinner_join(dist_site_ghcn, # join distances to values, retaining only records with bothby ="ghcn_id", # join by station idrelationship ="many-to-many"# allow multiple matches on both sides (multiple stations inform multiple sites) ) |>arrange(site_id, meteo_type, dt_meteo, dist) |># ascending sort by distance within each site, type, and dategroup_by(site_id, meteo_type, dt_meteo) |>mutate(neighbor =row_number()) |># label values within a site, type, date with increasing distanceungroup()
Calculate inverse distance squared weights, identify k-nearest neighbors, and compute simple and weighted means to interpolate values from stations to sampling sites for each study date:
df_meteo <- df_meteo_join |>mutate(wt =1/ (dist^2), # IDSW weightwt_val = meteo_val * wt, # weighted weather valuewt_nn =ifelse(neighbor ==1, meteo_val, NA), # identify nearest stationwt_5nn =ifelse(neighbor <=6, meteo_val, NA) # identify 5 nearest stations ) |>group_by(site_id, dt_meteo, meteo_type, meteo_unit) |># aggregate by site, date, and weather typesummarise(n_station =sum(!is.na(meteo_val)), # count stations with valueswt_tot =sum(wt, na.rm =TRUE), # sum of weightsdist_near = dist[neighbor ==1], # distance to nearest stationdist_5nn =ifelse(n_station <5, # distance to 5th nearest station dist[neighbor == n_station], dist[neighbor ==5]),dist_far =max(dist, na.rm =TRUE), # distance to furthest stationmeteo_nn =mean(wt_nn, na.rm =TRUE), # assign value from nearest stationmeteo_5nn =mean(wt_5nn, na.rm =TRUE), # assign mean of 5 nearest stations meteo_avg =mean(meteo_val, na.rm =TRUE), # simple mean of all stationsmeteo_idsw =sum(wt_val, na.rm =TRUE) / wt_tot # IDSW mean ) |>ungroup() |>select(-wt_tot) |>mutate(across(meteo_nn:meteo_idsw,~ifelse(is.nan(.x), NA, .x))) # replace NAN with NAsaveRDS(df_meteo, "../out/meteo_site_idsw.rds")
Rolling cumulative precipitation
Compute moving-window (“rolling”) sums of precipitation in the previous 2 days, 7 days, and 3 - 7 days:
df_roll_prcp <- df_meteo |>rename(phenom = meteo_type, unit = meteo_unit) |># change names for conveniencefilter(phenom =="prcp") |># just rainpivot_longer(starts_with("meteo_"), # long formatnames_to ="method", # weighting method (idsw, simple, nn, 5nn)values_to ="meteo0"# value on day 0 (sample date) ) |>mutate(method =str_remove(method, "meteo_")) |># drop prefixarrange(site_id, phenom, method, dt_meteo) |># order by date last (important for rolling)group_by(site_id, phenom, unit, method) |># by site, type, and weighting methodmutate(meteo1 =lag(meteo0, n =1L), # value on previous day (enforce temporality)meteo2 =rollapply(meteo1, width =2, # 2-day moving windowFUN = sum, # rolling sumfill =NA, # fill values outside window with NApartial =TRUE, # allow missing values in windowalign ="right"), # assign output to right-end of window (most recent)meteo7 =rollapply(meteo1,width =7, # 7-day moving windowFUN = sum,fill =NA,partial =TRUE,align ="right"),meteo3_7 =rollapply(lag(meteo0, 3), # start with 3 days priorwidth =5, # 5-day window starting day -3, ending day -7FUN = sum,fill =NA,partial =TRUE,align ="right"),rollfun ="sum"# record the rolling function applied ) |>ungroup()
Assign previous day’s temperatures to each study date. Calculate rolling mean temperatures over the same windows used for rainfall, for consistency in the data structure:
Code
df_roll_temp <- df_meteo |>rename(phenom = meteo_type, unit = meteo_unit) |># change names for conveniencefilter(phenom !="prcp") |># not rainpivot_longer(starts_with("meteo_"), # long formatnames_to ="method", # weighting method (idsw, simple, nn, 5nn)values_to ="meteo0"# value on day 0 (sample date) ) |>mutate(method =str_remove(method, "meteo_")) |># drop prefixarrange(site_id, phenom, method, dt_meteo) |># order by date last (important for rolling)group_by(site_id, phenom, unit, method) |># by site, type, and weighting methodmutate(meteo1 =lag(meteo0, n =1L), # value on previous day (enforce temporality)meteo2 =rollapply(meteo1, width =2, # 2-day moving windowFUN = mean, # rolling meanfill =NA, # fill values outside window with NApartial =TRUE, # allow missing values in windowalign ="right"), # assign output to right-end of window (most recent)meteo7 =rollapply(meteo1,width =7, # 7-day moving windowFUN = mean,fill =NA,partial =TRUE,align ="right"),meteo3_7 =rollapply(lag(meteo0, 3), # start with 3 days priorwidth =5, # 5-day window starting day -3, ending day -7FUN = mean,fill =NA,partial =TRUE,align ="right"),rollfun ="mean"# record the rolling function applied ) |>ungroup()
Treated wastewater and CSO inputs potentially impacting each site are represented by the aggregate maximum permitted load from all upstream outfalls, weighted by the inverse of the squared hydrological distance between the site and each upstream outfall along the stream network. While conceptually a similar IDSW operation, this approach requires measuring the distance along the stream network and restricting the calculations to outfalls that are both flow connected to, and upstream of, the sampling site. These operations are facilitated by constructing a spatial stream network (SSN) using the SSNbler and SSN2 packages, the later of which provides functions for identifying upstream, flow-connected points and calculating network distances.
Stream vector data were obtained from USGS National Hydrography Dataset Plus (NHDPlus) High Resolution and manually corrected in ArcGIS Pro to resolve topological inconsistencies, ensuring the entire network within the study watershed was flow-connected and discharged at a single terminal outlet located at the downstream-most sampling site (Fairburn Rd). The topology correction process created a small number of new, very short reaches to connect sections of the network. These new reaches are not in the NHD database and therefore lack cumulative drainage area values, which will later be used as weights as a proxy for flow contributions downstream and cannot be missing or zero. The smallest existing cummulative area (area_cum) value in the dataset is substituted for reaches with area_cum = 0 to overcome this limitation.
The pre-processed stream network vector data are provided as the shapefile stream_line.shp in the compressed folder stream.zip. After downloading and extracting the compressed file, read in the stream shapefile and substitute the minimum cumulative area value for reaches with zero cumulative area values and save the sf for future use:
Additive function values (AFV) are numerical values associated with each edge (stream reach) that represent the contribution of that edge to the edges downstream of it. For stream networks, the AFV can encode the flow contributed by each reach (which is an additive function of the flows from all reaches upstream). Because flow is variable and unknown for most reaches, the cumulative area drained by each reach is used as a proxy for its potential flow. The AFVs are not used directly for the network distance calculations but are required to create a valid SSN object, from which the distances are later calculated.
Assign AFVs to the stream network edges and compute the AFVs at each site:
Assemble an SNN from the LSN and various components generated above; while LSNs are network data objects and can be reused, SSNs are data models that are specific to the sites and data associated with them and are best assembled separately for each use case. For this application, record the sampling sites as the “observed sites” (obs_sites) and the outfall locations as the “prediction sites” (preds_list):
Two distance matrices are returned in a list; only the first, with the sites in rows and the outfalls in columns, provides the flow-connected, upstream network distances required. Extract the correct distance matrix and label with IDs extracted from the SSN elements:
Weighted total (cumulative) permitted loads are used to represent the influence of wastewater infrastructure on sites. Because weighted sums are not normalized by the total weights, as was the case for IDSW interpolation of weather data and IDSW land cover proportions. As such, the weighted output values are dependent on the distance unit used and cannot be interpreted in isolation or compared to other contexts. For convenience, distance is converted to kilometers, such that an outfall 1 km from the site receives a weight of 1, outfalls further than 1 km from the site receive smaller weights less than 1, and outfalls closer than 1 km to the site receive larger weights greater than 1.
Calculate weights and and IDSW cumulative load on the basis of permitted maximum discharge separately for WWTP treated effluent outfalls and for CSOs:
As point locations lacking additional scale information (like the permitted discharge for outfalls) and heterogeneously distributed within the study area, inverse distance weighting approaches are less appropriate for constructing site-level variables from septic system locations. Rather, the simple spatial density of all septic systems within a drainage basin is used to represent the potential influence of septic dystems on a sampling site.
Load the septic locations, tally up the systems in each drainage basin, and calculate the spatial density; because some basins contained no septic systems, restore the dropped basins and set septic density to zero:
septic_dens <-readRDS("../gis/sf/septic_pt.rds") |>st_drop_geometry() |>group_by(site_num, wshed_area_clp) |>summarise(septic_n =n()) |>ungroup() |>mutate(septic_dens = septic_n / wshed_area_clp) |># some basins had no septic so restore and set density to zerofull_join(tibble(site_num =3:13), by ="site_num") |>mutate(across(starts_with("septic_"), ~replace_na(.x, 0))) |>arrange(site_num)saveRDS(septic_dens, "../out/septic_density.rds")
Data set assembly
The generated variables are assembled into a full data set of E. coli observations and the environmental factors corresponding to each sampling site and date. Subsets of this full data set are extracted and assembled into specific SSNs to facilitate analysis by spatial generalized linear mixed models (GLMM) that account for flow-connected network distances between sites and and temporal clustering by time of sample collection.
Full data table
Read in the processed sample data (E. coli imputed concentrations, water quality; meropenem-selected and cefotaxime-selected at 4 μg/mL already excluded):
df_samp <-readRDS("../out/sample_imputed.rds")
Read in the rolling weather data, restrict to study time period and desired variables, pivot to wide format (each weather variable is a column), and scale the rain values to units of cm (10 mm):
Read in land cover proportions and mean imperviousness, select desired variables, and scale so all land cover variables are in units of 10 percentage points and basin area is log10-transformed:
Read in outfall loads, select desired variables, differentiate between WWTPs and CSOs, and restore sites without outfalls, substituting zero for missing load:
Read in septic density and join all the flat data together:
# septic densitydf_septic <-readRDS("../out/septic_density.rds") |>select(site_num, septic_n, septic_dens)# combine for analysisdf_analysis <- df_samp |>left_join(df_lulc, by =c("site_num", "site_id")) %>%left_join(df_outfall, by ="site_num") %>%left_join(df_septic, by ="site_num") %>%left_join(df_meteo, by =c("site_id", "dt_collect"="dt_meteo"))saveRDS(df_analysis, "../out/analysis_data.rds")
With exceptions for the outcome E. coli concentrations and raw counts, variables are standard deviation-scaled for analysis. Variables without a meaningful zero value are also centered around the mean prior to scaling. Calculate the mean and standard deviation of each variable by the selective drug used to inform centering and scaling:
The SSN objects are somewhat particular about the data used to assemble them and the opportunities for subsetting post-assembly. While inefficient, a reliable approach is to assemble separate SSNs for each set of outcome variables that are to be analyzed. The E. coli occurrence under each selective antibiotic is a separate outcome for the purposes of these analyses. Furthermore, the co-measured sample-level covariates (turbidity, optical brighteners, and generic E. coli for the antibiotic-selected observations) are missing for some of the E. coli observations, resulting in different sample sizes under complete-cases analyses. Separate SSNs are therefore assembled for each combination of E. coli outcome and covariate completeness. Fortunately, the majority of environmental factors covariates considered, which were derived from comprehensive spatial data sources, are complete for all observations, allowing a single SSN to be used for all analyses of a given outcome that do not include co-measured, sample-level covariates.
Landscape network
The LSN previously created in the process of calculating network distances to wastewater outfalls has those outfall locations associated with it as prediction locations. For simplicity, construct a new LSN with just the site data associated to use in assembling the new SSNs.
Generate new LSN stream network:
# stream vectorstream <-readRDS("../gis/sf/stream_line.rds")# create new LSN for full datasetstream_net <-lines_to_lsn(streams = stream,lsn_path ="../gis/ssn/lsn_full",snap_tolerance =0.1,check_topology =TRUE,topo_tolerance =20,overwrite =TRUE,verbose =TRUE,remove_ZM =TRUE)
Snap site locations to network:
# sample site locationssites <-readRDS("../gis/sf/site_pt.rds")# snap to stream LSNsites_lsn <-sites_to_lsn(sites = sites,edges = stream_net,lsn_path ="../gis/ssn/lsn_full",file_name ="sites",snap_tolerance =100,save_local =TRUE,overwrite =TRUE)
Calculate the upstream distance along each edge (stream segment) and for sample sites:
Join the full, scaled analysis data set with all observations (dropping any missing E. coli outcomes) to the LSN sites:
# full data set for all observationsobs_full <-readRDS("../out/analysis_data_scl.rds") |>filter(!is.na(detect))# merge with LSN sitessites_full <-readRDS("../gis/sf/lsn_sites_afv_pt.rds") |>left_join(obs_full, by ="site_num")saveRDS(sites_full, "../out/analysis_data_lsn.rds")
Cefotaxime-selected SSNs
Assemble SSNs for ESBL-Ec (cefotaxime-selected at 1 μg/mL), the main environmental AR outcome for the modeling analysis.
First assemble an SSN for all the ESBL-Ec observations and create the distance matrix:
Finally, assemble a comparable set of SSNs for the generic E. coli (no antibiotic selection) outcome, skipping the SSN for non-missing generic E. coli as a covariate:
Before developing statistical models to estimate associations between surface water E. coli and environmental factors, conduct descriptive analyses to better understand the contours of the prepared data.
E. coli time series
Explore the occurrence of generic E. coli, ESBL-Ec, CST-EC, and MEM-Ec at each site across the study period by plotting time series of their detection and concentration.
Read in the (non-imputed) sample data time series and arrange for plotting:
To allow for the nice arrangement of 11 different sites, plot the river and tributary sites separately and combine. First plot the river sites, adding LOESS smooths for each E. coli concentration outcome as a visual aid and differentiating between detect, non-detect, and TNTC observations by shape:
This plot reproduces Figure 2 from Holcomb et al. (2026):
Figure 2. Measured concentration (points) and trend (lines) of generic E. coli (orange, solid lines), ESBL-Ec (purple, dashed lines), and CST-Ec (green, dotted lines) by site. Dates on which MEM-Ec were detected are indicated by blue vertical lines. Non-detect (X) and too numerous to count (TNTC; square) samples are plotted at their limit of detection (LOD) or upper limit of quantification (ULOQ), respectively. Trend lines were generated using non-parametric locally estimated scatterplot smoothing (LOESS) to aid visual interpretation.
E. coli correlation
While overall associations between generic and resistant E. coli are estimated later using models to account for clustering in space and in time, compute non-parametric correlations to examine concordance between paired concentrations of generic and resistant E. coli measured in the same samples for each site and site type (river mainstem vs. tributary).
Define a function to compute Kendall’s \(\tau_{b}\) correlations between paired observations with 95% confidence intervals (CI) and export in tidy format:
corr_ci <-function(x, # vector of values y, # vector of paired values (same length as x)conf.level =0.95# CI width ){ nobs =length(x) # number of paired observations ci <-KendallTauB(x, y, # estimate Kendall's tau with CIconf.level = conf.level) pval <-cor.test(x, y, method ="kendall",exact =FALSE) out <-tibble(test ="Tau-b",nobs = nobs,corr =unname(ci[1]),corr_lo =unname(ci[2]),corr_hi =unname(ci[3]),pval = pval$p.value)return(out)}
Because Kendall’s \(\tau_{b}\) correlation is based on ranks and accounts for ties, it is robust to the values used to represent left- and right-censored observations. Prepare the generic and resistant E. coli concentration data by substituting LOD/2 for non-detects and ULOQ + 1 for TNTC observations:
df_corr_sub <-readRDS("../out/sample_data.rds") |>filter(site_num >2,!is.na(detect)) |>mutate(mpn =case_when(cens =="left"~ mpn_lod /2, cens =="right"~ mpn_uloq +1, cens =="no"~ mpn))
Filter for just the generic E. coli to join back to the resistant E. coli from the same samples:
This reproduces the results (if not quite the formatting) of Table S3 in the Supplementary Material of Holcomb et al. (2026):
ESBL-Ec vs Generic
CST-Ec vs Generic
Site
Type
N
τ
95% CI
p
N
τ
95% CI
p
river
river
78
0.75
0.67–0.83
<0.001
74
0.55
0.43–0.67
<0.001
Roswell Rd
river
7
0.84
0.63–1.00
0.011
7
0.75
0.40–1.00
0.028
Johnson Ferry
river
17
0.83
0.68–0.99
<0.001
17
0.75
0.61–0.89
<0.001
Powers Ferry
river
11
0.88
0.68–1.00
<0.001
11
0.60
0.28–0.92
0.021
Paces Ferry
river
12
0.76
0.53–1.00
<0.001
11
0.60
0.28–0.91
0.021
S Cobb Dr
river
10
0.72
0.35–1.00
0.006
9
0.33
-0.20–0.86
0.259
Fairburn Rd
river
21
0.71
0.57–0.85
<0.001
19
0.47
0.15–0.79
0.014
tributary
tributary
53
0.42
0.23–0.60
<0.001
53
-0.19
-0.40–0.02
0.073
Tanyard Creek
tributary
18
0.53
0.26–0.80
0.003
18
-0.01
-0.45–0.42
0.939
Clear Creek
tributary
4
0.67
0.01–1.00
0.174
4
-0.18
-1.00–0.88
0.718
Ansley Stream
tributary
4
0.24
-0.55–1.00
0.655
Proctor Creek
tributary
19
0.01
-0.37–0.38
0.972
19
0.09
-0.42–0.60
0.655
Utoy Creek
tributary
8
0.50
-0.05–1.00
0.083
8
0.50
0.11–0.89
0.127
Total
131
0.60
0.51–0.68
<0.001
127
0.25
0.13–0.36
<0.001
Regression modeling
The primary analysis approach in this study is the application of spatial generalized linear mixed models (GLMM) to estimate associations between environmental factors and the occurrence of generic and resistant E. coli in surface water. The spatial GLMMs are implemented as spatial stream network models using the SSN objects generated previously.
Linear SSN models are used to estimate the expected change in log10 concentration of generic E. coli and ESBL-Ec and logistic SSN models are used to estimate the odds ratio (OR) for CST-Ec detection associated with a unit-increase in each environmental factor. The SSN models all include a random effect for sampling event to account for temporal clustering among contemporaneous samples and a flow-weighted spatial random effect with exponential distance decay to account for spatial autocorrelation between flow-connected sampling locations (a “tail-up” covariance component in the SSN parlance).
Model-fitting function
Write a utility function to execute an SSN model and extract the model estimates in a convenient format with identifying labels:
fit_ssn <-function(lab, # model label to assign to output (string) set, # model set to label output (string) outcome, # outcome assessed (string) response, # response variable representing outcome (string) predictor, # explanatory variable set label (string) error, # covariance structure label (string) re, # random effect label (string) path, # path to SSN for model fit form, # model formularandom =NULL, # random effect formulafamily ="Gaussian", # likelihoodadditive ="afv_area", # AFV column name (string)tailup_type ="none", # correlation model for tail-up (flow-connected) covariance componenttaildown_type ="none", # correlation model for tail-down (flow-disconnected) covariance componenteuclid_type ="none", # correlation model for Euclidean distance covariance componentnugget_type ="nugget", # include unstructured error component (aka nugget)est_method ="reml"# estimation algorithm ){## import SSN net <-ssn_import(path,overwrite =TRUE)## check formula form <-as.formula(form)## fit model mod <-ssn_glm(formula = form,random = random,ssn.object = net,family = {{family}},tailup_type = tailup_type,taildown_type = taildown_type,euclid_type = euclid_type,nugget_type = nugget_type,additive = {{additive}},est_method = est_method)## extract estimates est_coef <-tidy(mod, conf.int =TRUE, effects ="fixed") |># regression coefs ("fixed effects")mutate(effect ="coef") est_var <-tidy(mod, conf.int =TRUE, effects ="ssn") |># variance componentsfilter(!is_known)if(!is.null(random)){ est_re <-tidy(mod, conf.int =TRUE, effects ="randcov") |># random effect variancemutate(effect ="re") est_raw <-bind_rows(est_coef, est_var, est_re) } else { est_raw <-bind_rows(est_coef, est_var) }# record estimation scale based on link (implied by likelihood) used scale <-case_when(family =="Gaussian"~"response", family %in%c("binomial", "beta") ~"logit", family %in%c("poisson", "nbinomial", "Gamma") ~"log")# combine estimates est <- est_raw |>mutate(scale = scale) |>select(effect, scale, term,est = estimate,se = std.error,lo = conf.low,hi = conf.high,pval = p.value)## extract fit fit <-glance(mod) |>bind_cols(loocv(mod)) |>mutate(method = est_method) |>select(n_obs = n,n_coef = p,n_varcomp = npar, method,aic = AIC,aicc = AICc,loglik = logLik,dev = deviance,dev_explain = pseudo.r.squared, bias,mse = MSPE,rmse = RMSPE,rav = RAV)## partition variance varpart <-varcomp(mod) |>mutate(effect ="varcomp",scale ="proportion") |>select(effect, scale,term = varcomp,est = proportion) |>mutate(term =ifelse(term =="Covariates (PR-sq)", "predictor", term))## output as dataframe out <-tibble_row(lab = lab, # model fit labelset = set, # model set labeloutcome = outcome, # outcome assessed labelresponse = response, # response variable representing outcome labelpredictor = predictor, # explanatory variable set labelerror = error, # covariance structure labelre = re, # random effect labelssn = path, # path to SSN used for modelformula =deparse1(form), # formula used in model as stringest =list(est), # model estimates as nested dataframefit =list(fit), # model fit metrics as nested dataframevarcomp =list(varpart) # model variance partitioning as nested dataframe )return(out)}
Define a formula for a random intercept by sampling event ID (as a categorical variable, nut numeric) to be reused across all models as the temporal random effect:
form_re <-formula(~1|as.factor(samp_event))
Univariable models
Univariable models are fit separately for each environmental factor as the explanatory variable to assess associations with each of the three outcome variables (generic E. coli log10 concentration, ESBL-Ec log10 concentration, and CST-Ec detection). While it would be possible to write a second wrapper function to execute the model fits and estimate extractions for each of the three outcomes, there is just enough idiosyncrasy across the covariates in naming convention and other nuances, and few enough separate covariates, to defend hard-coding the univariable model for each combination of outcome and covariate.
Generic E. coli
Start by estimating the univariable associations between the resistant outcomes (ESBL-Ec and CST-Ec) and generic E. coli. Apply the fit_ssn() function to the ESBL-Ec SSN assembled previously:
m_gen_esbl <-fit_ssn(lab ="gen", # label for generic Ec-as-predictor modelsset ="uni", # univariable model setoutcome ="esbl", # ESBL-Ec outcomeresponse ="mpn", # concentration (MPN) continuous response varpredictor ="mpn_gen", # generic EC concentration (MPN) explanatory varerror ="flow", # flow-connected (tail-up) covariance structurere ="event", # random effect by sampling eventpath ="../gis/ssn/esbl_gen.ssn", # path to ESBL SSN with complete generic Ec dataform =formula(mpn ~ mpn_gen), # model formularandom = form_re, # random effect formulafamily ="Gaussian", # linear model (identity link) (defaut)additive ="afv_area", # cumulative area drained as flow proxy (default)tailup_type ="exponential", # YES tailup covariance with exponential modeltaildown_type ="none", # NO taildown covariance (default)euclid_type ="none", # NO Euclidean covariance component (default)nugget_type ="nugget", # YES nugget (unstructured covariance) (default)est_method ="reml") # use restricted maximum likelihood (default)
Repeat for CST-Ec, specifying detection (rather than concentration) as the response, binomial family for logistic regression, and dropping default arguments for simplicity:
m_gen_cst <-fit_ssn(lab ="gen", # label for generic Ec-as-predictor modelsset ="uni", # univariable model setoutcome ="cst", # CST-Ec outcomeresponse ="detect", # detection binary response varpredictor ="mpn_gen", # generic EC concentration (MPN) explanatory varerror ="flow", # flow-connected (tail-up) covariance structurere ="event", # random effect by sampling eventpath ="../gis/ssn/cst_gen.ssn", # path to ESBL SSN with complete generic Ec dataform =formula(detect ~ mpn_gen), # model formulafamily ="binomial", # binomial likelihood (logistic link)random = form_re, # random effect formulatailup_type ="exponential") # YES tailup covariance with exponential model
Models using the co-measured water quality characteristics, turbidity and optical brighteners, also use the specific SSNs developed to account for missing covariate values, but are fit using generic E. coli as an outcome, as well as ESBL-Ec and CST-Ec.
Fit models with turbidity as the explanatory variable for all three outcomes:
Fit univariable models with basin features (land cover, imperviousness, and area) as the explanatory variable for each of the three outcomes. The full-observation SSNs are used because all basin feature covariate data are complete.
Although univariable models are not applicable for the interaction of septic density and rainfall (which also includes terms for the main effects of septic density and rainfall), fit a model for the interaction in the “univariable” set for consistency with the later multivariable estimates:
Multivariable models are also fit that include all environmental factors except for those measured in the water samples (i.e., excluding generic E. coli, turbidity, and optical brighteners) to explore the independent contribution of each external environmental factor conditional on the others. To account for potential subsurface transport of septic effluent mobilized by antecedent precipitation, include an interaction term in the multivariable models for the product of septic density and 3–7 day rainfall. Exclude the basin area and open land cover-class variables from the multivariable models to reduce collinearity among the basin-feature variable set.
Use the same fit_ssn() function to fit a multivariable model to the ESBL-Ec concentration data:
For convenience, a key mapping the shorthand variable names used in the analysis to their full display labels, along with grouping and ordering information, is provided in covar_labels.xlsx. Read in the labels from the key and the means and standard deviations previously used to scale the variables:
# covariate labelsvar_lab <-read_xlsx("../data/covar_labels.xlsx")# covariate scaling factorsvar_scl <-readRDS("../out/covar_summ_stats.rds") |>mutate(outcome =case_when(drug =="cefotaxime"~"esbl", drug =="colistin"~"cst", drug =="none"~"gen"))
Pull in the univariable model coefficient estimates, adding empty rows for the generic E. coli outcome with generic E. coli as the explanatory variable and setting the septic-rainfall interaction term to missing (and dropping the associated main effects) for display purposes when merged with the multivariable model estimates:
coef_uni <-readRDS("../est/ssn_est_uni.rds") |># empty row for generic Ec predicted by generic Ecadd_case(outcome ="gen",set ="uni",effect ="coef",term ="mpn_gen",predictor ="mpn_gen") |># set septic-rainfall interaction to NA and drop septic and rainfal main effects from interaction modelmutate(across(c(est, se, lo, hi, pval),~ifelse(term =="septic_dens:rain7", NA, .x))) |>filter(!(lab =="raindens"& term %in%c("septic_dens", "rain7")))
Obtain a list of unique effects estimated in univariable models, join the multivariable model estimates to match, stack with univariable estimates, and join to display labels and scaling factors:
coef_all <- coef_uni |># identify unique outcome-predictor combinations from univariable models select(outcome, response, effect, term) |>distinct() |>mutate(set ="multi") |># join matching multivariable model estimatesleft_join(readRDS("../est/ssn_est_multi.rds"),by =c("outcome", "response", "effect", "term", "set")) |># stack with univariable estimatesbind_rows(coef_uni) |># label with display names, grouping infoleft_join(var_lab, by =c("effect", "term")) |># join covariate scaling factorsleft_join(var_scl, by =c("outcome", "term"="var"))
Filter for just the main effect estimates (i.e., dropping intercepts and variance component hyperparameters), calculate descaled estimate values, exponentiate estimates from logistic models to obtain odds ratios (OR), and re-order for display:
coef_slct <- coef_all |># filter main effect estimates (identified in covar_label key)filter(effect =="coef", include =="yes") %>%mutate(# descale coef estimates to obtain on original covariate scaleest_descl =ifelse(scl_method =="none", est, est / sd),lo_descl =ifelse(scl_method =="none", lo, lo / sd),hi_descl =ifelse(scl_method =="none", hi, hi / sd),# exponentiate logistic model estimates for odds ratiosacross(c(est, lo, hi, ends_with("_descl")),~ifelse(scale =="logit", exp(.x), .x)),# order by covariate groupgroup =fct_relevel(group,c("wq", "out", "lc", "meteo")),set =fct_relevel(set, c("uni", "multi")) |>fct_recode("Univariable"="uni","Multivariable"="multi") |>fct_rev()) |>arrange(group, term_order) |>mutate(group_lab =fct_inorder(group_lab),term_lab =fct_inorder(term_lab),drug_lab =case_when(outcome =="esbl"~"ESBL-Ec<br>concentration", outcome =="cst"~"CST-Ec<br>detection", outcome =="gen"~"Generic *E. coli*<br>concentration"))
Display the univariable and multivariable coefficient estimates together as the point estimates and their 95% confidence intervals. Prepare a similar forest plot for each of the three E. coli outcomes and stack the plots together for display, beginning with the generic E. coli concentration outcome:
This reproduces Figure 3 from Holcomb et al. (2026). Note that the estimates are for the scaled covariates and thus correspond to the expected change in the outcome for a 1 standard deviation increase in the environmental factor:
Figure 3. Estimates (points) and 95% confidence intervals (lines) of the expected change in concentration of generic E. coli (left column) and extended-spectrum β-lactamase-producing E. coli (ESBL-Ec; center column) or odds of detecting colistin-resistant E. coli (CST-Ec; right column) associated with a standard deviation-increase in the value of environmental explanatory variables in univariable (circles, solid lines) and multivariable (diamonds, dotted lines) spatial generalized linear mixed models (GLMM).
Generate an alternative version of the forest plot with the descaled coefficient estimates that correspond to to the expected change in each outcome for a 1-unit increase in the environmental factor on its original measurement scale:
This alternative version of Figure 3 with coefficient estimates on the original scale of the covariate data makes comparisons between different environmental factors difficult, but does serve to highlight the large magnitude of the associations between resistant E. coli and generic E. coli concentrations (in units of log10 MPN/100 mL), as well as recent rainfall (in units of cm):
Figure 3-Alt. Estimates and 95% confidence intervals of the expected change in concentration of generic E. coli ESBL-Ec or odds of detecting CST-Ec associated with a one-unit increase in the value of environmental explanatory variables on the original scale of the environmental factor data.
Variance components
Examine the proportion of the total variance in the E. coli outcomes explained by each of the model structural components in the multivariable models, as well as the estimated spatial range (km) for the flow-connected spatial random effect.
Join the variable label key used for the coefficents to the previously extracted variance component estimates and filter for the variables of interest:
This essentially reproduces Table S9 from the Supplementary Material of Holcomb et al. (2026):
Variance explained (%)
Response
Nugget
Explanatory variables
Sampling event random effect
Flow-connected spatial random effect
Flow-connected spatial range (km)
Generic Ec concentration
42.9
6.1
9.0
42.0
9929
ESBL-Ec concentration
69.8
11.4
17.9
0.8
56
CST-Ec detection
0.2
31.2
68.4
0.2
33
Disclaimer:
The findings and conclusions of this report are those of the authors and do not necessarily represent the official position of the Centers for Disease Control and Prevention.
Source Code
---title: "Environmental Monitoring for Antimicrobial Resistance in a Metropolitan Watershed"date: 2026-08-28author: - name: David A Holcomb orcid: 0000-0003-4055-7164subtitle: Version 1.4csl: numbered-bracket.cslformat: html: theme: light: yeti fontsize: 0.9em page-layout: article toc: true toc-title: Contents toc-location: left toc-expand: 1 toc-depth: 3 self-contained: true fig-align: center fig-dpi: 300 fig-format: svg code-fold: false code-tools: true code-links: - text: Project repository icon: database href: https://github.com/CDCgov/EMEL/tree/main/manuscripts/river_ar - text: Water sample data icon: moisture href: https://cdcgov.github.io/EMEL/manuscripts/river_ar/sample_data.xlsx - text: Wastewater outfalls icon: nut href: https://cdcgov.github.io/EMEL/manuscripts/river_ar/outfall_info.xlsx - text: Drainage basins shapefile icon: layers href: https://cdcgov.github.io/EMEL/manuscripts/river_ar/basin.zip - text: Stream network shapefile icon: bezier2 href: https://cdcgov.github.io/EMEL/manuscripts/river_ar/stream.zip - text: Sequencing data icon: body-text href: http://www.ncbi.nlm.nih.gov/bioproject/1423195 - text: CDC Environmental Labs icon: droplet href: https://www.cdc.gov/healthy-water-data/laboratories/index.html cache: false---# DescriptionThis report provides **R** code to process the data and perform the analysis of antimicrobial-resistant *E. coli* in river water reported by:> Holcomb DA, Cantrell ME, Feistel DJ, Meyer M, Sterling J, Chen J, Krishnan K, Im S, Folster JP, Kim J, Burroughs M, Nobles S, Stahl A, Mattioli MC. (2026). Environmental monitoring for antimicrobial resistance in a metropolitan watershed. *Journal of Water and Health*. [https://doi.org/10.2166/wh.2026.219](https://doi.org/10.2166/wh.2026.219)This analysis code and the study data are available in the manuscript [repository](https://github.com/CDCgov/EMEL/tree/main/manuscripts/river_ar). The whole-genome sequencing data are available from NCBI under [BioProject PRJNA1423195](http://www.ncbi.nlm.nih.gov/bioproject/1423195) with accession numbers SRR55333857--SRR55334224.**Setup**Load packages and set options.```{r}#| label: setup#| message: false#| warning: false#| results: hide#| code-fold: true## Optionsoptions(tidyverse.quiet =TRUE)options(dplyr.summarise.inform =FALSE)options(knitr.kable.NA ='')## Packages# Utilitieslibrary(readxl) # I/Olibrary(writexl) # I/Olibrary(lubridate) # date formattinglibrary(prettyunits) # number formattinglibrary(zoo) # time series datalibrary(units) # handle units in dataframes# Spatiallibrary(sf) # spatial data & operationslibrary(nngeo) # spatial operationslibrary(rmapshaper) # spatial operationslibrary(terra) # raster data & operationslibrary(exactextractr) # summarize raster data within polygonslibrary(FedData) # download US federal data sourceslibrary(tigris) # download US census shapefileslibrary(SSNbler) # prepare spatial stream networkslibrary(SSN2) # fit spatial stream network models# Modelslibrary(broom) # extract model estimateslibrary(fitdistrplus) # fit distributions to datalibrary(truncnorm) # draw from truncated distributionlibrary(DescTools) # gives correlation CIs# Figureslibrary(patchwork) # combine plotslibrary(ragg) # more consistent graphics outputlibrary(ggtext) # plot label formattinglibrary(scales) # axis formatting# Tidyverselibrary(knitr) # table renderinglibrary(kableExtra) # table headerslibrary(tidyverse) # load last to avoid conflicts```# Data processing## Water samplesThe variables measured in surface water grab samples, including microbial concentrations, turbidity, and optical brightener fluorescence, along with sampling site information, are provided in the file [`sample_data.xlsx`](https://cdcgov.github.io/EMEL/manuscripts/river_ar/sample_data.xlsx). This file is structured with the data in the first sheet and a codebook in the second sheet. Each row in the data sheet corresponds to the measured *E. coli* concentration in a particular water sample from a single site and date, selected with a given antibiotic at a particular dose.Read in the sample data, fix the dates, and save as RDS for later use:```{r}#| label: sample_data# read sample datadf_samp <-read_xlsx("../data/sample_data.xlsx", sheet ="data") |>mutate(dt_collect =as_date(dt_collect))# save in RDS format for stabilitysaveRDS(df_samp, "../out/sample_data.rds")```## Sampling site locationsInformation about each of the sampling sites is contained in the waster sample data set. Extract the site information from the sample data, condense, and save for later:```{r}#| label: site_infodf_site <- df_samp |>select(starts_with("site_")) |>distinct()saveRDS(df_site, "../out/site_data.rds")```Extract the *coordinate reference system* (CRS) from the site data and use it to convert the lat/lon coordinates for each of the sites into `sf` format to facilitate future spatial operations. Project the site locations into the *NAD83(2011) / Georgia West (m)* projection (`EPSG:6446`) and save for later use:```{r}#| label: site_sf# extract existing CRScrs_site <- df_site |>select(site_crs) |>distinct() |>pull()# transform to simple featuressf_site <- df_site |>st_as_sf(coords =c("site_lon", "site_lat"),crs = crs_site) |>st_transform(crs =6446) |># project as NAD83(2011) / Georgia West (m)mutate(epsg ="EPSG:6446",proj ="NAD83(2011) / Georgia West")saveRDS(sf_site, "../gis/sf/site_pt.rds")```## Drainage basin boundariesPolygons representing the perimeter of the drainage basin contributing to each sampling location were obtained using the [USGS StreamStats Batch Processor](https://streamstats.usgs.gov/) v1.5.1. Sample site locations have already been manually adjusted in to coincide with the nearest feasible [USGS StreamGrid](https://streamstats.usgs.gov/streamgrids/ga/ga.zip) raster cell prior to basin delineation. In addition to the sampling sites, a drainage basin was delineated for all areas upstream of the USGS monitoring station at Buford Dam ([USGS02334430](https://waterdata.usgs.gov/monitoring-location/USGS-02334430/)) to encompass Lake Lanier and all the Upper Chattahoochee watershed. This basin upstream of Buford Dam is used to truncate the drainage basins for all the sampling sites to restrict them to downstream of Lake Lanier, as all the river mainstem sites would otherwise share fully overlapping drainage basins upstream of the dam.The raw drainage basin polygons output by USGS StreamStats are provided as the shapefile `basin_raw.shp` in the compressed folder [`basin.zip`](https://cdcgov.github.io/EMEL/manuscripts/river_ar/basin.zip). After downloading and extracting the compressed file, read in the raw basins shapefile and perform spatial quality control operations, including filling holes in the polygons and validating the resulting geometries:```{r}#| label: basin_shp# raw basin polygons from streamstats shapefilebasin_raw <-read_sf("../gis/shape/basin/basin_raw.shp")# join to site info, fill holes, validatebasin_full <-readRDS("../out/site_data.rds") |># start with site info# add a row corresponding to the Buford Dam basin, upstream of all the other sites.add_case(site_num =1, site_id ="BUD", site_lab ="Buford Dam", site_type ="lake") |>left_join(basin_raw, by =c("site_id"="Name")) |># join the basin polygonsselect(site_num, site_id, site_lab, site_type, geometry) |>st_as_sf(crs =st_crs(basin_raw)) |># convert back to sf format using original crsst_transform(crs =6446) |># project as NAD83(2011) / Georgia Westst_remove_holes() |># remove holes from the polygons st_make_valid() |># validate geometry arrange(site_num) |># order by site number (roughly upstream -> downstream)mutate(site_id =fct_inorder(site_id),epsg ="EPSG:6446",proj ="NAD83(2011) / Georgia West")saveRDS(basin_full, "../gis/sf/basin_full.rds")```Separate the Buford Dam drainage basin (upstream-most) from the rest:```{r}#| label: basin_sep# upper basin as geometry so the attributes don't conflictbasin_upper <- basin_full |>filter(site_id =="BUD")# everything elsebasin_lower <- basin_full |>filter(site_id !="BUD")```Use the upper basin to clip the rest of the basins so they are all truncated at Buford Dam, then re-validate the geometries:```{r}#| label: basin_clipbasin_clip <-ms_erase(basin_lower, basin_upper,remove_slivers =TRUE) |>st_make_valid()```Re-append the upper basin polygon to the other truncated ones, repeat geometry quality control measures, and save:```{r}#| label: basin_appendbasin <- basin_full |>filter(site_id =="BUD") |>bind_rows(basin_clip) |>st_remove_holes() |>st_make_valid() |>arrange(site_num) |>mutate(site_id =fct_inorder(site_id))saveRDS(basin, "../gis/sf/basin_poly.rds")```Take a quick look at the basin outlines:```{r}#| label: basin_peekplot(st_geometry(basin))```Merge (union) all the truncated basins below the dam to obtain the perimeter (scope) of the study area, repeating geometry quality control measures:```{r}#| label: basin_scopebasin_scope <- basin |>filter(site_num >2) |>st_union() |>st_remove_holes() |>st_make_valid()saveRDS(basin_scope, "../gis/sf/basin_scope.rds")```Plot the scope:```{r}#| label: basin_scope_peekplot(st_geometry(basin_scope))```## Land coverLand cover data were obtained from the USGS National Land Cover Database (NLCD) 2019 vintage via the `FedData` package. Both land cover classes and percent impervious surface data were obtained as separate raster files. For mapping purposes, the full study area scope (including the shared watershed upstream of Buford Dam) was used to select the raster coverage area.Merge all the basins, including the basin above the dam, to obtain the extended scope for which to download land cover rasters:```{r}#| label: basin_scope_full#| eval: falsebasin_scope_full <-readRDS("../gis/sf/basin_poly.rds") |>st_union() |>st_remove_holes() |>st_make_valid()```Add a small (2 km) buffer to the scope so the retrieved rasters fully cover the study area extent and transform to WGS84 to work with the `FedData` package:```{r}#| label: basin_scope_buff#| eval: falsebasin_scope_buff <- basin_scope_full |>st_buffer(dist =2000) |>st_transform(crs =4326) # transform to WGS84 to work with FedData```Download NLCD LULC 2019 and NLCD Impervious 2019 as separate rasters and save as GeoTIFF format:```{r}#| label: nlcd_down#| eval: false# LULCnlcd_lulc <-get_nlcd(template = basin_scope_buff,label ="river_ar",year =2019,dataset ="landcover",landmass ="L48",extraction.dir ="../gis/lulc/",force.redo =TRUE )writeRaster(nlcd_lulc, "../gis/lulc/nlcd_2019_lulc.tif",overwrite =TRUE)# imperviousnlcd_impv <-get_nlcd(template = basin_scope_buff,label ="river_ar",year =2019,dataset ="impervious",landmass ="L48",extraction.dir ="../gis/lulc/",force.redo =TRUE )writeRaster(nlcd_impv, "../gis/lulc/nlcd_2019_impv.tif",overwrite =TRUE)```Because re-projecting rasters requires a non-reversible transformation of the data, extract the existing CRS information from the downloaded rasters so it can be used to re-project the vector data for analyses:```{r}#| label: nlcd_crs# reload from saved TIFFnlcd_lulc <-rast("../gis/lulc/nlcd_2019_lulc.tif")# extract CRS as sfnlcd_crs <-st_crs(nlcd_lulc)saveRDS(nlcd_crs, "../gis/lulc/nlcd_crs_sf.rds")```## MeteorologyMeteorological data are obtained from the NOAA GHCN-Daily records using the `FedData` package. The following code demonstrates how the data were originally obtained, but because the data available from NOAA may change over time, re-executing this code could return a slightly different data set than was used for the canonical analysis reported in the manuscript.Set a search zone for weather stations within 20 km of the study area (defined as the truncated drainage basins), then transform to WGS84 to work with the `FedData` package:```{r}#| label: ghcn_template#| eval: falsetemplate_fd <-readRDS("../gis/sf/basin_scope.rds") |>st_buffer(dist =20000) |># search within 20 km buffer of study areast_transform(crs =4326) # transform to WGS84 to work with FedData```Different stations report different weather variables, so download the GHCN data separately for precipitation and temperature to ensure broadest swath of stations are available for each for both full years of the study period (2018 & 2019):```{r}#| label: ghcn_down#| message: false#| warning: false#| eval: false# stations with precipghcn_prcp <-get_ghcn_daily(template = template_fd,label ="river_ar",years =2018:2019,elements =c("prcp") )# stations with tempghcn_temp <-get_ghcn_daily(template = template_fd,label ="river_ar",years =2018:2019,elements =c("tmin", "tmax", "tavg") )``````{r}#| label: ghcn_save#| eval: false#| include: falsesaveRDS(ghcn_prcp, "../out/ghcn_prcp.rds")saveRDS(ghcn_temp, "../out/ghcn_temp.rds")``````{r}#| label: ghcn_load#| echo: false#| message: false#| warning: falseghcn_prcp <-readRDS("../out/ghcn_prcp.rds")ghcn_temp <-readRDS("../out/ghcn_temp.rds")```Extract station location information:```{r}#| label: ghcn_locghcn_loc <-bind_rows(ghcn_prcp$spatial, ghcn_temp$spatial) |>distinct() |>rename(ghcn_id = ID,ghcn_name = NAME) |>st_transform(crs =6446) |># project as NAD83(2011) / Georgia Westmutate(epsg ="EPSG:6446",proj ="NAD83(2011) / Georgia West")saveRDS(ghcn_loc, "../gis/sf/meteo_pt.rds")```Write a utility function to parse the GHCN tabluar data and wrangle into long format:```{r}#| label: ghcn_funcextract_ghcn <-function(x){tibble(val = x) |>unnest_longer(col = val,indices_to ="meteo_type") |>unnest(col = val) |>pivot_longer(D1:D31,names_to ="day",values_to ="meteo_val") |>mutate(day =str_remove(day, "D"),dt_meteo =ymd(str_c(YEAR, MONTH, day, sep ="-")),meteo_val = meteo_val /10) |># values for 1/10 of unit, divide by 10 to get desired unitsfilter(!is.na(dt_meteo)) |>select(ghcn_id = STATION, dt_meteo, meteo_type, meteo_val)}```Apply the utility function to the data downloads, combine the resulting long-format dataframes, and save:```{r}#| label: ghcn_parse#| warning: false# precipitation datadf_ghcn_prcp <- ghcn_prcp$tabular %>%extract_ghcn() %>%mutate(meteo_unit ="mm")# temperature datadf_ghcn_temp <- ghcn_temp$tabular %>%extract_ghcn() %>%mutate(meteo_unit ="degC")# combinedf_ghcn <-bind_rows(df_ghcn_prcp, df_ghcn_temp) %>%mutate(meteo_type =str_to_lower(meteo_type))saveRDS(df_ghcn, "../out/ghcn_data.rds")```## InfrastructureWastewater infrastructure data were obtained from two sources:- Wastewater treatment plant (WWTP) treated effluent outfall and combined sewer overflow (CSO) locations and permitted maximum discharge (million gallons per day, MGD) from [Georgia Environmental Protection Division online system](https://geos.epd.georgia.gov/ga/geos/public/govent/shared/pages/main/login.aspx), retrieve and collated by [Chattahoochee Riverkeeper](https://chattahoochee.org/)- Septic system locations from the [Metropolitan North Georgia Water Planning District](https://northgeorgiawater.org/) for 15 counties covering the study areaRequests for the raw infrastructure data should be directed to the respective organizations that provided them for the study.### Wastewater outfallsInformation on the outfalls in the study area is provided in the file [`outfall_info.xlsx'](https://cdcgov.github.io/EMEL/manuscripts/river_ar/outfall_info.xlsx). Read the outfall data:```{r}#| label: outfall_datadf_outfall <-read_xlsx("../data/outfall_info.xlsx", sheet ="data")```Convert coordinates to `sf` format and project for spatial operations:```{r}#| label: outfall_sfsf_outfall <- df_outfall |>st_as_sf(coords =c("lon", "lat"),crs =4269) |>st_transform(crs =6446) |># project as NAD83(2011) / Georgia West (m)mutate(epsg ="EPSG:6446",proj ="NAD83(2011) / Georgia West")```Perform a spatial join with the drainage basin polygons to associate the outfalls with the basins to which they contribute:```{r}#| label: outfall_join#| warning: false# reload basins and filter for study areabasin_out <-readRDS("../gis/sf/basin_poly.rds") |>filter(site_num >2)# spatial joinoutfall_loc <-st_intersection(sf_outfall, basin_out) |>select(starts_with("site"),starts_with("facil"),starts_with("out"), epsg, proj, geometry)saveRDS(outfall_loc, "../gis/sf/outfall_pt.rds")```### Septic systemsSeptic system locations are not shared for privacy reasons. Contact [Metropolitan North Georgia Water Planning District](https://northgeorgiawater.org/) to request access to septic system data.The data cover septic system locations for 15 counties in the metropolitan Atlanta, GA area. These 15 counties fully encompass the study area, with the exception of a small portion of the Fairburn Rd drainage basin (downstream-most site), which extends into Carroll County, not covered in the septic system data. In order to associate septic systems with drainage basins, the Carroll County portion of the Fairburn Rd basin must be masked and the truncated basin area used in density calculations.Download county boundary polygons from the US Census using the `tigris` package, project, and reformat:```{r}#| label: county_down#| eval: false# download Georgia Countiescnty_poly_raw <- tigris::counties(year =2020,state ="GA",cb =TRUE)# project and reformatcnty_poly <- cnty_poly_raw |>st_transform(crs =6446) %>%# project as NAD83(2011) / Georgia Westmutate(cnty_area =st_area(.) |>set_units(km^2),cnty_perim =st_perimeter(.) |>set_units(km),GEOID =as.character(GEOID),epsg ="EPSG:6446",proj ="NAD83(2011) / Georgia West") |>st_make_valid() |>select(cnty_id = GEOID,state_fips = STATEFP,state_id = STUSPS,state_name = STATE_NAME,cnty_fips = COUNTYFP,cnty_name = NAME, cnty_area, cnty_perim,cnty_land = ALAND,cnty_water = AWATER, epsg, proj, geometry) |>drop_units() |>select(cnty_id,starts_with("state"),starts_with("cnty"), epsg, proj, geometry) |>arrange(cnty_id)saveRDS(cnty_poly, "../gis/sf/cnty_poly_2020.rds")```Pull out Carroll County and validate geometry:```{r}#| label: county_carrollcnty_carroll <-readRDS("../gis/sf/cnty_poly_2020.rds") |>filter(cnty_name =="Carroll") |>st_remove_holes() |>st_make_valid()```Re-load the basin polygons, restrict to downstream of the dam, and calculate area:```{r}#| label: basin_septicbasin_septic <-readRDS("../gis/sf/basin_poly.rds") |>arrange(site_num) |>filter(site_num >2) %>%mutate(wshed_area =st_area(.) |> units::set_units("km^2"),wshed_perim =st_perimeter(.) |> units::set_units("km")) |>drop_units()```Clip the basins to exclude Carroll County, validate geometry, and re-calculate area:```{r}#| label: basin_septic_clipbasin_clip <-ms_erase(basin_septic, cnty_carroll,remove_slivers =TRUE) |>st_remove_holes() |>st_make_valid() %>%mutate(wshed_area_clp =st_area(.) |> units::set_units("km^2")) |>drop_units() %>%mutate(area_diff = wshed_area - wshed_area_clp,pct_reduction = (area_diff / wshed_area) *100)```Read and project septic system locations, then perform a spatial join with the basin polygons and save:```{r}#| label: septic_join#| eval: false#| warning: false# read and project septic locationssf_septic <-readRDS("../gis/sf/septic_loc.rds") |>st_transform(crs =6446) |># project as NAD83(2011) / Georgia Westmutate(epsg ="EPSG:6446",proj ="NAD83(2011) / Georgia West")# spatial join to basinsseptic_loc <-st_intersection(sf_septic, basin_clip) |>select(septic_id,starts_with("site"),starts_with("wshed"),starts_with("dt"), source,septic_lat = lat, septic_lon = lon, epsg, proj, geometry)saveRDS(septic_loc, "../gis/sf/septic_pt.rds")```# Variable generation```{r}#| label: refresh_vars#| include: falserm(list =ls())```The variables used to represent *E. coli* outcomes and environmental factors in the analyses were created from the data sets obtained above through additional processing.## *E. coli* concentration: imputation of censored observationsRaw concentration data for both generic and resistant *E. coli* include values that are problematic for statistical analysis of microbial concentrations. First, the concentration of the target in samples where it was not detected was not necessarily zero, just less than the concentration that corresponds to the lower limit of detection (LOD). Because zeros are also problematic for logarithmic transformation, which is typically appropriate for analyzing microbial concentration data that scale geometrically rather than arithmetically, it is useful to consider non-detects as censored below the LOD. Similarly, samples in which the target is so abundant that all the wells in the most probable number (MPN) test are positive (i.e., TNTC: "too numerous to count") are only known to have a concentration that is greater than that corresponding to the upper limit of quantification (ULOQ) of the test.Censoring below the LOD and above the ULOQ correspond to left and right censoring, respectively, and can be addressed using imputation techniques to substitute informed values for the left- and right-censored observations to facilitate log-transformation and downstream analysis. Imputation is performed separately for each target (i.e., each for each selective antibiotic and dose) by fitting a probability distribution to the censored data and drawing a random value for each censored observation from the fitted distribution truncated below the LOD or above the ULOQ for left- and right-censored observations, respectively.Prepare a data set of microbial concentrations for censored imputation processing. Because few meropenem-nonsusceptible *E. coli* (MEM-Ec) were detected and a relatively small number of samples were selected by cefotaxime at 4 μg/mL, restrict the data set to generic *E. coli*, cefotaxime-selected *E. coli* at 1 μg/mL (ESBL-Ec), and colistin-selected *E. coli* (CST-Ec). All concentrations will be log~10~-transformed for this analysis; for this reason, non-detect and TNTC observations have already been replaced with their respective LOD or ULOQ values, respectively, along with an indicator variables denoting whether they are left-, right-, or non-censored. These substituted values are simply for convenience and will be treated as missing in the imputation procedure.For compatibility with the function from the `fitdistrplus` package for fitting censored distributions to these data, the concentration variables (`mpn`, `mpn_lod`, and `mpn_uloq`) along with this censored indicator variable (`cens`) need to be converted into a pair of variables named `left` and `right` that both contain the observed concentration when not censored and contain a missing value in corresponding column and the censoring value in the opposite column for censored observations. I.e., a left-censored observation should be missing in `left` and contain the LOD value in `right`, while a right-censored observation will contain the ULOQ value in `left` and be missing in `right`:```{r}#| label: micro_censdf_micro_cens <-readRDS("../out/sample_data.rds") |>filter(drug %in%c("cefotaxime", "colistin", "none"), drug_dose !=4, site_num >2,!is.na(mpn)) |>group_by(drug, drug_dose) |>mutate(nobs =n(),npos =sum(cens =="no"),ncens =sum(cens =="left"| cens =="right")) |>ungroup() |>mutate(# log10-transform the concentration and limitsacross(c(mpn, mpn_lod, mpn_uloq), ~log10(.x)),# create left-censored values variableleft =case_when(cens =="left"~NA_real_, cens =="no"~ mpn, cens =="right"~ mpn_uloq),# create right-censored value variableright =case_when(cens =="right"~NA_real_, cens =="no"~ mpn, cens =="left"~ mpn_lod))```Extract the number of censored observations by outcome:```{r}#| label: ncensncens <- df_micro_cens |>select(drug, ncens) |>distinct()```Fit censored normal distributions to the concentrations separately for each outcome, starting with cefotaxime.First, filter the data set, extract just the `left` and `right` variables, and convert to basic dataframe for compatibility:```{r}#| label: cens_filter_ctxdf_cens_ctx <- df_micro_cens |>filter(drug =="cefotaxime") |>select(left, right) |>as.data.frame()```Fit the censored normal distribution and extracte the estimated mean and standard deviation with their standard errors:```{r}#| label: cens_fit_ctx# fit distributionfit_cens_ctx <-fitdistcens(censdata = df_cens_ctx,distr ="norm")# extract parameter estimatesest_cens_ctx <-tibble(drug ="cefotaxime",nobs = fit_cens_ctx$n, # number of observations used to fit the distributionmu = fit_cens_ctx$estimate[1], # estimated mean of the normal distributionmu_se = fit_cens_ctx$sd[1], # standard error of the mean estimatesigma = fit_cens_ctx$estimate[2], # estimated standard deviation of the normal distributionsigma_se = fit_cens_ctx$sd[2]) # standard error of the SD estimate```Repeat for colistin-selected and generic *E. coli*:```{r}#| label: cens_fit#| code-fold: true## Colistin-selected (CST)# dataframedf_cens_cst <- df_micro_cens |>filter(drug =="colistin") |>select(left, right) |>as.data.frame()# fit distributionfit_cens_cst <-fitdistcens(censdata = df_cens_cst,distr ="norm")# extract parameter estimatesest_cens_cst <-tibble(drug ="colistin",nobs = fit_cens_cst$n,mu = fit_cens_cst$estimate[1],mu_se = fit_cens_cst$sd[1],sigma = fit_cens_cst$estimate[2],sigma_se = fit_cens_cst$sd[2])## Generic# dataframedf_cens_gen <- df_micro_cens |>filter(drug =="none") |>select(left, right) |>as.data.frame()# fit distributionfit_cens_gen <-fitdistcens(censdata = df_cens_gen,distr ="norm")# extract parameter estimatesest_cens_gen <-tibble(drug ="none",nobs = fit_cens_gen$n,mu = fit_cens_gen$estimate[1],mu_se = fit_cens_gen$sd[1],sigma = fit_cens_gen$estimate[2],sigma_se = fit_cens_gen$sd[2])```Combine the estimates, join the counts of censored observations, and save:```{r}#| label: cens_est_comboest_cens <-bind_rows(est_cens_ctx, est_cens_cst, est_cens_gen) |>left_join(ncens, by ="drug") |>select(drug, nobs, ncens, mu, mu_se, sigma, sigma_se)saveRDS(est_cens, "../est/distr_cens_mpn.rds")```The MLE means and standard deviations from the fitted censored normal distributions are used to parameter truncated normal distributions, from which suitable values are drawn to represent the censored observations. A random draw for each censored observation is used for analysis, but the expected value of the distribution truncated at the sample-specific LOD/ULOQ is also calculated for each censored observation for additional context.Set a random seed for reproducibility, join the estimates to the censored observation data set, prepare the estimates for drawing from a truncated normal distribution. Define the lower bound `a` and the upper bound `b` for each sample to correspond with the censoring status of that observation. Imputes values for every observation (including the non-censored ones, which are then drawn from a normal distribution without truncation), retaining the imputed values only for censored observations:```{r}#| label: imp_draw# seed for reproducibilityset.seed(30333)# join estimates to data and imputeimp_cens <- df_micro_cens |>select(starts_with("site"),starts_with("drug"), wk_collect, dt_collect, samp_event, samp_id, test_id, lab_vol, uid, uid_gen, uid_wq, wq_turb, wq_optic, target, detect, mpn, cens, mpn_lod, mpn_uloq, left, right) |># merge distribution estimates left_join(est_cens, by ="drug") |># add bounds variablesmutate(# Define lower bound variablea =case_when(cens =="left"~-Inf, cens =="no"~-Inf, cens =="right"~ mpn_uloq),# Define upper bound variableb =case_when(cens =="left"~ mpn_lod, cens =="no"~Inf, cens =="right"~Inf)) |># perform imputation and only keep censoredmutate(imp_e =etruncnorm(a, b, mu, sigma), # expected value of truncated normal distimp_r =rtruncnorm(1, a, b, mu, sigma), # random draw from truncated normal distmpn_imp_e =ifelse(cens =="no", mpn, imp_e), # only keep expected value for censored obsmpn_imp_r =ifelse(cens =="no", mpn, imp_r)) # only keep random draw for censored obs```Separate out the generic *E. coli* values:```{r}#| label: imp_genimp_cens_gen <- imp_cens |>filter(drug =="none") |>select(uid,cens_gen = cens,mpn_lod_gen = mpn_lod,mpn_uloq_gen = mpn_uloq,mpn_e_gen = mpn_imp_e,mpn_r_gen = mpn_imp_r)```Clean up the data set and join the imputed values for generic *E. coli* to their corresponding resistance-selected (and non-selected [generic], for consistency) *E. coli* observations:```{r}#| label: imp_cleandf_imp <- imp_cens |>select(starts_with("site"),starts_with("drug"), wk_collect, dt_collect, samp_event, samp_id, test_id, lab_vol, uid, uid_gen, uid_wq, wq_turb, wq_optic, target, detect, cens, mpn_lod, mpn_uloq,mpn_e = mpn_imp_e,mpn_r = mpn_imp_r) |>left_join(imp_cens_gen,by =c("uid_gen"="uid"))saveRDS(df_imp, "../out/sample_imputed.rds")```## Basin land cover: distance-weighted, areal proportionVariables representing the influence of land cover classes on each sampling location are generated by weighting each raster cell of a given class contained within a site's drainage basin by the inverse of the squared Euclidian distance between that cell's centroid and the sampling location. Summing these cells by class provides inverse distance squared weighted (IDSW) cell counts, which are divided by the total IDSW cell count across all classes to obtain the IDSW proportion of the basin represented by that land cover class. Similarly, the IDSW mean imperviousness is obtained as a weighted average of the percent imperviousness of each cell in the drainage basin.### Calculate distance and weightsRead in the site locations and basin polygons and reproject to match the NLCD projection extracted previously:```{r}#| label: nlcd_reproj# CRSnlcd_crs <-readRDS("../gis/lulc/nlcd_crs_sf.rds")# basinbasin_nlcd <-readRDS("../gis/sf/basin_poly.rds") |>filter(site_num >2) |>st_transform(nlcd_crs) |>arrange(site_num, site_id)# sitessite_nlcd <-readRDS("../gis/sf/site_pt.rds") |>st_transform(nlcd_crs) |>arrange(site_num, site_id)```Write a function to calculate the distance between a sampling site and the centroid of each raster cell. The function should take as inputs a raster with the data to be extracted and an `sf` dataframe with the point locations to calculate distances from. The output will be a raster stack with a layer for each sampling location where the raster value in each cell records its distance from the sampling point. Include functionality to encode the resulting rasters with a weight based on the inverse of the distance (raised to an arbitrary power) rather than the raw distance itself.```{r}#| label: f_dist_rastsite_dist <-function(x_rast, # raster to calculate distances for y_sf, # sf tibble with point locations to measure distances fromidw =FALSE, # return inverse distance weights instead of the distance?power =1# distance raised to power for calculating weights ){ n_feature <-nrow(y_sf) # number of sites to calculate distances to# loop through each site dist_list <-list() # Prepare empty list object for(i in1:n_feature){ y <- y_sf[i, ] # select site x <-distance(x_rast, y) # calculate distances# apply inverse distance weightingif(idw){ x_idw <-1/ (x^power) # raise distances to specified power dist_list <-c(dist_list, x_idw) # insert weights into list object } else { dist_list <-c(dist_list, x) # insert distances into list object } dist_stack <-rast(dist_list) # convert list to raster stack }return(dist_stack)}```Apply function to calculate inverse distance squared weights by specifying `idw = TRUE` and `power = 2`. Because the LULC and impervious surface rasters contain identical cell locations and extents, and differ only in the values recorded for each cell, the distance function need only be applied once to obtain the raster stack of weights that can be used for both:```{r}#| label: nlcd_dist_wt#| eval: false# read in LULC raster (impervious would also work)nlcd_lulc <-rast("../gis/lulc/nlcd_2019_lulc.tif")# calculate inverse distance weightsnlcd_idsw <-site_dist(nlcd_lulc, site_nlcd,idw =TRUE, power =2)# name the raster layers by corresponding site idnames(nlcd_idsw) <- site_nlcd$site_id# save to drive as GeoTIFFwriteRaster(nlcd_idsw, "../gis/lulc/nlcd_site_idsw.tif",overwrite =TRUE)```### Extract weighted land cover class proportionsWrite a function for use with the `exactextractr` package to calculate the weighted counts of each class for a given basin, as well as basin-wide statistics. The `exactextractr::exact_extract()` function calculates fractional coverage for cells partially covered by the basin polygon and returns a dataframe with a row-per-cell with computed values given pre-specified names (e.g., `coverage_fraction`) that are accessed directly in the function before summarizing land cover class:```{r}#| label: f_extract_idwextract_idw <-function(df){ df |># calculate cell-level characteristicsmutate(cover_tot =sum(coverage_fraction), # total cells in basin, accounting for partial coveragefrac = coverage_fraction / cover_tot, # proportion of total in this cellcover_wt = coverage_fraction * weight, # cell's weighted coverage fractioncover_wt_tot =sum(coverage_fraction * weight), # total weighted cells in basinfrac_wt = cover_wt / cover_wt_tot, # proportion of weighted total in this cellwt_tot =sum(weight)) |># sum of all the weights across the basin# specify variables to retain after summarizing# only one site at a time so site and tot variables are static for all cells# `value` contains the class and is the only variable that separate summaries are calculated bygroup_by(site_num, site_id, value, cover_tot, cover_wt_tot, wt_tot) |># tally up cellssummarize(freq =sum(frac), # raw class proportion freq_wt =sum(frac_wt), # weighted class proportionwt =sum(weight), # sum of weights within classcount =sum(coverage_fraction), # raw class cell count in basincount_wt =sum(cover_wt), # weighted class cell count in basinarea =sum(area), # total area covered by classarea_wt =sum(area * weight)) # weighted area covered by class}```Write a second function to loop through each basin and its corresponding IDSW raster in the stack to apply the extraction function to each basin in turn. Include an option :```{r}#| label: f_loop_idwloop_idw <-function(x_rast, # raster with the values to extract y_poly, # sf with basin polygons (one per row) wt_rast, # raster stack with weights (must perfectly align with x_rast) include_cols =NULL# character vector of column names to retain in output ){ n_feature <-nrow(y_poly) # number of basins extract_df <-tibble() # prepare empty dataframe# loop through each basinfor(i in1:n_feature){ y <- y_poly[i, ] # select basin wt <- wt_rast[[i]] # select weights raster x <-exact_extract(x_rast, y,fun = extract_idw, # custom weighting functionweights = wt,summarize_df =TRUE, # produce summarized output dataframeinclude_cols = include_cols, # column names to retain (as character vector)include_area =TRUE, # calculate area within basinprogress =TRUE) # show progress indicator extract_df <-bind_rows(extract_df, x) # add basin summary to dataframe (inefficient) }return(extract_df)}```Apply the functions to the NLCD LULC raster to calculate IDSW land cover class proportions:```{r}#| label: lulc_idsw#| eval: falsebasin_lulc_idsw <-loop_idw(x_rast = nlcd_lulc,y_poly = basin_nlcd,wt_rast = nlcd_idsw,include_cols =c("site_num", "site_id"))saveRDS(basin_lulc_idsw, "../out/basin_lulc_idsw.rds")```### Extract weighted mean imperviousness Each cell in the land cover raster contains a single, discrete value corresponding to the land cover class assigned to that cell (cell value is a categorical variable). By contrast, the imperviousness raster contains an integer value between 0 and 100 that represents the proportion of that cell that is covered by impervious surfaces. Rather than a weighted sum of the cells assigned to each class, a weighted mean imperviousness across all cells in the basin in calculated.Weighted mean operations are built into `exactextractr::exact_extract()`, so revise the looping function to apply the built-in statistical functions (which require slightly different syntax in the arguments) to each basin in turn:```{r}#| label: f_loop_idw_impvloop_idw_impv <-function(x_rast, # raster with the values to extract y_poly, # sf with basin polygons (one per row) wt_rast, # raster stack with weights (must perfectly align with x_rast) include_cols =NULL# character vector of column names to retain in output ){ n_feature <-nrow(y_poly) # number of basins extract_df <-tibble() # prepare empty dataframe# loop through each basinfor(i in1:n_feature){ y <- y_poly[i, ] # select basin wt <- wt_rast[[i]] # select weights raster x <-exact_extract(x_rast, y,fun =c("count","mean", "weighted_mean","stdev", "weighted_stdev"),weights = wt,append_cols = include_cols,progress =TRUE) extract_df <-bind_rows(extract_df, x) # add basin summary to dataframe }return(extract_df)}```Apply the functions to the NLCD LULC raster to calculate IDSW land cover class proportions:```{r}#| label: impv_idsw#| eval: false# impervious rasternlcd_impv <-rast("../gis/lulc/nlcd_2019_impv.tif")# reload weights raster stacknlcd_idsw <-rast("../gis/lulc/nlcd_site_idsw.tif")# apply the loopbasin_impv_idsw <-loop_idw_impv(x_rast = nlcd_impv,y_poly = basin_nlcd,wt_rast = nlcd_idsw,include_cols =c("site_num", "site_id"))saveRDS(basin_impv_idsw, "../out/basin_impv_idsw.rds")```### Condense The NLCD land cover classes are narrowly defined and are collapsed into broader categories for analysis:- Developed: 22, 23, 24- Forest: 41, 42, 43, 90- Open: a combined category encompassing - Agriculture: 81, 82 - Barren: 31 - Grassland: 52, 71, 72, 95 - Open developed: 21For convenience, a key is provided at [`nlcd_lulc_key.xlsx`](https://cdcgov.github.io/EMEL/manuscripts/river_ar/nlcd_lulc_key.xlsx) to map the numerical class values contained in the NLCD LULC raster to the condensed categories used in the analysis. Read in the key, join to the IDSW proportion estimates, and aggregate the class proportions up to the broader categories: ```{r}#| label: nlcd_collapse# land use classesnlcd_class <-read_xlsx("../data/nlcd_lulc_key.xlsx")# join and aggregatebasin_lulc_long <-readRDS("../out/basin_lulc_idsw.rds") |>left_join(nlcd_class, by =c("value")) |>mutate(method ="idsw") |>group_by(site_num, site_id, cat, method) |>summarise(freq =sum(freq),freq_wt =sum(freq_wt),cells =sum(count)) |>group_by(site_num, site_id, method) |>mutate(freq_tot =sum(freq),freq_wt_tot =sum(freq_wt),cells_tot =sum(cells)) |>ungroup()```Pivot the aggregated LULC category values wider to obtain separate variables for each class:```{r}#| label: basin_lulc_widebasin_lulc <- basin_lulc_long |>select(site_num, site_id, method, cat, freq_wt) |>pivot_wider(names_from = cat,values_from = freq_wt,names_prefix ="lc_") |>mutate(across(starts_with("lc_"), ~replace_na(.x, 0))) |>rowwise() |>mutate(lc_tot =sum(c_across(starts_with("lc_")))) |>ungroup()```Trim the weighted mean impervious values:```{r}#| label: basin_impv_trimbasin_impv <-readRDS("../out/basin_impv_idsw.rds") |>mutate(method ="idsw") |>select(site_num, site_id, method, wshed_cells = count,lc_impv = weighted_mean,lc_impv_sd = weighted_stdev)```Calculate the basin area in km^2^, merge land cover variables, and save:```{r}#| label: nlcd_basin_areadf_basin_feat <-readRDS("../gis/sf/basin_poly.rds") |>arrange(site_num) |>filter(site_num >2) %>%mutate(wshed_area =st_area(.) |> units::set_units("km^2"),wshed_perim =st_perimeter(.) |> units::set_units("km")) |>drop_units() |>st_drop_geometry() |>left_join(basin_impv, by =c("site_num", "site_id")) |>left_join(basin_lulc,by =c("site_num", "site_id", "method"))saveRDS(df_basin_feat, "../out/basin_features.rds")```## Daily weather: distance-weighted, rolling aggregateWeather variables (precipitation in mm and maximum temperature in °C) are computed for each sampling location on each date of the study period (September 2018 - September 2019). Data from weather stations are weighted by the inverse of the squared distance of the station from the sample location to interpolate IDSW weather values for each site and date. Water samples are assigned the IDSW value for the day prior to that sample's collection to ensure the weather values used reflect antecedent meteorological conditions. For precipitation, rolling sums over the preceding 2 days ("2-day rainfall") and preceding 7 days, excluding the 2 most recent days ("3-7 day rainfall"), are computed to represent more- and less-recent cumulative antecedent precipitation, respectively.### Inverse distance squared interpolationLoad weather station and sampling site location data and re-project to match:```{r}#| label: meteo_loc# weather station locationsloc_ghcn <-readRDS("../gis/sf/meteo_pt.rds") |>st_transform(crs =6446) # project as NAD83(2011) / Georgia West# sample site locationsloc_site <-readRDS("../gis/sf/site_pt.rds") |>st_transform(crs =6446) # project as NAD83(2011) / Georgia West```Compute the Euclidean distance matrix between weather stations and sampling sites:```{r}#| label: meteo_dist# distance matrixdistmat_site_ghcn <-st_distance(loc_site, loc_ghcn)# add labelscolnames(distmat_site_ghcn) <- loc_ghcn$ghcn_id # weather stations in columnsrow.names(distmat_site_ghcn) <- loc_site$site_id # sample sites in rows```Convert distance matrix to a long-format dataframe:```{r}#| label: meteo_dist_longdist_site_ghcn <- distmat_site_ghcn |>as_tibble(rownames ="site_id") |>pivot_longer(-site_id,names_to ="ghcn_id",values_to ="dist") |>mutate(dist = units::drop_units(dist))```Read in weather data, join distances, drop missing values, and add an indicator ranking values for a given weather type on each day on the basis of distance from the sampling site to assist with identifying "neighboring" stations:```{r}#| label: meteo_dist_joindf_meteo_join <-readRDS("../out/ghcn_data.rds") |>filter(!is.na(meteo_val)) |># drop missing meteo valuesinner_join(dist_site_ghcn, # join distances to values, retaining only records with bothby ="ghcn_id", # join by station idrelationship ="many-to-many"# allow multiple matches on both sides (multiple stations inform multiple sites) ) |>arrange(site_id, meteo_type, dt_meteo, dist) |># ascending sort by distance within each site, type, and dategroup_by(site_id, meteo_type, dt_meteo) |>mutate(neighbor =row_number()) |># label values within a site, type, date with increasing distanceungroup()```Calculate inverse distance squared weights, identify k-nearest neighbors, and compute simple and weighted means to interpolate values from stations to sampling sites for each study date:```{r}#| label: meteo_idswdf_meteo <- df_meteo_join |>mutate(wt =1/ (dist^2), # IDSW weightwt_val = meteo_val * wt, # weighted weather valuewt_nn =ifelse(neighbor ==1, meteo_val, NA), # identify nearest stationwt_5nn =ifelse(neighbor <=6, meteo_val, NA) # identify 5 nearest stations ) |>group_by(site_id, dt_meteo, meteo_type, meteo_unit) |># aggregate by site, date, and weather typesummarise(n_station =sum(!is.na(meteo_val)), # count stations with valueswt_tot =sum(wt, na.rm =TRUE), # sum of weightsdist_near = dist[neighbor ==1], # distance to nearest stationdist_5nn =ifelse(n_station <5, # distance to 5th nearest station dist[neighbor == n_station], dist[neighbor ==5]),dist_far =max(dist, na.rm =TRUE), # distance to furthest stationmeteo_nn =mean(wt_nn, na.rm =TRUE), # assign value from nearest stationmeteo_5nn =mean(wt_5nn, na.rm =TRUE), # assign mean of 5 nearest stations meteo_avg =mean(meteo_val, na.rm =TRUE), # simple mean of all stationsmeteo_idsw =sum(wt_val, na.rm =TRUE) / wt_tot # IDSW mean ) |>ungroup() |>select(-wt_tot) |>mutate(across(meteo_nn:meteo_idsw,~ifelse(is.nan(.x), NA, .x))) # replace NAN with NAsaveRDS(df_meteo, "../out/meteo_site_idsw.rds")```### Rolling cumulative precipitationCompute moving-window ("rolling") sums of precipitation in the previous 2 days, 7 days, and 3 - 7 days:```{r}#| label: roll_raindf_roll_prcp <- df_meteo |>rename(phenom = meteo_type, unit = meteo_unit) |># change names for conveniencefilter(phenom =="prcp") |># just rainpivot_longer(starts_with("meteo_"), # long formatnames_to ="method", # weighting method (idsw, simple, nn, 5nn)values_to ="meteo0"# value on day 0 (sample date) ) |>mutate(method =str_remove(method, "meteo_")) |># drop prefixarrange(site_id, phenom, method, dt_meteo) |># order by date last (important for rolling)group_by(site_id, phenom, unit, method) |># by site, type, and weighting methodmutate(meteo1 =lag(meteo0, n =1L), # value on previous day (enforce temporality)meteo2 =rollapply(meteo1, width =2, # 2-day moving windowFUN = sum, # rolling sumfill =NA, # fill values outside window with NApartial =TRUE, # allow missing values in windowalign ="right"), # assign output to right-end of window (most recent)meteo7 =rollapply(meteo1,width =7, # 7-day moving windowFUN = sum,fill =NA,partial =TRUE,align ="right"),meteo3_7 =rollapply(lag(meteo0, 3), # start with 3 days priorwidth =5, # 5-day window starting day -3, ending day -7FUN = sum,fill =NA,partial =TRUE,align ="right"),rollfun ="sum"# record the rolling function applied ) |>ungroup()```Assign previous day's temperatures to each study date. Calculate rolling *mean* temperatures over the same windows used for rainfall, for consistency in the data structure:```{r}#| label: roll_temp#| code-fold: truedf_roll_temp <- df_meteo |>rename(phenom = meteo_type, unit = meteo_unit) |># change names for conveniencefilter(phenom !="prcp") |># not rainpivot_longer(starts_with("meteo_"), # long formatnames_to ="method", # weighting method (idsw, simple, nn, 5nn)values_to ="meteo0"# value on day 0 (sample date) ) |>mutate(method =str_remove(method, "meteo_")) |># drop prefixarrange(site_id, phenom, method, dt_meteo) |># order by date last (important for rolling)group_by(site_id, phenom, unit, method) |># by site, type, and weighting methodmutate(meteo1 =lag(meteo0, n =1L), # value on previous day (enforce temporality)meteo2 =rollapply(meteo1, width =2, # 2-day moving windowFUN = mean, # rolling meanfill =NA, # fill values outside window with NApartial =TRUE, # allow missing values in windowalign ="right"), # assign output to right-end of window (most recent)meteo7 =rollapply(meteo1,width =7, # 7-day moving windowFUN = mean,fill =NA,partial =TRUE,align ="right"),meteo3_7 =rollapply(lag(meteo0, 3), # start with 3 days priorwidth =5, # 5-day window starting day -3, ending day -7FUN = mean,fill =NA,partial =TRUE,align ="right"),rollfun ="mean"# record the rolling function applied ) |>ungroup()```Stack, sort, and save:```{r}#| label: roll_stackdf_meteo_roll <-bind_rows(df_roll_prcp, df_roll_temp) %>%arrange(site_id, phenom, method, dt_meteo)saveRDS(df_meteo_roll, "../out/meteo_site_roll.rds")```## Wastewater discharge: hydrological distance-weighted loadTreated wastewater and CSO inputs potentially impacting each site are represented by the aggregate maximum permitted load from all upstream outfalls, weighted by the inverse of the squared hydrological distance between the site and each upstream outfall along the stream network. While conceptually a similar IDSW operation, this approach requires measuring the distance along the stream network and restricting the calculations to outfalls that are both flow connected to, and upstream of, the sampling site. These operations are facilitated by constructing a *spatial stream network* (SSN) using the `SSNbler` and `SSN2` packages, the later of which provides functions for identifying upstream, flow-connected points and calculating network distances.Stream vector data were obtained from USGS National Hydrography Dataset Plus (NHDPlus) High Resolution and manually corrected in ArcGIS Pro to resolve topological inconsistencies, ensuring the entire network within the study watershed was flow-connected and discharged at a single terminal outlet located at the downstream-most sampling site (Fairburn Rd). The topology correction process created a small number of new, very short reaches to connect sections of the network. These new reaches are not in the NHD database and therefore lack cumulative drainage area values, which will later be used as weights as a proxy for flow contributions downstream and cannot be missing or zero. The smallest existing cummulative area (`area_cum`) value in the dataset is substituted for reaches with `area_cum = 0` to overcome this limitation.The pre-processed stream network vector data are provided as the shapefile `stream_line.shp` in the compressed folder [`stream.zip`](https://cdcgov.github.io/EMEL/manuscripts/river_ar/stream.zip). After downloading and extracting the compressed file, read in the stream shapefile and substitute the minimum cumulative area value for reaches with zero cumulative area values and save the `sf` for future use:```{r}#| label: stream_areastream <-read_sf("../gis/shape/stream/stream_line.shp") |>mutate(area_cum =ifelse(area_cum ==0,min(area_cum[area_cum >0], na.rm =TRUE), area_cum))saveRDS(stream, "../gis/sf/stream_line.rds")```Convert stream polyline vector to a *landscape network* (LSN) that represents the flow-directed network as nodes and edges using the `SSNbler` package:```{r}#| label: stream_lsn#| eval: falsestream_net <-lines_to_lsn(streams = stream,lsn_path ="../gis/ssn/lsn",snap_tolerance =0.1,check_topology =TRUE,topo_tolerance =20,overwrite =TRUE,verbose =TRUE,remove_ZM =TRUE)```Snap sample points to stream network: ```{r}#| label: out_site_snap#| eval: false# sample site locationssites <-readRDS("../gis/sf/site_pt.rds")# snap to stream LSNsites_lsn <-sites_to_lsn(sites = sites,edges = stream_net,lsn_path ="../gis/ssn/lsn",file_name ="sites",snap_tolerance =100,save_local =TRUE,overwrite =TRUE)```Identify unique outfall locations and snap to stream network```{r}#| label: out_falls_snap#| eval: false# outfall locationsoutfalls <-readRDS("../gis/sf/outfall_pt.rds") |>select(out_id, out_type, out_vol, geometry) |>distinct()# snap to stream LSNoutfalls_lsn <-sites_to_lsn(sites = outfalls,edges = stream_net,lsn_path ="../gis/ssn/lsn",file_name ="outfalls",snap_tolerance =100,save_local =TRUE,overwrite =TRUE)```Calculate the upstream distance along each edge (stream segment):```{r}#| label: out_updist_edge#| eval: falseedges <-updist_edges(edges = stream_net,save_local =TRUE,lsn_path ="../gis/ssn/lsn",calc_length =TRUE)```Calculate the upstream distance between sites and outfalls:```{r}#| label: out_updist_site#| eval: falsesite_list <-updist_sites(sites =list(sites = sites_lsn,outfalls = outfalls_lsn),edges = edges,length_col ="Length",save_local =TRUE,lsn_path ="../gis/ssn/lsn")```Additive function values (AFV) are numerical values associated with each edge (stream reach) that represent the contribution of that edge to the edges downstream of it. For stream networks, the AFV can encode the flow contributed by each reach (which is an additive function of the flows from all reaches upstream). Because flow is variable and unknown for most reaches, the cumulative area drained by each reach is used as a proxy for its potential flow. The AFVs are not used directly for the network distance calculations but are required to create a valid SSN object, from which the distances are later calculated.Assign AFVs to the stream network edges and compute the AFVs at each site:```{r}#| label: out_afv#| eval: false# edges (reaches)edges_afv <-afv_edges(edges = edges,lsn_path ="../gis/ssn/lsn",save_local =TRUE,overwrite =TRUE,infl_col ="area_cum",segpi_col ="segpi_area",afv_col ="afv_area")# sample sitessites_afv <-afv_sites(sites = site_list,edges = edges_afv,lsn_path ="../gis/ssn/lsn",save_local =TRUE,overwrite =TRUE,afv_col ="afv_area")```Assemble an SNN from the LSN and various components generated above; while LSNs are network data objects and can be reused, SSNs are data models that are specific to the sites and data associated with them and are best assembled separately for each use case. For this application, record the sampling sites as the "observed sites" (`obs_sites`) and the outfall locations as the "prediction sites" (`preds_list`):```{r}#| label: out_ssn#| eval: falseoutfall_ssn <-ssn_assemble(edges = edges_afv,lsn_path ="../gis/ssn/lsn",obs_sites = sites_afv$sites,preds_list = sites_afv[c("outfalls")],ssn_path ="../gis/ssn/outfalls.ssn",import =FALSE,check =TRUE,afv_col ="afv_area",overwrite =TRUE)```Calculate the network distance matrix between the sites and outfalls (writes to disk):```{r}#| label: out_distmat#| eval: false#| warning: false# import SSNoutfall_ssn <-ssn_import("../gis/ssn/outfalls.ssn",predpts ="outfalls",overwrite =TRUE)# calculate distance matrixssn_create_distmat(outfall_ssn,predpts ="outfalls",only_predpts =TRUE,among_predpts =FALSE,overwrite =TRUE)```Two distance matrices are returned in a list; only the first, with the sites in rows and the outfalls in columns, provides the flow-connected, upstream network distances required. Extract the correct distance matrix and label with IDs extracted from the SSN elements:```{r}#| label: out_ssn_import#| include: false# re-import SSNoutfall_ssn <-ssn_import("../gis/ssn/outfalls.ssn",predpts ="outfalls",overwrite =TRUE)``````{r}#| label: out_dist_extract# ssn elementsoutfall_ssn_stream <- outfall_ssn$edgesoutfall_ssn_site <- outfall_ssn$obsoutfall_ssn_out <- outfall_ssn$preds$outfalls# extract downstream distances (from outfalls to sites)outfall_distmat <-ssn_get_stream_distmat(outfall_ssn,name ="outfalls") |>pluck("dist.net1.a")colnames(outfall_distmat) <- outfall_ssn_out$out_idrownames(outfall_distmat) <- outfall_ssn_site$site_id```Re-format the matrix as a long dataframe and save:```{r}#| label: out_dist_longoutfall_dist <- outfall_distmat %>%as_tibble(rownames ="site_id") %>%pivot_longer(-site_id,names_to ="out_id",values_to ="dist")saveRDS(outfall_dist, "../out/outfall_dist_ssn.rds")```Merge the outfall data with the distance values:```{r}#| label: out_dist_joindf_out_dist <-readRDS("../gis/sf/outfall_pt.rds") |>st_drop_geometry() %>%inner_join(outfall_dist,by =c("site_id", "out_id"))saveRDS(df_out_dist, "../out/outfall_riverdist.rds")```Weighted total (cumulative) permitted loads are used to represent the influence of wastewater infrastructure on sites. Because weighted sums are not normalized by the total weights, as was the case for IDSW interpolation of weather data and IDSW land cover proportions. As such, the weighted output values are dependent on the distance unit used and cannot be interpreted in isolation or compared to other contexts. For convenience, distance is converted to kilometers, such that an outfall 1 km from the site receives a weight of 1, outfalls further than 1 km from the site receive smaller weights less than 1, and outfalls closer than 1 km to the site receive larger weights greater than 1.Calculate weights and and IDSW cumulative load on the basis of permitted maximum discharge separately for WWTP treated effluent outfalls and for CSOs:```{r}#| label: out_idswdf_out_idsw <-readRDS("../out/outfall_riverdist.rds") |>mutate(dist = dist /1000, # scale to kilometerswt =1/ (dist^2),wt_val = out_vol * wt) |>group_by(site_num, site_id, site_lab, site_type, out_type) |>summarise(n_outfall =sum(!is.na(out_vol)),wt_tot =sum(wt, na.rm =TRUE),out_tot =sum(out_vol, na.rm =TRUE),out_idsw =sum(wt_val, na.rm =TRUE),out_tot_avg = out_tot / n_outfall,out_idsw_avg = out_idsw / wt_tot) |>ungroup() |>select(-wt_tot) |>mutate(across(out_tot:out_idsw_avg, ~ifelse(is.nan(.x), NA, .x)))saveRDS(df_out_idsw, "../out/outfall_site_idsw.rds")```## Septic systems: basin spatial densityAs point locations lacking additional scale information (like the permitted discharge for outfalls) and heterogeneously distributed within the study area, inverse distance weighting approaches are less appropriate for constructing site-level variables from septic system locations. Rather, the simple spatial density of all septic systems within a drainage basin is used to represent the potential influence of septic dystems on a sampling site. Load the septic locations, tally up the systems in each drainage basin, and calculate the spatial density; because some basins contained no septic systems, restore the dropped basins and set septic density to zero:```{r}#| label: septic_densseptic_dens <-readRDS("../gis/sf/septic_pt.rds") |>st_drop_geometry() |>group_by(site_num, wshed_area_clp) |>summarise(septic_n =n()) |>ungroup() |>mutate(septic_dens = septic_n / wshed_area_clp) |># some basins had no septic so restore and set density to zerofull_join(tibble(site_num =3:13), by ="site_num") |>mutate(across(starts_with("septic_"), ~replace_na(.x, 0))) |>arrange(site_num)saveRDS(septic_dens, "../out/septic_density.rds")```# Data set assemblyThe generated variables are assembled into a full data set of *E. coli* observations and the environmental factors corresponding to each sampling site and date. Subsets of this full data set are extracted and assembled into specific SSNs to facilitate analysis by spatial generalized linear mixed models (GLMM) that account for flow-connected network distances between sites and and temporal clustering by time of sample collection.```{r}#| label: refresh_table#| include: falserm(list =ls())```## Full data tableRead in the processed sample data (*E. coli* imputed concentrations, water quality; meropenem-selected and cefotaxime-selected at 4 μg/mL already excluded): ```{r}#| label: flat_sampdf_samp <-readRDS("../out/sample_imputed.rds")```Read in the rolling weather data, restrict to study time period and desired variables, pivot to wide format (each weather variable is a column), and scale the rain values to units of cm (10 mm):```{r}#| label: flat_meteodf_meteo <-readRDS("../out/meteo_site_roll.rds") |>select(site_id, dt_meteo, phenom, method,`1`= meteo1, `2`= meteo2, `7`= meteo3_7) |>filter(method =="idsw", phenom %in%c("prcp", "tmax", "tavg"), dt_meteo >ymd("2018-09-01"), dt_meteo <ymd("2019-09-30")) |>mutate(phenom =ifelse(phenom =="prcp", "rain", phenom)) |>pivot_wider(names_from = phenom,values_from =c(`1`, `2`, `7`),names_glue ="{phenom}{.value}",names_vary ="slowest") |>select(-method) |>mutate(across(starts_with("rain"), ~.x/10))```Read in land cover proportions and mean imperviousness, select desired variables, and scale so all land cover variables are in units of 10 percentage points and basin area is log~10~-transformed: ```{r}#| label: flat_landdf_lulc <-readRDS("../out/basin_features.rds") |>filter(method =="idsw") |>select(site_num, site_id, wshed_area,starts_with("lc")) |>mutate(lc_impv = lc_impv /100,across(starts_with("lc_"), ~.x*10 ),wshed_area =log10(wshed_area))```Read in outfall loads, select desired variables, differentiate between WWTPs and CSOs, and restore sites without outfalls, substituting zero for missing load:```{r}#| label: flat_outdf_outfall <-readRDS("../out/outfall_site_idsw.rds") |>select(site_num, out_type,n = n_outfall,load = out_idsw) %>%mutate(out_type =ifelse(out_type =="treated", "wwtp", out_type)) %>%pivot_wider(names_from = out_type,values_from =c(n, load),names_glue ="out_{out_type}_{.value}",names_vary ="slowest") %>%full_join(tibble(site_num =3:13), by ="site_num") %>%mutate(across(starts_with("out_"), ~replace_na(.x, 0)))```Read in septic density and join all the flat data together:```{r}#| label: flat_merge# septic densitydf_septic <-readRDS("../out/septic_density.rds") |>select(site_num, septic_n, septic_dens)# combine for analysisdf_analysis <- df_samp |>left_join(df_lulc, by =c("site_num", "site_id")) %>%left_join(df_outfall, by ="site_num") %>%left_join(df_septic, by ="site_num") %>%left_join(df_meteo, by =c("site_id", "dt_collect"="dt_meteo"))saveRDS(df_analysis, "../out/analysis_data.rds")```With exceptions for the outcome *E. coli* concentrations and raw counts, variables are standard deviation-scaled for analysis. Variables without a meaningful zero value are also centered around the mean prior to scaling. Calculate the mean and standard deviation of each variable by the selective drug used to inform centering and scaling:```{r}#| label: flat_tabtab_analysis <- df_analysis |>select(drug, drug_dose,mpn = mpn_r,mpn_gen = mpn_r_gen,starts_with("wq_"), wshed_area,starts_with("lc_"),starts_with("out_"),starts_with("septic_"), rain2, rain7, tmax1, tavg1) |>pivot_longer(-drug,names_to ="var",values_to ="val") |>filter(!is.na(val)) |>group_by(drug, var) |>summarise(nobs =n(),avg =mean(val),sd =sd(val),min =min(val),max =max(val)) |>ungroup() |>mutate(scl_method =case_when(var %in%c("mpn") ~"none", # response unscaled (already log10) var %in%c("mpn_gen") ~"scl", # predictor EC scaled (already log10)str_detect(var, "lc_") ~"scl", # lulc proportions scaled (0 meaningful)str_detect(var, "wq_") ~"std", # water quality centered, scaledstr_detect(var, "load") ~"scl", # outfall loads scaled (0 meaningful)str_detect(var, "_n") ~"none", # outfall counts unscaled (0 meaningful)str_detect(var, "_dens") ~"scl", # septic density scaled (0 meaningful) var %in%c("wshed_area") ~"std", # area centered, scaled (already log10)str_detect(var, "rain") ~"scl", # rain scaled (0 meaningful) var %in%c("tmax1", "tavg1") ~"std")) # temp centered, scaledsaveRDS(tab_analysis, "../out/covar_summ_stats.rds")```Apply centering and scaling as determined:```{r}#| label: flat_scaledf_analysis_scl <- df_analysis |>select(site_num, drug, drug_dose, wk_collect, dt_collect, samp_event, samp_id, uid, uid_gen, uid_wq, detect,mpn = mpn_r,mpn_gen = mpn_r_gen,starts_with("wq_"), wshed_area,starts_with("lc_"),starts_with("out_"),starts_with("septic_"), rain2, rain7, tmax1, tavg1) |>select(-lc_tot) |>pivot_longer(-c(site_num, drug, drug_dose, wk_collect, dt_collect, samp_event, samp_id, uid, uid_gen, uid_wq, detect),names_to ="var",values_to ="val") |>left_join(tab_analysis, by =c("drug", "var")) |>mutate(val_scl =case_when(scl_method =="none"~ val, scl_method =="scl"~ val / sd, scl_method =="std"~ (val - avg) / sd) )|>select(-c(val, nobs, avg, sd, min, max, scl_method)) |>pivot_wider(names_from = var,values_from = val_scl)saveRDS(df_analysis_scl, "../out/analysis_data_scl.rds")```## Spatial stream networksThe SSN objects are somewhat particular about the data used to assemble them and the opportunities for subsetting post-assembly. While inefficient, a reliable approach is to assemble separate SSNs for each set of *outcome* variables that are to be analyzed. The *E. coli* occurrence under each selective antibiotic is a separate outcome for the purposes of these analyses. Furthermore, the co-measured sample-level covariates (turbidity, optical brighteners, and generic *E. coli* for the antibiotic-selected observations) are missing for some of the *E. coli* observations, resulting in different sample sizes under complete-cases analyses. Separate SSNs are therefore assembled for each combination of *E. coli* outcome and covariate completeness. Fortunately, the majority of environmental factors covariates considered, which were derived from comprehensive spatial data sources, are complete for all observations, allowing a single SSN to be used for all analyses of a given outcome that do not include co-measured, sample-level covariates.```{r}#| label: refresh_ssn#| include: falserm(list =ls())```### Landscape networkThe LSN previously created in the process of calculating network distances to wastewater outfalls has those outfall locations associated with it as prediction locations. For simplicity, construct a new LSN with just the site data associated to use in assembling the new SSNs.Generate new LSN stream network:```{r}#| label: lsn_full_stream#| eval: false# stream vectorstream <-readRDS("../gis/sf/stream_line.rds")# create new LSN for full datasetstream_net <-lines_to_lsn(streams = stream,lsn_path ="../gis/ssn/lsn_full",snap_tolerance =0.1,check_topology =TRUE,topo_tolerance =20,overwrite =TRUE,verbose =TRUE,remove_ZM =TRUE)```Snap site locations to network:```{r}#| label: lsn_full_site#| eval: false# sample site locationssites <-readRDS("../gis/sf/site_pt.rds")# snap to stream LSNsites_lsn <-sites_to_lsn(sites = sites,edges = stream_net,lsn_path ="../gis/ssn/lsn_full",file_name ="sites",snap_tolerance =100,save_local =TRUE,overwrite =TRUE)```Calculate the upstream distance along each edge (stream segment) and for sample sites:```{r}#| label: lsn_full_updist#| eval: false# edges (stream reaches)edges <-updist_edges(edges = stream_net,save_local =TRUE,lsn_path ="../gis/ssn/lsn_full",calc_length =TRUE)# sample sitessite_list <-updist_sites(sites =list(sites = sites_lsn),edges = edges,length_col ="Length",save_local =TRUE,lsn_path ="../gis/ssn/lsn_full")```Calculate additive function values for edges and sites, saving as `sf` for easy access when generating SSNs:```{r}#| label: lsn_full_afv#| eval: false# edges (stream reaches)edges_afv <-afv_edges(edges = edges,lsn_path ="../gis/ssn/lsn_full",save_local =TRUE,overwrite =TRUE,infl_col ="area_cum",segpi_col ="segpi_area",afv_col ="afv_area")saveRDS(edges_afv, "../gis/sf/lsn_edges_afv_line.rds")# sample sitessites_afv <-afv_sites(sites = site_list,edges = edges_afv,lsn_path ="../gis/ssn/lsn_full",save_local =TRUE,overwrite =TRUE,afv_col ="afv_area")saveRDS(sites_afv$sites, "../gis/sf/lsn_sites_afv_pt.rds")```Join the full, scaled analysis data set with all observations (dropping any missing *E. coli* outcomes) to the LSN sites:```{r}#| label: lsn_full_obs# full data set for all observationsobs_full <-readRDS("../out/analysis_data_scl.rds") |>filter(!is.na(detect))# merge with LSN sitessites_full <-readRDS("../gis/sf/lsn_sites_afv_pt.rds") |>left_join(obs_full, by ="site_num")saveRDS(sites_full, "../out/analysis_data_lsn.rds")```### Cefotaxime-selected SSNsAssemble SSNs for ESBL-Ec (cefotaxime-selected at 1 μg/mL), the main environmental AR outcome for the modeling analysis.First assemble an SSN for all the ESBL-Ec observations and create the distance matrix:```{r}#| label: ssn_esbl_full#| eval: false# LSN edgesedges_afv <-readRDS("../gis/sf/lsn_edges_afv_line.rds")# filter for all ESBL obssites_esbl <-readRDS("../out/analysis_data_lsn.rds") |>filter(drug =="cefotaxime")# assemble full ESBL SSNesbl_ssn <-ssn_assemble(edges = edges_afv,lsn_path ="../gis/ssn/lsn_full",obs_sites = sites_esbl,preds_list =NULL,ssn_path ="../gis/ssn/esbl.ssn",import =TRUE,check =TRUE,afv_col ="afv_area",overwrite =TRUE)# generate distance matrixssn_create_distmat(esbl_ssn,overwrite =TRUE)```Assemble an SSN and distance matrix for only the ESBL-Ec observations with non-missing generic *E. coli* concentrations:```{r}#| label: ssn_esbl_gen#| eval: false# non-missing generic E. coli datasites_esbl_gen <- sites_esbl %>%filter(!is.na(mpn_gen))# assemble ssnesbl_ssn_gen <-ssn_assemble(edges = edges_afv,lsn_path ="../gis/ssn/lsn_full",obs_sites = sites_esbl_gen,preds_list =NULL,ssn_path ="../gis/ssn/esbl_gen.ssn",import =TRUE,check =TRUE,afv_col ="afv_area",overwrite =TRUE)# distance matrixssn_create_distmat(esbl_ssn_gen,overwrite =TRUE)```Assemble SSNs for only the ESBL-Ec observations with complete water quality covariates:```{r}#| label: ssn_esbl_wq#| eval: false## Turbidity# non-missing turbiditysites_esbl_turb <- sites_esbl %>%filter(!is.na(wq_turb))# assemble ssnesbl_ssn_turb <-ssn_assemble(edges = edges_afv,lsn_path ="../gis/ssn/lsn_full",obs_sites = sites_esbl_turb,preds_list =NULL,ssn_path ="../gis/ssn/esbl_turb.ssn",import =TRUE,check =TRUE,afv_col ="afv_area",overwrite =TRUE)# distance matrixssn_create_distmat(esbl_ssn_turb,overwrite =TRUE)## Optical brighteners# non-missing optical brightenerssites_esbl_optic <- sites_esbl %>%filter(!is.na(wq_optic))# assemble ssnesbl_ssn_optic <-ssn_assemble(edges = edges_afv,lsn_path ="../gis/ssn/lsn_full",obs_sites = sites_esbl_optic,preds_list =NULL,ssn_path ="../gis/ssn/esbl_optic.ssn",import =TRUE,check =TRUE,afv_col ="afv_area",overwrite =TRUE)# distance matrixssn_create_distmat(esbl_ssn_optic,overwrite =TRUE)```### Colistin-selected SSNsAssemble the same set of SSNs for the CST-Ec (colistin-selected at 2 μg/mL) outcome:```{r}#| label: ssn_cst#| eval: false#| code-fold: true## All observations# filter for all ESBL obssites_cst <-readRDS("../out/analysis_data_lsn.rds") |>filter(drug =="colistin")# assemble full ESBL SSNcst_ssn <-ssn_assemble(edges = edges_afv,lsn_path ="../gis/ssn/lsn_full",obs_sites = sites_cst,preds_list =NULL,ssn_path ="../gis/ssn/cst.ssn",import =TRUE,check =TRUE,afv_col ="afv_area",overwrite =TRUE)# generate distance matrixssn_create_distmat(cst_ssn,overwrite =TRUE)## Generic E. coli# non-missing generic E. coli datasites_cst_gen <- sites_cst %>%filter(!is.na(mpn_gen))# assemble ssncst_ssn_gen <-ssn_assemble(edges = edges_afv,lsn_path ="../gis/ssn/lsn_full",obs_sites = sites_cst_gen,preds_list =NULL,ssn_path ="../gis/ssn/cst_gen.ssn",import =TRUE,check =TRUE,afv_col ="afv_area",overwrite =TRUE)# distance matrixssn_create_distmat(cst_ssn_gen,overwrite =TRUE)## Turbidity# non-missing turbiditysites_cst_turb <- sites_cst %>%filter(!is.na(wq_turb))# assemble ssncst_ssn_turb <-ssn_assemble(edges = edges_afv,lsn_path ="../gis/ssn/lsn_full",obs_sites = sites_cst_turb,preds_list =NULL,ssn_path ="../gis/ssn/cst_turb.ssn",import =TRUE,check =TRUE,afv_col ="afv_area",overwrite =TRUE)# distance matrixssn_create_distmat(cst_ssn_turb,overwrite =TRUE)## Optical brighteners# non-missing optical brightenerssites_cst_optic <- sites_cst %>%filter(!is.na(wq_optic))# assemble ssncst_ssn_optic <-ssn_assemble(edges = edges_afv,lsn_path ="../gis/ssn/lsn_full",obs_sites = sites_cst_optic,preds_list =NULL,ssn_path ="../gis/ssn/cst_optic.ssn",import =TRUE,check =TRUE,afv_col ="afv_area",overwrite =TRUE)# distance matrixssn_create_distmat(cst_ssn_optic,overwrite =TRUE)```### Generic *E. coli* SSNsFinally, assemble a comparable set of SSNs for the generic *E. coli* (no antibiotic selection) outcome, skipping the SSN for non-missing generic *E. coli* as a covariate:```{r}#| label: ssn_gen#| eval: false#| code-fold: true## All observations# filter for all ESBL obssites_gen <-readRDS("../out/analysis_data_lsn.rds") |>filter(drug =="none")# assemble full ESBL SSNgen_ssn <-ssn_assemble(edges = edges_afv,lsn_path ="../gis/ssn/lsn_full",obs_sites = sites_gen,preds_list =NULL,ssn_path ="../gis/ssn/gen.ssn",import =TRUE,check =TRUE,afv_col ="afv_area",overwrite =TRUE)# generate distance matrixssn_create_distmat(gen_ssn,overwrite =TRUE)## Turbidity# non-missing turbiditysites_gen_turb <- sites_gen %>%filter(!is.na(wq_turb))# assemble ssngen_ssn_turb <-ssn_assemble(edges = edges_afv,lsn_path ="../gis/ssn/lsn_full",obs_sites = sites_gen_turb,preds_list =NULL,ssn_path ="../gis/ssn/gen_turb.ssn",import =TRUE,check =TRUE,afv_col ="afv_area",overwrite =TRUE)# distance matrixssn_create_distmat(gen_ssn_turb,overwrite =TRUE)## Optical brighteners# non-missing optical brightenerssites_gen_optic <- sites_gen %>%filter(!is.na(wq_optic))# assemble ssngen_ssn_optic <-ssn_assemble(edges = edges_afv,lsn_path ="../gis/ssn/lsn_full",obs_sites = sites_gen_optic,preds_list =NULL,ssn_path ="../gis/ssn/gen_optic.ssn",import =TRUE,check =TRUE,afv_col ="afv_area",overwrite =TRUE)# distance matrixssn_create_distmat(gen_ssn_optic,overwrite =TRUE)```# Descriptive analysisBefore developing statistical models to estimate associations between surface water *E. coli* and environmental factors, conduct descriptive analyses to better understand the contours of the prepared data.## *E. coli* time seriesExplore the occurrence of generic *E. coli*, ESBL-Ec, CST-EC, and MEM-Ec at each site across the study period by plotting time series of their detection and concentration.Read in the (non-imputed) sample data time series and arrange for plotting:```{r}#| label: df_plot_ecdf_plot_ec <-readRDS("../out/sample_data.rds") |>filter(dt_collect >ymd("2018-09-05"), dt_collect <ymd("2019-09-06")) |>filter(!is.na(mpn), drug %in%c("cefotaxime", "colistin", "none"), drug_dose <4) |>mutate(mpn_log10 =log10(mpn),cens_lab =fct_relevel(cens, c("no", "left", "right")) |>fct_recode("detect"="no","non-detect"="left","TNTC"="right"),site_type =fct_relevel(site_type, c("river", "tributary"))) |>arrange(site_type, site_num) |>mutate(site_id =fct_inorder(site_id),site_lab =fct_inorder(site_lab),type =ifelse(drug =="none", "generic", "resistant"),drug_lab =fct_recode(drug,"ESBL"="cefotaxime","CST"="colistin","Generic"="none") |>fct_relevel(c("Generic", "ESBL", "CST")))```Separately, read in the MEM-Ec detections and arrange for plotting:```{r}#| label: df_plot_mem#| code-fold: truedf_plot_mem <-readRDS("../out/sample_data.rds") |>filter(dt_collect >ymd("2018-09-05"), dt_collect <ymd("2019-09-06")) |>filter(site_num >2,!is.na(mpn), drug =="meropenem", detect ==1) |>mutate(site_type =fct_relevel(site_type, c("river", "tributary"))) |>arrange(site_type, site_num) |>mutate(site_id =fct_inorder(site_id),site_lab =fct_inorder(site_lab)) |>select(site_type, site_lab, dt_collect, detect)```Extract the range of dates and concentration values to set axis limits:```{r}#| label: df_rangedf_range <- df_plot_ec %>%ungroup() %>%summarise(dt_min =min(dt_collect),dt_max =max(dt_collect),mpn_max =max(mpn_log10))```Set color and line-type palettes for plotting:```{r}#| label: p_ts_palval_pal <-c("Generic"="#e66101", "ESBL"="#5e3c99", "CST"="seagreen")val_line <-c("Generic"="solid", "ESBL"="32", "CST"="11")```To allow for the nice arrangement of 11 different sites, plot the river and tributary sites separately and combine. First plot the river sites, adding LOESS smooths for each *E. coli* concentration outcome as a visual aid and differentiating between detect, non-detect, and TNTC observations by shape:```{r}#| label: p_ts_riverp_river <- df_plot_ec |>filter(site_type =="river") %>%ggplot(data = .,aes(x = dt_collect, y = mpn_log10)) +facet_grid(rows =vars(site_lab)) +geom_vline(aes(xintercept = dt_collect),data = df_plot_mem %>%filter(site_type =="river"),color ="cornflowerblue",lwd =1.3, alpha =0.85) +geom_smooth(aes(color = drug_lab, fill = drug_lab, linetype = drug_lab),se =FALSE) +geom_point(aes(shape = cens_lab, fill = drug_lab, colour = drug_lab),show.legend =TRUE,size =2.5, stroke =1.5,alpha =0.85) +scale_x_date(breaks =date_breaks(width ="2 month"),limits =c(df_range$dt_min, df_range$dt_max),date_labels ="%b '%y") +annotation_logticks(side ="l") +scale_y_continuous(labels =math_format(10^.x)) +coord_cartesian(ylim =c(0, 5)) +scale_fill_manual(values = val_pal,drop =FALSE,name ="*E. coli* phenotype") +scale_color_manual(values = val_pal,drop =FALSE,name ="*E. coli* phenotype") +scale_linetype_manual(values = val_line,drop =FALSE,name ="*E. coli* phenotype") +scale_shape_manual(values =c("detect"=16, "non-detect"=4, "TNTC"=0),drop =FALSE,name ="Status") +labs(x =NULL,y ="MPN/100 mL",subtitle ="River sites")```Repeat for tributary sites:```{r}#| label: p_ts_trib#| code-fold: truep_trib <- df_plot_ec |>filter(site_type =="tributary") %>%ggplot(data = .,aes(x = dt_collect, y = mpn_log10)) +facet_grid(rows =vars(site_lab)) +geom_vline(aes(xintercept = dt_collect),data = df_plot_mem %>%filter(site_type =="tributary"),color ="cornflowerblue",lwd =1.3, alpha =0.85) +geom_smooth(aes(color = drug_lab, fill = drug_lab, linetype = drug_lab),se =FALSE) +geom_point(aes(shape = cens_lab, fill = drug_lab, colour = drug_lab),show.legend =TRUE,size =2.5, stroke =1.5,alpha =0.85) +scale_x_date(breaks =date_breaks(width ="2 month"),limits =c(df_range$dt_min, df_range$dt_max),date_labels ="%b '%y") +annotation_logticks(side ="l") +scale_y_continuous(labels =math_format(10^.x)) +coord_cartesian(ylim =c(0, 5)) +scale_fill_manual(values = val_pal,drop =FALSE,name ="*E. coli* phenotype") +scale_color_manual(values = val_pal,drop =FALSE,name ="*E. coli* phenotype") +scale_linetype_manual(values = val_line,drop =FALSE,name ="*E. coli* phenotype") +scale_shape_manual(values =c("detect"=16, "non-detect"=4, "TNTC"=0),drop =FALSE,name ="Status") +labs(x =NULL,y =NULL,subtitle ="Tributary sites")```Arrange the two plots and save the image as a .png with specific dimensions to ensure proper display:```{r}#| label: p_ts#| warning: falsep_ts <- p_river + p_trib +plot_layout(ncol =2,design =" 12 12 12 12 12 1# ",heights =c(6, 4.6),guides ="collect") &guides(colour =guide_legend(order =1),fill =guide_legend(order =1),linetype =guide_legend(order =1),shape =guide_legend(order =2)) &theme_bw() &theme(axis.title.y =element_markdown(),plot.title =element_markdown(),panel.grid.minor.y =element_blank(),legend.title =element_markdown(),legend.position ="bottom",legend.justification ="center",legend.key.width =unit(.4, "in"),legend.direction ="horizontal")ggsave(filename ="../fig/ecoli_timeseries.png",plot = p_ts,device = agg_png, dpi =300,height =11, width =8.5, units ="in")```This plot reproduces **Figure 2** from Holcomb et al. (2026):## *E. coli* correlationWhile overall associations between generic and resistant *E. coli* are estimated later using models to account for clustering in space and in time, compute non-parametric correlations to examine concordance between paired concentrations of generic and resistant *E. coli* measured in the same samples for each site and site type (river mainstem vs. tributary).Define a function to compute Kendall's $\tau_{b}$ correlations between paired observations with 95% confidence intervals (CI) and export in tidy format:```{r}#| label: f_corr_cicorr_ci <-function(x, # vector of values y, # vector of paired values (same length as x)conf.level =0.95# CI width ){ nobs =length(x) # number of paired observations ci <-KendallTauB(x, y, # estimate Kendall's tau with CIconf.level = conf.level) pval <-cor.test(x, y, method ="kendall",exact =FALSE) out <-tibble(test ="Tau-b",nobs = nobs,corr =unname(ci[1]),corr_lo =unname(ci[2]),corr_hi =unname(ci[3]),pval = pval$p.value)return(out)}```Because Kendall's $\tau_{b}$ correlation is based on ranks and accounts for ties, it is robust to the values used to represent left- and right-censored observations. Prepare the generic and resistant *E. coli* concentration data by substituting LOD/2 for non-detects and ULOQ + 1 for TNTC observations:```{r}#| label: df_corr_subdf_corr_sub <-readRDS("../out/sample_data.rds") |>filter(site_num >2,!is.na(detect)) |>mutate(mpn =case_when(cens =="left"~ mpn_lod /2, cens =="right"~ mpn_uloq +1, cens =="no"~ mpn))```Filter for just the generic *E. coli* to join back to the resistant *E. coli* from the same samples:```{r}#| label: df_corr_gendf_corr_gen <- df_corr_sub |>filter(drug =="none") |>select(uid, mpn, mpn_lod, mpn_uloq)```Join generic *E. coli* to resistant *E. coli* and retain only non-missing pairs:```{r}#| label: df_corrdf_corr <- df_corr_sub |>filter(drug %in%c("cefotaxime", "colistin"), drug_dose <4) |>left_join(df_corr_gen,by =c("uid_gen"="uid"),suffix =c("", "_gen")) |>filter(!is.na(mpn),!is.na(mpn_gen))```Nest the data by selective antibiiotic and site and apply the `corr_ci()` by group to calculate site-specific correlations:```{r}#| label: tab_corr_sitetab_corr_site <- df_corr |>group_by(site_num, site_id, site_lab, site_type, drug) |>mutate(n_detect =sum(detect)) |>filter(n_detect >0) |>nest() |>ungroup() |>mutate(cor =map(data, ~corr_ci(.x$mpn, .x$mpn_gen))) |>select(-data)```Repeat, nesting by site type:```{r}#| label: tab_corr_type#| code-fold: truetab_corr_type <- df_corr |>group_by(site_type, drug) |>mutate(n_detect =sum(detect)) |>filter(n_detect >0) |>nest() |>ungroup() |>mutate(cor =map(data, ~corr_ci(.x$mpn, .x$mpn_gen))) |>select(-data) |>mutate(site_lab = site_type,site_num =0)```Calculate across all sites and types, then arrange to output as table:```{r}#| label: tab_corrtab_corr <- df_corr |>group_by(drug) |>mutate(n_detect =sum(detect)) |>filter(n_detect >0) |>nest() |>ungroup() |>mutate(site_lab ="Total",cor =map(data, ~corr_ci(.x$mpn, .x$mpn_gen))) |>select(-data) |>bind_rows(tab_corr_type, tab_corr_site) |>unnest(cor) |>mutate(ci_val =str_c(pretty_round(corr_lo, 2), "–",pretty_round(corr_hi, 2)), pval_val =ifelse(pval <0.001, "<0.001", pretty_round(pval, 3)),drug =case_when(drug =="cefotaxime"~"CTX", drug =="colistin"~"CST")) |>select(starts_with("site_"), drug, nobs, corr,ci = ci_val,pval = pval_val) |>pivot_wider(names_from = drug,values_from =c(nobs, corr, ci, pval),names_glue ="{drug}_{.value}",names_vary ="slowest") |>arrange(site_type, site_num) |>select(-site_num, -site_id)```This reproduces the results (if not quite the formatting) of **Table S3** in the Supplementary Material of Holcomb et al. (2026):```{r}#| label: tab_corr_kable#| echo: falsetab_corr |>kable(col.names =c("Site", "Type","N", "τ", "95% CI", "p","N", "τ", "95% CI", "p"),digits =2) |>add_header_above(header =c(" "=2,"ESBL-Ec vs Generic"=4,"CST-Ec vs Generic"=4))```# Regression modelingThe primary analysis approach in this study is the application of spatial generalized linear mixed models (GLMM) to estimate associations between environmental factors and the occurrence of generic and resistant *E. coli* in surface water. The spatial GLMMs are implemented as spatial stream network models using the SSN objects generated previously.Linear SSN models are used to estimate the expected change in log~10~ concentration of generic *E. coli* and ESBL-Ec and logistic SSN models are used to estimate the odds ratio (OR) for CST-Ec detection associated with a unit-increase in each environmental factor. The SSN models all include a random effect for sampling event to account for temporal clustering among contemporaneous samples and a flow-weighted spatial random effect with exponential distance decay to account for spatial autocorrelation between flow-connected sampling locations (a "tail-up" covariance component in the SSN parlance).```{r}#| label: refresh_model#| include: falserm(list =ls())```## Model-fitting functionWrite a utility function to execute an SSN model and extract the model estimates in a convenient format with identifying labels:```{r}#| label: f_fit_ssnfit_ssn <-function(lab, # model label to assign to output (string) set, # model set to label output (string) outcome, # outcome assessed (string) response, # response variable representing outcome (string) predictor, # explanatory variable set label (string) error, # covariance structure label (string) re, # random effect label (string) path, # path to SSN for model fit form, # model formularandom =NULL, # random effect formulafamily ="Gaussian", # likelihoodadditive ="afv_area", # AFV column name (string)tailup_type ="none", # correlation model for tail-up (flow-connected) covariance componenttaildown_type ="none", # correlation model for tail-down (flow-disconnected) covariance componenteuclid_type ="none", # correlation model for Euclidean distance covariance componentnugget_type ="nugget", # include unstructured error component (aka nugget)est_method ="reml"# estimation algorithm ){## import SSN net <-ssn_import(path,overwrite =TRUE)## check formula form <-as.formula(form)## fit model mod <-ssn_glm(formula = form,random = random,ssn.object = net,family = {{family}},tailup_type = tailup_type,taildown_type = taildown_type,euclid_type = euclid_type,nugget_type = nugget_type,additive = {{additive}},est_method = est_method)## extract estimates est_coef <-tidy(mod, conf.int =TRUE, effects ="fixed") |># regression coefs ("fixed effects")mutate(effect ="coef") est_var <-tidy(mod, conf.int =TRUE, effects ="ssn") |># variance componentsfilter(!is_known)if(!is.null(random)){ est_re <-tidy(mod, conf.int =TRUE, effects ="randcov") |># random effect variancemutate(effect ="re") est_raw <-bind_rows(est_coef, est_var, est_re) } else { est_raw <-bind_rows(est_coef, est_var) }# record estimation scale based on link (implied by likelihood) used scale <-case_when(family =="Gaussian"~"response", family %in%c("binomial", "beta") ~"logit", family %in%c("poisson", "nbinomial", "Gamma") ~"log")# combine estimates est <- est_raw |>mutate(scale = scale) |>select(effect, scale, term,est = estimate,se = std.error,lo = conf.low,hi = conf.high,pval = p.value)## extract fit fit <-glance(mod) |>bind_cols(loocv(mod)) |>mutate(method = est_method) |>select(n_obs = n,n_coef = p,n_varcomp = npar, method,aic = AIC,aicc = AICc,loglik = logLik,dev = deviance,dev_explain = pseudo.r.squared, bias,mse = MSPE,rmse = RMSPE,rav = RAV)## partition variance varpart <-varcomp(mod) |>mutate(effect ="varcomp",scale ="proportion") |>select(effect, scale,term = varcomp,est = proportion) |>mutate(term =ifelse(term =="Covariates (PR-sq)", "predictor", term))## output as dataframe out <-tibble_row(lab = lab, # model fit labelset = set, # model set labeloutcome = outcome, # outcome assessed labelresponse = response, # response variable representing outcome labelpredictor = predictor, # explanatory variable set labelerror = error, # covariance structure labelre = re, # random effect labelssn = path, # path to SSN used for modelformula =deparse1(form), # formula used in model as stringest =list(est), # model estimates as nested dataframefit =list(fit), # model fit metrics as nested dataframevarcomp =list(varpart) # model variance partitioning as nested dataframe )return(out)}```Define a formula for a random intercept by sampling event ID (as a categorical variable, nut numeric) to be reused across all models as the temporal random effect:```{r}#| label: form_reform_re <-formula(~1|as.factor(samp_event))```## Univariable modelsUnivariable models are fit separately for each environmental factor as the explanatory variable to assess associations with each of the three outcome variables (generic E. coli log10 concentration, ESBL-Ec log10 concentration, and CST-Ec detection). While it would be possible to write a second wrapper function to execute the model fits and estimate extractions for each of the three outcomes, there is just enough idiosyncrasy across the covariates in naming convention and other nuances, and few enough separate covariates, to defend hard-coding the univariable model for each combination of outcome and covariate.### Generic *E. coli*Start by estimating the univariable associations between the resistant outcomes (ESBL-Ec and CST-Ec) and generic *E. coli*. Apply the `fit_ssn()` function to the ESBL-Ec SSN assembled previously:```{r}#| label: ssn_uni_gen_esbl#| eval: falsem_gen_esbl <-fit_ssn(lab ="gen", # label for generic Ec-as-predictor modelsset ="uni", # univariable model setoutcome ="esbl", # ESBL-Ec outcomeresponse ="mpn", # concentration (MPN) continuous response varpredictor ="mpn_gen", # generic EC concentration (MPN) explanatory varerror ="flow", # flow-connected (tail-up) covariance structurere ="event", # random effect by sampling eventpath ="../gis/ssn/esbl_gen.ssn", # path to ESBL SSN with complete generic Ec dataform =formula(mpn ~ mpn_gen), # model formularandom = form_re, # random effect formulafamily ="Gaussian", # linear model (identity link) (defaut)additive ="afv_area", # cumulative area drained as flow proxy (default)tailup_type ="exponential", # YES tailup covariance with exponential modeltaildown_type ="none", # NO taildown covariance (default)euclid_type ="none", # NO Euclidean covariance component (default)nugget_type ="nugget", # YES nugget (unstructured covariance) (default)est_method ="reml") # use restricted maximum likelihood (default)```Repeat for CST-Ec, specifying detection (rather than concentration) as the response, binomial family for logistic regression, and dropping default arguments for simplicity:```{r}#| label: ssn_uni_gen_cst#| eval: falsem_gen_cst <-fit_ssn(lab ="gen", # label for generic Ec-as-predictor modelsset ="uni", # univariable model setoutcome ="cst", # CST-Ec outcomeresponse ="detect", # detection binary response varpredictor ="mpn_gen", # generic EC concentration (MPN) explanatory varerror ="flow", # flow-connected (tail-up) covariance structurere ="event", # random effect by sampling eventpath ="../gis/ssn/cst_gen.ssn", # path to ESBL SSN with complete generic Ec dataform =formula(detect ~ mpn_gen), # model formulafamily ="binomial", # binomial likelihood (logistic link)random = form_re, # random effect formulatailup_type ="exponential") # YES tailup covariance with exponential model```Stack the model outputs and save:```{r}#| label: ssn_uni_gen#| eval: falsem_gen <-bind_rows(m_gen_esbl, m_gen_cst)saveRDS(m_gen, "../est/ssn_uni_gen.rds")```### Water qualityModels using the co-measured water quality characteristics, turbidity and optical brighteners, also use the specific SSNs developed to account for missing covariate values, but are fit using generic *E. coli* as an outcome, as well as ESBL-Ec and CST-Ec.Fit models with turbidity as the explanatory variable for all three outcomes:```{r}#| label: ssn_uni_turb#| code-fold: true#| eval: false# ESBL-Ec concentration outcomem_turb_esbl <-fit_ssn(lab ="turb",set ="uni",outcome ="esbl",response ="mpn",predictor ="wq_turb",error ="flow",re ="event",path ="../gis/ssn/esbl_turb.ssn",form =formula(mpn ~ wq_turb),random = form_re,tailup_type ="exponential")# generic Ec concentration outcomem_turb_gen <-fit_ssn(lab ="turb",set ="uni",outcome ="gen",response ="mpn",predictor ="wq_turb",error ="flow",re ="event",path ="../gis/ssn/gen_turb.ssn",form =formula(mpn ~ wq_turb),random = form_re,tailup_type ="exponential")# CST-Ec detection outcomem_turb_cst <-fit_ssn(lab ="turb",set ="uni",outcome ="cst",response ="detect",predictor ="wq_turb",error ="flow",re ="event",path ="../gis/ssn/cst_turb.ssn",form =formula(detect ~ wq_turb),random = form_re,family ="binomial",tailup_type ="exponential")# combinem_turb <-bind_rows(m_turb_esbl, m_turb_gen, m_turb_cst)saveRDS(m_turb, "../est/ssn_uni_turb.rds")```Repeat for optical brighteners, using the appropriate prepared SSNs:```{r}#| label: ssn_uni_optic#| code-fold: true#| eval: false# ESBL-Ec concentration outcomem_optic_esbl <-fit_ssn(lab ="optic",set ="uni",outcome ="esbl",response ="mpn",predictor ="wq_optic",error ="flow",re ="event",path ="../gis/ssn/esbl_optic.ssn",form =formula(mpn ~ wq_optic),random = form_re,tailup_type ="exponential")# generic Ec concentration outcomem_optic_gen <-fit_ssn(lab ="optic",set ="uni",outcome ="gen",response ="mpn",predictor ="wq_optic",error ="flow",re ="event",path ="../gis/ssn/gen_optic.ssn",form =formula(mpn ~ wq_optic),random = form_re,tailup_type ="exponential")# CST-Ec detection outcomem_optic_cst <-fit_ssn(lab ="optic",set ="uni",outcome ="cst",response ="detect",predictor ="wq_optic",error ="flow",re ="event",path ="../gis/ssn/cst_optic.ssn",form =formula(detect ~ wq_optic),random = form_re,family ="binomial",tailup_type ="exponential")# combinem_optic <-bind_rows(m_optic_esbl, m_optic_gen, m_optic_cst)saveRDS(m_optic, "../est/ssn_uni_optic.rds")```### Land coverFit univariable models with basin features (land cover, imperviousness, and area) as the explanatory variable for each of the three outcomes. The full-observation SSNs are used because all basin feature covariate data are complete.Fit models for basin area:```{r}#| label: ssn_uni_area#| code-fold: true#| eval: false# ESBL-Ec concentration outcomem_area_esbl <-fit_ssn(lab ="area",set ="uni",outcome ="esbl",response ="mpn",predictor ="wshed_area",error ="flow",re ="event",path ="../gis/ssn/esbl.ssn",form =formula(mpn ~ wshed_area),random = form_re,tailup_type ="exponential")# generic Ec concentration outcomem_area_gen <-fit_ssn(lab ="area",set ="uni",outcome ="gen",response ="mpn",predictor ="wshed_area",error ="flow",re ="event",path ="../gis/ssn/gen.ssn",form =formula(mpn ~ wshed_area),random = form_re,tailup_type ="exponential")# CST-Ec detection outcomem_area_cst <-fit_ssn(lab ="area",set ="uni",outcome ="cst",response ="detect",predictor ="wshed_area",error ="flow",re ="event",path ="../gis/ssn/cst.ssn",form =formula(detect ~ wshed_area),random = form_re,family ="binomial",tailup_type ="exponential")# combinem_area <-bind_rows(m_area_esbl, m_area_gen, m_area_cst)saveRDS(m_area, "../est/ssn_uni_area.rds")```Fit models for imperviousness:```{r}#| label: ssn_uni_impv#| code-fold: true#| eval: false# ESBL-Ec concentration outcomem_impv_esbl <-fit_ssn(lab ="impv",set ="uni",outcome ="esbl",response ="mpn",predictor ="lc_impv",error ="flow",re ="event",path ="../gis/ssn/esbl.ssn",form =formula(mpn ~ lc_impv),random = form_re,tailup_type ="exponential")# generic Ec concentration outcomem_impv_gen <-fit_ssn(lab ="impv",set ="uni",outcome ="gen",response ="mpn",predictor ="lc_impv",error ="flow",re ="event",path ="../gis/ssn/gen.ssn",form =formula(mpn ~ lc_impv),random = form_re,tailup_type ="exponential")# CST-Ec detection outcomem_impv_cst <-fit_ssn(lab ="impv",set ="uni",outcome ="cst",response ="detect",predictor ="lc_impv",error ="flow",re ="event",path ="../gis/ssn/cst.ssn",form =formula(detect ~ lc_impv),random = form_re,family ="binomial",tailup_type ="exponential")# combinem_impv <-bind_rows(m_impv_esbl, m_impv_gen, m_impv_cst)saveRDS(m_impv, "../est/ssn_uni_impv.rds")```Fit models for forested land cover:```{r}#| label: ssn_uni_forest#| code-fold: true#| eval: false# ESBL-Ec concentration outcomem_forest_esbl <-fit_ssn(lab ="forest",set ="uni",outcome ="esbl",response ="mpn",predictor ="lc_forest",error ="flow",re ="event",path ="../gis/ssn/esbl.ssn",form =formula(mpn ~ lc_forest),random = form_re,tailup_type ="exponential")# generic Ec concentration outcomem_forest_gen <-fit_ssn(lab ="forest",set ="uni",outcome ="gen",response ="mpn",predictor ="lc_forest",error ="flow",re ="event",path ="../gis/ssn/gen.ssn",form =formula(mpn ~ lc_forest),random = form_re,tailup_type ="exponential")# CST-Ec detection outcomem_forest_cst <-fit_ssn(lab ="forest",set ="uni",outcome ="cst",response ="detect",predictor ="lc_forest",error ="flow",re ="event",path ="../gis/ssn/cst.ssn",form =formula(detect ~ lc_forest),random = form_re,family ="binomial",tailup_type ="exponential")# combinem_forest <-bind_rows(m_forest_esbl, m_forest_gen, m_forest_cst)saveRDS(m_forest, "../est/ssn_uni_forest.rds")```Fit models for open land cover:```{r}#| label: ssn_uni_open#| code-fold: true#| eval: false# ESBL-Ec concentration outcomem_open_esbl <-fit_ssn(lab ="open",set ="uni",outcome ="esbl",response ="mpn",predictor ="lc_open",error ="flow",re ="event",path ="../gis/ssn/esbl.ssn",form =formula(mpn ~ lc_open),random = form_re,tailup_type ="exponential")# generic Ec concentration outcomem_open_gen <-fit_ssn(lab ="open",set ="uni",outcome ="gen",response ="mpn",predictor ="lc_open",error ="flow",re ="event",path ="../gis/ssn/gen.ssn",form =formula(mpn ~ lc_open),random = form_re,tailup_type ="exponential")# CST-Ec detection outcomem_open_cst <-fit_ssn(lab ="open",set ="uni",outcome ="cst",response ="detect",predictor ="lc_open",error ="flow",re ="event",path ="../gis/ssn/cst.ssn",form =formula(detect ~ lc_open),random = form_re,family ="binomial",tailup_type ="exponential")# combinem_open <-bind_rows(m_open_esbl, m_open_gen, m_open_cst)saveRDS(m_open, "../est/ssn_uni_open.rds")```### MeteorologyFit models for rain in the previous two days:```{r}#| label: ssn_uni_rain2#| code-fold: true#| eval: false# ESBL-Ec concentration outcomem_rain2_esbl <-fit_ssn(lab ="rain2",set ="uni",outcome ="esbl",response ="mpn",predictor ="rain2",error ="flow",re ="event",path ="../gis/ssn/esbl.ssn",form =formula(mpn ~ rain2),random = form_re,tailup_type ="exponential")# generic Ec concentration outcomem_rain2_gen <-fit_ssn(lab ="rain2",set ="uni",outcome ="gen",response ="mpn",predictor ="rain2",error ="flow",re ="event",path ="../gis/ssn/gen.ssn",form =formula(mpn ~ rain2),random = form_re,tailup_type ="exponential")# CST-Ec detection outcomem_rain2_cst <-fit_ssn(lab ="rain2",set ="uni",outcome ="cst",response ="detect",predictor ="rain2",error ="flow",re ="event",path ="../gis/ssn/cst.ssn",form =formula(detect ~ rain2),random = form_re,family ="binomial",tailup_type ="exponential")# combinem_rain2 <-bind_rows(m_rain2_esbl, m_rain2_gen, m_rain2_cst)saveRDS(m_rain2, "../est/ssn_uni_rain2.rds")```Fit models for rain in the previous 3--7 days:```{r}#| label: ssn_uni_rain7#| code-fold: true#| eval: false# ESBL-Ec concentration outcomem_rain7_esbl <-fit_ssn(lab ="rain7",set ="uni",outcome ="esbl",response ="mpn",predictor ="rain7",error ="flow",re ="event",path ="../gis/ssn/esbl.ssn",form =formula(mpn ~ rain7),random = form_re,tailup_type ="exponential")# generic Ec concentration outcomem_rain7_gen <-fit_ssn(lab ="rain7",set ="uni",outcome ="gen",response ="mpn",predictor ="rain7",error ="flow",re ="event",path ="../gis/ssn/gen.ssn",form =formula(mpn ~ rain7),random = form_re,tailup_type ="exponential")# CST-Ec detection outcomem_rain7_cst <-fit_ssn(lab ="rain7",set ="uni",outcome ="cst",response ="detect",predictor ="rain7",error ="flow",re ="event",path ="../gis/ssn/cst.ssn",form =formula(detect ~ rain7),random = form_re,family ="binomial",tailup_type ="exponential")# combinem_rain7 <-bind_rows(m_rain7_esbl, m_rain7_gen, m_rain7_cst)saveRDS(m_rain7, "../est/ssn_uni_rain7.rds")```Fit models for previous day maximum temperature:```{r}#| label: ssn_uni_tmax1#| code-fold: true#| eval: false# ESBL-Ec concentration outcomem_tmax1_esbl <-fit_ssn(lab ="tmax1",set ="uni",outcome ="esbl",response ="mpn",predictor ="tmax1",error ="flow",re ="event",path ="../gis/ssn/esbl.ssn",form =formula(mpn ~ tmax1),random = form_re,tailup_type ="exponential")# generic Ec concentration outcomem_tmax1_gen <-fit_ssn(lab ="tmax1",set ="uni",outcome ="gen",response ="mpn",predictor ="tmax1",error ="flow",re ="event",path ="../gis/ssn/gen.ssn",form =formula(mpn ~ tmax1),random = form_re,tailup_type ="exponential")# CST-Ec detection outcomem_tmax1_cst <-fit_ssn(lab ="tmax1",set ="uni",outcome ="cst",response ="detect",predictor ="tmax1",error ="flow",re ="event",path ="../gis/ssn/cst.ssn",form =formula(detect ~ tmax1),random = form_re,family ="binomial",tailup_type ="exponential")# combinem_tmax1 <-bind_rows(m_tmax1_esbl, m_tmax1_gen, m_tmax1_cst)saveRDS(m_tmax1, "../est/ssn_uni_tmax1.rds")```### WastewaterFit models for WWTP treated effluent loads:```{r}#| label: ssn_uni_wwtp#| code-fold: true#| eval: false# ESBL-Ec concentration outcomem_wwtp_esbl <-fit_ssn(lab ="wwtp",set ="uni",outcome ="esbl",response ="mpn",predictor ="out_wwtp_load",error ="flow",re ="event",path ="../gis/ssn/esbl.ssn",form =formula(mpn ~ out_wwtp_load),random = form_re,tailup_type ="exponential")# generic Ec concentration outcomem_wwtp_gen <-fit_ssn(lab ="wwtp",set ="uni",outcome ="gen",response ="mpn",predictor ="out_wwtp_load",error ="flow",re ="event",path ="../gis/ssn/gen.ssn",form =formula(mpn ~ out_wwtp_load),random = form_re,tailup_type ="exponential")# CST-Ec detection outcomem_wwtp_cst <-fit_ssn(lab ="wwtp",set ="uni",outcome ="cst",response ="detect",predictor ="out_wwtp_load",error ="flow",re ="event",path ="../gis/ssn/cst.ssn",form =formula(detect ~ out_wwtp_load),random = form_re,family ="binomial",tailup_type ="exponential")# combinem_wwtp <-bind_rows(m_wwtp_esbl, m_wwtp_gen, m_wwtp_cst)saveRDS(m_wwtp, "../est/ssn_uni_wwtp.rds")```Fit models for CSO effluent loads:```{r}#| label: ssn_uni_cso#| code-fold: true#| eval: false# ESBL-Ec concentration outcomem_cso_esbl <-fit_ssn(lab ="cso",set ="uni",outcome ="esbl",response ="mpn",predictor ="out_cso_load",error ="flow",re ="event",path ="../gis/ssn/esbl.ssn",form =formula(mpn ~ out_cso_load),random = form_re,tailup_type ="exponential")# generic Ec concentration outcomem_cso_gen <-fit_ssn(lab ="cso",set ="uni",outcome ="gen",response ="mpn",predictor ="out_cso_load",error ="flow",re ="event",path ="../gis/ssn/gen.ssn",form =formula(mpn ~ out_cso_load),random = form_re,tailup_type ="exponential")# CST-Ec detection outcomem_cso_cst <-fit_ssn(lab ="cso",set ="uni",outcome ="cst",response ="detect",predictor ="out_cso_load",error ="flow",re ="event",path ="../gis/ssn/cst.ssn",form =formula(detect ~ out_cso_load),random = form_re,family ="binomial",tailup_type ="exponential")# combinem_cso <-bind_rows(m_cso_esbl, m_cso_gen, m_cso_cst)saveRDS(m_cso, "../est/ssn_uni_cso.rds")```Fit models for septic density:```{r}#| label: ssn_uni_dens#| code-fold: true#| eval: false# ESBL-Ec concentration outcomem_dens_esbl <-fit_ssn(lab ="dens",set ="uni",outcome ="esbl",response ="mpn",predictor ="septic_dens",error ="flow",re ="event",path ="../gis/ssn/esbl.ssn",form =formula(mpn ~ septic_dens),random = form_re,tailup_type ="exponential")# generic Ec concentration outcomem_dens_gen <-fit_ssn(lab ="dens",set ="uni",outcome ="gen",response ="mpn",predictor ="septic_dens",error ="flow",re ="event",path ="../gis/ssn/gen.ssn",form =formula(mpn ~ septic_dens),random = form_re,tailup_type ="exponential")# CST-Ec detection outcomem_dens_cst <-fit_ssn(lab ="dens",set ="uni",outcome ="cst",response ="detect",predictor ="septic_dens",error ="flow",re ="event",path ="../gis/ssn/cst.ssn",form =formula(detect ~ septic_dens),random = form_re,family ="binomial",tailup_type ="exponential")# combinem_dens <-bind_rows(m_dens_esbl, m_dens_gen, m_dens_cst)saveRDS(m_dens, "../est/ssn_uni_dens.rds")```Although univariable models are not applicable for the interaction of septic density and rainfall (which also includes terms for the main effects of septic density and rainfall), fit a model for the interaction in the "univariable" set for consistency with the later multivariable estimates:```{r}#| label: ssn_uni_raindens#| code-fold: true#| eval: false# ESBL-Ec concentration outcomem_raindens_esbl <-fit_ssn(lab ="raindens",set ="uni",outcome ="esbl",response ="mpn",predictor ="interaction",error ="flow",re ="event",path ="../gis/ssn/esbl.ssn",form =formula(mpn ~ septic_dens*rain7),random = form_re,tailup_type ="exponential")# generic Ec concentration outcomem_raindens_gen <-fit_ssn(lab ="raindens",set ="uni",outcome ="gen",response ="mpn",predictor ="interaction",error ="flow",re ="event",path ="../gis/ssn/gen.ssn",form =formula(mpn ~ septic_dens*rain7),random = form_re,tailup_type ="exponential")# CST-Ec detection outcomem_raindens_cst <-fit_ssn(lab ="raindens",set ="uni",outcome ="cst",response ="detect",predictor ="interaction",error ="flow",re ="event",path ="../gis/ssn/cst.ssn",form =formula(detect ~ septic_dens*rain7),random = form_re,family ="binomial",tailup_type ="exponential")# combinem_raindens <-bind_rows(m_raindens_esbl, m_raindens_gen, m_raindens_cst)saveRDS(m_raindens, "../est/ssn_uni_raindens.rds")```## Multivariable modelsMultivariable models are also fit that include all environmental factors except for those measured in the water samples (i.e., excluding generic *E. coli*, turbidity, and optical brighteners) to explore the independent contribution of each external environmental factor conditional on the others. To account for potential subsurface transport of septic effluent mobilized by antecedent precipitation, include an interaction term in the multivariable models for the product of septic density and 3–7 day rainfall. Exclude the basin area and open land cover-class variables from the multivariable models to reduce collinearity among the basin-feature variable set.Use the same `fit_ssn()` function to fit a multivariable model to the ESBL-Ec concentration data:```{r}#| label: ssn_multi_esbl#| eval: false# ESBL-Ec concentration outcomem_multi_esbl <-fit_ssn(lab ="multi",set ="multi",outcome ="esbl",response ="mpn",predictor ="external",error ="flow",re ="event",path ="../gis/ssn/esbl.ssn",form =formula(mpn ~ out_wwtp_load + out_cso_load + lc_impv + lc_forest + rain2 + tmax1 + septic_dens*rain7),random = form_re,tailup_type ="exponential")```Repeat for generic *E. coli* concentration and CST-Ec detection, then stack the estimates and save:```{r}#| label: ssn_multi#| code-fold: true#| eval: false# generic Ec concentration outcomem_multi_gen <-fit_ssn(lab ="multi",set ="multi",outcome ="gen",response ="mpn",predictor ="external",error ="flow",re ="event",path ="../gis/ssn/gen.ssn",form =formula(mpn ~ out_wwtp_load + out_cso_load + lc_impv + lc_forest + rain2 + tmax1 + septic_dens*rain7),random = form_re,tailup_type ="exponential")# CST-Ec detection outcomem_multi_cst <-fit_ssn(lab ="multi",set ="multi",outcome ="cst",response ="detect",predictor ="external",error ="flow",re ="event",path ="../gis/ssn/cst.ssn",form =formula(detect ~ out_wwtp_load + out_cso_load + lc_impv + lc_forest + rain2 + tmax1 + septic_dens*rain7),random = form_re,family ="binomial",tailup_type ="exponential")# combinem_multi <-bind_rows(m_multi_esbl, m_multi_gen, m_multi_cst)saveRDS(m_multi, "../est/ssn_multi.rds")```## Model estimatesRead and stack the univariable model output files:```{r}#| label: mod_stack#| eval: false#| error: false#| warning: falsem_uni <-bind_rows(readRDS("../est/ssn_uni_gen.rds"),readRDS("../est/ssn_uni_turb.rds"),readRDS("../est/ssn_uni_optic.rds"),readRDS("../est/ssn_uni_area.rds"),readRDS("../est/ssn_uni_impv.rds"),readRDS("../est/ssn_uni_forest.rds"),readRDS("../est/ssn_uni_open.rds"),readRDS("../est/ssn_uni_rain2.rds"),readRDS("../est/ssn_uni_rain7.rds"),readRDS("../est/ssn_uni_tmax1.rds"),readRDS("../est/ssn_uni_wwtp.rds"),readRDS("../est/ssn_uni_cso.rds"),readRDS("../est/ssn_uni_dens.rds"),readRDS("../est/ssn_uni_raindens.rds"))saveRDS(m_uni, "../est/ssn_mod_uni.rds")```Drop the nested model fit and variance component dataframes and unnest the univariable model coefficient estimates:```{r}#| label: est_uni#| eval: falseest_uni <- m_uni |>select(-fit, -varcomp) |>unnest(est)saveRDS(est_uni, "../est/ssn_est_uni.rds")```Separately unnest the model fit and variance component estimates:```{r}#| label: unnest_uni#| code-fold: true#| eval: false# model fit metricsfit_uni <- m_uni |>select(-est, -varcomp) |>unnest(fit)saveRDS(fit_uni, "../est/ssn_fit_uni.rds")# variance partitioning estimatesvar_uni <- m_uni |>select(-est, - fit) |>unnest(varcomp)saveRDS(var_uni, "../est/ssn_var_uni.rds")```Likewise, extract the estimates, model fits, and variance components for the multivariable models:```{r}#| label: unnest_multi#| eval: false#| code-fold: true# multivariable model outputm_multi <-readRDS("../est/ssn_multi.rds")# coefficient estimatesest_multi <- m_multi |>select(-fit, -varcomp) |>unnest(est)saveRDS(est_multi, "../est/ssn_est_multi.rds")# model fit metricsfit_multi <- m_multi |>select(-est, -varcomp) |>unnest(fit)saveRDS(fit_multi, "../est/ssn_fit_multi.rds")# variance partitioning estimatesvar_multi <- m_multi |>select(-est, - fit) |>unnest(varcomp)saveRDS(var_multi, "../est/ssn_var_multi.rds")```### Environmental factor associationsFor convenience, a key mapping the shorthand variable names used in the analysis to their full display labels, along with grouping and ordering information, is provided in [`covar_labels.xlsx`](https://cdcgov.github.io/EMEL/manuscripts/river_ar/covar_labels.xlsx). Read in the labels from the key and the means and standard deviations previously used to scale the variables:```{r}#| label: coef_key# covariate labelsvar_lab <-read_xlsx("../data/covar_labels.xlsx")# covariate scaling factorsvar_scl <-readRDS("../out/covar_summ_stats.rds") |>mutate(outcome =case_when(drug =="cefotaxime"~"esbl", drug =="colistin"~"cst", drug =="none"~"gen"))```Pull in the univariable model coefficient estimates, adding empty rows for the generic *E. coli* outcome with generic *E. coli* as the explanatory variable and setting the septic-rainfall interaction term to missing (and dropping the associated main effects) for display purposes when merged with the multivariable model estimates:```{r}#| label: coef_unicoef_uni <-readRDS("../est/ssn_est_uni.rds") |># empty row for generic Ec predicted by generic Ecadd_case(outcome ="gen",set ="uni",effect ="coef",term ="mpn_gen",predictor ="mpn_gen") |># set septic-rainfall interaction to NA and drop septic and rainfal main effects from interaction modelmutate(across(c(est, se, lo, hi, pval),~ifelse(term =="septic_dens:rain7", NA, .x))) |>filter(!(lab =="raindens"& term %in%c("septic_dens", "rain7")))```Obtain a list of unique effects estimated in univariable models, join the multivariable model estimates to match, stack with univariable estimates, and join to display labels and scaling factors:```{r}#| label: coef_allcoef_all <- coef_uni |># identify unique outcome-predictor combinations from univariable models select(outcome, response, effect, term) |>distinct() |>mutate(set ="multi") |># join matching multivariable model estimatesleft_join(readRDS("../est/ssn_est_multi.rds"),by =c("outcome", "response", "effect", "term", "set")) |># stack with univariable estimatesbind_rows(coef_uni) |># label with display names, grouping infoleft_join(var_lab, by =c("effect", "term")) |># join covariate scaling factorsleft_join(var_scl, by =c("outcome", "term"="var"))```Filter for just the main effect estimates (i.e., dropping intercepts and variance component hyperparameters), calculate descaled estimate values, exponentiate estimates from logistic models to obtain odds ratios (OR), and re-order for display:```{r}#| label: coef_slctcoef_slct <- coef_all |># filter main effect estimates (identified in covar_label key)filter(effect =="coef", include =="yes") %>%mutate(# descale coef estimates to obtain on original covariate scaleest_descl =ifelse(scl_method =="none", est, est / sd),lo_descl =ifelse(scl_method =="none", lo, lo / sd),hi_descl =ifelse(scl_method =="none", hi, hi / sd),# exponentiate logistic model estimates for odds ratiosacross(c(est, lo, hi, ends_with("_descl")),~ifelse(scale =="logit", exp(.x), .x)),# order by covariate groupgroup =fct_relevel(group,c("wq", "out", "lc", "meteo")),set =fct_relevel(set, c("uni", "multi")) |>fct_recode("Univariable"="uni","Multivariable"="multi") |>fct_rev()) |>arrange(group, term_order) |>mutate(group_lab =fct_inorder(group_lab),term_lab =fct_inorder(term_lab),drug_lab =case_when(outcome =="esbl"~"ESBL-Ec<br>concentration", outcome =="cst"~"CST-Ec<br>detection", outcome =="gen"~"Generic *E. coli*<br>concentration"))```Display the univariable and multivariable coefficient estimates together as the point estimates and their 95% confidence intervals. Prepare a similar forest plot for each of the three *E. coli* outcomes and stack the plots together for display, beginning with the generic *E. coli* concentration outcome:```{r}#| label: p_coef_genp_coef_gen <- coef_slct |>filter(outcome =="gen") %>%ggplot(data = .) +facet_grid(cols =vars(drug_lab),rows =vars(group_lab),drop =TRUE,scales ="free_y",space ="free_y") +geom_vline(aes(xintercept =0),linetype ="dashed") +geom_pointrange(aes(x = est, y = term_lab,xmin = lo, xmax = hi,color = group_lab, fill = group_lab,shape = set, linetype = set),size =0.6, linewidth =0.9,position =position_dodge(width =0.75)) +coord_cartesian(xlim =c(-0.75, 1)) +scale_y_discrete(limits = rev) +scale_shape_manual(values =c("Univariable"=21, "Multivariable"=23),name =NULL) +scale_linetype_manual(values =c("Univariable"="solid", "Multivariable"="11"),name =NULL) +scale_color_brewer(palette ="Set2") +scale_fill_brewer(palette ="Set2") +guides(color ="none", fill ="none",shape =guide_legend(reverse =TRUE,override.aes =list(fill ="black")),linetype =guide_legend(reverse =TRUE)) +labs(x ="Δ log~10~ MPN/100 mL",y =NULL) +theme_bw() +theme(axis.text =element_markdown(),axis.title =element_markdown(),strip.text.y =element_blank())```Repeat for the ESBL-Ec concentration outcome, dropping axis ticks and related decorations:```{r}#| label: p_coef_esbl#| code-fold: truep_coef_esbl <- coef_slct |>filter(outcome =="esbl") %>%ggplot(data = .) +facet_grid(cols =vars(drug_lab),rows =vars(group_lab),drop =TRUE,scales ="free_y",space ="free_y") +geom_vline(aes(xintercept =0),linetype ="dashed") +geom_pointrange(aes(x = est, y = term_lab,xmin = lo, xmax = hi,color = group_lab, fill = group_lab,shape = set, linetype = set),size =0.6, linewidth =0.9,position =position_dodge(width =0.75)) +coord_cartesian(xlim =c(-0.75, 1)) +scale_y_discrete(limits = rev) +scale_shape_manual(values =c("Univariable"=21, "Multivariable"=23),name =NULL) +scale_linetype_manual(values =c("Univariable"="solid", "Multivariable"="11"),name =NULL) +scale_color_brewer(palette ="Set2") +scale_fill_brewer(palette ="Set2") +guides(color ="none", fill ="none",shape =guide_legend(reverse =TRUE,override.aes =list(fill ="black")),linetype =guide_legend(reverse =TRUE)) +labs(x ="Δ log~10~ MPN/100 mL",y =NULL) +theme_bw() +theme(axis.text.y =element_blank(),axis.ticks.y =element_blank(),axis.title =element_markdown(),strip.text.y =element_blank())```Finally, plot odds ratio estimates for CST-Ec detection, adjusting the axes and labels as appropriate:```{r}#| label: p_coef_cst#| code-fold: truep_coef_cst <- coef_slct |>filter(outcome =="cst") %>%ggplot(data = .) +facet_grid(cols =vars(drug_lab),rows =vars(group_lab),drop =TRUE,scales ="free_y",space ="free_y") +geom_vline(aes(xintercept =1),linetype ="dashed") +geom_pointrange(aes(x = est, y = term_lab,xmin = lo, xmax = hi,color = group_lab, fill = group_lab,shape = set, linetype = set),size =0.6, linewidth =0.9,position =position_dodge(width =0.75)) +scale_x_log10() +scale_y_discrete(limits = rev) +scale_shape_manual(values =c("Univariable"=21, "Multivariable"=23),name =NULL) +scale_linetype_manual(values =c("Univariable"="solid", "Multivariable"="11"),name =NULL) +scale_color_brewer(palette ="Set2") +scale_fill_brewer(palette ="Set2") +guides(color ="none", fill ="none",shape =guide_legend(reverse =TRUE,override.aes =list(fill ="black")),linetype =guide_legend(reverse =TRUE)) +labs(x ="odds ratio",y =NULL) +theme_bw() +theme(axis.text.y =element_blank(),axis.ticks.y =element_blank(),axis.title =element_markdown())```Arrange the three forest plots together and save the image as a .png with specific dimensions to ensure proper display:```{r}#| label: p_coef#| warning: false#| message: false# arrange plotsp_coef <- p_coef_gen + p_coef_esbl + p_coef_cst +plot_layout(ncol =3,guides ="collect") &theme(legend.key.width =unit(1.5, "lines"),legend.key =element_rect(fill ="white", colour =NA),legend.position ="bottom",legend.direction ="horizontal",strip.text =element_markdown(),panel.border =element_rect(colour ="black", fill=NA),panel.spacing =unit(0, "lines"),,legend.margin =margin(c(1, 1, 1, 1)),plot.margin =margin(c(2, 1, 1, 1)))# save PNG imageggsave(filename ="../fig/coef_multi.png",plot = p_coef,device = agg_png, dpi =300,width =9, height =8, units ="in")```This reproduces **Figure 3** from Holcomb et al. (2026). Note that the estimates are for the *scaled* covariates and thus correspond to the expected change in the outcome for a **1 standard deviation increase** in the environmental factor:Generate an alternative version of the forest plot with the descaled coefficient estimates that correspond to to the expected change in each outcome for a 1-unit increase in the environmental factor on its original measurement scale:```{r}#| label: p_coef_descl#| code-fold: true#| warning: false#| message: false# generic E. colip_descl_gen <- coef_slct |>filter(outcome =="gen") %>%ggplot(data = .) +facet_grid(cols =vars(drug_lab),rows =vars(group_lab),drop =TRUE,scales ="free_y",space ="free_y") +geom_vline(aes(xintercept =0),linetype ="dashed") +geom_pointrange(aes(x = est_descl, y = term_lab,xmin = lo_descl, xmax = hi_descl,color = group_lab, fill = group_lab,shape = set, linetype = set),size =0.6, linewidth =0.9,position =position_dodge(width =0.75)) +scale_y_discrete(limits = rev) +scale_shape_manual(values =c("Univariable"=21, "Multivariable"=23),name =NULL) +scale_linetype_manual(values =c("Univariable"="solid", "Multivariable"="11"),name =NULL) +scale_color_brewer(palette ="Set2") +scale_fill_brewer(palette ="Set2") +guides(color ="none", fill ="none",shape =guide_legend(reverse =TRUE,override.aes =list(fill ="black")),linetype =guide_legend(reverse =TRUE)) +labs(x ="Δ log~10~ MPN/100 mL",y =NULL) +theme_bw() +theme(axis.text =element_markdown(),axis.title =element_markdown(),strip.text.y =element_blank())# ESBL-Ecp_descl_esbl <- coef_slct |>filter(outcome =="esbl") %>%ggplot(data = .) +facet_grid(cols =vars(drug_lab),rows =vars(group_lab),drop =TRUE,scales ="free_y",space ="free_y") +geom_vline(aes(xintercept =0),linetype ="dashed") +geom_pointrange(aes(x = est_descl, y = term_lab,xmin = lo_descl, xmax = hi_descl,color = group_lab, fill = group_lab,shape = set, linetype = set),size =0.6, linewidth =0.9,position =position_dodge(width =0.75)) +scale_y_discrete(limits = rev) +scale_shape_manual(values =c("Univariable"=21, "Multivariable"=23),name =NULL) +scale_linetype_manual(values =c("Univariable"="solid", "Multivariable"="11"),name =NULL) +scale_color_brewer(palette ="Set2") +scale_fill_brewer(palette ="Set2") +guides(color ="none", fill ="none",shape =guide_legend(reverse =TRUE,override.aes =list(fill ="black")),linetype =guide_legend(reverse =TRUE)) +labs(x ="Δ log~10~ MPN/100 mL",y =NULL) +theme_bw() +theme(axis.text.y =element_blank(),axis.ticks.y =element_blank(),axis.title =element_markdown(),strip.text.y =element_blank())# CST-Ecp_descl_cst <- coef_slct |>filter(outcome =="cst") %>%ggplot(data = .) +facet_grid(cols =vars(drug_lab),rows =vars(group_lab),drop =TRUE,scales ="free_y",space ="free_y") +geom_vline(aes(xintercept =1),linetype ="dashed") +geom_pointrange(aes(x = est_descl, y = term_lab,xmin = lo_descl, xmax = hi_descl,color = group_lab, fill = group_lab,shape = set, linetype = set),size =0.6, linewidth =0.9,position =position_dodge(width =0.75)) +scale_x_log10() +scale_y_discrete(limits = rev) +scale_shape_manual(values =c("Univariable"=21, "Multivariable"=23),name =NULL) +scale_linetype_manual(values =c("Univariable"="solid", "Multivariable"="11"),name =NULL) +scale_color_brewer(palette ="Set2") +scale_fill_brewer(palette ="Set2") +guides(color ="none", fill ="none",shape =guide_legend(reverse =TRUE,override.aes =list(fill ="black")),linetype =guide_legend(reverse =TRUE)) +labs(x ="odds ratio",y =NULL) +theme_bw() +theme(axis.text.y =element_blank(),axis.ticks.y =element_blank(),axis.title =element_markdown())## combinep_descl <- p_descl_gen + p_descl_esbl + p_descl_cst +plot_layout(ncol =3,guides ="collect") &theme(legend.key.width =unit(1.5, "lines"),legend.key =element_rect(fill ="white", colour =NA),legend.position ="bottom",legend.direction ="horizontal",strip.text =element_markdown(),panel.border =element_rect(colour ="black", fill=NA),panel.spacing =unit(0, "lines"),,legend.margin =margin(c(1, 1, 1, 1)),plot.margin =margin(c(2, 1, 1, 1)))# save PNG imageggsave(filename ="../fig/coef_descl.png",plot = p_descl,device = agg_png, dpi =300,width =9, height =8, units ="in")```This alternative version of Figure 3 with coefficient estimates on the original scale of the covariate data makes comparisons between different environmental factors difficult, but does serve to highlight the large magnitude of the associations between resistant *E. coli* and generic *E. coli* concentrations (in units of log~10~ MPN/100 mL), as well as recent rainfall (in units of cm):### Variance componentsExamine the proportion of the total variance in the *E. coli* outcomes explained by each of the model structural components in the multivariable models, as well as the estimated spatial range (km) for the flow-connected spatial random effect.Join the variable label key used for the coefficents to the previously extracted variance component estimates and filter for the variables of interest:```{r}#| label: varcompvarcomp <-readRDS("../est/ssn_var_multi.rds") |>left_join(var_lab, by =c("effect", "term")) |>filter(include =="yes") |>select(outcome, set, effect, term_order, term, term_lab, est)```Extract the spatial range estimates:```{r}#| label: rangerange <-readRDS("../est/ssn_est_multi.rds") |>filter(effect %in%c("tailup", "nugget", "re")) |>left_join(var_lab, by =c("effect", "term")) |>filter(include =="yes") |>mutate(effect ="est") |>select(outcome, set, effect, term_order, term, term_lab, est)```Stack and reformat for presentation:```{r}#| label: tab_varcomptab_varcomp <- range |>bind_rows(varcomp) |>mutate(outcome =fct_relevel(outcome, c("gen", "esbl", "cst")) |>fct_recode("Generic Ec concentration"="gen","ESBL-Ec concentration"="esbl","CST-Ec detection"="cst"),term_lab =fct_relevel(term_lab,c("Residual variance","Explanatory variables","Sampling event variance","Flow-connected variance","Flow-connected range"))) |>arrange(outcome, term_lab, set) |>filter(effect =="varcomp"| term =="range") |>mutate(est =ifelse(term =="range", est /1000, est *100)) |>select(outcome, term_lab, est) |>pivot_wider(names_from = term_lab, values_from = est)```This essentially reproduces **Table S9** from the Supplementary Material of Holcomb et al. (2026):```{r}#| label: kable_varcomp#| echo: falsetab_varcomp |>kable(col.names =c("Response","Nugget","Explanatory variables","Sampling event random effect","Flow-connected spatial random effect","Flow-connected spatial range (km)"),digits =c(0, 1, 1, 1, 1, 0)) |>add_header_above(header =c(" "=1,"Variance explained (%)"=4," "=1))```<br /><br />------------------------------------------------------------------------**Disclaimer:**The findings and conclusions of this report are those of the authors and do not necessarily represent the official position of the Centers for Disease Control and Prevention.
