Example function for calculating Working-Hotelling and Bonferroni confidence intervals at a 95% confidence level

working_hotelling_bonferroni_intervals.R

# Example function for calculating Working-Hotelling and Bonferroni confidence intervals at a 95% level.
# Function takes two arguments:
#   x: predictor variable. Must be a vector.
#   y: response variable. Must be a vector.
# Plots are built using ggplot2
# gridExtra package, https://cran.r-project.org/web/packages/gridExtra/index.html, is used to plot multiple ggplots
# Function used to demonstrate how to construct simultaneous confidence intervals in post here: http://wp.me/p4aZEo-5Mg


working.hotelling.bonferroni.intervals <- function(x, y) {
  require(ggplot2)
  require(gridExtra)
  
  y <- as.matrix(y)
  if (ncol(y) > 1) {
    stop('y must be a vector')
  }
  x <- as.matrix(x)
  if (ncol(x) > 1) {
    stop('x must be a vector')
  }
  n <- length(y)
  if (n != length(x)) {
    stop('x and y must be the same length')
  }
  
  # Get the fitted values of the linear model
  fit <- lm(y ~ x)
  fit <- fit$fitted.values
  
  # Find standard error as defined above
  se <- sqrt(sum((y - fit)^2) / (n - 2)) * 
    sqrt(1 / n + (x - mean(x))^2 / 
           sum((x - mean(x))^2))
  
  # Calculate B and W statistics for both procedures
  W <- sqrt(2 * qf(p = 0.95, df1 = 2, df2 = n - 2))
  B <- 1-qt(.95/(2 * 3), n - 1)
  
  # Compute the simulatenous confidence intervals
  
  # Working-Hotelling
  wh.upper <- fit + W * se
  wh.lower <- fit - W * se
  
  # Bonferroni
  bon.upper <- fit + B * se
  bon.lower <- fit - B * se
  
  xy <- data.frame(cbind(x,y))
  
  # Plot the Working-Hotelling intervals
  wh <- ggplot(xy, aes(x=x, y=y)) + 
    geom_point(size=2.5) + 
    geom_line(aes(y=fit, x=x), size=1) + 
    geom_line(aes(x=x, y=wh.upper), colour='blue', linetype='dashed', size=1) + 
    geom_line(aes(x=x, wh.lower), colour='blue', linetype='dashed', size=1) +
    labs(title='Working-Hotelling')
  
  # Plot the Bonferroni intervals
  bonn <- ggplot(xy, aes(x=x, y=y)) + 
    geom_point(size=2.5) + 
    geom_line(aes(y=fit, x=x), size=1) + 
    geom_line(aes(x=x, y=bon.upper), colour='blue', linetype='dashed', size=1) + 
    geom_line(aes(x=x, bon.lower), colour='blue', linetype='dashed', size=1) +
    labs(title='Bonferroni')
  
  grid.arrange(wh, bonn, ncol = 2)
  
  # Collect results of procedures into a data.frame and return
  res <- data.frame(round(cbind(W, B), 3), row.names = c('Result'))
  colnames(res) <- c('W', 'B')
  
  return(res)
}