Simultaneous Inference"

In api2lm: Functions and Data Sets for the Book "A Progressive Introduction to Linear Models"

The api2lm package makes it easier to perform simultaneous inference for confidence intervals of regression coefficients and the response mean, as well as the prediction intervals for a new response.

Let's fit a basic basic linear model using the mtcars available in the datasets package. We consider the model fit in the confint.lm function.

fit <- lm(100/mpg ~ disp + hp + wt + am, data = mtcars)

Confidence for regression coefficients

We construct typical $t$-based confidence intervals using the confint function, as shown below. We use the default individual confidence level of 0.95.

confint(fit)

The family-wise confidence level is guaranteed to be at least $1-5(0.05)\geq 0.75$ based on Boole's inequality. We can use the Bonferroni correction or the Working-Hotelling procedure (the latter applies to all linear combinations of the regression coefficients) to control the family-wise confidence level. These adjusted intervals are available using the confint_adjust function in the api2lm package. By default, the function makes no adjustments, but indicates the lower bound of the family-wise confidence level.

library(api2lm)
confint_adjust(fit)

We get the same intervals as before, but the print method for the object returned by confint_adjust provides the family-wise confidence level.

To use a Bonferroni adjustment, we set method = "bonferroni" for confint_adjust, as shown below.

(ci_b <- confint_adjust(fit, method = "bonferroni"))

Naturally, these intervals are wider than the unadjusted intervals but control the family-wise confidence level at 0.95.

To construct the Working-Hotelling-adjusted intervals, adjustment, we set method = "wh" for confint_adjust, as shown below.

confint_adjust(fit, method = "wh")

The Working-Hotelling adjusted intervals are even wider than the Bonferroni-adjusted intervals, so the Bonferroni correction is preferred here.

knitr::opts_chunk$set(echo = TRUE)

We can easily plot our confidence intervals using plot or autoplot (if the ggplot2 package is available. We plot the Bonferroni-adjusted intervals stored in ci_b using the code below.

plot(ci_b)

The intervals for hp and disp are difficult to see using because of their scale relative to the other intervals, so we use the parm function to look at them specifically in the code below.

plot(ci_b, parm = c("hp", "disp"))

The intervals are now easier to see. We can use the mar argument to reduce the margin along the y-axis (which is modified internally by default so that all interval labels are shown).

plot(ci_b, parm = c("hp", "disp"), mar = c(4.1, 4.1, 2.1, 2.1))

Alternatively, we can use the autoplot function, which makes these adjustments automatically.

library(ggplot2)
autoplot(ci_b, parm = c("hp", "disp"))

Confidence intervals for response mean

Interval procedures, whether for the response mean or new observations, suffer from the same type of multiple comparisons problems that intervals for regression coefficients have.

The predict_adjust function can be used to create adjusted intervals that control the family-wise confidence level for the mean response. The function can be used to produce unadjusted intervals (method = "none"), Bonferroni-adjusted intervals (method = "bonferroni"), and Working-Hotelling-adjusted intervals (method = "wh"). The Working-Hotelling intervals will tend to be narrower the more intervals considered because the Working-Hotelling procedure is valid for ALL linear combinations of the regression coefficients and not only the ones being produced. We produce unadjusted, Bonferroni-adjusted, and Working-Hotelling-adjusted intervals for two combinations of predictors in the code below.

# observations for which to predict the mean response
newdata <- as.data.frame(rbind(
               apply(mtcars, 2, mean),
               apply(mtcars, 2, median)))
# unadjusted intervals
predict_adjust(fit, newdata = newdata,
               interval = "confidence",
               method = "none")
# bonferroni-adjusted intervals
predict_adjust(fit, newdata = newdata,
               interval = "confidence",
               method = "bonferroni")
# working-hotelling-adjusted intervals
predict_adjust(fit, newdata = newdata,
               interval = "confidence",
               method = "wh")

Prediction intervals for a new response

The predict_adjust function can be used to create adjusted intervals that control the family-wise confidence level for new responses. The function can be used to produce unadjusted intervals (method = "none"), Bonferroni-adjusted intervals (method = "bonferroni"), and Scheffe-adjusted intervals (method = "scheffe").We produce unadjusted, Bonferroni-adjusted, and Scheffe-adjusted predictions intervals for four combinations of predictors in the code below.

# observations for which to predict the mean response
newdata <- as.data.frame(rbind(
               apply(mtcars, 2, mean),
               apply(mtcars, 2, median),
               apply(mtcars, 2, quantile, prob = 0.25),
               apply(mtcars, 2, quantile, prob = 0.75)))
# unadjusted intervals
predict_adjust(fit, newdata = newdata,
               interval = "prediction",
               method = "none")
# bonferroni-adjusted intervals
predict_adjust(fit, newdata = newdata,
               interval = "prediction",
               method = "bonferroni")
# scheffe-adjusted intervals
predict_adjust(fit, newdata = newdata,
               interval = "prediction",
               method = "scheffe")

In this case, the Bonferroni and Scheffe-adjusted intervals produce the same results.

api2lm documentation built on July 9, 2023, 5:52 p.m.