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.6105023 10.649480
#> GPQ_sigma2b   0.5249176  2.477168
#> GPQ_sigma2w   0.7128439  1.328449
#> GPQ_sigma2tot 1.4472506  3.473547

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.194718 13.037688

One-sided generalized tolerance intervals

To get the bound of a one-sided generalized tolerance interval with tolerance level p and confidence level 1 − α, generate the simulations of the generalized pivotal quantity associated to the 100p%-quantile of the distribution of the response, then take the 100α%-quantile of these simulations for the right-sided tolerance interval and the 100(1 − α)%-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.492507 12.769368