R/predIntNparSimultaneous.R
predIntNparSimultaneous.RdConstruct a nonparametric simultaneous prediction interval for the next \(r\) sampling “occasions” based on one of three possible rules: k-of-m, California, or Modified California. The simultaneous prediction interval assumes the observations from from a continuous distribution.
Extended methods and derivations are available in
Prediction intervals for environmental monitoring: predIntNparSimultaneous.
Open the vignette locally with
vignette("prediction-intervals", package = "EnvStats").
a numeric vector of observations. Missing (NA), undefined (NaN), and
infinite (Inf, -Inf) values are allowed but will be removed.
positive odd integer specifying the sample size associated with the future medians.
The default value is n.median=1 (i.e., individual observations). Note that
all future medians must be based on the same sample size.
for the \(k\)-of-\(m\) rule (rule="k.of.m"), a positive integer
specifying the minimum number of observations (or medians) out of \(m\)
observations (or medians) (all obtained on one future sampling “occassion”)
the prediction interval should contain.
The default value is k=1. This argument is ignored when the argument
rule is not equal to "k.of.m".
positive integer specifying the maximum number of future observations (or
medians) on one future sampling “occasion”.
The default value is m=2, except when rule="Modified.CA", in which
case this argument is ignored and m is automatically set equal to 4.
positive integer specifying the number of future sampling “occasions”.
The default value is r=1.
character string specifying which rule to use. The possible values are
"k.of.m" (\(k\)-of-\(m\) rule; the default), "CA" (California rule),
and "Modified.CA" (modified California rule).
See the DETAILS section below for more information.
positive integer indicating the rank of the order statistic to use for the lower
bound of the prediction interval. When pi.type="lower", the default value
is lpl.rank=1 (implying the minimum value of x is used as the lower
bound of the prediction interval). When pi.type="upper", the argument
lpl.rank is set equal to 0 and the value of lb is used as the
lower bound of the tolerance interval.
positive integer related to the rank of the order statistic to use for the upper
bound of the prediction interval. A value of n.plus.one.minus.upl.rank=1
(the default) means use the first largest value, and in general a value of n.plus.one.minus.upl.rank=\(i\) means use the \(i\)'th largest value.
When pi.type="lower", the argument n.plus.one.minus.upl.rank is set
equal to 0 and the value of ub is used as the upper bound of the
prediction interval.
scalars indicating lower and upper bounds on the distribution. By default,
lb=-Inf and ub=Inf. If you are constructing a prediction interval for
a distribution that you know has a lower bound other than -Inf
(e.g., 0), set lb to this value. Similarly, if you know the
distribution has an upper bound other than Inf, set ub to this value.
The argument lb is ignored if pi.type="two-sided" or
pi.type="lower". The argument ub is ignored if
pi.type="two-sided" or pi.type="upper".
character string indicating what kind of prediction interval to compute.
The possible values are "upper" (the default) and "lower".
a list of arguments to supply to the integrate function. The
default value is integrate.args.list=NULL which means that the
default values of integrate are used.
a list of class "estimate" containing the simultaneous prediction interval
and other information. See the help file for estimate.object for
details.
Cameron, Kirk. (2011). Personal communication, February 16, 2011. MacStat Consulting, Ltd., Colorado Springs, Colorado.
Chew, V. (1968). Simultaneous Prediction Intervals. Technometrics 10(2), 323–331.
Danziger, L., and S. Davis. (1964). Tables of Distribution-Free Tolerance Limits. Annals of Mathematical Statistics 35(5), 1361–1365.
Davis, C.B. (1994). Environmental Regulatory Statistics. In Patil, G.P., and C.R. Rao, eds., Handbook of Statistics, Vol. 12: Environmental Statistics. North-Holland, Amsterdam, a division of Elsevier, New York, NY, Chapter 26, 817–865.
Davis, C.B., and R.J. McNichols. (1987). One-sided Intervals for at Least p of m Observations from a Normal Population on Each of r Future Occasions. Technometrics 29, 359–370.
Davis, C.B., and R.J. McNichols. (1994a). Ground Water Monitoring Statistics Update: Part I: Progress Since 1988. Ground Water Monitoring and Remediation 14(4), 148–158.
Davis, C.B., and R.J. McNichols. (1994b). Ground Water Monitoring Statistics Update: Part II: Nonparametric Prediction Limits. Ground Water Monitoring and Remediation 14(4), 159–175.
Davis, C.B., and R.J. McNichols. (1999). Simultaneous Nonparametric Prediction Limits (with Discusson). Technometrics 41(2), 89–112.
Gibbons, R.D. (1987a). Statistical Prediction Intervals for the Evaluation of Ground-Water Quality. Ground Water 25, 455–465.
Gibbons, R.D. (1991b). Statistical Tolerance Limits for Ground-Water Monitoring. Ground Water 29, 563–570.
Gibbons, R.D., and J. Baker. (1991). The Properties of Various Statistical Prediction Intervals for Ground-Water Detection Monitoring. Journal of Environmental Science and Health A26(4), 535–553.
Gibbons, R.D., D.K. Bhaumik, and S. Aryal. (2009). Statistical Methods for Groundwater Monitoring, Second Edition. John Wiley & Sons, Hoboken.
Hahn, G.J., and W.Q. Meeker. (1991). Statistical Intervals: A Guide for Practitioners. John Wiley and Sons, New York, 392pp.
Hahn, G., and W. Nelson. (1973). A Survey of Prediction Intervals and Their Applications. Journal of Quality Technology 5, 178–188.
Hall, I.J., R.R. Prairie, and C.K. Motlagh. (1975). Non-Parametric Prediction Intervals. Journal of Quality Technology 7(3), 109–114.
Millard, S.P., and Neerchal, N.K. (2001). Environmental Statistics with S-PLUS. CRC Press, Boca Raton, Florida.
USEPA. (2009). Statistical Analysis of Groundwater Monitoring Data at RCRA Facilities, Unified Guidance. EPA 530/R-09-007, March 2009. Office of Resource Conservation and Recovery Program Implementation and Information Division. U.S. Environmental Protection Agency, Washington, D.C.
USEPA. (2010). Errata Sheet - March 2009 Unified Guidance. EPA 530/R-09-007a, August 9, 2010. Office of Resource Conservation and Recovery, Program Information and Implementation Division. U.S. Environmental Protection Agency, Washington, D.C.
Prediction and tolerance intervals have long been applied to quality control and life testing problems (Hahn, 1970b,c; Hahn and Nelson, 1973; Krishnamoorthy and Mathew, 2009). In the context of environmental statistics, prediction intervals are useful for analyzing data from groundwater detection monitoring programs at hazardous and solid waste facilities (e.g., Gibbons et al., 2009; Millard and Neerchal, 2001; USEPA, 2009).
# Generate 20 observations from a lognormal mixture distribution with
# parameters mean1=1, cv1=0.5, mean2=5, cv2=1, and p.mix=0.1. Use
# predIntNparSimultaneous to construct an upper one-sided prediction interval
# using the maximum observed value using the 1-of-3 rule.
# (Note: the call to set.seed simply allows you to reproduce this example.)
set.seed(250)
dat <- rlnormMixAlt(n = 20, mean1 = 1, cv1 = 0.5,
mean2 = 5, cv2 = 1, p.mix = 0.1)
predIntNparSimultaneous(dat, k = 1, m = 3, lb = 0)
#>
#> Results of Distribution Parameter Estimation
#> --------------------------------------------
#>
#> Assumed Distribution: None
#>
#> Data: dat
#>
#> Sample Size: 20
#>
#> Prediction Interval Method: exact
#>
#> Prediction Interval Type: upper
#>
#> Confidence Level: 99.94353%
#>
#> Prediction Limit Rank(s): 20
#>
#> Minimum Number of
#> Future Observations
#> Interval Should Contain: 1
#>
#> Total Number of
#> Future Observations: 3
#>
#> Prediction Interval: LPL = 0.000000
#> UPL = 1.817311
#>
#Results of Distribution Parameter Estimation
#--------------------------------------------
#
#Assumed Distribution: None
#
#Data: dat
#
#Sample Size: 20
#
#Prediction Interval Method: exact
#
#Prediction Interval Type: upper
#
#Confidence Level: 99.94353%
#
#Prediction Limit Rank(s): 20
#
#Minimum Number of
#Future Observations
#Interval Should Contain: 1
#
#Total Number of
#Future Observations: 3
#
#Prediction Interval: LPL = 0.000000
# UPL = 1.817311
#----------
# Compare the confidence levels for the 1-of-3 rule, California Rule, and
# Modified California Rule.
predIntNparSimultaneous(dat, k = 1, m = 3, lb = 0)$interval$conf.level
#> [1] 0.9994353
#[1] 0.9994353
predIntNparSimultaneous(dat, m = 3, rule = "CA", lb = 0)$interval$conf.level
#> [1] 0.9919066
#[1] 0.9919066
predIntNparSimultaneous(dat, rule = "Modified.CA", lb = 0)$interval$conf.level
#> [1] 0.9984943
#[1] 0.9984943
#=========
# Repeat the above example, but create the baseline data using just
# n=8 observations and set r to 4 future sampling occasions
set.seed(598)
dat <- rlnormMixAlt(n = 8, mean1 = 1, cv1 = 0.5,
mean2 = 5, cv2 = 1, p.mix = 0.1)
predIntNparSimultaneous(dat, k = 1, m = 3, r = 4, lb = 0)
#>
#> Results of Distribution Parameter Estimation
#> --------------------------------------------
#>
#> Assumed Distribution: None
#>
#> Data: dat
#>
#> Sample Size: 8
#>
#> Prediction Interval Method: exact
#>
#> Prediction Interval Type: upper
#>
#> Confidence Level: 97.7599%
#>
#> Prediction Limit Rank(s): 8
#>
#> Minimum Number of
#> Future Observations
#> Interval Should Contain
#> (per Sampling Occasion): 1
#>
#> Total Number of
#> Future Observations
#> (per Sampling Occasion): 3
#>
#> Number of Future
#> Sampling Occasions: 4
#>
#> Prediction Interval: LPL = 0.000000
#> UPL = 5.683453
#>
#Results of Distribution Parameter Estimation
#--------------------------------------------
#
#Assumed Distribution: None
#
#Data: dat
#
#Sample Size: 8
#
#Prediction Interval Method: exact
#
#Prediction Interval Type: upper
#
#Confidence Level: 97.7599%
#
#Prediction Limit Rank(s): 8
#
#Minimum Number of
#Future Observations
#Interval Should Contain
#(per Sampling Occasion): 1
#
#Total Number of
#Future Observations
#(per Sampling Occasion): 3
#
#Number of Future
#Sampling Occasions: 4
#
#Prediction Interval: LPL = 0.000000
# UPL = 5.683453
#----------
# Compare the confidence levels for the 1-of-3 rule, California Rule, and
# Modified California Rule.
predIntNparSimultaneous(dat, k = 1, m = 3, r = 4, lb = 0)$interval$conf.level
#> [1] 0.977599
#[1] 0.977599
predIntNparSimultaneous(dat, m = 3, r = 4, rule = "CA", lb = 0)$interval$conf.level
#> [1] 0.8737798
#[1] 0.8737798
predIntNparSimultaneous(dat, r = 4, rule = "Modified.CA", lb = 0)$interval$conf.level
#> [1] 0.9510178
#[1] 0.9510178
#==========
# Example 19-5 of USEPA (2009, p. 19-33) shows how to compute nonparametric upper
# simultaneous prediction limits for various rules based on trace mercury data (ppb)
# collected in the past year from a site with four background wells and 10 compliance
# wells (data for two of the compliance wells are shown in the guidance document).
# The facility must monitor the 10 compliance wells for five constituents
# (including mercury) annually.
# Here we will compute the confidence level associated with two different sampling plans:
# 1) the 1-of-2 retesting plan for a median of order 3 using the background maximum and
# 2) the 1-of-4 plan on individual observations using the 3rd highest background value.
# The data for this example are stored in EPA.09.Ex.19.5.mercury.df.
# We will pool data from 4 background wells that were sampled on
# a number of different occasions, giving us a sample size of
# n = 20 to use to construct the prediction limit.
# There are 10 compliance wells and we will monitor 5 different
# constituents at each well annually. For this example, USEPA (2009)
# recommends setting r to the product of the number of compliance wells and
# the number of evaluations per year.
# To determine the minimum confidence level we require for
# the simultaneous prediction interval, USEPA (2009) recommends
# setting the maximum allowed individual Type I Error level per constituent to:
# 1 - (1 - SWFPR)^(1 / Number of Constituents)
# which translates to setting the confidence limit to
# (1 - SWFPR)^(1 / Number of Constituents)
# where SWFPR = site-wide false positive rate. For this example, we
# will set SWFPR = 0.1. Thus, the required individual Type I Error level
# and confidence level per constituent are given as follows:
# n = 20 based on 4 Background Wells
# nw = 10 Compliance Wells
# nc = 5 Constituents
# ne = 1 Evaluation per year
n <- 20
nw <- 10
nc <- 5
ne <- 1
# Set number of future sampling occasions r to
# Number Compliance Wells x Number Evaluations per Year
r <- nw * ne
conf.level <- (1 - 0.1)^(1 / nc)
conf.level
#> [1] 0.9791484
#[1] 0.9791484
alpha <- 1 - conf.level
alpha
#> [1] 0.02085164
#[1] 0.02085164
#----------
# Look at the data:
head(EPA.09.Ex.19.5.mercury.df)
#> Event Well Well.type Mercury.ppb.orig Mercury.ppb Censored
#> 1 1 BG-1 Background 0.21 0.21 FALSE
#> 2 2 BG-1 Background <.2 0.20 TRUE
#> 3 3 BG-1 Background <.2 0.20 TRUE
#> 4 4 BG-1 Background <.2 0.20 TRUE
#> 5 5 BG-1 Background <.2 0.20 TRUE
#> 6 6 BG-1 Background NA FALSE
# Event Well Well.type Mercury.ppb.orig Mercury.ppb Censored
#1 1 BG-1 Background 0.21 0.21 FALSE
#2 2 BG-1 Background <.2 0.20 TRUE
#3 3 BG-1 Background <.2 0.20 TRUE
#4 4 BG-1 Background <.2 0.20 TRUE
#5 5 BG-1 Background <.2 0.20 TRUE
#6 6 BG-1 Background NA FALSE
longToWide(EPA.09.Ex.19.5.mercury.df, "Mercury.ppb.orig",
"Event", "Well", paste.row.name = TRUE)
#> BG-1 BG-2 BG-3 BG-4 CW-1 CW-2
#> Event.1 0.21 <.2 <.2 <.2 0.22 0.36
#> Event.2 <.2 <.2 0.23 0.25 0.2 0.41
#> Event.3 <.2 <.2 <.2 0.28 <.2 0.28
#> Event.4 <.2 0.21 0.23 <.2 0.25 0.45
#> Event.5 <.2 <.2 0.24 <.2 0.24 0.43
#> Event.6 <.2 0.54
# BG-1 BG-2 BG-3 BG-4 CW-1 CW-2
#Event.1 0.21 <.2 <.2 <.2 0.22 0.36
#Event.2 <.2 <.2 0.23 0.25 0.2 0.41
#Event.3 <.2 <.2 <.2 0.28 <.2 0.28
#Event.4 <.2 0.21 0.23 <.2 0.25 0.45
#Event.5 <.2 <.2 0.24 <.2 0.24 0.43
#Event.6 <.2 0.54
# Construct the upper simultaneous prediction limit using the 1-of-2
# retesting plan for a median of order 3 based on the background maximum
Hg.Back <- with(EPA.09.Ex.19.5.mercury.df,
Mercury.ppb[Well.type == "Background"])
pred.int.1.of.2.med.3 <- predIntNparSimultaneous(Hg.Back, n.median = 3,
k = 1, m = 2, r = r, lb = 0)
#> Warning: There were 4 nonfinite values in x : 4 NA's
#> Warning: 4 observations with NA/NaN/Inf in 'x' removed.
pred.int.1.of.2.med.3
#>
#> Results of Distribution Parameter Estimation
#> --------------------------------------------
#>
#> Assumed Distribution: None
#>
#> Data: Hg.Back
#>
#> Sample Size: 20
#>
#> Number NA/NaN/Inf's: 4
#>
#> Prediction Interval Method: exact
#>
#> Prediction Interval Type: upper
#>
#> Confidence Level: 99.40354%
#>
#> Prediction Limit Rank(s): 20
#>
#> Minimum Number of
#> Future Medians
#> Interval Should Contain
#> (per Sampling Occasion): 1
#>
#> Total Number of
#> Future Medians
#> (per Sampling Occasion): 2
#>
#> Number of Future
#> Sampling Occasions: 10
#>
#> Sample Size for Medians: 3
#>
#> Prediction Interval: LPL = 0.00
#> UPL = 0.28
#>
#Results of Distribution Parameter Estimation
#--------------------------------------------
#
#Assumed Distribution: None
#
#Data: Hg.Back
#
#Sample Size: 20
#
#Number NA/NaN/Inf's: 4
#
#Prediction Interval Method: exact
#
#Prediction Interval Type: upper
#
#Confidence Level: 99.40354%
#
#Prediction Limit Rank(s): 20
#
#Minimum Number of
#Future Medians
#Interval Should Contain
#(per Sampling Occasion): 1
#
#Total Number of
#Future Medians
#(per Sampling Occasion): 2
#
#Number of Future
#Sampling Occasions: 10
#
#Sample Size for Medians: 3
#
#Prediction Interval: LPL = 0.00
# UPL = 0.28
# Note that the achieved confidence level of 99.4% is greater than the
# required confidence level of 97.9%.
# Now determine whether either compliance well indicates evidence of
# Mercury contamination.
# Compliance Well 1
#------------------
Hg.CW.1 <- with(EPA.09.Ex.19.5.mercury.df, Mercury.ppb.orig[Well == "CW-1"])
Hg.CW.1
#> [1] "0.22" "0.2" "<.2" "0.25" "0.24" "<.2"
#[1] "0.22" "0.2" "<.2" "0.25" "0.24" "<.2"
# The median of the first 3 observations is 0.2, which is less than
# the UPL of 0.28, so there is no evidence of contamination.
# Compliance Well 2
#------------------
Hg.CW.2 <- with(EPA.09.Ex.19.5.mercury.df, Mercury.ppb.orig[Well == "CW-2"])
Hg.CW.2
#> [1] "0.36" "0.41" "0.28" "0.45" "0.43" "0.54"
#[1] "0.36" "0.41" "0.28" "0.45" "0.43" "0.54"
# The median of the first 3 observations is 0.36, so 3 more observations have to
# be looked at. The median of the second 3 observations is 0.45, which is
# larger than the UPL of 0.28, so there is evidence of contamination.
#----------
# Now create the upper simultaneous prediction limit using the 1-of-4 plan
# on individual observations using the 3rd highest background value.
pred.int.1.of.4.3rd <- predIntNparSimultaneous(Hg.Back, k = 1, m = 4,
r = r, lb = 0, n.plus.one.minus.upl.rank = 3)
#> Warning: There were 4 nonfinite values in x : 4 NA's
#> Warning: 4 observations with NA/NaN/Inf in 'x' removed.
pred.int.1.of.4.3rd
#>
#> Results of Distribution Parameter Estimation
#> --------------------------------------------
#>
#> Assumed Distribution: None
#>
#> Data: Hg.Back
#>
#> Sample Size: 20
#>
#> Number NA/NaN/Inf's: 4
#>
#> Prediction Interval Method: exact
#>
#> Prediction Interval Type: upper
#>
#> Confidence Level: 98.64909%
#>
#> Prediction Limit Rank(s): 18
#>
#> Minimum Number of
#> Future Observations
#> Interval Should Contain
#> (per Sampling Occasion): 1
#>
#> Total Number of
#> Future Observations
#> (per Sampling Occasion): 4
#>
#> Number of Future
#> Sampling Occasions: 10
#>
#> Prediction Interval: LPL = 0.00
#> UPL = 0.24
#>
#Results of Distribution Parameter Estimation
#--------------------------------------------
#
#Assumed Distribution: None
#
#Data: Hg.Back
#
#Sample Size: 20
#
#Number NA/NaN/Inf's: 4
#
#Prediction Interval Method: exact
#
#Prediction Interval Type: upper
#
#Confidence Level: 98.64909%
#
#Prediction Limit Rank(s): 18
#
#Minimum Number of
#Future Observations
#Interval Should Contain
#(per Sampling Occasion): 1
#
#Total Number of
#Future Observations
#(per Sampling Occasion): 4
#
#Number of Future
#Sampling Occasions: 10
#
#Prediction Interval: LPL = 0.00
# UPL = 0.24
# Note that the achieved confidence level of 98.6% is greater than the
# required confidence level of 97.9%.
# Now determine whether either compliance well indicates evidence of
# Mercury contamination.
# Compliance Well 1
#------------------
Hg.CW.1 <- with(EPA.09.Ex.19.5.mercury.df, Mercury.ppb.orig[Well == "CW-1"])
Hg.CW.1
#> [1] "0.22" "0.2" "<.2" "0.25" "0.24" "<.2"
#[1] "0.22" "0.2" "<.2" "0.25" "0.24" "<.2"
# The first observation is less than the UPL of 0.24, which is less than
# the UPL of 0.28, so there is no evidence of contamination.
# Compliance Well 2
#------------------
Hg.CW.2 <- with(EPA.09.Ex.19.5.mercury.df, Mercury.ppb.orig[Well == "CW-2"])
Hg.CW.2
#> [1] "0.36" "0.41" "0.28" "0.45" "0.43" "0.54"
#[1] "0.36" "0.41" "0.28" "0.45" "0.43" "0.54"
# All of the first 4 observations are greater than the UPL of 0.24, so there
# is evidence of contamination.
#==========
# Cleanup
#--------
rm(dat, n, nw, nc, ne, r, conf.level, alpha, Hg.Back, pred.int.1.of.2.med.3,
pred.int.1.of.4.3rd, Hg.CW.1, Hg.CW.2)