Two-sided simultaneous prediction intervals are invalid

[SteveMillard]

Hello Alex, Justin, and everyone else who is maintaining EnvStats, the package that I originally authored. I am so grateful for all of your work! I am pretty much retired now, but plan to periodically help with improving the code. When EnvStats was updated from Version 2.3.1 to Version 2.4.0, Justin added the option for two-sided confidence intervals in the functions predIntNormSimultaneous() and predIntNormSimultaneousK(). Justin assumed that the way to derive K for a two-sided simultaneous prediction interval is the same as deriving it for a simple prediction intervals, i.e., he assumed that Equation (1) in the help files for predIntNormSimultaneous() and predIntNormSimultaneousK() applied to simultaneous prediction intervals. The way that the text was worded in these help files could lead the reader to assume this, however, this is not correct, and there is no theoretical foundation for this assumption; the math is much more complicated to derive a simultaneous prediction interval compared to a simple prediction interval. The references cited in the help file for predIntNormSimultaneous() only discuss one-sided prediction intervals. This is why in the help files there is the statement:

“Note: For simultaneous prediction intervals, only lower (pi.type="lower") and upper (pi.type="upper") prediction intervals are available.”

(Note that Justin forgot to omit this statement in the help files when he updated the function.)

Below my signature is code for an example showing that the two-sided prediction intervals are not valid.

So Justin, it would great if you could update the code for all of the functions that were changed to allow pi.type = "two-sided". I had suggested to Alex that we revert the code back to how it was in version 2.3.1, but Alex had a better suggestion of keeping "two-sided" as a possible value for pi.type, and if that is selected, then the functions should return an error indicating that two-sided simultaneous prediction intervals are not available, so that if someone runs a script that they had used for versions 2.4.0 - 2.8.1 they will see the error. We need to fix code in the following functions:

• predIntNormSimultaneous o Update it to throw an error when pi.type = "two-sided"

• predIntNormSimultaneousK o Update it to throw an error when pi.type = "two-sided"
o It does not need to supply pi.type in the call to any of the following functions:  pred.int.norm.k.of.m.on.r.K  pred.int.norm.CA.on.r.K  pred.int.norm.Modified.CA.on.r.K

• Omit the argument pi.type for the following functions: o pred.int.norm.k.of.m.on.r.K o pred.int.norm.CA.on.r.K o pred.int.norm.Modified.CA.on.r.K

• predIntNormSimultaneousTestPower o Update it to throw an error when pi.type = "two-sided".

• Update the help files for predIntNormSimultaneous(), predIntNormSimultaneousK(), and predIntNormSimultaneousTestPower.

Thanks again to everyone for keeping EnvStats going!!

Sincerely, Steve Millard

=================================================================================

#####
# File:    Example Invalid Values 2-Sided Simultaneous Pred Int.R
#
# Purpose: This R script gives an example that shows that the 
#          implementation of a two-sided simultaneous prediction 
#          interval in the function predIntNormSimultaneous() is 
#          invalid.
#
# Author: Steve Millard
#
# Last Updated: 2024.01.24
#####

library(EnvStats)

#-------------------------------------------------------------------
# Generate 8 observations from a normal distribution with parameters
# mean=10 and sd=2, then use predIntNormSimultaneous to estimate the
# mean and standard deviation of the true distribution and construct
# an *upper* 90% simultaneous prediction interval to contain 
# at least 1 out of the next 3 observations.
#
# Then generate 3 new observations from a N(10, 2) distribution 
# and record how many are less than or equal to the 
# upper prediciton limit.
#
# Do this 1000 times and estimate the probability that the prediction 
# intervals contain at least 1 of 3 future observations.  The 
# underlying probability should be 90%, so the estimate should be 
# close to 0.9, and the 95% confidence interval for the probability 
# should contain 0.9.
#
# (The call to set.seed allows you to reproduce this example.)
#--------------------------------------------------------------------

set.seed(479)

N <- 1000

In_interval <- rep(as.logical(NA), N)

for (i in 1:N) {
     dat <- rnorm(8, mean = 10, sd = 2)
     UPL <- predIntNormSimultaneous(dat, k = 1, m = 3, 
          pi.type = "upper", conf.level = 0.9)$interval$limits["UPL"]
     new_dat <- rnorm(3, mean = 10, sd = 2)
     In_interval[i] <- sum(new_dat <= UPL) >= 1
}

mean(In_interval)
#[1] 0.909

ebinom(In_interval, ci = TRUE)
#Results of Distribution Parameter Estimation
#--------------------------------------------
#
#Assumed Distribution:            Binomial
#
#Estimated Parameter(s):          size = 1000.000
#                                 prob =    0.909
#
#Estimation Method:               mle/mme/mvue for 'prob'
#
#Data:                            In_interval
#
#Sample Size:                     1000
#
#Confidence Interval for:         prob
#
#Confidence Interval Method:      Score normal approximation
#                                 (With continuity correction)
#
#Confidence Interval Type:        two-sided
#
#Confidence Level:                95%
#
#Confidence Interval:             LCL = 0.8890328
#                                 UCL = 0.9257497

