CryptoAllAtOnce
Anika Larson
5/30/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 <- 100
nsu <- 100
ndvar(nsv) # Variability## [1] 100ndunc(nsu) # Uncertainty## [1] 100rdv <- (nsu/3)
cwv <- (nsu/4)
div.nsv <- (nsv/3)
div5.nsv <- (nsv/5)
div5.nsu <- (nsu/5)
div7.nsu <- (nsu/7)
ns <- nsu * nsvThe 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("output_dupont_crypto/alphalist.csv", sep = ',')
n50 <- fread("output_dupont_crypto/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 NSULeaving that ^ for convenience, but since what we really want is a k-value list made of samples from the three different dose-response estimates (Dupont, Tamu, and UCP):
k_dup <- fread("output_dupont_crypto/klist_dupont.csv", sep = ',')
gdrkD <- k_dup[sample(nrow(k_dup), div.nsv, replace = TRUE), list(k)] #sampling rows in the dose-response bootstrapped table equal to 1/3 the number of simulations in the uncertainty dimension, aka NSU
k_tamu <- fread("output_tamu_crypto/klist.csv", sep = ',')
gdrkT <- k_tamu[sample(nrow(k_tamu), div.nsv, replace = TRUE), list(k)]
k_ucp <- fread("output_UCP_crypto/klist.csv", sep = ',')
gdrkU <- k_ucp[sample(nrow(k_ucp), div.nsv, replace = TRUE), list(k)]
kvalues <- rbind(k_dup,k_tamu,k_ucp)
num.kvals <- as.numeric(kvalues$k)
what <- kvalues[sample(nrow(kvalues), nsv, replace = TRUE), list(k)]
kval <- as.numeric(what$k)Sampling from the Cryptosporidium custom distribution estimates:
#Using the original custom distribution of estimates from aggregated cistern estimates:
cw.vlg1 <- fread("conc_data_aggreg.csv", sep = ',')
cw.vl.samp1 <- sample(cw.vlg1[,c(oocysts)], nsu, replace = TRUE)
#Using a combination of log-normal distributions from each of the the 5 different cisterns:
cw.vl.ln <- fread("concdata_percis.csv", sep = ',')
v1 <- cw.vl.ln[sample(nrow(cw.vl.ln), div7.nsu, replace = TRUE), list(c_v1_tmed)]
colnames(v1) <- c("cw.vals")
v11 <- cw.vl.ln[sample(nrow(cw.vl.ln), div7.nsu, replace = TRUE), list(c_v11_tmed)]
colnames(v11) <- c("cw.vals")
v2 <- cw.vl.ln[sample(nrow(cw.vl.ln), div7.nsu, replace = TRUE), list(c_v2_tmed)]
colnames(v2) <- c("cw.vals")
v5 <- cw.vl.ln[sample(nrow(cw.vl.ln), div7.nsu, replace = TRUE), list(c_v5_tmed)]
colnames(v5) <- c("cw.vals")
v8 <- cw.vl.ln[sample(nrow(cw.vl.ln), div7.nsu, replace = TRUE), list(c_v8a_tmed)]
colnames(v8) <- c("cw.vals")
v9 <- cw.vl.ln[sample(nrow(cw.vl.ln), div7.nsu, replace = TRUE), list(c_v9a_tmed)]
colnames(v9) <- c("cw.vals")
v10 <- cw.vl.ln[sample(nrow(cw.vl.ln), div7.nsu, replace = TRUE), list(c_v10_tmed)]
colnames(v10) <- c("cw.vals")
cw.vlg2 <- unlist(c(v1$cw.vals,v11$cw.vals,v2$cw.vals,v5$cw.vals,v8$cw.vals,v9$cw.vals,v10$cw.vals))
cw.vl.samp2 <- sample(cw.vlg2, nsu, replace = TRUE)
#Combine the estimates prior to sampling
cw.vl.ln <- fread("concdata_percis.csv", sep = ',')
cw.vlg3 <- unlist(c(cw.vl.ln$c_v1_tmed, cw.vl.ln$c_v2_tmed, cw.vl.ln$c_v5_tmed, cw.vl.ln$c_v8a_tmed, cw.vl.ln$c_v9a_tmed, cw.vl.ln$c_v10_tmed))
#Determine which estimate to use: the original custom distribution (cw.vl.samp1) or the combined lognormal ditributions for each cistern (cw.vl.samp2).
cw.vl.samp <- cw.vlg3```
Sampling from the distribution of the estimates of Crypto reduction by porous ceramic filtration, and converting them from log reduction to percent reduction:
#pfc_reduc set using the distribution
gcfilt <- fread("cer_filts/CerFiltCrypto.csv", sep = ',')
colnames(gcfilt) <- c("trial","logr")
gcfilt$reduc <- 10^(-gcfilt$logr)
pfc_reduc2 <- gcfilt[sample(nrow(gcfilt), nsu, replace = TRUE), list(reduc)]Sampling from the daily hosing exposure estimates distribution:
#Creating a distribution for hosing
hose <- fread("exp_params/Hosing.csv", sep = ',')
colnames(hose) <- c("trial","num", "mL")
hose$hose_exp1 <- (hose$num * hose$mL)
hose_exp <- hose[sample(nrow(hose), nsu, replace = TRUE), list(hose_exp1)]Sampling from the daily showering exposure estimates distribution:
#Creating a distribution for showering
shower <- fread("exp_params/Showers.csv", sep = ',')
colnames(shower) <- c("trial","num", "mL")
shower$shower_exp1 <- (shower$num * shower$mL)
shower_exp <- shower[sample(nrow(shower), nsu, replace = TRUE), list(shower_exp1)]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.vl2 <- mcstoc(rempiricalD, type = "VU", values = cw.vl.samp)
cw.vl <- cw.vl2
# Water consumption (L/day):
dw.ing1 <- mcdata(1.466, type = "U")
dw.ing2 <- mcstoc(rpois, type = "V", lambda = 1.466)
dw.ing <- dw.ing1
# Hosing per event ingestion rate (L/event): worst case estimate = 10mL or .01L per day
hose.perevent <- mcdata(.01, type = "U")
hose.perevent.dis <- mcdata(hose_exp$hose_exp1, type = "U")
#Daily hosing exposure sampling from distribution
daily_hose_exp <- mcstoc(rempiricalD, type = "V", values = hose_exp$hose_exp1)
hosing <- daily_hose_exp
# Showering per event ingestion rate (L/event): worst case estimate = 10mL or .01L per day
shower.perevent <- mcdata(.01, type = "U")
shower.perevent.dis <- mcdata(shower_exp$shower_exp1, type = "U")
#Daily showering sampling from distribution:
daily_shower_exp <- mcstoc(rempiricalD, type = "V", values = shower_exp$shower_exp1)
shower <- daily_shower_exp
# Estimate the exposure based on drinking water, hosing and ingestion, assuming no particular distribution:
exporaw.mc1 <- (cw.vl * (dw.ing + hosing + shower))
#Reduction by one of the three estimates in Table 3, assuming all are equally likely:
exporaw.freechl.tog <- mcstoc(rempiricalD, values = c(.53,.5096,.01), prob = c(.33,33,33), type = "V")
exporaw_free_chl <- (exporaw.mc1 * exporaw.freechl.tog)
#Reduction by porous ceramic filtration based on the distribution entered earlier:
pfc_reduc <- mcstoc(rempiricalD, type = "V", values = pfc_reduc2$reduc)
exporaw.porous_ceram <- (exporaw.mc1 * pfc_reduc)
#Reduction by porous ceramic filtration, assuming 99% reduction:
exporaw.por_ceram_singval <- (exporaw.mc1 * .01)
#Reduction by UV filtration:
exporaw.uv <- (exporaw.mc1 * .0001)
#Final exposure - can be assigned to whichever filtration system is in use:
expo.mc1 <- exporaw.porous_ceram
dose1 <- mcstoc(rpois, type = "VU", lambda = expo.mc1)
dose2 <- expo.mc1
dose <- dose2
# Create mcnodes from previously sampled dose-response bootstrap estimates - a and n50 just left here for convenience.
a <- mcdata(scdr$alpha, type = "U")
n50 <- mcdata(scdr$N50, type = "U")
k1 <- mcdata(kval, type = "V")
k2 <- mcstoc(rempiricalD, type = "V", values = num.kvals)
k <- k2
#estimating risk with exponential equation:
risk1 <- (1 - exp(-dose*k))
risk2 <- mcstoc(rpois, type = "VU", lambda = risk1)
risk <- risk1
#yearlyrisk
yearrisk <- 1- ((1 - risk)^365)
#RAW VERSION:
#Either directly consider the exposure the dose, or add uncertainty by considering it a Poisson distribution centered at the previous estimate (? - this showed up in an earlier version) :
dose1.raw <- mcstoc(rpois, type = "VU", lambda = exporaw.mc1)
dose2.raw <- exporaw.mc1
dose.raw <- dose2.raw
#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. BUT it's beta-poisson, and just left here for convenience.
riskbp <- 1 - (1 + dose*((2^(1/a)-1)/n50))^(-a)
#Exponential version of the equation (also from QMRA Wiki and please lmk if it looks wrong):
risk1.raw <- (1 - exp(-dose.raw*k))
risk2 <- mcstoc(rpois, type = "VU", lambda = risk1)
risk.raw <- risk1.raw
#yearlyrisk
yearrisk.raw <- 1- ((1 - risk.raw)^365)
#FREE CHLORINE VERSION:
#Either directly consider the exposure the dose, or add uncertainty by considering it a Poisson distribution centered at the previous estimate (? - this showed up in an earlier version) :
dose1 <- mcstoc(rpois, type = "VU", lambda = exporaw_free_chl)
dose2 <- exporaw_free_chl
dose.chl <- dose2
#Exponential version of the equation (also from QMRA Wiki and please lmk if it looks wrong):
risk1.chl <- (1 - exp(-dose.chl*k))
risk2 <- mcstoc(rpois, type = "VU", lambda = risk1)
risk.chl <- risk1.chl
#yearlyrisk
yearrisk.chl <- 1- ((1 - risk.chl)^365)
#UV VERSION:
#Either directly consider the exposure the dose, or add uncertainty by considering it a Poisson distribution centered at the previous estimate (? - this showed up in an earlier version) :
dose1.uv <- mcstoc(rpois, type = "VU", lambda = exporaw.uv)
dose2.uv <- exporaw.uv
dose.uv <- dose2.uv
#Exponential version of the equation (also from QMRA Wiki and please lmk if it looks wrong):
risk1.uv <- (1 - exp(-dose.uv*k))
risk2.uv <- mcstoc(rpois, type = "VU", lambda = risk1)
risk.uv <- risk1.uv
#yearlyrisk
yearrisk.uv <- 1- ((1 - risk.uv)^365)
#POROUS CERAMIC FILTRATION VERSION:
#Either directly consider the exposure the dose, or add uncertainty by considering it a Poisson distribution centered at the previous estimate (? - this showed up in an earlier version) :
dose1.pcf <- mcstoc(rpois, type = "VU", lambda = exporaw.porous_ceram)
dose2.pcf <- exporaw.porous_ceram
dose.pcf <- dose2.pcf
#Exponential version of the equation (also from QMRA Wiki and please lmk if it looks wrong):
risk1.pcf <- (1 - exp(-dose.pcf*k))
risk2.pcf <- mcstoc(rpois, type = "VU", lambda = risk1.pcf)
risk.pcf <- risk1.pcf
#yearlyrisk
yearrisk.pcf <- 1- ((1 - risk.pcf)^365)
# Build a mc model from all of the mcnode objects.
mc(cw.vl2, cw.vl, dw.ing, hose.perevent, shower.perevent, exporaw_free_chl, exporaw.uv, exporaw.porous_ceram, pfc_reduc, a, n50, k1, k2, k, expo.mc1, dose, riskbp, risk1, risk, yearrisk, risk.raw, dose.raw, yearrisk.raw,risk.chl, dose.chl, yearrisk.chl,risk.uv, dose.uv, yearrisk.uv,risk.pcf, dose.pcf, yearrisk.pcf)
})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)
print(expo.ev1, digits = digits)## node mode nsv nsu nva variate min mean
## 1 cw.vl2 numeric 100 100 1 1 1.0000e-08 5.2345e-02
## 2 cw.vl numeric 100 100 1 1 1.0000e-08 5.2345e-02
## 3 dw.ing numeric 1 100 1 1 1.4660e+00 1.4660e+00
## 4 hose.perevent numeric 1 100 1 1 1.0000e-02 1.0000e-02
## 5 shower.perevent numeric 1 100 1 1 1.0000e-02 1.0000e-02
## 6 exporaw_free_chl numeric 100 100 1 1 1.4660e-10 2.1726e-01
## 7 exporaw.uv numeric 100 100 1 1 1.4660e-12 9.8770e-05
## 8 exporaw.porous_ceram numeric 100 100 1 1 4.7541e-11 1.4439e-02
## 9 pfc_reduc numeric 100 1 1 1 2.8184e-03 1.4929e-02
## 10 a numeric 1 100 1 1 1.7844e-01 3.6898e+03
## 11 n50 numeric 1 100 1 1 8.8966e+00 1.4795e+02
## 12 k1 numeric 100 1 1 1 1.0759e-04 2.6559e-03
## 13 k2 numeric 100 1 1 1 8.3790e-05 2.4924e-03
## 14 k numeric 100 1 1 1 8.3790e-05 2.4924e-03
## 15 expo.mc1 numeric 100 100 1 1 4.7541e-11 1.4439e-02
## 16 dose numeric 100 100 1 1 4.7541e-11 1.4439e-02
## 17 riskbp numeric 100 100 1 1 0.0000e+00 2.2345e-04
## 18 risk1 numeric 100 100 1 1 6.3283e-15 3.0010e-05
## 19 risk numeric 100 100 1 1 6.3283e-15 3.0010e-05
## 20 yearrisk numeric 100 100 1 1 2.3098e-12 9.8796e-03
## 21 risk.raw numeric 100 100 1 1 1.7099e-12 2.1536e-03
## 22 dose.raw numeric 100 100 1 1 1.4660e-08 9.8770e-01
## 23 yearrisk.raw numeric 100 100 1 1 6.2410e-10 2.1693e-01
## 24 risk.chl numeric 100 100 1 1 1.7097e-14 4.3336e-04
## 25 dose.chl numeric 100 100 1 1 1.4660e-10 2.1726e-01
## 26 yearrisk.chl numeric 100 100 1 1 6.2406e-12 6.9903e-02
## 27 risk.uv numeric 100 100 1 1 2.2204e-16 2.1926e-07
## 28 dose.uv numeric 100 100 1 1 1.4660e-12 9.8770e-05
## 29 yearrisk.uv numeric 100 100 1 1 8.1046e-14 7.9977e-05
## 30 risk.pcf numeric 100 100 1 1 6.3283e-15 3.0010e-05
## 31 dose.pcf numeric 100 100 1 1 4.7541e-11 1.4439e-02
## 32 yearrisk.pcf numeric 100 100 1 1 2.3098e-12 9.8796e-03
## median max Nas type outm
## 1 2.0000e-02 6.8000e-01 0 VU each
## 2 2.0000e-02 6.8000e-01 0 VU each
## 3 1.4660e+00 1.4660e+00 0 U each
## 4 1.0000e-02 1.0000e-02 0 U each
## 5 1.0000e-02 1.0000e-02 0 U each
## 6 4.0891e-03 2.2106e+01 0 VU each
## 7 1.9743e-05 6.0401e-03 0 VU each
## 8 2.0132e-03 1.3575e+00 0 VU each
## 9 1.2303e-02 5.8884e-02 0 V each
## 10 5.5055e+00 2.2114e+04 0 U each
## 11 1.4059e+02 3.1242e+02 0 U each
## 12 2.5617e-03 9.1829e-03 0 V each
## 13 2.2169e-03 1.8999e-02 0 V each
## 14 2.2169e-03 1.8999e-02 0 V each
## 15 2.0132e-03 1.3575e+00 0 VU each
## 16 2.0132e-03 1.3575e+00 0 VU each
## 17 1.1307e-05 8.8880e-02 0 VU each
## 18 1.2657e-06 4.0127e-03 0 VU each
## 19 1.2657e-06 4.0127e-03 0 VU each
## 20 4.6187e-04 7.6952e-01 0 VU each
## 21 9.5105e-05 2.0141e-01 0 VU each
## 22 1.9743e-01 6.0401e+01 0 VU each
## 23 3.4119e-02 1.0000e+00 0 VU each
## 24 4.9748e-06 5.2686e-02 0 VU each
## 25 4.0891e-03 2.2106e+01 0 VU each
## 26 1.8142e-03 1.0000e+00 0 VU each
## 27 9.5110e-09 2.2491e-05 0 VU each
## 28 1.9743e-05 6.0401e-03 0 VU each
## 29 3.4715e-06 8.1755e-03 0 VU each
## 30 1.2657e-06 4.0127e-03 0 VU each
## 31 2.0132e-03 1.3575e+00 0 VU each
## 32 4.6187e-04 7.6952e-01 0 VU eachSummarize results
Print a summary and a plot of the evaluation results (expo.ev1).
# Summarize the results.
summary(expo.ev1)## cw.vl2 :
## mean sd Min 2.5% 25% 50% 75% 97.5% Max nsv
## median 0.0523 0.1051 1e-08 1e-08 1.00e-08 0.0198 0.0377 0.403 0.560 100
## mean 0.0523 0.1049 1e-08 1e-08 2.73e-04 0.0198 0.0387 0.397 0.552 100
## 2.5% 0.0325 0.0647 1e-08 1e-08 1.00e-08 0.0147 0.0307 0.210 0.348 100
## 97.5% 0.0740 0.1408 1e-08 1e-08 4.48e-03 0.0246 0.0513 0.538 0.677 100
## Na's
## median 0
## mean 0
## 2.5% 0
## 97.5% 0
##
## cw.vl :
## mean sd Min 2.5% 25% 50% 75% 97.5% Max nsv
## median 0.0523 0.1051 1e-08 1e-08 1.00e-08 0.0198 0.0377 0.403 0.560 100
## mean 0.0523 0.1049 1e-08 1e-08 2.73e-04 0.0198 0.0387 0.397 0.552 100
## 2.5% 0.0325 0.0647 1e-08 1e-08 1.00e-08 0.0147 0.0307 0.210 0.348 100
## 97.5% 0.0740 0.1408 1e-08 1e-08 4.48e-03 0.0246 0.0513 0.538 0.677 100
## Na's
## median 0
## mean 0
## 2.5% 0
## 97.5% 0
##
## dw.ing :
## NoVar
## median 1.47
## mean 1.47
## 2.5% 1.47
## 97.5% 1.47
##
## hose.perevent :
## NoVar
## median 0.01
## mean 0.01
## 2.5% 0.01
## 97.5% 0.01
##
## shower.perevent :
## NoVar
## median 0.01
## mean 0.01
## 2.5% 0.01
## 97.5% 0.01
##
## exporaw_free_chl :
## mean sd Min 2.5% 25% 50% 75% 97.5% Max
## median 0.2022 0.688 1.53e-10 4.41e-10 7.59e-08 0.00429 0.0912 1.736 5.10
## mean 0.2173 0.771 1.69e-10 4.57e-10 1.58e-05 0.00424 0.0957 1.982 5.96
## 2.5% 0.0972 0.269 1.47e-10 1.83e-10 1.77e-08 0.00190 0.0433 0.615 1.56
## 97.5% 0.3821 1.871 2.89e-10 8.77e-10 2.51e-04 0.00715 0.1682 4.882 16.12
## nsv Na's
## median 100 0
## mean 100 0
## 2.5% 100 0
## 97.5% 100 0
##
## exporaw.uv :
## mean sd Min 2.5% 25% 50% 75%
## median 9.50e-05 0.000249 1.47e-12 2.69e-12 2.18e-11 2.00e-05 6.38e-05
## mean 9.88e-05 0.000269 1.58e-12 3.20e-12 1.18e-07 2.04e-05 6.74e-05
## 2.5% 5.58e-05 0.000113 1.47e-12 1.60e-12 1.69e-11 1.04e-05 4.86e-05
## 97.5% 1.57e-04 0.000518 2.25e-12 6.44e-12 1.90e-06 3.22e-05 9.66e-05
## 97.5% Max nsv Na's
## median 0.000733 0.001618 100 0
## mean 0.000765 0.001975 100 0
## 2.5% 0.000363 0.000622 100 0
## 97.5% 0.001309 0.004748 100 0
##
## exporaw.porous_ceram :
## mean sd Min 2.5% 25% 50% 75% 97.5%
## median 0.01438 0.0422 8.83e-11 3.27e-10 2.84e-09 0.00205 0.00800 0.0980
## mean 0.01444 0.0468 1.04e-10 3.48e-10 1.06e-05 0.00206 0.00845 0.1042
## 2.5% 0.00766 0.0155 4.75e-11 1.69e-10 1.66e-09 0.00111 0.00571 0.0465
## 97.5% 0.02449 0.1136 2.69e-10 6.53e-10 1.95e-04 0.00331 0.01344 0.1891
## Max nsv Na's
## median 0.3036 100 0
## mean 0.3809 100 0
## 2.5% 0.0942 100 0
## 97.5% 1.0866 100 0
##
## pfc_reduc :
## mean sd Min 2.5% 25% 50% 75% 97.5% Max nsv
## NoUnc 0.0149 0.0123 0.00282 0.00295 0.00794 0.0123 0.0165 0.055 0.0589 100
## Na's
## NoUnc 0
##
## a :
## NoVar
## median 5.51e+00
## mean 3.69e+03
## 2.5% 3.47e-01
## 97.5% 2.20e+04
##
## n50 :
## NoVar
## median 140.6
## mean 147.9
## 2.5% 31.5
## 97.5% 283.0
##
## k1 :
## mean sd Min 2.5% 25% 50% 75% 97.5%
## NoUnc 0.00266 0.00231 0.000108 0.000135 0.00021 0.00256 0.00419 0.00783
## Max nsv Na's
## NoUnc 0.00918 100 0
##
## k2 :
## mean sd Min 2.5% 25% 50% 75% 97.5%
## NoUnc 0.00249 0.00276 8.38e-05 0.000108 0.00017 0.00222 0.00378 0.00819
## Max nsv Na's
## NoUnc 0.019 100 0
##
## k :
## mean sd Min 2.5% 25% 50% 75% 97.5%
## NoUnc 0.00249 0.00276 8.38e-05 0.000108 0.00017 0.00222 0.00378 0.00819
## Max nsv Na's
## NoUnc 0.019 100 0
##
## expo.mc1 :
## mean sd Min 2.5% 25% 50% 75% 97.5%
## median 0.01438 0.0422 8.83e-11 3.27e-10 2.84e-09 0.00205 0.00800 0.0980
## mean 0.01444 0.0468 1.04e-10 3.48e-10 1.06e-05 0.00206 0.00845 0.1042
## 2.5% 0.00766 0.0155 4.75e-11 1.69e-10 1.66e-09 0.00111 0.00571 0.0465
## 97.5% 0.02449 0.1136 2.69e-10 6.53e-10 1.95e-04 0.00331 0.01344 0.1891
## Max nsv Na's
## median 0.3036 100 0
## mean 0.3809 100 0
## 2.5% 0.0942 100 0
## 97.5% 1.0866 100 0
##
## dose :
## mean sd Min 2.5% 25% 50% 75% 97.5%
## median 0.01438 0.0422 8.83e-11 3.27e-10 2.84e-09 0.00205 0.00800 0.0980
## mean 0.01444 0.0468 1.04e-10 3.48e-10 1.06e-05 0.00206 0.00845 0.1042
## 2.5% 0.00766 0.0155 4.75e-11 1.69e-10 1.66e-09 0.00111 0.00571 0.0465
## 97.5% 0.02449 0.1136 2.69e-10 6.53e-10 1.95e-04 0.00331 0.01344 0.1891
## Max nsv Na's
## median 0.3036 100 0
## mean 0.3809 100 0
## 2.5% 0.0942 100 0
## 97.5% 1.0866 100 0
##
## riskbp :
## mean sd Min 2.5% 25% 50% 75%
## median 0.000074 0.000217 3.68e-13 1.86e-12 1.70e-11 1.15e-05 4.69e-05
## mean 0.000223 0.000586 2.95e-12 8.96e-12 5.77e-08 4.21e-05 1.63e-04
## 2.5% 0.000034 0.000084 0.00e+00 0.00e+00 7.15e-12 4.20e-06 2.06e-05
## 97.5% 0.001462 0.003654 2.79e-11 4.30e-11 8.33e-07 2.78e-04 8.97e-04
## 97.5% Max nsv Na's
## median 0.000550 0.00165 100 0
## mean 0.001555 0.00429 100 0
## 2.5% 0.000218 0.00056 100 0
## 97.5% 0.013253 0.02321 100 0
##
## risk1 :
## mean sd Min 2.5% 25% 50% 75%
## median 2.72e-05 9.30e-05 1.01e-14 9.10e-14 5.39e-12 1.29e-06 1.36e-05
## mean 3.00e-05 1.13e-04 2.68e-14 1.10e-13 2.53e-09 1.36e-06 1.36e-05
## 2.5% 1.46e-05 3.23e-05 6.33e-15 5.44e-14 2.71e-12 4.41e-07 7.99e-06
## 97.5% 5.93e-05 3.12e-04 8.79e-14 2.21e-13 3.80e-08 2.70e-06 2.00e-05
## 97.5% Max nsv Na's
## median 0.000187 0.000722 100 0
## mean 0.000223 0.000939 100 0
## 2.5% 0.000101 0.000167 100 0
## 97.5% 0.000496 0.002803 100 0
##
## risk :
## mean sd Min 2.5% 25% 50% 75%
## median 2.72e-05 9.30e-05 1.01e-14 9.10e-14 5.39e-12 1.29e-06 1.36e-05
## mean 3.00e-05 1.13e-04 2.68e-14 1.10e-13 2.53e-09 1.36e-06 1.36e-05
## 2.5% 1.46e-05 3.23e-05 6.33e-15 5.44e-14 2.71e-12 4.41e-07 7.99e-06
## 97.5% 5.93e-05 3.12e-04 8.79e-14 2.21e-13 3.80e-08 2.70e-06 2.00e-05
## 97.5% Max nsv Na's
## median 0.000187 0.000722 100 0
## mean 0.000223 0.000939 100 0
## 2.5% 0.000101 0.000167 100 0
## 97.5% 0.000496 0.002803 100 0
##
## yearrisk :
## mean sd Min 2.5% 25% 50% 75% 97.5%
## median 0.00926 0.0305 3.69e-12 3.32e-11 1.97e-09 0.000472 0.00494 0.0659
## mean 0.00988 0.0337 9.79e-12 4.01e-11 9.25e-07 0.000496 0.00497 0.0772
## 2.5% 0.00524 0.0115 2.31e-12 1.99e-11 9.90e-10 0.000161 0.00291 0.0362
## 97.5% 0.01685 0.0747 3.21e-11 8.07e-11 1.39e-05 0.000986 0.00729 0.1654
## Max nsv Na's
## median 0.2318 100 0
## mean 0.2671 100 0
## 2.5% 0.0592 100 0
## 97.5% 0.6391 100 0
##
## risk.raw :
## mean sd Min 2.5% 25% 50% 75%
## median 0.00203 0.00634 2.13e-12 1.07e-11 5.74e-10 9.34e-05 0.001177
## mean 0.00215 0.00722 3.79e-12 1.14e-11 4.16e-07 1.03e-04 0.001206
## 2.5% 0.00113 0.00243 1.71e-12 5.60e-12 2.43e-10 4.37e-05 0.000697
## 97.5% 0.00401 0.01738 1.08e-11 1.97e-11 5.29e-06 1.98e-04 0.001876
## 97.5% Max nsv Na's
## median 0.01589 0.0454 100 0
## mean 0.01811 0.0568 100 0
## 2.5% 0.00718 0.0143 100 0
## 97.5% 0.03917 0.1620 100 0
##
## dose.raw :
## mean sd Min 2.5% 25% 50% 75% 97.5% Max nsv
## median 0.950 2.49 1.47e-08 2.69e-08 2.18e-07 0.200 0.638 7.33 16.18 100
## mean 0.988 2.69 1.58e-08 3.20e-08 1.18e-03 0.204 0.674 7.65 19.75 100
## 2.5% 0.558 1.13 1.47e-08 1.60e-08 1.69e-07 0.104 0.486 3.63 6.22 100
## 97.5% 1.566 5.18 2.25e-08 6.44e-08 1.90e-02 0.322 0.966 13.09 47.48 100
## Na's
## median 0
## mean 0
## 2.5% 0
## 97.5% 0
##
## yearrisk.raw :
## mean sd Min 2.5% 25% 50% 75% 97.5% Max nsv
## median 0.218 0.306 7.76e-10 3.89e-09 2.10e-07 0.0335 0.349 0.995 1.000 100
## mean 0.217 0.308 1.38e-09 4.15e-09 1.51e-04 0.0369 0.353 0.986 1.000 100
## 2.5% 0.179 0.273 6.24e-10 2.05e-09 8.88e-08 0.0158 0.224 0.913 0.994 100
## 97.5% 0.259 0.344 3.94e-09 7.17e-09 1.93e-03 0.0698 0.496 1.000 1.000 100
## Na's
## median 0
## mean 0
## 2.5% 0
## 97.5% 0
##
## risk.chl :
## mean sd Min 2.5% 25% 50% 75%
## median 0.000416 0.001558 2.13e-14 2.08e-13 2.71e-11 4.98e-06 6.47e-05
## mean 0.000433 0.001891 6.90e-14 2.99e-13 6.31e-09 5.32e-06 7.01e-05
## 2.5% 0.000179 0.000512 1.71e-14 1.32e-13 8.58e-12 1.30e-06 2.96e-05
## 97.5% 0.000934 0.004464 2.06e-13 6.54e-13 5.88e-08 1.16e-05 1.70e-04
## 97.5% Max nsv Na's
## median 0.00308 0.01216 100 0
## mean 0.00368 0.01570 100 0
## 2.5% 0.00150 0.00401 100 0
## 97.5% 0.00917 0.03914 100 0
##
## dose.chl :
## mean sd Min 2.5% 25% 50% 75% 97.5% Max
## median 0.2022 0.688 1.53e-10 4.41e-10 7.59e-08 0.00429 0.0912 1.736 5.10
## mean 0.2173 0.771 1.69e-10 4.57e-10 1.58e-05 0.00424 0.0957 1.982 5.96
## 2.5% 0.0972 0.269 1.47e-10 1.83e-10 1.77e-08 0.00190 0.0433 0.615 1.56
## 97.5% 0.3821 1.871 2.89e-10 8.77e-10 2.51e-04 0.00715 0.1682 4.882 16.12
## nsv Na's
## median 100 0
## mean 100 0
## 2.5% 100 0
## 97.5% 100 0
##
## yearrisk.chl :
## mean sd Min 2.5% 25% 50% 75% 97.5% Max
## median 0.0699 0.176 7.78e-12 7.60e-11 9.90e-09 0.001816 0.0233 0.666 0.988
## mean 0.0699 0.177 2.52e-11 1.09e-10 2.30e-06 0.001941 0.0252 0.660 0.948
## 2.5% 0.0453 0.125 6.24e-12 4.82e-11 3.13e-09 0.000475 0.0107 0.415 0.769
## 97.5% 0.0963 0.237 7.50e-11 2.39e-10 2.15e-05 0.004234 0.0601 0.955 1.000
## nsv Na's
## median 100 0
## mean 100 0
## 2.5% 100 0
## 97.5% 100 0
##
## risk.uv :
## mean sd Min 2.5% 25% 50% 75%
## median 2.05e-07 6.45e-07 2.22e-16 1.05e-15 5.74e-14 9.34e-09 1.18e-07
## mean 2.19e-07 7.48e-07 4.04e-16 1.13e-15 4.16e-11 1.03e-08 1.21e-07
## 2.5% 1.13e-07 2.44e-07 2.22e-16 5.69e-16 2.43e-14 4.37e-09 6.97e-08
## 97.5% 4.22e-07 1.88e-06 1.06e-15 1.94e-15 5.29e-10 1.98e-08 1.88e-07
## 97.5% Max nsv Na's
## median 1.60e-06 4.65e-06 100 0
## mean 1.83e-06 5.94e-06 100 0
## 2.5% 7.21e-07 1.44e-06 100 0
## 97.5% 4.00e-06 1.77e-05 100 0
##
## dose.uv :
## mean sd Min 2.5% 25% 50% 75%
## median 9.50e-05 0.000249 1.47e-12 2.69e-12 2.18e-11 2.00e-05 6.38e-05
## mean 9.88e-05 0.000269 1.58e-12 3.20e-12 1.18e-07 2.04e-05 6.74e-05
## 2.5% 5.58e-05 0.000113 1.47e-12 1.60e-12 1.69e-11 1.04e-05 4.86e-05
## 97.5% 1.57e-04 0.000518 2.25e-12 6.44e-12 1.90e-06 3.22e-05 9.66e-05
## 97.5% Max nsv Na's
## median 0.000733 0.001618 100 0
## mean 0.000765 0.001975 100 0
## 2.5% 0.000363 0.000622 100 0
## 97.5% 0.001309 0.004748 100 0
##
## yearrisk.uv :
## mean sd Min 2.5% 25% 50% 75%
## median 7.49e-05 0.000235 8.10e-14 3.82e-13 2.10e-11 3.41e-06 4.30e-05
## mean 8.00e-05 0.000273 1.48e-13 4.13e-13 1.52e-08 3.77e-06 4.40e-05
## 2.5% 4.13e-05 0.000089 8.10e-14 2.08e-13 8.86e-12 1.60e-06 2.54e-05
## 97.5% 1.54e-04 0.000684 3.88e-13 7.09e-13 1.93e-07 7.24e-06 6.85e-05
## 97.5% Max nsv Na's
## median 0.000584 0.001695 100 0
## mean 0.000668 0.002166 100 0
## 2.5% 0.000263 0.000524 100 0
## 97.5% 0.001458 0.006432 100 0
##
## risk.pcf :
## mean sd Min 2.5% 25% 50% 75%
## median 2.72e-05 9.30e-05 1.01e-14 9.10e-14 5.39e-12 1.29e-06 1.36e-05
## mean 3.00e-05 1.13e-04 2.68e-14 1.10e-13 2.53e-09 1.36e-06 1.36e-05
## 2.5% 1.46e-05 3.23e-05 6.33e-15 5.44e-14 2.71e-12 4.41e-07 7.99e-06
## 97.5% 5.93e-05 3.12e-04 8.79e-14 2.21e-13 3.80e-08 2.70e-06 2.00e-05
## 97.5% Max nsv Na's
## median 0.000187 0.000722 100 0
## mean 0.000223 0.000939 100 0
## 2.5% 0.000101 0.000167 100 0
## 97.5% 0.000496 0.002803 100 0
##
## dose.pcf :
## mean sd Min 2.5% 25% 50% 75% 97.5%
## median 0.01438 0.0422 8.83e-11 3.27e-10 2.84e-09 0.00205 0.00800 0.0980
## mean 0.01444 0.0468 1.04e-10 3.48e-10 1.06e-05 0.00206 0.00845 0.1042
## 2.5% 0.00766 0.0155 4.75e-11 1.69e-10 1.66e-09 0.00111 0.00571 0.0465
## 97.5% 0.02449 0.1136 2.69e-10 6.53e-10 1.95e-04 0.00331 0.01344 0.1891
## Max nsv Na's
## median 0.3036 100 0
## mean 0.3809 100 0
## 2.5% 0.0942 100 0
## 97.5% 1.0866 100 0
##
## yearrisk.pcf :
## mean sd Min 2.5% 25% 50% 75% 97.5%
## median 0.00926 0.0305 3.69e-12 3.32e-11 1.97e-09 0.000472 0.00494 0.0659
## mean 0.00988 0.0337 9.79e-12 4.01e-11 9.25e-07 0.000496 0.00497 0.0772
## 2.5% 0.00524 0.0115 2.31e-12 1.99e-11 9.90e-10 0.000161 0.00291 0.0362
## 97.5% 0.01685 0.0747 3.21e-11 8.07e-11 1.39e-05 0.000986 0.00729 0.1654
## Max nsv Na's
## median 0.2318 100 0
## mean 0.2671 100 0
## 2.5% 0.0592 100 0
## 97.5% 0.6391 100 0# Plot the results.
plot(expo.ev1)
Report the median and median of the exposure medians with a 95% confidence interval (CI95).
median.expo <- sapply(1:ndunc(), function(j) median(expo.ev1$expo.mc1[, j, ]))
median(median.expo)## [1] 0.00205099quantile(median.expo, probs = seq(0, 1, 0.025))[c("2.5%", "50%", "97.5%")]## 2.5% 50% 97.5%
## 0.001105943 0.002050990 0.003307286Report the median and median of the DAILY risk probability medians with a 95% confidence interval (CI95).
median.risk <- sapply(1:ndunc(), function(j) median(expo.ev1$risk[, j, ]))
median(median.risk)## [1] 1.294016e-06quantile(median.risk, probs = seq(0, 1, 0.025))[c("2.5%", "50%", "97.5%")]## 2.5% 50% 97.5%
## 4.412978e-07 1.294016e-06 2.702145e-06Report the median and median of the YEARLY risk probability medians with a 95% confidence interval (CI95).
median.yearrisk <- sapply(1:ndunc(), function(j) median(expo.ev1$yearrisk[, j, ]))
median(median.yearrisk)## [1] 0.0004722045quantile(median.yearrisk, probs = seq(0, 1, 0.025))[c("2.5%", "50%", "97.5%")]## 2.5% 50% 97.5%
## 0.0001610607 0.0004722045 0.0009857955Plot 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)
Label the plot:
# 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)
RAW:
Report the median and median of the DAILY risk probability medians with a 95% confidence interval (CI95).
median.risk.raw <- sapply(1:ndunc(), function(j) median(expo.ev1$risk.raw[, j, ]))
median(median.risk.raw)## [1] 9.339484e-05quantile(median.risk.raw, probs = seq(0, 1, 0.025))[c("2.5%", "50%", "97.5%")]## 2.5% 50% 97.5%
## 4.369814e-05 9.339484e-05 1.984299e-04Report the median and median of the YEARLY risk probability medians with a 95% confidence interval (CI95).
median.yearrisk.raw <- sapply(1:ndunc(), function(j) median(expo.ev1$yearrisk.raw[, j, ]))
median(median.yearrisk.raw)## [1] 0.03351579quantile(median.yearrisk.raw, probs = seq(0, 1, 0.025))[c("2.5%", "50%", "97.5%")]## 2.5% 50% 97.5%
## 0.01582241 0.03351579 0.06980533Plot 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.raw)
Label the plot:
# Generate an "ecdf" plot. This actually calls plot.mcnode().
plot(expo.ev1$risk.raw, main = "", ylab = "Proportion of Estimates", xlab = "Daily Probability of Infection")
title("Daily Infec. Risk with Cryptosporidium by RCRW
from Raw Water", cex.main = 1, font.main= 1, line = 7)
# Generate an "ecdf" plot. This actually calls plot.mcnode().
plot(expo.ev1$yearrisk.raw, main = "", ylab = "Proportion of Estimates", xlab = "Yearly Probability of Infection")
title("Yearly Infec. Risk with Cryptosporidium by RCRW
from Raw Water", cex.main = 1, font.main= 1, line = 7)
CHLORINATED:
Report the median and median of the DAILY risk probability medians with a 95% confidence interval (CI95).
median.risk.chl <- sapply(1:ndunc(), function(j) median(expo.ev1$risk.chl[, j, ]))
median(median.risk.chl)## [1] 4.980547e-06quantile(median.risk.chl, probs = seq(0, 1, 0.025))[c("2.5%", "50%", "97.5%")]## 2.5% 50% 97.5%
## 1.300651e-06 4.980547e-06 1.162402e-05Report the median and median of the YEARLY risk probability medians with a 95% confidence interval (CI95).
median.yearrisk.chl <- sapply(1:ndunc(), function(j) median(expo.ev1$yearrisk.chl[, j, ]))
median(median.yearrisk.chl)## [1] 0.001816241quantile(median.yearrisk.chl, probs = seq(0, 1, 0.025))[c("2.5%", "50%", "97.5%")]## 2.5% 50% 97.5%
## 0.0004746239 0.0018162412 0.0042337751Plot 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.chl)
Label the plot:
# Generate an "ecdf" plot. This actually calls plot.mcnode().
plot(expo.ev1$risk.chl, main = "", ylab = "Proportion of Estimates", xlab = "Daily Probability of Infection")
title("Daily Infec. Risk with Cryptosporidium by RCRW
from Chlorinated Water", cex.main = 1, font.main= 1, line = 7)
# Generate an "ecdf" plot. This actually calls plot.mcnode().
plot(expo.ev1$yearrisk.chl, main = "", ylab = "Proportion of Estimates", xlab = "Yearly Probability of Infection")
title("Yearly Infec. Risk with Cryptosporidium by RCRW
from Chlorinated Water", cex.main = 1, font.main= 1, line = 7)
UV TREATED:
Report the median and median of the DAILY risk probability medians with a 95% confidence interval (CI95).
median.risk.uv <- sapply(1:ndunc(), function(j) median(expo.ev1$risk.uv[, j, ]))
median(median.risk.uv)## [1] 9.33992e-09quantile(median.risk.uv, probs = seq(0, 1, 0.025))[c("2.5%", "50%", "97.5%")]## 2.5% 50% 97.5%
## 4.369910e-09 9.339920e-09 1.984501e-08Report the median and median of the YEARLY risk probability medians with a 95% confidence interval (CI95).
median.yearrisk.uv <- sapply(1:ndunc(), function(j) median(expo.ev1$yearrisk.uv[, j, ]))
median(median.yearrisk.uv)## [1] 3.409065e-06quantile(median.yearrisk.uv, probs = seq(0, 1, 0.025))[c("2.5%", "50%", "97.5%")]## 2.5% 50% 97.5%
## 1.595016e-06 3.409065e-06 7.243402e-06Plot 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.uv)
Label the plot:
# Generate an "ecdf" plot. This actually calls plot.mcnode().
plot(expo.ev1$risk.uv, main = "", ylab = "Proportion of Estimates", xlab = "Daily Probability of Infection")
title("Daily Infec. Risk with Cryptosporidium by RCRW
from UV Treated Water", cex.main = 1, font.main= 1, line = 7)
# Generate an "ecdf" plot. This actually calls plot.mcnode().
plot(expo.ev1$yearrisk.uv, main = "", ylab = "Proportion of Estimates", xlab = "Yearly Probability of Infection")
title("Yearly Infec. Risk with Cryptosporidium by RCRW
from UV Treated Water", cex.main = 1, font.main= 1, line = 7)
POROUS CERAMIC FILTRATION:
Report the median and median of the DAILY risk probability medians with a 95% confidence interval (CI95).
median.risk.pcf <- sapply(1:ndunc(), function(j) median(expo.ev1$risk.pcf[, j, ]))
median(median.risk.pcf)## [1] 1.294016e-06quantile(median.risk.pcf, probs = seq(0, 1, 0.025))[c("2.5%", "50%", "97.5%")]## 2.5% 50% 97.5%
## 4.412978e-07 1.294016e-06 2.702145e-06Report the median and median of the YEARLY risk probability medians with a 95% confidence interval (CI95).
median.yearrisk.pcf <- sapply(1:ndunc(), function(j) median(expo.ev1$yearrisk.pcf[, j, ]))
median(median.yearrisk.pcf)## [1] 0.0004722045quantile(median.yearrisk.pcf, probs = seq(0, 1, 0.025))[c("2.5%", "50%", "97.5%")]## 2.5% 50% 97.5%
## 0.0001610607 0.0004722045 0.0009857955Plot 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.pcf)
Label the plot:
# Generate an "ecdf" plot. This actually calls plot.mcnode().
plot(expo.ev1$risk.pcf, main = "", ylab = "Proportion of Estimates", xlab = "Daily Probability of Infection")
title("Daily Infec. Risk with Cryptosporidium by RCRW
from Water Treated with PCF", cex.main = 1, font.main= 1, line = 7)
# Generate an "ecdf" plot. This actually calls plot.mcnode().
plot(expo.ev1$yearrisk.pcf, main = "", ylab = "Proportion of Estimates", xlab = "Yearly Probability of Infection")
title("Yearly Infec. Risk with Cryptosporidium by RCRW
from Water Treated with PCF", cex.main = 1, font.main= 1, line = 7)
Manually remake the graph, this time using yearrisk, which should just be another mcnode object with all values between 0 and 1:
The first one plots yearrisk without a log x scale:
library(reshape)##
## Attaching package: 'reshape'## The following object is masked from 'package:data.table':
##
## meltlibrary(ggplot2)
Risk.mat.1 <- expo.ev1$yearrisk[,,]
probs1 <- c(0.025, 0.25, 0.50, 0.75, 0.975)
quant1 <- as.data.frame(t(apply(Risk.mat.1, 1, quantile, probs = probs1)))
ecdfs1 <- sapply(names(quant1), function(q) ecdf(quant1[[q]]))
plot(ecdfs1[['50%']], main = '')
grays <- c('gray75', 'gray35', 'black', 'gray35', 'gray75')
lines <- mapply(function(e, g) lines(e, col = g), ecdfs1, grays)
Plotting yearrisk ECDF plot with ggplot2:
quant.melt1 <- suppressMessages(melt(quant1))
names(quant.melt1) <- c('q', 'x')
grays <- c('gray75', 'gray35', 'black', 'gray35', 'gray75')
ggplot(quant.melt1, aes(x = x)) + theme_bw() + theme(legend.position = 'none') + geom_hline(yintercept = c(0, 1), linetype = "dashed", color = 'gray') + stat_ecdf(aes(group = q, color = q)) + xlab('x') + ylab('Fn(x)') + scale_colour_manual(values = grays)
Plot Version 2 with ggplot2, and a log-scale x axis:
quant.melt1 <- suppressMessages(melt(quant1))
names(quant.melt1) <- c('q', 'x')
grays <- c('gray75', 'gray35', 'black', 'gray35', 'gray75')
ggplot(quant.melt1, aes(x = x)) + theme_bw() + theme(legend.position = 'none') + geom_hline(yintercept = c(0, 1), linetype = "dashed", color = 'gray') + stat_ecdf(aes(group = q, color = q)) + xlab('x') + ylab('Fn(x)') + scale_colour_manual(values = grays) + scale_x_log10()