For this risk estimate, we will use mcmodel and evalmcmod from the mc2d package.

The following describes a distribution of estimates of of daily probability of Cryptosporidium infection for all households relying solely on rooftop-collected rainwater for drinking and hosing, where the household water was treated with UV filtration.

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

Define variables

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

nsv <- 401
nsu <- 100
ndvar(nsv)  # Variability
## [1] 401
ndunc(nsu)   # Uncertainty
## [1] 100

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

Introducing the bootstrapped dose-response data, fread command and sampling from Brian:

alpha <- fread("~/u_drive/envh543/output/alphalist.csv", sep = ',')
n50 <- fread("~/u_drive/envh543/output/N50list.csv", sep = ',')
cdr <- merge(alpha,n50, by = "V1")
scdr <- cdr[sample(nrow(cdr), nsu, replace = TRUE), list(alpha, N50)] #sampling rows in the dose-response bootstrapped table equal to the number of simulations in the uncertainty dimension, aka NSU

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.

# 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, rtrunc = TRUE, linf = 0)
  
  # Water consumption (L/day):
  dw.ing <- mcdata(1.466, type = "0")

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

  # UV filtration at UV dose < 10 mQ-s/cm2:
  UV.use <- mcdata(.001, type = "0")
  
  # Estimate the exposure based on drinking water ingestion using the 0, V and U nodes to create a VU node.
  expo.mc1 <- (cw.vl * (dw.ing + hose.perevent) * UV.use) 
  
  #Estimate the risk using a dose-response model pulling from a bootstrapped model of the dose-response
  dose <- mcstoc(rpois, type = "VU", lambda = expo.mc1)
  
  # Create mcnodes from previously sampled dose-response bootstrap estimates. - Brian
  
  a <- mcdata(scdr$alpha, type = "U")
  n50 <- mcdata(scdr$N50, type = "U")
  
  risk <- 1 - (1 + dose*((2^(1/a)-1)/n50))^(-a) #this equation is my version of what came off QMRA Wiki <http://qmrawiki.canr.msu.edu/index.php?title=Dose_response_assessment> and someone should check it -A
  
  lrisk <- ifelse(risk==0, c(log10(.0000000000000000000000001)), c(log10(risk)))
  
  # Build a mc model from all of the mcnode objects.
  mc(cw.vl, dw.ing, hose.perevent, UV.use, a, n50, expo.mc1, risk, lrisk) 
})

Evaluate the model

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

expo.ev1 <- evalmcmod(expo.mod1, seed = seed)
print(expo.ev1, digits = digits)
##            node    mode nsv nsu nva variate         min        mean
## 1         cw.vl numeric 401 100   1       1  3.2242e-04  2.0792e+01
## 2        dw.ing numeric   1   1   1       1  1.4660e+00  1.4660e+00
## 3 hose.perevent numeric   1   1   1       1  1.0000e-02  1.0000e-02
## 4        UV.use numeric   1   1   1       1  1.0000e-03  1.0000e-03
## 5             a numeric   1 100   1       1  2.5232e-01  4.2369e+03
## 6           n50 numeric   1 100   1       1  3.0449e+01  1.4968e+02
## 7      expo.mc1 numeric 401 100   1       1  4.7590e-07  3.0689e-02
## 8          risk numeric 401 100   1       1  0.0000e+00  2.4537e-04
## 9         lrisk numeric 401 100   1       1 -2.5000e+01 -2.4294e+01
##       median         max Nas type outm
## 1  18.195494  9.3749e+01   0   VU each
## 2   1.466000  1.4660e+00   0    0 each
## 3   0.010000  1.0000e-02   0    0 each
## 4   0.001000  1.0000e-03   0    0 each
## 5   5.750011  3.1706e+04   0    U each
## 6 137.570101  3.9119e+02   0    U each
## 7   0.026857  1.3837e-01   0   VU each
## 8   0.000000  6.3853e-02   0   VU each
## 9 -25.000000 -1.1948e+00   0   VU each

Summarize results

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