#===================================================================

#-------------------------------------------------------------------
# Now do the same thing, but use Justin Jent's two-sided 
# simultaneous prediction interval.  For this example, note that 
# the two-sided prediction intervals contain at least 1 of the 
# next 3 observations only about 38% of the time, not 90%.
#-------------------------------------------------------------------

In_interval <- rep(as.logical(NA), N)

for (i in 1:N) {
     dat <- rnorm(8, mean = 10, sd = 2)
     limits <- predIntNormSimultaneous(dat, k = 1, m = 3, 
          pi.type = "two-sided", conf.level = 0.9)$interval$limits
     new_dat <- rnorm(3, mean = 10, sd = 2)
     In_interval[i] <- sum(new_dat >= limits["LPL"] & 
          new_dat <= limits["UPL"]) >= 1
}

mean(In_interval)
#[1] 0.379

ebinom(In_interval, ci = TRUE)
#Results of Distribution Parameter Estimation
#--------------------------------------------
#
#Assumed Distribution:            Binomial
#
#Estimated Parameter(s):          size = 1000.000
#                                 prob =    0.379
#
#Estimation Method:               mle/mme/mvue for 'prob'
#
#Data:                            In_interval
#
#Sample Size:                     1000
#
#Confidence Interval for:         prob
#
#Confidence Interval Method:      Score normal approximation
#                                 (With continuity correction)
#
#Confidence Interval Type:        two-sided
#
#Confidence Level:                95%
#
#Confidence Interval:             LCL = 0.3489580
#                                 UCL = 0.4099834

rm(N, In_interval, i, dat, new_dat, UPL, limits)

[alexkowa]

@jentjr this is related to your pull request https://github.com/alexkowa/EnvStats/pull/1

[jentjr]

Thanks for pointing this out, and sorry for causing the issue. I can be assigned, or help out, with implementing Alex's suggestion as described.

[SteveMillard]

No worries, Justin. Alex and Justin: I am happy to update the code for the functions and associated help files (including making the help files clearer as to why two-sided simultaneous prediction intervals are not available). Just let me know.

[alexkowa]

Steve, yes please. Implement it.

[SteveMillard]

Thanks Alex, will do. It turns out there is literature describing how to create two-sided simultaneous prediction intervals for a Normal distribution, but at this point I'm just going to update the functions in EnvStats to only compute Upper or Lower intervals. It may take me a while to implement the algorithms in the publications listed below.

Robert E. Odeh (1989) Simultaneous Two-Sided Prediction Intervals to Contain at Least l Out of k Future Means from a Normal Distribution, Communications in Statistics - Simulation and Computation, 18:2, 429-457, DOI: 10.1080/03610918908812769

Yimeng Xie, Yili Hong, Luis A. Escobar & William Q. Meeker (2017) A general algorithm for computing simultaneous prediction intervals for the (log)-location-scale family of distributions, Journal of Statistical Computation and Simulation, 87:8, 1559-1576, DOI: 10.1080/00949655.2016.1277426

[SteveMillard]

@alexkowa and @jentjr , do either of you have a way to get the PDF of this article?

Robert E. Odeh (1989) Simultaneous Two-Sided Prediction Intervals to Contain at Least l Out of k Future Means from a Normal Distribution, Communications in Statistics - Simulation and Computation, 18:2, 429-457, DOI: 10.1080/03610918908812769

[jentjr]

Unfortunately, I have not found a way to access a PDF of the paper.

[SteveMillard]

No worries. Thanks for trying. I'll see if I can get a physical copy from a local library.

[SteveMillard]

@alexkowa and @jentjr : FYI, I will be getting a hard copy of this article within a week.

Robert E. Odeh (1989) Simultaneous Two-Sided Prediction Intervals to Contain at Least l Out of k Future Means from a Normal Distribution, Communications in Statistics - Simulation and Computation, 18:2, 429-457, DOI: 10.1080/03610918908812769

I'll look at it and see if I can easily use the described algorithm to include in predIntNormSimultaneous and predIntNormSimultaneousK. If it will take too long I will simply update the current functions to return an error with information as we previously discussed, and work on including two-sided simultaneous prediction intervals in another future version.

[alexkowa]

Great. Thank you.

[SteveMillard]

@alexkowa I created a fix for this issue. There are 4 pull requests under my fork, but I can't figure out how to link them to this issue.

[alexkowa]

You can simply write in the description of the pull request closes #X where X is in this case 28. See https://docs.github.com/en/issues/tracking-your-work-with-issues/linking-a-pull-request-to-an-issue

[SteveMillard]

Thanks!!Sincerely,—Steve M.On Apr 8, 2024, at 9:42 PM, Alexander Kowarik @.***> wrote: You can simply write in the description of the pull request closes #X where X is in this case 28. See https://docs.github.com/en/issues/tracking-your-work-with-issues/linking-a-pull-request-to-an-issue

—Reply to this email directly, view it on GitHub, or unsubscribe.You are receiving this because you were assigned.Message ID: @.***>