Written by: Paul Rubin

Primary Source: OR in an OB World

I just wrapped up (knock on wood!) a coding project using R and Shiny. (Shiny, while *way* cool, is incidental to this post.) It was a favor for a friend, something she intends to use teaching an online course. Two of the tasks, while fairly mundane, generated code that was just barely obscure enough to be possibly worth sharing. It’s straight R code, so you need not use (or have installed) Shiny to use it.

The first task was to output, in tabular form, the coefficients of a linear regression model, along with their respective confidence intervals. The second task was to output, again in tabular form, the fitted values, confidence intervals and prediction intervals for the same model. Here is the function I wrote to do the first task (with Roxygen comments):

#' #' Summarize a fitted linear model, displaying both coefficient significance #' and confidence intervals. #' #' @param model an instance of class lm #' @param level the confidence level (default 0.95) #' #' @return a matrix combining the coefficient summary and confidence intervals #' model.ctable <- function(model, level = 0.95) { cbind(summary(model)$coefficients, confint(model, level = level)) }

To demonstrate its operation, I’ll generate a small sample with random data and run a linear regression on it.

I’ll generate the coefficient table using confidence level 0.9, rather than the default 0.95, for the coefficients.

model.ctable(m, level = 0.9)

The output is as follows:

Estimate Std. Error t value Pr(>|t|) 5 % 95 % (Intercept) 6.039951 0.2285568 26.42648 3.022477e-15 5.642352 6.437550 x 3.615331 0.2532292 14.27691 6.763279e-11 3.174812 4.055850 y -5.442428 0.3072587 -17.71285 2.156161e-12 -5.976937 -4.907918

The code for the second table (fits, confidence intervals and prediction intervals) is a bit longer:

#' #' Compute a table of fitted values, confidence intervals and #' prediction intervals from a regression model. #' #' @param model a fitted regression model #' @param level the desired confidence level (default 0.95) #' @param names the names to assign to the columns (after #' resequencing if necessry) #' @param order the order in which to list the columns #' (1 = fitted, 2 = lower c.i. limit, 3 = upper c.i. limit, #' 4 = lower p.i. limit, 5 = upper p.i. limit) #' #' @return a matrix with one row per observation and five #' columns (fitted value, lower/upper c.i. bounds, lower/upper #' p.i. bounds) in the order specified by the user #' intervals <- function(model, level = 0.95, names = c("Fitted", "CI Low", "CI High", "PI Low", "PI High"), order = c(4, 2, 1, 3, 5)) { # generate fits and confidence intervals temp <- predict(model, interval = "confidence", level = level) # generate fits and prediciton intervals (suppressing # the warning about predicting past values) temp2 <- suppressWarnings( predict(model, interval = "prediction", level = level) ) # drop the redundant fit column temp2 <- temp2[,2:3] # merge the tables and reorder the columns temp <- cbind(temp, temp2)[, order] # rename the columns colnames(temp) <- names[order] temp }

Here is the call with default arguments (using head() to limit the amount of output):

The output is this:

PI Low CI Low Fitted CI High PI High 1 -0.7928115 0.65769280 1.5196870 2.381681 3.832185 2 7.9056270 9.40123642 10.1928094 10.984382 12.479992 3 4.9125024 6.61897662 7.1149000 7.610823 9.317298 4 7.3386447 8.66123993 9.7406923 10.820145 12.142740 5 -1.4295587 0.05464529 0.8637503 1.672855 3.157059 6 4.1962493 5.84893725 6.4156619 6.982387 8.635074

Finally, I’ll run it again, changing the confidence level to 0.9, tweaking the column headings a bit, and reordering them:

The output is:

Fit CI_l CI_u PI_l PI_u 1 1.5196870 0.8089467 2.230427 -0.3870379 3.426412 2 10.1928094 9.5401335 10.845485 8.3069584 12.078660 3 7.1149000 6.7059962 7.523804 5.2989566 8.930843 4 9.7406923 8.8506512 10.630733 7.7601314 11.721253 5 0.8637503 0.1966188 1.530882 -1.0271523 2.754653 6 6.4156619 5.9483803 6.882943 4.5856891 8.245635

By the way, all syntax highlighting was Created by Pretty R at inside-R.org.

#### Paul Rubin

#### Latest posts by Paul Rubin (see all)

- CPLEX Callbacks: ThreadLocal v. Clone - April 19, 2019
- Pseudocode in LyX Revisited - March 13, 2019
- Guessing Pareto Solutions: A Test - February 26, 2019