ex0618prob.Rmd

---
title: "1-D Monte Carlo simulation of microbial exposure"
author: "Jane Pouzou and Brian High"
date: '![CC BY-SA 4.0](cc_by-sa_4.png)'
output:
  html_document:
    toc: true
    theme: united
    keep_md: yes
    self-contained: yes
---

## Introduction

This document offers a 1-D 
[Monte Carlo](https://en.wikipedia.org/wiki/Monte_Carlo_method) probabilistic 
solution in R for the daily microbial exposure from drinking water consumption, 
swimming in surface water and shellfish consumption for 
[Example 6.18](images/ex0618.png) from pages 
[215-216](https://onlinelibrary.wiley.com/doi/pdf/10.1002/9781118910030.ch6#page=57) of:

[Quantitative Microbial Risk Assessment, 2nd Edition](http://www.wiley.com/WileyCDA/WileyTitle/productCd-1118145291,subjectCd-CH20.html) 
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]( https://creativecommons.org/licenses/by-sa/4.0/). 
Details may be found at the end of this document.

## Set global options

```{r}
# Set knitr options for use when rendering this document.
library("knitr")
opts_chunk$set(cache=TRUE, message=FALSE)
```

## Define variables

Define variables provided in the example for three exposure types.

```{r}
# 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 year
```

## Sample from probablity distributions

Sample from the probablity distributions for drinking water and swim duration. 
Also plot from these distributions as a quick visual check.

### Drinking water

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.

```{r}
set.seed(1)
dw.cons.L <- rlnorm(5000, meanlog = 7.49, sdlog = 0.407) / 1000
```

Plot the kernal density curve of the generated values just as a check.

```{r dw-dens-plot}
plot(density(dw.cons.L))
```

### Swim duration

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.

```{r}
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.

```{r swim-hist-plot}
hist(swim.duration)
```

## Estimate daily dose

Calculate estimated daily dose using a probabilistic simulation model.

### Define risk function

Define a function to calculate microbial exposure risk.

```{r}
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 the simulation

Compute 5000 simulated daily dose results and store as a vector.

```{r}
# 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)
```

## Summarize results

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.

```{r}
# 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))
```

### Calculate kernel density estimates

Calculate and print the kernel density estimates using the `density()` function 
from the *stats* package.

```{r}
dens <- density(daily.dose)
dens
```

### Calculate measures of central tendency

Calculate the mean, geometric mean, median, and mode.

```{r}
# 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])
```

### Plot kernel density estimates

Plot the kernel density estimates with measures of central tendency.

```{r kern-dens-plot}
# 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 empirical cumulative distribution

Plot the empirical cumulative distribution with measures of central tendency.

```{r ecdf-plot}
# 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)
```

## License

Except where noted otherwise, this work is licensed under a 
[Creative Commons Attribution-ShareAlike 4.0 International License]( https://creativecommons.org/licenses/by-sa/4.0/). 

### You are free to:

* 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.

### Under the following 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.

### Notices:

* 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.