Using the generalized pivotal quantities
In AOV1R: Inference in the Balanced One-Way ANOVA Model with Random Factor

knitr::opts_chunk$set(
  collapse = FALSE,
  comment = "#>"
)

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

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)

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