This is a 2-dimensional Monte Carlo simulation using the mc2d package.

We will use mcmodel and evalmcmod from the mc2d package. This is a simple way to implement the 2-D simulation because we just need to define the model with mcmodel() and evaluate it with evalmcmod().

Load packages

We did this earlier in the document, but will add it again here as a reminder of what packages are needed for this example.

# Load packages.
library(mc2d)
## Loading required package: mvtnorm
## 
## Attaching package: 'mc2d'
## The following objects are masked from 'package:base':
## 
##     pmax, pmin

Define variables

Set the number of simulations for the variability and uncertainty dimensions.

ndvar(401)  # Variability
## [1] 401
ndunc(59)   # Uncertainty
## [1] 59

The mc2d functions will set nsv to ndvar() and nsu to ndunc() by default. Alternatively, we could supply these values to the mc2d model functions when we call them.

Define the variables we will use to set the seed for random sampling and the number of digits for print() statements.

seed <- 1
digits <- 5

We will use this variable to explicitly set the seed with the various mc2d functions. Another approach would be to only set it through the model evaluation functions, or not at all. Since we want to do our best to provide the most reproducible results, we set the seed explicitly.

Define exposure model

Within the mcmodel function, use mcstoc to define “mc nodes” for each component of the model.

For each stochastic variable, use the mcstoc() function to create a “mcnode”, supply a probability function, the node type as “V” for variability or “U” for uncertainty, and any additional parameters to be passed to the probability function.

For deterministic variables, create nodes with the mcdata() function using the “0” type.

We will model the deterministic factors as uniform probablity distributions.

The last statement makes an “mc” object using the mc function.

Here is where you first create all your variables. Look for what the variables mean, and then for the section after “Estimate the exposure using the 0, V and U nodes to create a VU node,” because that is the equation that connects them.

The program will then run through a bunch of loops for that equation, taking into account uncertainty for whichver variables you called “U” or “UV”.

If we wanted to, I think we could throw in a dose-response equation at the end here and generate different plots for the exposure and for the likelihood of infection. TBD.

# Define an exposure model for evaluation by evalmcmod().
expo.mod1 <- mcmodel({
  # Values from Example 6.18 from Quantitative Microbial Risk Assessment, 
  # 2nd Edition by Charles N. Haas, Joan B. Rose, and Charles P. Gerba. 
  # (Wiley, 2014), pp. 215-216. Other fictitious values are noted below.
  
  # Cistern water viral loading (viruses/L):
  cw.vl <- cw.vl <- mcstoc(rnorm, type = "VU", seed = seed, mean = 8.39, sd = 21.8, linf = 0)
  
  # Water consumption (L/day):
  dw.ing <- mcdata(1.466, type = "U")

  # Hosing per event ingestion rate (L/event): worst case estimate = 10mL or .01L
  hose.perevent <- mcdata(.01, type = "0")

  # Hosing frequency (hosing/day), highly conservative estimate:
  hosing.frequency <- mcdata(1, type = "0")

  # Estimate the exposure using the 0, V and U nodes to create a VU node.
  expo.mc1 <- (cw.vl * dw.ing)

  dose <- mcstoc(rpois, type = "VU", lambda = expo.mc1)

  risk <- 1 - (1 - .00419)^dose

  # Build a mc model from all of the mcnode objects.
  mc(cw.vl, dw.ing, expo.mc1, dose, risk) 
})

Evaluate the model

Evaluate the model with 401 iterations in the variability dimension and 59 iterations in the uncertainty dimension, as set previously.

expo.ev1 <- evalmcmod(expo.mod1, seed = seed)
## Warning in (function (n, lambda) : NAs produced
print(expo.ev1, digits = digits)
##       node    mode nsv nsu nva variate      min     mean   median
## 1    cw.vl numeric 401  59   1       1  -85.411  8.36581  8.21693
## 2   dw.ing numeric   1  59   1       1    1.466  1.46600  1.46600
## 3 expo.mc1 numeric 401  59   1       1 -125.212 12.26428 12.04602
## 4     dose numeric 401  59   1       1    0.000 30.54859 26.00000
## 5     risk numeric 401  59   1       1    0.000  0.11656  0.10342
##         max  Nas type outm
## 1  91.45403    0   VU each
## 2   1.46600    0    U each
## 3 134.07161    0   VU each
## 4 162.00000 8287   VU each
## 5   0.49349 8287   VU each

Summarize results

Print a summary and a plot of the evaluation results (expo.ev1).

# Summarize the results.
summary(expo.ev1)
## Warning: NA's value in mc object

