Jane Pouzou and Brian High
This document offers a 1-D Monte Carlo probabilistic solution in R for the daily microbial exposure from drinking water consumption, swimming in surface water and shellfish consumption 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 under a Creative Commons Attribution-ShareAlike 4.0 International License. Details may be found at the end of this document.
# Set knitr options for use when rendering this document.
library("knitr")
opts_chunk$set(cache=TRUE, message=FALSE)Define variables provided in the example for three exposure types.
# Shellfish consumption
sf.viral.load <- 1
sf.cons.g <- 9e-4 * 150 # 9e-4 days/year * 150 g/day
# Drinking water consumption
dw.viral.load <- 0.001
# Surface water consumption while swimming
sw.viral.load <- 0.1
sw.daily.IR <- 50 # Ingestion rate in mL of surface water
sw.frequency <- 7 # Exposure frequency of 7 swims per yearSample from the probablity distributions for drinking water and swim duration. Also plot from these distributions as a quick visual check.
Generate 5000 random values from a log-normal distribution to estimate exposure from consumption of drinking water (ml/day). Divide by 1000 mL/L to get consumption in liters/day. Values for meanlog and sdlog are from the QMRA textbook (Haas, 2014), page 216, Table 6.30.
set.seed(1)
dw.cons.L <- rlnorm(5000, meanlog = 7.49, sdlog = 0.407) / 1000Plot the kernal density curve of the generated values just as a check.
plot(density(dw.cons.L))
Sample 5000 times from a discrete distribution of swim duration with assigned probabilities of each outcome. These values are hypothetical and are not found in the text, but are defined here to provide an example of sampling from a discrete distribution.
set.seed(1)
swim.duration <- sample(x = c(0.5, 1, 2, 2.6), 5000, replace = TRUE,
prob = c(0.1, 0.1, 0.2, 0.6))Create a simple histogram of our distribution as a check.
hist(swim.duration)
Calculate estimated daily dose using a probabilistic simulation model.
Define a function to calculate microbial exposure risk.
Risk.fcn <- function(sf.vl, sf.cons.g, dw.cons.L, dw.vl, sw.vl,
sw.daily.IR, sw.duration, sw.frequency) {
return(((sf.vl * sf.cons.g) + (dw.cons.L * dw.vl) +
((sw.vl * (sw.daily.IR * sw.duration * sw.frequency)) / 365 / 1000)))
}Compute 5000 simulated daily dose results and store as a vector.
# First compute the simulation using sapply().
daily.dose <- sapply(1:5000,
function(i) Risk.fcn(dw.cons.L = dw.cons.L[i],
sw.duration = swim.duration[i],
sf.vl = sf.viral.load,
dw.vl = dw.viral.load,
sf.cons.g = sf.cons.g,
sw.vl = sw.viral.load,
sw.daily.IR = sw.daily.IR,
sw.frequency = sw.frequency))
# Second, just for comparison, compute the simulation using a for loop.
daily.dose2 <- as.vector(NULL)
for (i in 1:5000) {
daily.dose2[i] <- Risk.fcn(dw.cons.L = dw.cons.L[i],
sw.duration = swim.duration[i],
sf.vl = sf.viral.load,
dw.vl = dw.viral.load,
sf.cons.g = sf.cons.g,
sw.vl = sw.viral.load,
sw.daily.IR = sw.daily.IR,
sw.frequency = sw.frequency)
}
# Are the results the same?
identical(daily.dose, daily.dose2)## [1] TRUE
For the vector of simulated daily dose results, first report the geometric mean, then plot the kernel density estimates and finally the empirical cumulative distribution.
# Set display options for use with the print() function.
options(digits=3)
# Print the geometric mean of the vector of simulated daily dose results.
print(format(exp(mean(log(daily.dose))), scientific = TRUE))## [1] "1.37e-01"
Calculate and print the kernel density estimates using the density() function
from the stats package.
dens <- density(daily.dose)
dens##
## Call:
## density.default(x = daily.dose)
##
## Data: daily.dose (5000 obs.); Bandwidth 'bw' = 0.0001264
##
## x y
## Min. :0.135 Min. : 0
## 1st Qu.:0.137 1st Qu.: 1
## Median :0.140 Median : 11
## Mean :0.140 Mean :112
## 3rd Qu.:0.142 3rd Qu.:158
## Max. :0.144 Max. :577
Calculate the mean, geometric mean, median, and mode.
# Calculate measures of central tendency.
meas <- data.frame(
measure = c("mean", "g. mean", "median", "mode"),
value = round(c(
mean(daily.dose), exp(mean(log(daily.dose))),
median(daily.dose), dens$x[which.max(dens$y)]
), 6),
color = c("red", "orange", "green", "blue"),
stringsAsFactors = FALSE
)
# Set display options for use with the print() function.
options(digits=6)
# Print measures of central tendency.
print(meas[1:2])## measure value
## 1 mean 0.137151
## 2 g. mean 0.137149
## 3 median 0.136990
## 4 mode 0.136810
Plot the kernel density estimates with measures of central tendency.
# Contruct text labels by combining each measure with its value.
meas$label <- sapply(1:nrow(meas), function(x)
paste(meas$measure[x], as.character(meas$value[x]), sep = ' = '))
# Add lines for measures of central tendency and a legend to a plot.
add_lines_and_legend <- function(meas, x.pos = 0, y.pos = 0, cex = 1) {
n <- nrow(meas)
# Plot measures of central tendency as vertical lines.
res <- sapply(1:n, function(x)
abline(v = meas$value[x], col = meas$color[x]))
# Add a legend to the plot.
legend(x.pos, y.pos, meas$label, col = meas$color,
cex = cex, lty = rep(1, n), lwd = rep(2, n))
}
# Plot the kernel density estimates.
plot(dens)
# Add lines for measures of central tendency and a legend.
add_lines_and_legend(meas, 0.139, 550)
Plot the empirical cumulative distribution with measures of central tendency.
# Plot the empirical cumulative distribution for the exposure estimates.
plot(ecdf(daily.dose))
# Add lines for measures of central tendency and a legend.
add_lines_and_legend(meas, 0.139, 0.8)
Except where noted otherwise, this work is licensed under a Creative Commons Attribution-ShareAlike 4.0 International License.
- Share — copy and redistribute the material in any medium or format
- Adapt — remix, transform, and build upon the material for any purpose, even commercially.
The licensor cannot revoke these freedoms as long as you follow the license terms.
-
Attribution — You must give appropriate credit, provide a link to the license, and indicate if changes were made. You may do so in any reasonable manner, but not in any way that suggests the licensor endorses you or your use.
-
ShareAlike — If you remix, transform, or build upon the material, you must distribute your contributions under the same license as the original.
No additional restrictions — You may not apply legal terms or technological measures that legally restrict others from doing anything the license permits.
- You do not have to comply with the license for elements of the material in the public domain or where your use is permitted by an applicable exception or limitation.
- No warranties are given. The license may not give you all of the permissions necessary for your intended use. For example, other rights such as publicity, privacy, or moral rights may limit how you use the material.