Skip to content

Latest commit

 

History

History
365 lines (317 loc) · 11.7 KB

ex0618.md

File metadata and controls

365 lines (317 loc) · 11.7 KB

QMRA, 2nd ed., Example 6.18

Brian High and John Kissel

Introduction

This document offers a solution in R for Example 6.18 from pages 215-216 of:

Quantitative Microbial Risk Assessment, 2nd Edition by Charles N. Haas, Joan B. Rose, and Charles P. Gerba. (Wiley, 2014).

This is the copyright statement for the book:

© Haas, Charles N.; Rose, Joan B.; Gerba, Charles P., Jun 02, 2014, Quantitative Microbial Risk Assessment Wiley, Somerset, ISBN: 9781118910528

The data for this example comes from this book, but the R code presented below is an original work, released into the public domain.

Set display options

# Set display options for use with the print() function.
options(digits=3)

Enter input data

Store the data provided in the textbook into one vector per variable.

# Create a data frame with input data.
exposure.source <- c("drinking water", "swimming", "shellfish")
viral.load <- c(0.001, 0.1, 1)
load.units <- c("viruses/L", "viruses/L", "viruses/g")
IR <- c(NA, 50/1000, NA)            # units: L/hr; "IR" = "ingestion rate"
duration <- c(NA, 2.6, NA)          # units: hrs
cons <- c(NA, NA, 150)              # units: L or g
freq <- c(NA, 7, .0009*365)         # units: days/yr
avg.daily.IR <- c(1.4, NA, NA)      # units: L/day or g/day; see: load.units

Combine the vectors into a data frame, where each vector becomes a column.

df <- data.frame(exposure.source, viral.load, load.units, IR, duration, 
                 cons, freq, avg.daily.IR)

View the input data

# Display the dataset, using scientific notation for the viral load.
df$v.load.sn <- format(df$viral.load, scientific = TRUE)
knitr::kable(df[, c(1, 9, 3:8)], format = "html", 
             table.attr = "style='width:100%;'")
exposure.source v.load.sn load.units IR duration cons freq avg.daily.IR
drinking water 1e-03 viruses/L NA NA NA NA 1.4
swimming 1e-01 viruses/L 0.05 2.6 NA 7.000 NA
shellfish 1e+00 viruses/g NA NA 150 0.328 NA

Estimate amount consumed

Before estimating the average daily ingestion rate, estimate the amount of medium which has been consumed per occurrence. Calculate the estimated consumption for those exposure sources which are missing both the ingestion rate and the consumption.

# Function: Calculate the consumption (amount of medium consumed) per occurrence.
calc.cons.per.occur <- function(IR, duration) {
    # Multiply the IR in L/hr by the duration in hrs.
    IR * duration
}

# Calculate the consumption for those exposure sources missing this value.
missing.cons <- is.na(df$avg.daily.IR) & is.na(df$cons)
df[missing.cons, "cons"] <- with(df[missing.cons, ], 
                                 calc.cons.per.occur(IR, duration))
knitr::kable(df[, c(1, 9, 3:8)], format = "html", 
             table.attr = "style='width:100%;'")
exposure.source v.load.sn load.units IR duration cons freq avg.daily.IR
drinking water 1e-03 viruses/L NA NA NA NA 1.4
swimming 1e-01 viruses/L 0.05 2.6 0.13 7.000 NA
shellfish 1e+00 viruses/g NA NA 150.00 0.328 NA

Estimate average daily ingestion rate

Calculate the estimated average daily ingestion rate for those exposure sources which are missing it.

# Function: Calculate the avg. daily IR given the consumption and frequency.
calc.avg.daily.IR <- function(consumption, frequency) {
    # Multiply the consumption (L or g) and frequency in days/yr and 
    # then divide by 365 days/year.
    consumption * frequency / 365
}

# Calculate the average daily IR for those exposure sources missing this value.
missing.IR <- is.na(df$avg.daily.IR)
df[missing.IR, "avg.daily.IR"] <- with(df[missing.IR, ], 
                                       calc.avg.daily.IR(cons, freq))
knitr::kable(df[, c(1, 9, 3:8)], format = "html", 
             table.attr = "style='width:100%;'")
exposure.source v.load.sn load.units IR duration cons freq avg.daily.IR
drinking water 1e-03 viruses/L NA NA NA NA 1.400
swimming 1e-01 viruses/L 0.05 2.6 0.13 7.000 0.002
shellfish 1e+00 viruses/g NA NA 150.00 0.328 0.135

Estimate exposures

Estimate the daily exposure for each source by multiplying the viral load by the average daily ingestion rate.

# Function: Calculate the dose given the viral load and average daily IR.
calc.dose <- function(viral.load, avg.daily.IR) {
    # Multiply the viral load in viruses/L or viruses/g by the 
    # average daily IR in L/day or g/day, respectively.
    viral.load * avg.daily.IR
}

# Calculate the average daily dose of viruses per day by exposure source.
df$dose <- with(df, calc.dose(viral.load, avg.daily.IR))
df$dose.sn <- format(df$dose, scientific = TRUE)  # Use scientific notation.
df$dose.units <- "viruses/day"
knitr::kable(df[, c(1, 9, 3, 8, 11, 12)], format = "html", 
             table.attr = "style='width:80%;'")
exposure.source v.load.sn load.units avg.daily.IR dose.sn dose.units
drinking water 1e-03 viruses/L 1.400 1.40e-03 viruses/day
swimming 1e-01 viruses/L 0.002 2.49e-04 viruses/day
shellfish 1e+00 viruses/g 0.135 1.35e-01 viruses/day

Estimate total exposure

Estimate the total daily exposure by calculating the sum of the estimated exposures for all of the sources.

# Calculate the total average daily dose for all exposure sources combined.
daily.dose <- sum(df$dose)
knitr::kable(format(daily.dose, scientific = TRUE), format = "html")
1.37e-01

Estimate exposures as fraction of total

Estimate the fraction of the total daily exposure due to each source by dividing the estimated exposure per source by the total estimated daily exposure.

# Calculate the average daily dose by source as a fraction of the total.
df$fraction <- round(df$dose / daily.dose, digits = 3)
knitr::kable(df[, c(1, 13)], format = "html", table.attr = "style='width:30%;'")
exposure.source fraction
drinking water 0.010
swimming 0.002
shellfish 0.988