Using the generalized pivotal quantities

The generalized pivotal quantities were introduced by Weerahandi. These are random variables, which are simulated by the function rGPQ.

Statistical inference based on the generalized pivotal quantities is similar to Bayesian posterior inference. For example, a generalized confidence interval of a parameter is obtained by taking the quantiles of the generalized pivotal quantity associated to this parameter.

Generalized confidence interval

Below is an example. We derive generalized confidence intervals for the three parameters defining the ANOVA model as well as for the total variance.

library(AOV1R)
dat <- simAOV1R(I = 20, J = 5, mu = 10, sigmab = 1, sigmaw = 1)
fit <- aov1r(y ~ group, data = dat)
nsims <- 50000L
gpq <- rGPQ(fit, nsims)
gpq[["GPQ_sigma2tot"]] <- with(gpq, GPQ_sigma2b + GPQ_sigma2w)
# Generalized confidence intervals:
t(vapply(gpq, quantile, numeric(2L), probs = c(2.5, 97.5)/100))

#>                    2.5%     97.5%
#> GPQ_mu        9.8270838 10.674387
#> GPQ_sigma2b   0.2433838  1.526248
#> GPQ_sigma2w   0.8020371  1.503742
#> GPQ_sigma2tot 1.2547967  2.658914

Generalized prediction interval

Here we generate simulations of the generalized predictive distribution:

ypred <- with(gpq, rnorm(nsims, GPQ_mu, sqrt(GPQ_sigma2tot)))

And then we get the generalized prediction interval by taking quantiles:

quantile(ypred, probs = c(2.5, 97.5)/100)

#>      2.5%     97.5% 
#>  7.586075 12.972425

One-sided generalized tolerance intervals

To get the bound of a one-sided generalized tolerance interval with tolerance level \(p\) and confidence level \(1-\alpha\), generate the simulations of the generalized pivotal quantity associated to the \(100p\%\)-quantile of the distribution of the response, then take the \(100\alpha\%\)-quantile of these simulations for the right-sided tolerance interval and the \(100(1-\alpha)\%\)-quantile for the left-sided tolerance interval:

p <- 90/100
alpha <- 2.5/100
z <- qnorm(p)
GPQ_lowerQuantile <- with(gpq, GPQ_mu - z*sqrt(GPQ_sigma2tot))
GPQ_upperQuantile <- with(gpq, GPQ_mu + z*sqrt(GPQ_sigma2tot))
c(
  quantile(GPQ_lowerQuantile, probs = alpha),
  quantile(GPQ_upperQuantile, probs = 1-alpha)
)

#>      2.5%     97.5% 
#>  7.943014 12.557365