predIntNparSimultaneousTestPower.RdCompute the probability that at least one set of future observations violates the given rule based on a nonparametric simultaneous prediction interval for the next \(r\) future sampling occasions. The three possible rules are: \(k\)-of-\(m\), California, or Modified California. The probability is based on assuming the true distribution of the observations is normal.
predIntNparSimultaneousTestPower(n, n.median = 1, k = 1, m = 2, r = 1,
rule = "k.of.m", lpl.rank = ifelse(pi.type == "upper", 0, 1),
n.plus.one.minus.upl.rank = ifelse(pi.type == "lower", 0, 1),
delta.over.sigma = 0, pi.type = "upper", r.shifted = r,
method = "approx", NMC = 100, ci = FALSE, ci.conf.level = 0.95,
integrate.args.list = NULL, evNormOrdStats.method = "royston")vector of positive integers specifying the sample sizes.
Missing (NA), undefined (NaN), and infinite (Inf, -Inf)
values are not allowed.
vector of positive odd integers 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 vector of positive integers
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".
vector of positive integers 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.
vector of positive integers 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).
vector of non-negative integers 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.
vector of non-negative integers 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.
numeric vector indicating the ratio \(\Delta/\sigma\). The quantity
\(\Delta\) (delta) denotes the difference between the mean of the population
that was sampled to construct the prediction interval, and the mean of the
population that will be sampled to produce the future observations. The quantity
\(\sigma\) (sigma) denotes the population standard deviation for both populations.
The default value is delta.over.sigma=0.
character string indicating what kind of prediction interval to compute.
The possible values are "two.sided" (the default), "lower", and
"upper".
vector of positive integers specifying the number of future sampling occasions for
which the scaled mean is shifted by \(\Delta/\sigma\). All values must be
integeters between 1 and the corresponding element of r.
The default value is r.shifted=r.
character string indicating what method to use to compute the power. The possible
values are "approx" (approximation based on predIntNormSimultaneousTestPower; the default) and
"simulate" (Monte Carlo simulation).
positive integer indicating the number of Monte Carlo trials to run when method="simulate". The default value is NMC=100.
logical scalar indicating whether to compute a confidence interval for the power
when method="simulate". The default value is ci=FALSE.
numeric scalar between 0 and 1 indicating the confidence level associated with the
confidence interval for the power. The argument is ignored if ci=FALSE
or method="approx".
list of arguments to supply to the integrate function. The default
value is NULL.
character string indicating which method to use in the call to
evNormOrdStatsScalar when method="approx". The default
value is evNormOrdStats.method="royston". See the DETAILS section for
more information.
What is a Nonparametric Simultaneous Prediction Interval?
A nonparametric prediction interval for some population is an interval on the real line
constructed so that it will contain at least \(k\) of \(m\) future observations from
that population with some specified probability \((1-\alpha)100\%\), where
\(0 < \alpha < 1\) and \(k\) and \(m\) are some pre-specified positive integers
and \(k \le m\). The quantity \((1-\alpha)100\%\) is called
the confidence coefficient or confidence level associated with the prediction
interval. The function predIntNpar computes a standard
nonparametric prediction interval.
The function predIntNparSimultaneous computes a nonparametric simultaneous
prediction interval that will contain a certain number of future observations
with probability \((1-\alpha)100\%\) for each of \(r\) future sampling
“occasions”,
where \(r\) is some pre-specified positive integer. The quantity \(r\) may
refer to \(r\) distinct future sampling occasions in time, or it may for example
refer to sampling at \(r\) distinct locations on one future sampling occasion,
assuming that the population standard deviation is the same at all of the \(r\)
distinct locations.
The function predIntNparSimultaneous computes a nonparametric simultaneous
prediction interval based on one of three possible rules:
For the \(k\)-of-\(m\) rule (rule="k.of.m"), at least \(k\) of
the next \(m\) future observations will fall in the prediction
interval with probability \((1-\alpha)100\%\) on each of the \(r\) future
sampling occasions. If obserations are being taken sequentially, for a particular
sampling occasion, up to \(m\) observations may be taken, but once
\(k\) of the observations fall within the prediction interval, sampling can stop.
Note: For this rule, when \(r=1\), the results of predIntNparSimultaneous
are equivalent to the results of predIntNpar.
For the California rule (rule="CA"), with probability
\((1-\alpha)100\%\), for each of the \(r\) future sampling occasions, either
the first observation will fall in the prediction interval, or else all of the next
\(m-1\) observations will fall in the prediction interval. That is, if the first
observation falls in the prediction interval then sampling can stop. Otherwise,
\(m-1\) more observations must be taken.
For the Modified California rule (rule="Modified.CA"), with probability
\((1-\alpha)100\%\), for each of the \(r\) future sampling occasions, either the
first observation will fall in the prediction interval, or else at least 2 out of
the next 3 observations will fall in the prediction interval. That is, if the first
observation falls in the prediction interval then sampling can stop. Otherwise, up
to 3 more observations must be taken.
Nonparametric simultaneous prediction intervals can be extended to using medians
in place of single observations (USEPA, 2009, Chapter 19). That is, you can
create a nonparametric simultaneous prediction interval that will contain a
specified number of medians (based on which rule you choose) on each of \(r\)
future sampling occassions, where each each median is based on \(b\) individual
observations. For the function predIntNparSimultaneous, the argument
n.median corresponds to \(b\).
The Form of a Nonparametric Prediction Interval
Let \(\underline{x} = x_1, x_2, \ldots, x_n\) denote a vector of \(n\)
independent observations from some continuous distribution, and let
\(x_{(i)}\) denote the the \(i\)'th order statistics in \(\underline{x}\).
A two-sided nonparametric prediction interval is constructed as:
$$[x_{(u)}, x_{(v)}] \;\;\;\;\;\; (1)$$
where \(u\) and \(v\) are positive integers between 1 and \(n\), and
\(u < v\). That is, \(u\) denotes the rank of the lower prediction limit, and
\(v\) denotes the rank of the upper prediction limit. To make it easier to write
some equations later on, we can also write the prediction interval (1) in a slightly
different way as:
$$[x_{(u)}, x_{(n + 1 - w)}] \;\;\;\;\;\; (2)$$
where
$$w = n + 1 - v \;\;\;\;\;\; (3)$$
so that \(w\) is a positive integer between 1 and \(n-1\), and
\(u < n+1-w\). In terms of the arguments to the function predIntNparSimultaneous,
the argument lpl.rank corresponds to \(u\), and the argument
n.plus.one.minus.upl.rank corresponds to \(w\).
If we allow \(u=0\) and \(w=0\) and define lower and upper bounds as: $$x_{(0)} = lb \;\;\;\;\;\; (4)$$ $$x_{(n+1)} = ub \;\;\;\;\;\; (5)$$ then Equation (2) above can also represent a one-sided lower or one-sided upper prediction interval as well. That is, a one-sided lower nonparametric prediction interval is constructed as: $$[x_{(u)}, x_{(n + 1)}] = [x_{(u)}, ub] \;\;\;\;\;\; (6)$$ and a one-sided upper nonparametric prediction interval is constructed as: $$[x_{(0)}, x_{(n + 1 - w)}] = [lb, x_{(n + 1 - w)}] \;\;\;\;\;\; (7)$$ Usually, \(lb = -\infty\) or \(lb = 0\) and \(ub = \infty\).
Note: For nonparametric simultaneous prediction intervals, only lower
(pi.type="lower") and upper (pi.type="upper") prediction
intervals are available.
Computing Power
The "power" of the prediction interval is defined as the probability that
at least one set of future observations violates the given rule based on a
simultaneous prediction interval for the next \(r\) future sampling occasions,
where the population for the future observations is allowed to differ from
the population for the observations used to construct the prediction interval.
For the function predIntNparSimultaneousTestPower, power is computed assuming
both the background and future the observations come from normal distributions with
the same standard deviation, but the means of the distributions are allowed to differ.
The quantity \(\Delta\) (upper case delta) denotes the difference between the
mean of the population that was sampled to construct the prediction interval, and
the mean of the population that will be sampled to produce the future observations.
The quantity \(\sigma\) (sigma) denotes the population standard deviation of both
of these populations. The argument delta.over.sigma corresponds to the
quantity \(\Delta/\sigma\).
Approximate Power (method="approx")
Based on Gansecki (2009), the power of a nonparametric simultaneous prediction
interval when the underlying observations come from a nomral distribution
can be approximated by the power of a normal simultaneous prediction
interval (see predIntNormSimultaneousTestPower) where the multiplier
\(K\) is replaced with the expected value of the normal order statistic that
corresponds to the rank of the order statistic used for the upper or lower bound
of the prediction interval. Gansecki (2009) uses the approximation:
$$K = \Phi^{-1}(\frac{i - 0.5}{n}) \;\;\;\;\;\; (8)$$
where \(\Phi\) denotes the cumulative distribution function of the standard
normal distribution and \(i\) denotes the rank of the order statistic used
as the prediction limit. By default, the value of the argument evNormOrdStats.method="royston", so the function
predIntNparSimultaneousTestPower uses the exact value of the
expected value of the normal order statistic in the call to
evNormOrdStatsScalar. You can change the
method of computing the expected value of the normal order statistic by
changing the value of the argument evNormOrdStats.method.
Power Based on Monte Carlo Simulation (method="simulate")
When method="simulate", the power of the nonparametric simultaneous
prediction interval is estimated based on a Monte Carlo simulation. The argument
NMC determines the number of Monte Carlo trials. If ci=TRUE, a
confidence interval for the power is created based on the NMC Monte Carlo
estimates of power.
vector of values between 0 and 1 equal to the probability that the rule will be violated.
See the help file for predIntNparSimultaneous.
Gansecki, M. (2009). Using the Optimal Rank Values Calculator. US Environmental Protection Agency, Region 8, March 10, 2009.
See the help file for predIntNparSimultaneous.
In the course of designing a sampling program, an environmental scientist may wish
to determine the relationship between sample size, significance level, power, and
scaled difference if one of the objectives of the sampling program is to determine
whether two distributions differ from each other. The functions
predIntNparSimultaneousTestPower and plotPredIntNparSimultaneousTestPowerCurve can be
used to investigate these relationships for the case of normally-distributed
observations.
# 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 levels and powers 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.
# Power will be computed assuming a normal distribution and setting
# delta.over.sigma equal to 2, 3, and 4.
# 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
# So the required confidence level is 0.98, or 98%.
# Now determine the confidence level associated with each plan.
# Note that both plans achieve the required confidence level.
# 1) the 1-of-2 retesting plan for a median of order 3 using the
# background maximum
predIntNparSimultaneousConfLevel(n = 20, n.median = 3, k = 1, m = 2, r = r)
#> [1] 0.9940354
#[1] 0.9940354
# 2) the 1-of-4 plan based on individual observations using the 3rd highest
# background value.
predIntNparSimultaneousConfLevel(n = 20, k = 1, m = 4, r = r,
n.plus.one.minus.upl.rank = 3)
#> [1] 0.9864909
#[1] 0.9864909
#------------------------------------------------------------------------------
# Compute approximate power of each plan to detect contamination at just 1 well
# assuming true underying distribution of Hg is Normal at all wells and
# using delta.over.sigma equal to 2, 3, and 4.
#------------------------------------------------------------------------------
# Computer aproximate power for
# 1) the 1-of-2 retesting plan for a median of order 3 using the
# background maximum
predIntNparSimultaneousTestPower(n = 20, n.median = 3, k = 1, m = 2, r = r,
delta.over.sigma = 2:4, r.shifted = 1)
#> [1] 0.3953712 0.9129671 0.9983054
#[1] 0.3953712 0.9129671 0.9983054
# Compute approximate power for
# 2) the 1-of-4 plan based on individual observations using the 3rd highest
# background value.
predIntNparSimultaneousTestPower(n = 20, k = 1, m = 4, r = r,
n.plus.one.minus.upl.rank = 3, delta.over.sigma = 2:4, r.shifted = 1)
#> [1] 0.4367972 0.8694664 0.9888779
#[1] 0.4367972 0.8694664 0.9888779
#----------
if (FALSE) { # \dontrun{
# Compare estimated power using approximation method with estimated power
# using Monte Carlo simulation for the 1-of-4 plan based on individual
# observations using the 3rd highest background value.
predIntNparSimultaneousTestPower(n = 20, k = 1, m = 4, r = r,
n.plus.one.minus.upl.rank = 3, delta.over.sigma = 2:4, r.shifted = 1,
method = "simulate", ci = TRUE, NMC = 1000)
#[1] 0.437 0.863 0.989
#attr(,"conf.int")
# [,1] [,2] [,3]
#LCL 0.4111999 0.8451148 0.9835747
#UCL 0.4628001 0.8808852 0.9944253
} # }
#==========
# Cleanup
#--------
rm(n, nw, nc, ne, r, conf.level)