# plotPredIntLnormAltSimultaneousTestPowerCurve: Power Curves for Sampling Design for Test Based on... In EnvStats: Package for Environmental Statistics, Including US EPA Guidance

## Description

Plot power vs. θ_1/θ_2 (ratio of means) for a sampling design for a test based on a simultaneous prediction interval for a lognormal distribution.

## Usage

 ```1 2 3 4 5 6 7``` ``` plotPredIntLnormAltSimultaneousTestPowerCurve(n = 8, df = n - 1, n.geomean = 1, k = 1, m = 2, r = 1, rule = "k.of.m", cv = 1, range.ratio.of.means = c(1, 5), pi.type = "upper", conf.level = 0.95, r.shifted = r, K.tol = .Machine\$double.eps^(1/2), integrate.args.list = NULL, plot.it = TRUE, add = FALSE, n.points = 20, plot.col = "black", plot.lwd = 3 * par("cex"), plot.lty = 1, digits = .Options\$digits, cex.main = par("cex"), ..., main = NULL, xlab = NULL, ylab = NULL, type = "l") ```

## Arguments

 `n` positive integer greater than 2 indicating the sample size upon which the prediction interval is based. The default is value is `n=8`. `df` positive integer indicating the degrees of freedom associated with the sample size. The default value is `df=n-1`. `n.geomean` positive integer specifying the sample size associated with the future geometric mean(s). The default value is `n.geomean=1` (i.e., individual observations). Note that all future geometric means must be based on the same sample size. `k` for the k-of-m rule (`rule="k.of.m"`), positive integer specifying the minimum number of observations (or averages) out of m observations (or averages) (all obtained on one future sampling “occassion”) the prediction interval should contain with confidence level `conf.level`. The default value is `k=1`. This argument is ignored when the argument `rule` is not equal to `"k.of.m"`. `m` positive integer specifying the maximum number of future observations (or averages) 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` positive integer 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). `cv` positive value specifying the coefficient of variation for both the population that was sampled to construct the prediction interval and the population that will be sampled to produce the future observations. The default value is `cv=1`. `range.ratio.of.means` numeric vector of length 2 indicating the range of the x-variable to use for the plot. The default value is `range.ratio.of.means=c(1,5)`. `pi.type` character string indicating what kind of prediction interval to compute. The possible values are `pi.type="upper"` (the default), and `pi.type="lower"`. `conf.level` numeric scalar between 0 and 1 indicating the confidence level of the prediction interval. The default value is `conf.level=0.95`. `r.shifted` positive integer between `1` and coder specifying the number of future sampling occasions for which the mean is shifted. The default value is `r.shifted=r`. `K.tol` numeric scalar indicating the tolerance to use in the nonlinear search algorithm to compute K. The default value is `K.tol=.Machine\$double.eps^(1/2)`. For many applications, the value of K needs to be known only to the second decimal place, in which case setting `K.tol=1e-4` will speed up computation a bit. `integrate.args.list` 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. `plot.it` a logical scalar indicating whether to create a plot or add to the existing plot (see explanation of the argument `add` below) on the current graphics device. If `plot.it=FALSE`, no plot is produced, but a list of (x,y) values is returned (see the section VALUE). The default value is `plot.it=TRUE`. `add` a logical scalar indicating whether to add the design plot to the existing plot (`add=TRUE`), or to create a plot from scratch (`add=FALSE`). The default value is `add=FALSE`. This argument is ignored if `plot.it=FALSE`. `n.points` a numeric scalar specifying how many (x,y) pairs to use to produce the plot. There are `n.points` x-values evenly spaced between `range.x.var[1]` and `range.x.var[2]`. The default value is `n.points=100`. `plot.col` a numeric scalar or character string determining the color of the plotted line or points. The default value is `plot.col="black"`. See the entry for `col` in the help file for `par` for more information. `plot.lwd` a numeric scalar determining the width of the plotted line. The default value is `3*par("cex")`. See the entry for `lwd` in the help file for `par` for more information. `plot.lty` a numeric scalar determining the line type of the plotted line. The default value is `plot.lty=1`. See the entry for `lty` in the help file for `par` for more information. `digits` a scalar indicating how many significant digits to print out on the plot. The default value is the current setting of `options("digits")`. `cex.main, main, xlab, ylab, type, ...` additional graphical parameters (see `par`).

## Details

See the help file for `predIntLnormAltSimultaneousTestPower` for information on how to compute the power of a hypothesis test for the difference between two means of lognormal distributions based on a simultaneous prediction interval for a lognormal distribution.

## Value

`plotPredIntLnormAltSimultaneousTestPowerCurve` invisibly returns a list with components:

 `x.var` x-coordinates of points that have been or would have been plotted. `y.var` y-coordinates of points that have been or would have been plotted.

## Note

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 `predIntLnormAltSimultaneousTestPower` and
`plotPredIntLnormAltSimultaneousTestPowerCurve` can be used to investigate these relationships for the case of normally-distributed observations.

## Author(s)

Steven P. Millard ([email protected])

## References

See the help file for `predIntNormSimultaneous`.

`predIntLnormAltSimultaneousTestPower`, `predIntLnormAltSimultaneous`, `predIntLnormAlt`, `predIntLnormAltTestPower`, Prediction Intervals, LognormalAlt.
 ``` 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127``` ``` # USEPA (2009) contains an example on page 19-23 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 (1-of-m, 1-of-m 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 site-wide false positive rate # (SWFPR) of 10% per year over all 4,000 comparisons # (100 wells x 20 chemicals x semi-annual 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 semi-annual 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-(1-alpha)^(1/(nw x nc)). The individual # confidence level is simply 1-alpha.test. Plugging in 0.1 for alpha, # 100 for nw, and 20 for nc yields an individual test confidence level of # 1-alpha.test = 0.9999473. nc <- 20 nw <- 100 conf.level <- (1 - 0.1)^(1 / (nc * nw)) conf.level #[1] 0.9999473 # The help file for predIntNormSimultaneousTestPower shows how to # create the results below for various sampling plans: # 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 K-multipliers 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. # Here we will create a variation of this example based on # using a lognormal distribution and plotting power versus ratio of the # means assuming cv=1. # Here is the power curve for the 1-of-4 sampling strategy: dev.new() plotPredIntLnormAltSimultaneousTestPowerCurve(n = 25, k = 1, m = 4, r = 2, rule="k.of.m", range.ratio.of.means = c(1, 10), pi.type = "upper", conf.level = conf.level, ylim = c(0, 1), main = "") title(main = paste("Power Curves for 1-of-4 Sampling Strategy Based on 25 Background", "Samples, SWFPR=10%, and 2 Future Sampling Periods", sep = "\n")) mtext("Assuming Lognormal Data with CV=1", line = 0) #---------- # Here are the power curves for the first four sampling strategies. # Because this takes several seconds to run, here we have commented out # the R commands. To run this example, just remove the pound signs (#) # from in front of the R commands. #dev.new() #plotPredIntLnormAltSimultaneousTestPowerCurve(n = 25, k = 1, m = 4, r = 2, # rule="k.of.m", range.ratio.of.means = c(1, 10), pi.type = "upper", # conf.level = conf.level, ylim = c(0, 1), main = "") #plotPredIntLnormAltSimultaneousTestPowerCurve(n = 25, k = 1, m = 3, r = 2, # rule="k.of.m", range.ratio.of.means = c(1, 10), pi.type = "upper", # conf.level = conf.level, add = TRUE, plot.col = "red", plot.lty = 2) #plotPredIntLnormAltSimultaneousTestPowerCurve(n = 25, k = 1, m = 2, r = 2, # rule="k.of.m", range.ratio.of.means = c(1, 10), pi.type = "upper", # conf.level = conf.level, add = TRUE, plot.col = "blue", plot.lty = 3) #plotPredIntLnormAltSimultaneousTestPowerCurve(n = 25, r = 2, rule="Modified.CA", # range.ratio.of.means = c(1, 10), pi.type = "upper", conf.level = conf.level, # add = TRUE, plot.col = "green3", plot.lty = 4) #legend("topleft", c("1-of-4", "Modified CA", "1-of-3", "1-of-2"), # col = c("black", "green3", "red", "blue"), lty = c(1, 4, 2, 3), # lwd = 3 * par("cex"), bty = "n") #title(main = paste("Power Curves for 4 Sampling Strategies Based on 25 Background", # "Samples, SWFPR=10%, and 2 Future Sampling Periods", sep = "\n")) #mtext("Assuming Lognormal Data with CV=1", line = 0) #---------- # Here are the power curves for the last 3 sampling strategies: # Because this takes several seconds to run, here we have commented out # the R commands. To run this example, just remove the pound signs (#) # from in front of the R commands. #dev.new() #plotPredIntLnormAltSimultaneousTestPowerCurve(n = 25, k = 1, m = 2, n.geomean = 2, # r = 2, rule="k.of.m", range.ratio.of.means = c(1, 10), pi.type = "upper", # conf.level = conf.level, ylim = c(0, 1), main = "") #plotPredIntLnormAltSimultaneousTestPowerCurve(n = 25, k = 1, m = 1, n.geomean = 2, # r = 2, rule="k.of.m", range.ratio.of.means = c(1, 10), pi.type = "upper", # conf.level = conf.level, add = TRUE, plot.col = "red", plot.lty = 2) #plotPredIntLnormAltSimultaneousTestPowerCurve(n = 25, k = 1, m = 1, n.geomean = 3, # r = 2, rule="k.of.m", range.ratio.of.means = c(1, 10), pi.type = "upper", # conf.level = conf.level, add = TRUE, plot.col = "blue", plot.lty = 3) #legend("topleft", c("1-of-2, Order 2", "1-of-1, Order 3", "1-of-1, Order 2"), # col = c("black", "blue", "red"), lty = c(1, 3, 2), lwd = 3 * par("cex"), # bty="n") #title(main = paste("Power Curves for 3 Sampling Strategies Based on 25 Background", # "Samples, SWFPR=10%, and 2 Future Sampling Periods", sep = "\n")) #mtext("Assuming Lognormal Data with CV=1", line = 0) #========== # Clean up #--------- rm(nc, nw, conf.level) graphics.off() ```