predIntNormSimultaneousTestPower {EnvStats}  R Documentation 
Compute 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 for a normal distribution. The three possible rules are:
k
ofm
, California, or Modified California.
predIntNormSimultaneousTestPower(n, df = n  1, n.mean = 1, k = 1, m = 2, r = 1,
rule = "k.of.m", delta.over.sigma = 0, pi.type = "upper", conf.level = 0.95,
r.shifted = r, K.tol = .Machine$double.eps^0.5, integrate.args.list = NULL)
n 
vector of positive integers greater than 2 indicating the sample size upon which the prediction interval is based. 
df 
vector of positive integers indicating the degrees of freedom associated with
the sample size. The default value is 
n.mean 
positive integer specifying the sample size associated with the future averages.
The default value is 
k 
for the 
m 
vector of positive integers specifying the maximum number of future observations (or
averages) on one future sampling “occasion”.
The default value is 
r 
vector of positive integers specifying the number of future sampling “occasions”.
The default value is 
rule 
character string specifying which rule to use. The possible values are

delta.over.sigma 
numeric vector indicating the ratio 
pi.type 
character string indicating what kind of prediction interval to compute.
The possible values are 
conf.level 
vector of values between 0 and 1 indicating the confidence level of the prediction interval.
The default value is 
r.shifted 
vector of positive integers specifying the number of future sampling occasions for
which the scaled mean is shifted by 
K.tol 
numeric scalar indicating the tolerance to use in the nonlinear search algorithm to
compute 
integrate.args.list 
a list of arguments to supply to the 
What is a Simultaneous Prediction Interval?
A prediction interval for some population is an interval on the real line constructed
so that it will contain k
future observations from that population
with some specified probability (1\alpha)100\%
, where
0 < \alpha < 1
and k
is some prespecified positive integer.
The quantity (1\alpha)100\%
is called
the confidence coefficient or confidence level associated with the prediction
interval. The function predIntNorm
computes a standard prediction
interval based on a sample from a normal distribution.
The function predIntNormSimultaneous
computes a 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 prespecified 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 predIntNormSimultaneous
computes a simultaneous prediction
interval based on one of three possible rules:
For the k
ofm
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: When k=m
and r=1
, the results of
predIntNormSimultaneous
are equivalent to the results of
predIntNorm
.
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
m1
observations will fall in the prediction interval. That is, if the first
observation falls in the prediction interval then sampling can stop. Otherwise,
m1
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.
Simultaneous prediction intervals can be extended to using averages (means) in place
of single observations (USEPA, 2009, Chapter 19). That is, you can create a
simultaneous prediction interval
that will contain a specified number of averages (based on which rule you choose) on
each of r
future sampling occassions, where each each average is based on
w
individual observations. For the function predIntNormSimultaneous
,
the argument n.mean
corresponds to w
.
The Form of a Prediction Interval
Let \underline{x} = x_1, x_2, \ldots, x_n
denote a vector of n
observations from a normal distribution with parameters
mean=
\mu
and sd=
\sigma
. Also, let w
denote the
sample size associated with the future averages (i.e., n.mean=
w
).
When w=1
, each average is really just a single observation, so in the rest of
this help file the term “averages” will replace the phrase
“observations or averages”.
For a normal distribution, the form of a twosided (1\alpha)100\%
prediction interval is:
[\bar{x}  Ks, \bar{x} + Ks] \;\;\;\;\;\; (1)
where \bar{x}
denotes the sample mean:
\bar{x} = \frac{1}{n} \sum_{i=1}^n x_i \;\;\;\;\;\; (2)
s
denotes the sample standard deviation:
s^2 = \frac{1}{n1} \sum_{i=1}^n (x_i  \bar{x})^2 \;\;\;\;\;\; (3)
and K
denotes a constant that depends on the sample size n
, the
confidence level, the number of future sampling occassions r
, and the
sample size associated with the future averages, w
. Do not confuse the
constant K
(uppercase K) with the number of future averages k
(lowercase k) in the k
ofm
rule. The symbol K
is used here
to be consistent with the notation used for tolerance intervals
(see tolIntNorm
).
Similarly, the form of a onesided lower prediction interval is:
[\bar{x}  Ks, \infty] \;\;\;\;\;\; (4)
and the form of a onesided upper prediction interval is:
[\infty, \bar{x} + Ks] \;\;\;\;\;\; (5)
Note: For simultaneous prediction intervals, only lower
(pi.type="lower"
) and upper
(pi.type="upper"
) prediction
intervals are available.
The derivation of the constant K
is explained in the help file for
predIntNormSimultaneousK
.
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 mean for the future observations is allowed to differ from
the population mean for the observations used to construct the prediction interval.
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
.
Power Based on the kofm Rule (rule="k.of.m"
)
For the k
ofm
rule (rule="k.of.m"
) with w=1
(i.e., n.mean=1
), 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 observations 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: When k=m
and r=1
, this kind of simultaneous prediction
interval becomes the same as a standard prediction interval for the next
k
observations (see predIntNorm
).
Davis and McNichols (1987) show that for a onesided upper prediction interval
(pi.type="upper"
), the probability p
that at least k
of the
next m
future observations will be contained in the interval given in
Equation (5) above, for each of r
future sampling occasions, is given by:
p = \int_0^1 T(\sqrt{n}K; n1, \sqrt{n}[\Phi^{1}(v) + \frac{\Delta}{\sigma}]) r[I(v; k, m+1k)]^{r1} [\frac{v^{k1}(1v)^{mk}}{B(k, m+1k)}] dv \;\;\;\;\;\; (6)
where T(x; \nu, \delta)
denotes the cdf of the
noncentral Student's tdistribution with parameters
df=
\nu
and ncp=
\delta
evaluated at x
;
\Phi(x)
denotes the cdf of the standard normal distribution
evaluated at x
; I(x; \nu, \omega)
denotes the cdf of the
beta distribution with parameters shape1=
\nu
and
shape2=
\omega
; and B(\nu, \omega)
denotes the value of the
beta function with parameters a=
\nu
and
b=
\omega
.
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. Usually you assume \Delta=0
unless you are interested
in computing the power of the rule to detect a change in means between the
populations, as we are here.
If we are interested in using averages instead of single observations, with
w \ge 1
(i.e., n.mean
\ge 1
), the first
term in the integral in Equation (6) that involves the cdf of the
noncentral Student's tdistribution becomes:
T(\sqrt{n}K; n1, \frac{\sqrt{n}}{\sqrt{w}}[\Phi^{1}(v) + \frac{\sqrt{w}\Delta}{\sigma}]) \;\;\;\;\;\; (7)
For a given confidence level (1\alpha)100\%
, the power of the rule to detect
a change in means is simply given by:
Power = 1  p \;\;\;\;\;\; (8)
where p
is defined in Equation (6) above using the value of K
that
corresponds to \Delta/\sigma = 0
. Thus, when the argument
delta.over.sigma=0
, the value of p
is 1\alpha
and the power is
simply \alpha 100\%
. As delta.over.sigma
increases above 0, the power
increases.
When pi.type="lower"
, the same value of K
is used as when
pi.type="upper"
, but Equation (4) is used to construct the prediction
interval. Thus, the power increases as delta.over.sigma
decreases below 0.
Power Based on the California Rule (rule="CA"
)
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 m1
observations will
fall in the prediction interval. That is, if the first observation falls in the
prediction interval then sampling can stop. Otherwise, m1
more observations
must be taken.
The derivation of the power is the same as for the k
ofm
rule, except
that Equation (6) becomes the following (Davis, 1998b):
p = \int_0^1 T(\sqrt{n}K; n1, \sqrt{n}[\Phi^{1}(v) + \frac{\Delta}{\sigma}]) r\{v[1+v^{m2}(1v)]\}^{r1} [1+v^{m2}(m1mv)] dv \;\;\;\;\;\; (9)
Power Based on the Modified California Rule (rule="Modified.CA"
)
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.
The derivation of the power is the same as for the k
ofm
rule, except
that Equation (6) becomes the following (Davis, 1998b):
p = \int_0^1 T(\sqrt{n}K; n1, \sqrt{n}[\Phi^{1}(v) + \frac{\Delta}{\sigma}]) r\{v[1+v(3v[52v])]\}^{r1} \{1+v[6v(158v)]\} dv \;\;\;\;\;\; (10)
vector of values between 0 and 1 equal to the probability that the rule will be violated.
See the help file for predIntNormSimultaneous
.
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
predIntNormSimultaneousTestPower
and
plotPredIntNormSimultaneousTestPowerCurve
can be
used to investigate these relationships for the case of normallydistributed
observations.
Steven P. Millard (EnvStats@ProbStatInfo.com)
See the help file for predIntNormSimultaneous
.
predIntNormSimultaneous
, predIntNormSimultaneousK
,
plotPredIntNormSimultaneousTestPowerCurve
,
predIntNorm
, predIntNormK
,
predIntNormTestPower
, Prediction Intervals,
Normal.
# For the kofm rule with n=4, k=1, m=3, and r=1, show how the power increases
# as delta.over.sigma increases. Assume a 95% upper prediction interval.
predIntNormSimultaneousTestPower(n = 4, m = 3, delta.over.sigma = 0:2)
#[1] 0.0500000 0.2954156 0.7008558
#
# Look at how the power increases with sample size for an upper onesided
# prediction interval using the kofm rule with k=1, m=3, r=20,
# delta.over.sigma=2, and a confidence level of 95%.
predIntNormSimultaneousTestPower(n = c(4, 8), m = 3, r = 20, delta.over.sigma = 2)
#[1] 0.6075972 0.9240924
#
# Compare the power for the 1of3 rule with the power for the California and
# Modified California rules, based on a 95% upper prediction interval and
# delta.over.sigma=2. Assume a sample size of n=8. Note that in this case the
# power for the Modified California rule is greater than the power for the
# 1of3 rule and California rule.
predIntNormSimultaneousTestPower(n = 8, k = 1, m = 3, delta.over.sigma = 2)
#[1] 0.788171
predIntNormSimultaneousTestPower(n = 8, m = 3, rule = "CA", delta.over.sigma = 2)
#[1] 0.7160434
predIntNormSimultaneousTestPower(n = 8, rule = "Modified.CA", delta.over.sigma = 2)
#[1] 0.8143687
#
# Show how the power for an upper 95% simultaneous prediction limit increases
# as the number of future sampling occasions r increases. Here, we'll use the
# 1of3 rule with n=8 and delta.over.sigma=1.
predIntNormSimultaneousTestPower(n = 8, k = 1, m = 3, r=c(1, 2, 5, 10),
delta.over.sigma = 1)
#[1] 0.3492512 0.4032111 0.4503603 0.4633773
#==========
# USEPA (2009) contains an example on page 1923 that involves monitoring
# nw=100 compliance wells at a large facility with minimal natural spatial
# variation every 6 months for nc=20 separate chemicals.
# There are n=25 background measurements for each chemical to use to create
# simultaneous prediction intervals. We would like to determine which kind of
# resampling plan based on normal distribution simultaneous prediction intervals to
# use (1ofm, 1ofm based on means, or Modified California) in order to have
# adequate power of detecting an increase in chemical concentration at any of the
# 100 wells while at the same time maintaining a sitewide false positive rate
# (SWFPR) of 10% per year over all 4,000 comparisons
# (100 wells x 20 chemicals x semiannual sampling).
# The function predIntNormSimultaneousTestPower includes the argument "r"
# that is the number of future sampling occasions (r=2 in this case because
# we are performing semiannual sampling), so to compute the individual test
# Type I error level alpha.test (and thus the individual test confidence level),
# we only need to worry about the number of wells (100) and the number of
# constituents (20): alpha.test = 1(1alpha)^(1/(nw x nc)). The individual
# confidence level is simply 1alpha.test. Plugging in 0.1 for alpha,
# 100 for nw, and 20 for nc yields an individual test confidence level of
# 1alpha.test = 0.9999473.
nc < 20
nw < 100
conf.level < (1  0.1)^(1 / (nc * nw))
conf.level
#[1] 0.9999473
# Now we can compute the power of any particular sampling strategy using
# predIntNormSimultaneousTestPower. For example, here is the power of
# detecting an increase of three standard deviations in concentration using
# the prediction interval based on the "1of2" resampling rule:
predIntNormSimultaneousTestPower(n = 25, k = 1, m = 2, r = 2, rule = "k.of.m",
delta.over.sigma = 3, pi.type = "upper", conf.level = conf.level)
#[1] 0.3900202
# The following commands will reproduce the table shown in Step 2 on page
# 1923 of USEPA (2009). Because these commands can take more than a few
# seconds to execute, we have commented them out here. To run this example,
# just remove the pound signs (#) that are in front of R commands.
#rule.vec < c(rep("k.of.m", 3), "Modified.CA", rep("k.of.m", 3))
#m.vec < c(2, 3, 4, 4, 1, 2, 1)
#n.mean.vec < c(rep(1, 4), 2, 2, 3)
#n.scenarios < length(rule.vec)
#K.vec < numeric(n.scenarios)
#Power.vec < numeric(n.scenarios)
#K.vec < predIntNormSimultaneousK(n = 25, k = 1, m = m.vec, n.mean = n.mean.vec,
# r = 2, rule = rule.vec, pi.type = "upper", conf.level = conf.level)
#Power.vec < predIntNormSimultaneousTestPower(n = 25, k = 1, m = m.vec,
# n.mean = n.mean.vec, r = 2, rule = rule.vec, delta.over.sigma = 3,
# pi.type = "upper", conf.level = conf.level)
#Power.df < data.frame(Rule = rule.vec, k = rep(1, n.scenarios), m = m.vec,
# N.Mean = n.mean.vec, K = round(K.vec, 2), Power = round(Power.vec, 2),
# Total.Samples = m.vec * n.mean.vec)
#Power.df
# Rule k m N.Mean K Power Total.Samples
#1 k.of.m 1 2 1 3.16 0.39 2
#2 k.of.m 1 3 1 2.33 0.65 3
#3 k.of.m 1 4 1 1.83 0.81 4
#4 Modified.CA 1 4 1 2.57 0.71 4
#5 k.of.m 1 1 2 3.62 0.41 2
#6 k.of.m 1 2 2 2.33 0.85 4
#7 k.of.m 1 1 3 2.99 0.71 3
# The above table shows the Kmultipliers for each prediction interval, along with
# the power of detecting a change in concentration of three standard deviations at
# any of the 100 wells during the course of a year, for each of the sampling
# strategies considered. The last three rows of the table correspond to sampling
# strategies that involve using the mean of two or three observations.
#==========
# Clean up
#
rm(nc, nw, conf.level, rule.vec, m.vec, n.mean.vec, n.scenarios, K.vec,
Power.vec, Power.df)