Exponential Crypto
Anika Larson
5/20/2016
For this risk estimate, we will use mcmodel and evalmcmod from the mc2d package.
# Load packages.
library(mc2d)## Loading required package: mvtnorm##
## Attaching package: 'mc2d'## The following objects are masked from 'package:base':
##
## pmax, pminlibrary(data.table)Define variables
Set the number of simulations for the variability and uncertainty dimensions.
nsv <- 401
nsu <- 100
ndvar(nsv) # Variability## [1] 401ndunc(nsu) # Uncertainty## [1] 100The 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 <- 5Introducing the bootstrapped dose-response data, fread command and sampling from Brian:
k_dup <- fread("~/u_drive/envh543/output_dupont_crypto/klist_dupont.csv", sep = ',')
gdrk <- k_dup[sample(nrow(k_dup), nsu, replace = TRUE), list(k)] #sampling rows in the dose-response bootstrapped table equal to the number of simulations in the uncertainty dimension, aka NSUWe 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")
# Hosing frequency (hosing/day), fictitious discrete distribution (this basically says that there is a 10% chance someone never hoses, a 70% chance they use the hose every other day, and a 10% chance they use the hose every day) - will be IGNORED for the moment:
hosing.frequency <- mcstoc(rempiricalD, type = "V", seed = seed,
values = c(0, 0.5, 1),
prob = c(0.1, 0.7, 0.2))
# Showering per event ingestion rate (L/event): worst case estimate = 10mL or .01L per day
shower.perevent <- mcdata(.01, type = "0")
# Estimate the exposure based ONLY on drinking water ingestion using the 0, V and U nodes to create a VU node.
exporaw.mc1 <- (cw.vl * (dw.ing + hose.perevent + shower.perevent))
#Reduction by porous ceramic filtration
expo.mc1 <- (exporaw.mc1 * .001)
#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
k <- mcdata(gdrk$k, type = "U")
risk <- (1 - exp(-dose*k)) #this equation is my version of a crypto exponential DR equation from what came off QMRA Wiki <http://qmrawiki.canr.msu.edu/index.php?title=Dose_response_assessment> and someone should check it -A
#yearlyrisk
yearrisk <- (1 - risk)^365
# Build a mc model from all of the mcnode objects.
mc(cw.vl, dw.ing, hose.perevent, shower.perevent, k, expo.mc1, risk, yearrisk)
})Evaluate the model
Evaluate the model with 5000 iterations in the variability dimension and 250 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 shower.perevent numeric 1 1 1 1 1.0000e-02 1.0000e-02
## 5 k numeric 1 100 1 1 2.0796e-03 4.6348e-03
## 6 expo.mc1 numeric 401 100 1 1 4.7912e-07 3.0897e-02
## 7 risk numeric 401 100 1 1 0.0000e+00 1.4144e-04
## 8 yearrisk numeric 401 100 1 1 1.1000e-03 9.7628e-01
## median max Nas type outm
## 1 18.1954938 93.748620 0 VU each
## 2 1.4660000 1.466000 0 0 each
## 3 0.0100000 0.010000 0 0 each
## 4 0.0100000 0.010000 0 0 each
## 5 0.0042515 0.011239 0 U each
## 6 0.0270385 0.139310 0 VU each
## 7 0.0000000 0.018491 0 VU each
## 8 1.0000000 1.000000 0 VU eachSummarize 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
##
## shower.perevent :
## NoUnc
## NoVar 0.01
##
## k :
## NoVar
## median 0.00425
## mean 0.00463
## 2.5% 0.00245
## 97.5% 0.00933
##
## expo.mc1 :
## mean sd Min 2.5% 25% 50% 75% 97.5% Max
## median 0.0310 0.0218 1.03e-04 0.001512 0.0136 0.0269 0.0446 0.0808 0.1099
## mean 0.0309 0.0218 1.46e-04 0.001581 0.0136 0.0270 0.0444 0.0814 0.1112
## 2.5% 0.0292 0.0201 6.70e-06 0.000926 0.0116 0.0243 0.0411 0.0750 0.0933
## 97.5% 0.0327 0.0236 5.85e-04 0.002872 0.0154 0.0298 0.0472 0.0886 0.1341
## 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.27e-04 0.000746 0 0 0 0 0 0.00332 0.00505 401 0
## mean 1.41e-04 0.000806 0 0 0 0 0 0.00302 0.00573 401 0
## 2.5% 4.88e-05 0.000370 0 0 0 0 0 0.00000 0.00248 401 0
## 97.5% 2.73e-04 0.001582 0 0 0 0 0 0.00701 0.01417 401 0
##
## yearrisk :
## mean sd Min 2.5% 25% 50% 75% 97.5% Max nsv Na's
## median 0.978 0.131 0.15772 0.2970 1 1 1 1 1 401 0
## mean 0.976 0.133 0.17786 0.4741 1 1 1 1 1 401 0
## 2.5% 0.962 0.084 0.00552 0.0767 1 1 1 1 1 401 0
## 97.5% 0.989 0.178 0.40427 1.0000 1 1 1 1 1 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.03089672quantile(mean.expo, probs = seq(0, 1, 0.025))[c("2.5%", "50%", "97.5%")]## 2.5% 50% 97.5%
## 0.02921265 0.03101782 0.03271462Report 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.0001414401quantile(mean.risk, probs = seq(0, 1, 0.025))[c("2.5%", "50%", "97.5%")]## 2.5% 50% 97.5%
## 4.880896e-05 1.267994e-04 2.725295e-04Plot 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$risk, main = "", ylab = "Proportion of Estimates", xlab = "Daily Probability of Infection")
title("Daily Infec. Risk with Cryptosporidium by RCRW
in the USVI", cex.main = 1, font.main= 1, line = 7)
# Generate an "ecdf" plot. This actually calls plot.mcnode().
plot(expo.ev1$yearrisk, main = "", ylab = "Proportion of Estimates", xlab = "Yearly Probability of Infection")
title("Yearly Infec. Risk with Cryptosporidium by RCRW
in the USVI", cex.main = 1, font.main= 1, line = 7)