## Warning: NA's value in mc object
## cw.vl :
##         mean   sd   Min  2.5%   25%   50%  75% 97.5%  Max nsv Na's
## median  8.34 21.8 -56.7 -33.3 -6.34  8.22 23.3  50.5 71.4 401    0
## mean    8.37 21.8 -56.7 -33.9 -6.16  8.24 22.9  51.0 72.7 401    0
## 2.5%    6.52 20.4 -71.8 -39.3 -8.83  5.93 19.8  43.1 60.2 401    0
## 97.5%  10.85 23.6 -43.1 -29.5 -3.55 11.00 25.7  57.5 88.9 401    0
## 
## dw.ing :
##        NoVar
## median  1.47
## mean    1.47
## 2.5%    1.47
## 97.5%   1.47
## 
## expo.mc1 :
##         mean   sd    Min  2.5%    25%  50%  75% 97.5%   Max nsv Na's
## median 12.23 31.9  -83.1 -48.8  -9.30 12.0 34.1  74.0 104.6 401    0
## mean   12.26 32.0  -83.1 -49.7  -9.03 12.1 33.6  74.7 106.5 401    0
## 2.5%    9.56 29.9 -105.2 -57.6 -12.95  8.7 29.0  63.2  88.3 401    0
## 97.5%  15.90 34.5  -63.2 -43.2  -5.20 16.1 37.6  84.3 130.3 401    0
## 
## dose :
##        mean   sd Min  2.5%  25%  50%  75% 97.5%   Max nsv Na's
## median 30.5 22.5   0 1.000 13.0 26.0 44.0  83.0 111.0 401  141
## mean   30.5 22.5   0 0.938 12.7 26.5 44.0  82.8 112.2 401  140
## 2.5%   28.3 20.2   0 0.000 10.0 24.0 39.5  71.5  88.3 401  121
## 97.5%  33.3 24.3   0 2.000 15.0 29.0 48.8  95.5 142.2 401  157
## 
## risk :
##         mean     sd Min    2.5%    25%    50%   75% 97.5%   Max nsv Na's
## median 0.116 0.0797   0 0.00419 0.0531 0.1034 0.169 0.294 0.373 401  141
## mean   0.117 0.0801   0 0.00393 0.0520 0.1053 0.169 0.294 0.375 401  140
## 2.5%   0.109 0.0726   0 0.00000 0.0411 0.0959 0.153 0.259 0.310 401  121
## 97.5%  0.126 0.0862   0 0.00836 0.0610 0.1146 0.185 0.330 0.449 401  157
# Plot the results.
plot(expo.ev1)

Report the mean and median of the means with a 95% confidence interval (CI95).

mean.exp <- sapply(1:ndunc(), function(j) mean(expo.ev1$expo.mc1[, j, ]))
mean(mean.exp)
## [1] 12.26428
quantile(mean.exp, probs = seq(0, 1, 0.025))[c("2.5%", "50%", "97.5%")]
##      2.5%       50%     97.5% 
##  9.564655 12.231495 15.898801

Plot the empirical cumulative distribution function (ecdf) of the exposure model (expo.mc1) estimates.

I added that title and those X and Y labels based on what I think it probably means! Tell me if you think it’s wrong.

# Generate an "ecdf" plot. This actually calls plot.mcnode().
plot(expo.ev1$expo.mc1, main = "Exposure Assessment", xlab = "Daily Ingested Crypto Oocysts",ylab = "Proportion of Population")

Now, this should print the risk assessment curve for Campylobacter using the data we’ve given it so far:

Plot the empirical cumulative distribution function (ecdf) of the exposure model (risk) estimates.

Again, I added that title and those X and Y labels based on what I think it probably means! Tell me if you think it’s wrong.

# Generate an "ecdf" plot. This actually calls plot.mcnode().
plot(expo.ev1$risk, main = "Risk Assessment", xlab = "Risk of Crypto Infection",ylab = "Proportion of Population")

To Dos:

  1. add actual values to the model, like not fictitious guesses about how often people are mowing their lawns.

-decide which variables have variation, which are uncertain, and which are both.

  1. do this for pathogens other than Crypto.

  2. figure out a way to integrate data on the same pathogen from multiple sources. Not at all sure how to do this so far.

  3. add dose-response component

  4. also, figure out how to tell the model that NO ONE GETS AN EXPOSURE LESS THAN 0. Actually, WHY it is generating exposures less than 0 in the first place. THAT SHOULD NOT BE POSSIBLE.

-Update: it has something to do with the variability and uncertainty in the cistern water pathogen estimates, because if you just set that to be a value of 8.39, you get a nice narrow range of dose estimates around 12, with just a little variation based on how often you use a hose. So how to fix that????????