withCladded
Anika Larson
5/10/2016
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 chlorine.
# 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:
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 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")
# Chlorine use, discrete distribution weighting viability reduction after 120 minutes or more given 2ppm, 10ppm, or 15ppm Cl in the water according to estimates from Pereira et al. (2008), Keegan et al. (2003), and LeChavallier and Au (2004) respectively:
cl.use <- mcstoc(rempiricalD, type = "V", seed = seed,
values = c(0.49, 0.99, 0.47),
prob = c(0.3, 0.3, 0.3))
# Estimate the exposure based ONLY on drinking water ingestion using the 0, V and U nodes to create a VU node.
expo.mc1 <- (cw.vl * (dw.ing + hose.perevent) * cl.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
# Build a mc model from all of the mcnode objects.
mc(cw.vl, dw.ing, hose.perevent, cl.use, a, n50, expo.mc1, risk)
})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 20.79187
## 2 dw.ing numeric 1 1 1 1 1.4660e+00 1.46600
## 3 hose.perevent numeric 1 1 1 1 1.0000e-02 0.01000
## 4 cl.use numeric 401 1 1 1 4.7000e-01 0.64990
## 5 a numeric 1 100 1 1 3.2323e-01 4828.28223
## 6 n50 numeric 1 100 1 1 1.2129e+01 149.96985
## 7 expo.mc1 numeric 401 100 1 1 4.7114e-04 19.86528
## 8 risk numeric 401 100 1 1 0.0000e+00 0.11716
## median max Nas type outm
## 1 18.19549 9.3749e+01 0 VU each
## 2 1.46600 1.4660e+00 0 0 each
## 3 0.01000 1.0000e-02 0 0 each
## 4 0.49000 9.9000e-01 0 V each
## 5 276.46798 2.7950e+04 0 U each
## 6 136.84941 4.1148e+02 0 U each
## 7 15.81103 1.3443e+02 0 VU each
## 8 0.08293 7.4138e-01 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
##
## cl.use :
## mean sd Min 2.5% 25% 50% 75% 97.5% Max nsv Na's
## NoUnc 0.65 0.241 0.47 0.47 0.47 0.49 0.99 0.99 0.99 401 0
##
## a :
## NoVar
## median 2.76e+02
## mean 4.83e+03
## 2.5% 3.66e-01
## 97.5% 2.33e+04
##
## n50 :
## NoVar
## median 136.8
## mean 150.0
## 2.5% 46.5
## 97.5% 325.1
##
## expo.mc1 :
## mean sd Min 2.5% 25% 50% 75% 97.5% Max nsv Na's
## median 19.9 16.7 0.05457 0.832 7.64 15.9 26.9 63.5 96.5 401 0
## mean 19.9 16.6 0.08920 0.924 7.81 15.8 27.0 63.4 96.8 401 0
## 2.5% 18.6 14.7 0.00313 0.506 6.73 14.5 24.9 55.5 76.7 401 0
## 97.5% 21.2 18.4 0.34650 1.692 9.18 17.1 29.6 71.6 122.8 401 0
##
## risk :
## mean sd Min 2.5% 25% 50% 75% 97.5% Max nsv Na's
## median 0.0998 0.0773 0 0.00000 0.0389 0.0808 0.1419 0.297 0.395 401 0
## mean 0.1172 0.0820 0 0.00255 0.0542 0.1045 0.1645 0.309 0.411 401 0
## 2.5% 0.0415 0.0328 0 0.00000 0.0160 0.0323 0.0561 0.119 0.175 401 0
## 97.5% 0.2589 0.1617 0 0.01918 0.1566 0.2591 0.3708 0.582 0.679 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] 19.86528quantile(mean.expo, probs = seq(0, 1, 0.025))[c("2.5%", "50%", "97.5%")]## 2.5% 50% 97.5%
## 18.57567 19.93447 21.24414Report 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.1171586quantile(mean.risk, probs = seq(0, 1, 0.025))[c("2.5%", "50%", "97.5%")]## 2.5% 50% 97.5%
## 0.04154890 0.09981471 0.25886257Plot 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 = 6)
Questions: