Probability That at Least One Set of Future Observations Violates the Given Rule Based on a Nonparametric Simultaneous Prediction Interval

Source: R/predIntNparSimultaneousTestPower.R

Compute 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")

Details

Extended methods and derivations are available in Prediction intervals for environmental monitoring: predIntNparSimultaneousTestPower. Open the vignette locally with vignette("prediction-intervals", package = "EnvStats").

Arguments

n

vector of positive integers specifying the sample sizes. Missing (NA), undefined (NaN), and infinite (Inf, -Inf) values are not allowed.

n.median

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.

k

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".

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.

r

vector of positive integers specifying the number of future sampling “occasions”. The default value is r=1.

rule

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).

lpl.rank

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.

n.plus.one.minus.upl.rank

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.

delta.over.sigma

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.

pi.type

character string indicating what kind of prediction interval to compute. The possible values are "two.sided" (the default), "lower", and "upper".

r.shifted

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.

method

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).

NMC

positive integer indicating the number of Monte Carlo trials to run when
method="simulate". The default value is NMC=100.

ci

logical scalar indicating whether to compute a confidence interval for the power when method="simulate". The default value is ci=FALSE.

ci.conf.level

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".

integrate.args.list

list of arguments to supply to the integrate function. The default value is NULL.

evNormOrdStats.method

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.

Value

vector of values between 0 and 1 equal to the probability that the rule will be violated.

References

See the help file for predIntNparSimultaneous.

Gansecki, M. (2009). Using the Optimal Rank Values Calculator. US Environmental Protection Agency, Region 8, March 10, 2009.

Author

Steven P. Millard (EnvStats@ProbStatInfo.com)

Note

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.

Examples

  # 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)