# Summarize the results.
summary(expo.ev1)
## cw.vl :
##        mean   sd     Min  2.5%   25%  50%  75% 97.5%  Max nsv Na's
## median 20.9 14.6 0.06900 1.018  9.13 18.1 30.0  54.4 74.0 401    0
## mean   20.8 14.7 0.09840 1.064  9.15 18.2 29.8  54.8 74.9 401    0
## 2.5%   19.7 13.5 0.00451 0.623  7.83 16.3 27.6  50.5 62.8 401    0
## 97.5%  22.0 15.9 0.39339 1.933 10.39 20.1 31.8  59.6 90.2 401    0
## 
## dw.ing :
##       NoUnc
## NoVar  1.47
## 
## hose.perevent :
##       NoUnc
## NoVar  0.01
## 
## UV.use :
##       NoUnc
## NoVar 0.001
## 
## a :
##           NoVar
## median 5.75e+00
## mean   4.24e+03
## 2.5%   4.23e-01
## 97.5%  2.20e+04
## 
## n50 :
##        NoVar
## median 137.6
## mean   149.7
## 2.5%    55.7
## 97.5%  322.8
## 
## expo.mc1 :
##          mean     sd      Min     2.5%    25%    50%    75%  97.5%    Max
## median 0.0308 0.0216 1.02e-04 0.001502 0.0135 0.0268 0.0443 0.0803 0.1092
## mean   0.0307 0.0217 1.45e-04 0.001570 0.0135 0.0268 0.0441 0.0808 0.1105
## 2.5%   0.0290 0.0200 6.65e-06 0.000919 0.0116 0.0241 0.0408 0.0745 0.0927
## 97.5%  0.0325 0.0234 5.81e-04 0.002852 0.0153 0.0296 0.0469 0.0880 0.1332
##        nsv Na's
## median 401    0
## mean   401    0
## 2.5%   401    0
## 97.5%  401    0
## 
## risk :
##            mean       sd Min 2.5% 25% 50% 75%   97.5%     Max nsv Na's
## median 1.81e-04 0.000955   0    0   0   0   0 0.00365 0.00614 401    0
## mean   2.45e-04 0.001385   0    0   0   0   0 0.00525 0.00966 401    0
## 2.5%   6.15e-05 0.000422   0    0   0   0   0 0.00000 0.00259 401    0
## 97.5%  7.94e-04 0.003831   0    0   0   0   0 0.01927 0.03329 401    0
## 
## lrisk :
##         mean   sd Min 2.5% 25% 50% 75%  97.5%   Max nsv Na's
## median -24.3 3.95 -25  -25 -25 -25 -25  -2.44 -2.21 401    0
## mean   -24.3 3.90 -25  -25 -25 -25 -25 -10.42 -2.14 401    0
## 2.5%   -24.6 2.91 -25  -25 -25 -25 -25 -25.00 -2.59 401    0
## 97.5%  -23.8 5.00 -25  -25 -25 -25 -25  -1.72 -1.48 401    0
# Plot the results.
plot(expo.ev1)

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

mean.expo <- sapply(1:ndunc(), function(j) mean(expo.ev1$expo.mc1[, j, ]))
mean(mean.expo)
## [1] 0.0306888
quantile(mean.expo, probs = seq(0, 1, 0.025))[c("2.5%", "50%", "97.5%")]
##       2.5%        50%      97.5% 
## 0.02901606 0.03080909 0.03249447

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

mean.risk <- sapply(1:ndunc(), function(j) mean(expo.ev1$risk[, j, ]))
mean(mean.risk)
## [1] 0.0002453741
quantile(mean.risk, probs = seq(0, 1, 0.025))[c("2.5%", "50%", "97.5%")]
##         2.5%          50%        97.5% 
## 0.0000614522 0.0001813395 0.0007937976

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

# Generate an "ecdf" plot. This actually calls plot.mcnode().
plot(expo.ev1$expo.mc1)

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

# Generate an "ecdf" plot. This actually calls plot.mcnode().
plot(expo.ev1$risk)

So if we did this right, could the plot above can be correctly labeled as follows?

# Generate an "ecdf" plot. This actually calls plot.mcnode().
plot(expo.ev1$lrisk, main = "", ylab = "Proportion of Estimates", xlab = "Daily Probability of Infection", log = x)
title("Daily Infec. Risk with Cryptosporidium by RCRW
in the USVI for individual, UV filtration", cex.main = 1, font.main= 1, line = 6)

Questions: