ex0618.Rmd

---
title: "QMRA, 2nd ed., Example 6.18"
author: "Brian High and John Kissel"
output:
  html_document:
    keep_md: yes
---

## Introduction

This document offers a solution in R 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 into the 
[public domain](https://creativecommons.org/publicdomain/zero/1.0/).

## Set display options

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

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

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

## View the input data

```{r}
# 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%;'")
```

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

```{r}
# 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%;'")
```

## Estimate average daily ingestion rate

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

```{r}
# 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%;'")
```

## Estimate exposures

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

```{r}
# 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%;'")
```

## Estimate total exposure

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

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

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

```{r}
# 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%;'")
```