# NOT RUN {
dat <- simAOV1R(I=20, J=5, mu=10, sigmab=1, sigmaw=1)
fit <- aov1r(y ~ group, data=dat)
nsims <- 20000
pivsims <- rGPQ(fit, nsims)
pivsims$GPQ_sigma2tot <- pivsims$GPQ_sigma2b + pivsims$GPQ_sigma2w
# Generalized confidence intervals:
lapply(pivsims, quantile, probs = c(0.025, 0.975))
# compare with the frequentist confidence intervals:
confint(fit, SDs = FALSE)
# Generalized prediction interval:
with(
pivsims,
quantile(rnorm(nsims, GPQ_mu, sqrt(GPQ_sigma2tot)),
probs = c(0.025, 0.975))
)
# compare with the frequentist prediction interval:
predict(fit)
# }
Run the code above in your browser using DataLab