Linear Regression: interactive

PLS 206 — Applied Multivariate Modeling

Apps for Week 2, in the order the lecture covers them. They run in the browser, so no R installation is needed. Tick “Show R code” in any app to get code that reproduces the current view in RStudio.

---

Least Squares by Hand

Move the slope and watch the residuals and their squares grow or shrink. The right panel scores every slope by its residual sum of squares (RSS): the bottom of that bowl is the line lm() picks. Find the bottom yourself first, then tick “Show lm() line” or press “Jump to lm()” to check.

#| '!! shinylive warning !!': |
#|   shinylive does not work in self-contained HTML documents.
#|   Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 820

library(shiny)

# the 40 wheat plants from the lecture (data/wheat_plants.csv)
height <- c(86.9, 71.6, 73.1, 75.9, 91, 94.4, 71.5, 71.8, 77.8, 72.8, 88.7,
            86.6, 58.6, 101.8, 67.8, 76.8, 75.1, 60.3, 62.3, 86.5, 77.4,
            55.7, 78.9, 51.1, 83.6, 45.6, 105, 89.6, 64.7, 76, 70.7, 70.5,
            92.8, 62.4, 90.6, 78.5, 76.4, 89.1, 75.9, 66.1)
biomass <- c(18.1, 15.5, 14.8, 15.6, 21.6, 17.2, 16, 17.4, 18.2, 20.1, 19.9,
             12.8, 11.6, 24.4, 10.5, 16.6, 14.3, 15, 9.4, 13.9, 16.6, 7.5,
             23.3, 8.4, 23.5, 5.8, 24.5, 18.1, 11.5, 12.6, 12.6, 18.1, 13.8,
             9.8, 20.7, 12.9, 12.4, 23.1, 17, 10.7)

fit_lm  <- lm(biomass ~ height)
b_lm    <- coef(fit_lm)
rss_lm  <- sum(residuals(fit_lm)^2)

# RSS of a line through the mean point, for a given slope (as in the lecture)
model_line <- function(slope) {
  b0 <- mean(biomass) - slope * mean(height)
  predicted <- b0 + slope * height
  sum((biomass - predicted)^2)
}
slopes   <- seq(0, 0.6, by = 0.005)
rss_bowl <- sapply(slopes, model_line)

# RSS over intercept x slope, for the two-knob view
rss2 <- function(b0, b1) sum((biomass - (b0 + b1 * height))^2)
b0s  <- seq(-30, 30, length.out = 60)
b1s  <- seq(0, 0.6, length.out = 60)
bowl2 <- outer(b0s, b1s, Vectorize(rss2))

fmt <- function(x) format(x, digits = 10)

ui <- fluidPage(
  sidebarLayout(
    sidebarPanel(
      sliderInput("slope", "Slope (g per cm)", min = 0, max = 0.6, value = 0.12, step = 0.005),
      checkboxInput("through_mean", "Line goes through the mean point (x-bar, y-bar)", TRUE),
      conditionalPanel("input.through_mean == false",
        sliderInput("intercept", "Intercept", min = -30, max = 30, value = -5, step = 0.1)
      ),
      checkboxInput("show_resid", "Show residuals", TRUE),
      checkboxInput("show_sq", "Show squares", FALSE),
      checkboxInput("show_lm", "Show lm() line", FALSE),
      actionButton("jump", "Jump to lm()"),
      hr(),
      checkboxInput("show_code", "Show R code", FALSE),
      hr(),
      helpText("Line A in the lecture has slope 0.12, line B 0.279, line C 0.45."),
      helpText("Untick the mean-point box to move the intercept too: the bowl becomes a valley over two knobs.")
    ),
    mainPanel(
      plotOutput("lsPlot", height = "380px"),
      verbatimTextOutput("lsStats"),
      conditionalPanel("input.show_code == true",
        verbatimTextOutput("ls_code")
      )
    )
  )
)

server <- function(input, output, session) {

  observeEvent(input$jump, {
    updateSliderInput(session, "slope", value = unname(b_lm[2]))
    updateSliderInput(session, "intercept", value = unname(b_lm[1]))
  })

  line_now <- reactive({
    slope <- input$slope
    b0 <- if (isTRUE(input$through_mean)) mean(biomass) - slope * mean(height) else input$intercept
    predicted <- b0 + slope * height
    residual  <- biomass - predicted
    list(b0 = b0, slope = slope, predicted = predicted, residual = residual,
         rss = sum(residual^2))
  })

  output$lsPlot <- renderPlot({
    L <- line_now()
    par(mfrow = c(1, 2), mar = c(4, 4, 2, 1))

    # left: the data, the line, the residuals
    if (isTRUE(input$show_sq)) {
      plot(height, biomass, xlab = "Height (cm)", ylab = "Biomass (g)", asp = 1,
           pch = 16, col = "steelblue", main = sprintf("RSS = %.1f", L$rss))
      rect(height, pmin(biomass, L$predicted), height + abs(L$residual),
           pmax(biomass, L$predicted), col = adjustcolor("#0070C0", 0.2), border = "#0070C0")
    } else {
      plot(height, biomass, xlab = "Height (cm)", ylab = "Biomass (g)",
           pch = 16, col = "steelblue", main = sprintf("RSS = %.1f", L$rss))
    }
    if (isTRUE(input$show_resid))
      segments(height, biomass, height, L$predicted, col = "#0070C0", lty = 2)
    abline(L$b0, L$slope, col = "#7030A0", lwd = 2)
    if (isTRUE(input$show_lm)) abline(fit_lm, col = "#002855", lwd = 2, lty = 2)
    points(mean(height), mean(biomass), pch = 3, cex = 1.5)
    legend("topleft", bty = "n",
           legend = c("Your line", if (isTRUE(input$show_lm)) "lm() line"),
           col = c("#7030A0", if (isTRUE(input$show_lm)) "#002855"),
           lty = c(1, if (isTRUE(input$show_lm)) 2), lwd = 2)

    # right: the RSS bowl
    if (isTRUE(input$through_mean)) {
      plot(slopes, rss_bowl, type = "l", xlab = "Slope (g per cm)",
           ylab = "Residual sum of squares", main = "RSS for every slope")
      points(b_lm[2], rss_lm, pch = 16, col = "#002855", cex = 1.6)
      points(L$slope, L$rss, pch = 16, col = "#7030A0", cex = 1.6)
    } else {
      contour(b0s, b1s, bowl2, levels = round(rss_lm * c(1.1, 1.5, 2, 3, 5, 10, 20, 50, 100)),
              xlab = "Intercept", ylab = "Slope (g per cm)", main = "RSS over intercept and slope")
      points(b_lm[1], b_lm[2], pch = 16, col = "#002855", cex = 1.6)
      points(L$b0, L$slope, pch = 16, col = "#7030A0", cex = 1.6)
    }
    legend("top", bty = "n", pch = 16, col = c("#7030A0", "#002855"),
           legend = c("Your line", "lm() (the bottom)"))
  })

  output$lsStats <- renderPrint({
    L <- line_now()
    b1_hand <- cov(height, biomass) / var(height)
    b0_hand <- mean(biomass) - b1_hand * mean(height)
    cat(sprintf("Your line:  biomass = %.3f + %.3f * height\n", L$b0, L$slope))
    cat(sprintf("sum(residual):     %.3f\n", sum(L$residual)))
    cat(sprintf("RSS, your line:    %.2f\n", L$rss))
    cat(sprintf("RSS, lm() line:    %.2f\n", rss_lm))
    cat("\nBy hand:\n")
    cat(sprintf("  slope     = cov(height, biomass) / var(height) = %.4f\n", b1_hand))
    cat(sprintf("  intercept = mean(biomass) - slope * mean(height) = %.4f\n", b0_hand))
    cat("\ncoef(lm(biomass ~ height)):\n")
    print(round(b_lm, 4))
    if (abs(L$rss - rss_lm) < 1e-6) cat("\nYou are at the bottom of the bowl.\n")
    else if (abs(L$slope - b_lm[2]) < 0.003 && isTRUE(input$through_mean))
      cat("\nAs close as the slider steps allow: the slider rounds to 0.005.\n")
  })

  output$ls_code <- renderPrint({
    cat("# R code to reproduce this view\n\n")
    cat("wheat <- read.csv(\"https://greymonroe.github.io/PLS_206/data/wheat_plants.csv\")\n\n")
    cat(sprintf("slope <- %s\n", fmt(input$slope)))
    if (isTRUE(input$through_mean)) {
      cat("b0 <- mean(wheat$biomass) - slope * mean(wheat$height)   # line through the mean point\n")
    } else {
      cat(sprintf("b0 <- %s\n", fmt(input$intercept)))
    }
    cat("predicted <- b0 + slope * wheat$height\n")
    cat("residual  <- wheat$biomass - predicted\n")
    cat("sum(residual)\n")
    cat("sum(residual^2)    # RSS of this line\n\n")
    cat("par(mfrow = c(1, 2))\n")
    if (isTRUE(input$show_sq)) {
      cat("plot(wheat$height, wheat$biomass, xlab = \"Height (cm)\", ylab = \"Biomass (g)\",\n")
      cat("     pch = 16, col = \"steelblue\", asp = 1)\n")
      cat("rect(wheat$height, pmin(wheat$biomass, predicted), wheat$height + abs(residual),\n")
      cat("     pmax(wheat$biomass, predicted), col = adjustcolor(\"#0070C0\", 0.2), border = \"#0070C0\")\n")
    } else {
      cat("plot(wheat$height, wheat$biomass, xlab = \"Height (cm)\", ylab = \"Biomass (g)\",\n")
      cat("     pch = 16, col = \"steelblue\")\n")
    }
    if (isTRUE(input$show_resid))
      cat("segments(wheat$height, wheat$biomass, wheat$height, predicted, col = \"#0070C0\", lty = 2)\n")
    cat("abline(b0, slope, col = \"#7030A0\", lwd = 2)\n")
    if (isTRUE(input$show_lm))
      cat("abline(lm(biomass ~ height, data = wheat), col = \"#002855\", lwd = 2, lty = 2)\n")
    cat("\n")
    if (isTRUE(input$through_mean)) {
      cat("# the RSS for every slope (lines through the mean point) is a bowl\n")
      cat("model_line <- function(slope) {\n")
      cat("  b0 <- mean(wheat$biomass) - slope * mean(wheat$height)\n")
      cat("  predicted <- b0 + slope * wheat$height\n")
      cat("  sum((wheat$biomass - predicted)^2)\n")
      cat("}\n")
      cat("slopes <- seq(0, 0.6, by = 0.005)\n")
      cat("rss <- sapply(slopes, model_line)\n")
      cat("plot(slopes, rss, type = \"l\", xlab = \"Slope (g per cm)\", ylab = \"Residual sum of squares\")\n")
      cat("points(slope, model_line(slope), pch = 16, col = \"#7030A0\", cex = 1.6)\n")
      cat("slopes[which.min(rss)]\n\n")
    } else {
      cat("# both knobs at once: the RSS over intercept and slope\n")
      cat("rss2 <- function(b0, b1) sum((wheat$biomass - (b0 + b1 * wheat$height))^2)\n")
      cat("b0s <- seq(-30, 30, length.out = 60); b1s <- seq(0, 0.6, length.out = 60)\n")
      cat("bowl <- outer(b0s, b1s, Vectorize(rss2))\n")
      cat("rss_lm <- sum(residuals(lm(biomass ~ height, data = wheat))^2)\n")
      cat("contour(b0s, b1s, bowl, levels = round(rss_lm * c(1.1, 1.5, 2, 3, 5, 10, 20, 50, 100)),\n")
      cat("        xlab = \"Intercept\", ylab = \"Slope (g per cm)\")\n")
      cat("points(b0, slope, pch = 16, col = \"#7030A0\", cex = 1.6)\n\n")
    }
    cat("# the least-squares slope and intercept by hand\n")
    cat("b1 <- cov(wheat$height, wheat$biomass) / var(wheat$height)\n")
    cat("b1\n")
    cat("mean(wheat$biomass) - b1 * mean(wheat$height)\n\n")
    cat("# lm() finds the same line\n")
    cat("model <- lm(biomass ~ height, data = wheat)\n")
    cat("coef(model)\n")
    cat("sum(residuals(model)^2)    # RSS of the lm() line\n")
  })
}

shinyApp(ui, server)

---

Regression Explorer

Set the true intercept and slope, add noise, shift the predictor mean, and watch the estimated line track the true one. With the same sample, changing the true slope does not move the estimated intercept: the change adds (change in slope) * x to every y, and the fitted slope absorbs all of it. Shifting the mean of x away from 0 makes the intercept an extrapolation to x = 0, so its SE grows.

#| '!! shinylive warning !!': |
#|   shinylive does not work in self-contained HTML documents.
#|   Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 720

library(shiny)

ui <- fluidPage(
  sidebarLayout(
    sidebarPanel(
      sliderInput("n",         "Sample size (n)",       min = 5,   max = 300, value = 30,  step = 5),
      sliderInput("intercept", "True intercept (\u03b2\u2080)", min = -10, max = 10,  value = 0,   step = 0.5),
      sliderInput("slope",     "True slope (\u03b2\u2081)",     min = -2,  max = 2,   value = 0.5, step = 0.1),
      sliderInput("noise",     "Noise (SD)",             min = 0.1, max = 5,   value = 1,   step = 0.1),
      sliderInput("mean_x",    "Mean of x (predictor)", min = -10, max = 10,  value = 0,   step = 0.5),
      actionButton("resample", "New sample"),
      hr(),
      checkboxInput("show_code", "Show R code", FALSE),
      hr(),
      helpText("Change the true slope: the estimated intercept stays put, because the fitted slope takes up the whole change."),
      helpText("Move the mean of x away from 0: the intercept (the line at x = 0) is now far from the data, and its SE grows.")
    ),
    mainPanel(
      plotOutput("regPlot", height = "300px"),
      verbatimTextOutput("regStats"),
      conditionalPanel("input.show_code == true",
        verbatimTextOutput("reg_code")
      )
    )
  )
)

server <- function(input, output, session) {
  seed <- reactive(42 + input$resample)

  dat <- reactive({
    set.seed(seed())
    x <- rnorm(input$n, mean = input$mean_x)
    y <- input$intercept + input$slope * x + rnorm(input$n, sd = input$noise)
    data.frame(x = x, y = y)
  })

  output$regPlot <- renderPlot({
    d   <- dat()
    fit <- lm(y ~ x, data = d)

    par(mar = c(4, 4, 2, 1))
    plot(d$x, d$y,
         col  = adjustcolor("steelblue", 0.6),
         pch  = 16, cex = 1.2,
         xlab = "x", ylab = "y",
         main = paste0("n=", input$n,
                       "  \u03b2\u2080=", input$intercept,
                       "  \u03b2\u2081=", input$slope,
                       "  noise=", input$noise,
                       "  mean(x)=", input$mean_x))
    abline(fit, col = "#002855", lwd = 2)
    abline(a = input$intercept, b = input$slope,
           col = "orange", lwd = 1.5, lty = 2)
    abline(v = input$mean_x, col = "gray60", lty = 3, lwd = 1.2)
    legend("topleft",
           legend = c("Estimated", "True", "Mean of x"),
           col    = c("#002855", "orange", "gray60"),
           lty    = c(1, 2, 3), lwd = c(2, 1.5, 1.2), bty = "n")
  })

  output$regStats <- renderPrint({
    d     <- dat()
    fit   <- lm(y ~ x, data = d)
    s     <- summary(fit)
    coefs <- coef(s)

    cat("--- Intercept ---\n")
    cat(sprintf("  Estimated: %+.3f  (true: %+.1f)\n",
                coefs[1,1], input$intercept))
    cat(sprintf("  p-value:   %.4f  %s\n", coefs[1,4],
                ifelse(coefs[1,4] < 0.05, "[*]", "[ ]")))

    cat("\n--- Slope ---\n")
    cat(sprintf("  Estimated: %+.3f  (true: %+.1f)\n",
                coefs[2,1], input$slope))
    cat(sprintf("  p-value:   %.4f  %s\n", coefs[2,4],
                ifelse(coefs[2,4] < 0.05, "[*]", "[ ]")))

    cat(sprintf("\nR\u00b2: %.3f\n", s$r.squared))
  })

  output$reg_code <- renderPrint({
    cat("# Minimal R code to reproduce this view:\n\n")
    cat(sprintf("set.seed(%d)\n", seed()))
    cat(sprintf("x   <- rnorm(%d, mean = %.1f)\n", input$n, input$mean_x))
    cat(sprintf("y   <- %.1f + %.1f * x + rnorm(%d, sd = %.1f)\n",
                input$intercept, input$slope, input$n, input$noise))
    cat("fit <- lm(y ~ x)\n")
    cat("summary(fit)\n\n")
    cat("plot(x, y, pch = 16, col = \"steelblue\")\n")
    cat("abline(fit, col = \"#002855\", lwd = 2)         # estimated\n")
    cat(sprintf("abline(a = %.1f, b = %.1f, col = \"orange\", lty = 2)  # true\n",
                input$intercept, input$slope))
    cat(sprintf("abline(v = %.1f, col = \"gray60\", lty = 3)               # mean of x\n",
                input$mean_x))
  })
}

shinyApp(ui, server)

---

Confidence vs Prediction Intervals

Two questions about plants at one height: where is the average plant (confidence interval), and where will one new plant land (prediction interval)? Raise n and watch which band shrinks, then slide the height past the grey data range to see why we interpolate and don’t extrapolate.

#| '!! shinylive warning !!': |
#|   shinylive does not work in self-contained HTML documents.
#|   Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 900

library(shiny)

ui <- fluidPage(
  sidebarLayout(
    sidebarPanel(
      sliderInput("n", "Plants (n)", 5, 300, 40, 1),
      sliderInput("sigma", "Scatter around the line (sigma, g)", 0.5, 8, 3.1, 0.1),
      sliderInput("h", "Predict at height (cm)", 0, 200, 100, 1),
      selectInput("level", "Confidence level",
                  choices = c("0.90" = "0.9", "0.95" = "0.95", "0.99" = "0.99"),
                  selected = "0.95"),
      checkboxInput("show_true", "Show the true line", TRUE),
      actionButton("resample", "New sample"),
      hr(),
      checkboxInput("show_code", "Show R code", FALSE),
      helpText("Raise n to 300: the confidence band shrinks toward the line, the prediction band does not."),
      helpText("True model: biomass = -5.65 + 0.279 * height + noise, heights between 45 and 105 cm.")
    ),
    mainPanel(
      plotOutput("plot", height = "500px"),
      verbatimTextOutput("stats"),
      conditionalPanel("input.show_code == true",
        verbatimTextOutput("code_out")
      )
    )
  )
)

server <- function(input, output, session) {

  seed <- reactiveVal(206)
  observeEvent(input$resample, seed(seed() + 1), ignoreInit = TRUE)

  # the calculations, written as the code students will copy
  calc_code <- reactive({
    paste0(
'# Confidence vs prediction intervals: simulated wheat plants
set.seed(', seed(), ')
n <- ', input$n, '; sigma <- ', input$sigma, '; h <- ', input$h, '; L <- ', input$level, '
height  <- runif(n, 45, 105)
biomass <- -5.65 + 0.279 * height + rnorm(n, 0, sigma)
plants  <- data.frame(height, biomass)
model   <- lm(biomass ~ height, data = plants)

# two intervals at height h: the line (average plant) or one new plant
new  <- data.frame(height = h)
conf <- predict(model, new, interval = "confidence", level = L)
pred <- predict(model, new, interval = "prediction", level = L)
conf
pred
conf[, "upr"] - conf[, "lwr"]   # width of the confidence interval
pred[, "upr"] - pred[, "lwr"]   # width of the prediction interval
sigma(model)                    # residual standard error: the typical miss
range(plants$height)            # outside this range we are extrapolating

# check: simulate 1000 new plants at height h from the true model
new_biomass <- -5.65 + 0.279 * h + rnorm(1000, 0, sigma)
mean(new_biomass >= pred[, "lwr"] & new_biomass <= pred[, "upr"])   # share inside
')
  })

  plot_code <- reactive({
    paste0(
'# both bands over heights 0 to 200
grid <- data.frame(height = 0:200)
ci <- predict(model, grid, interval = "confidence", level = L)
pi <- predict(model, grid, interval = "prediction", level = L)
plot(plants$height, plants$biomass, type = "n", xlim = c(0, 200),
     ylim = range(pi, plants$biomass, -5.65 + 0.279 * c(0, 200)),
     xlab = "Height (cm)", ylab = "Biomass (g)")
rect(min(plants$height), par("usr")[3], max(plants$height), par("usr")[4],
     col = "grey92", border = NA)                     # where the data are
points(plants$height, plants$biomass, pch = 16, col = "steelblue")
abline(model, col = "#002855", lwd = 2)                # fitted line
lines(grid$height, ci[, "lwr"], col = "#002855", lty = 2)   # confidence band
lines(grid$height, ci[, "upr"], col = "#002855", lty = 2)
lines(grid$height, pi[, "lwr"], col = "#0070C0", lty = 3, lwd = 2)  # prediction band
lines(grid$height, pi[, "upr"], col = "#0070C0", lty = 3, lwd = 2)
', if (input$show_true) 'abline(-5.65, 0.279, col = "orange", lwd = 2)            # true line
' else '',
'
# at height h: the two intervals as bars, and the normal curve of biomass
abline(v = h, col = "grey60")
segments(h - 2, pred[, "lwr"], h - 2, pred[, "upr"], col = "#0070C0", lwd = 4)
segments(h + 2, conf[, "lwr"], h + 2, conf[, "upr"], col = "#002855", lwd = 7)
y <- seq(pred[, "fit"] - 4 * sigma(model), pred[, "fit"] + 4 * sigma(model), length.out = 100)
lines(h + 4 + 20 * dnorm(y, pred[, "fit"], sigma(model)) / dnorm(0, 0, sigma(model)), y,
      col = "#0070C0")
points(h, -5.65 + 0.279 * h, pch = 21, bg = "orange", cex = 1.6)   # true mean at h
legend("topleft", bty = "n", cex = 0.85,
       legend = c("fitted line", "confidence band", "prediction band", "true line"),
       col = c("#002855", "#002855", "#0070C0", "orange"),
       lty = c(1, 2, 3, 1), lwd = 2)
')
  })

  res <- reactive({
    env <- new.env()
    eval(parse(text = calc_code()), envir = env)
    env
  })

  output$plot <- renderPlot({
    env <- res()
    par(mar = c(4.5, 4.5, 1, 1))
    eval(parse(text = plot_code()), envir = env)
  })

  output$stats <- renderPrint({
    e <- res()
    rng <- range(e$plants$height)
    inside <- e$h >= rng[1] && e$h <= rng[2]
    cat(sprintf("Seed %d,  n = %d,  height h = %d cm,  level = %s\n\n",
                seed(), e$n, e$h, input$level))
    cat("Confidence interval (where the average plant at h is):\n")
    print(round(e$conf, 3))
    cat("Prediction interval (where one new plant at h lands):\n")
    print(round(e$pred, 3))
    cat(sprintf("\nWidths: confidence = %.3f g,  prediction = %.3f g\n",
                e$conf[, "upr"] - e$conf[, "lwr"], e$pred[, "upr"] - e$pred[, "lwr"]))
    cat(sprintf("sigma(model) = %.3f g  (true sigma = %s)\n", sigma(e$model), e$sigma))
    cat(sprintf("Data heights: %.1f to %.1f cm  ->  %s\n", rng[1], rng[2],
                if (inside) "interpolating" else "extrapolating"))
    cat(sprintf("Of 1000 new plants simulated at h = %d cm, %.1f%% fall inside the prediction interval.\n",
                e$h, 100 * mean(e$new_biomass >= e$pred[, "lwr"] & e$new_biomass <= e$pred[, "upr"])))
  })

  output$code_out <- renderText({
    paste0(calc_code(), "\n", plot_code())
  })
}

shinyApp(ui, server)

---

Residual Diagnostics Gallery

Each broken assumption (L I N E) leaves its own signature in plot(model), and each has its own fix. Pick a scenario from the lecture (Part C, sections 9 to 11d), find the signature in the four diagnostic panels, then tick “Apply the fix” and watch it go away.

#| '!! shinylive warning !!': |
#|   shinylive does not work in self-contained HTML documents.
#|   Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 1000

library(shiny)

# ---- the 40 wheat plants (data/wheat_plants.csv: plant, height, biomass) ----
wheat_df <- data.frame(
  plant = 1:40,
  height = c(86.9, 71.6, 73.1, 75.9, 91, 94.4, 71.5, 71.8, 77.8, 72.8, 88.7,
             86.6, 58.6, 101.8, 67.8, 76.8, 75.1, 60.3, 62.3, 86.5, 77.4,
             55.7, 78.9, 51.1, 83.6, 45.6, 105, 89.6, 64.7, 76, 70.7, 70.5,
             92.8, 62.4, 90.6, 78.5, 76.4, 89.1, 75.9, 66.1),
  biomass = c(18.1, 15.5, 14.8, 15.6, 21.6, 17.2, 16, 17.4, 18.2, 20.1, 19.9,
              12.8, 11.6, 24.4, 10.5, 16.6, 14.3, 15, 9.4, 13.9, 16.6, 7.5,
              23.3, 8.4, 23.5, 5.8, 24.5, 18.1, 11.5, 12.6, 12.6, 18.1, 13.8,
              9.8, 20.7, 12.9, 12.4, 23.1, 17, 10.7)
)
wheat_url <- "https://greymonroe.github.io/PLS_206/data/wheat_plants.csv"

base_seed <- c(good = 3, curve = 1, fan = 2, counts = 4, typo = NA, pots = 22)

# ---- the code for each scenario, as text: the app runs this same text ----
# returns list(data = ..., model = ..., plot = ..., print = ...)
scenario_code <- function(scn, fix, seed, typo_val) {
  tail_plots <- c(
    "plot(model, which = c(1, 2, 3, 5))",
    "hist(resid(model), col = \"steelblue\", main = \"Histogram of residuals\", xlab = \"Residual\")",
    "par(mfrow = c(1, 1))")
  r2_line <- "summary(model)$r.squared"
  common_print <- c(
    "round(coef(summary(model)), 3)",
    r2_line,
    "shapiro.test(resid(model))$p.value",
    "which.max(cooks.distance(model))",
    "max(cooks.distance(model))")

  if (scn == "good") {
    data <- c(sprintf("set.seed(%d)", seed),
              "x <- runif(100, 0, 10)",
              "y <- 2 + 0.5 * x + rnorm(100)",
              "dat <- data.frame(x, y)")
    model <- "model <- lm(y ~ x, data = dat)"
    plot <- c("par(mfrow = c(2, 3))",
              "plot(dat$x, dat$y, col = \"steelblue\", pch = 16, xlab = \"x\", ylab = \"y\")",
              "abline(model, col = \"#002855\", lwd = 2)",
              tail_plots)
    prnt <- common_print
  }

  if (scn == "curve") {
    data <- c(sprintf("set.seed(%d)", seed),
              "water  <- runif(100, 0, 60)   # mm per week",
              "growth <- 20 + 2.2 * water - 0.03 * water^2 + rnorm(100, sd = 4)",
              "pots   <- data.frame(water, growth)")
    if (!fix) {
      model <- "model <- lm(growth ~ water, data = pots)"
      line <- "abline(model, col = \"#002855\", lwd = 2)"
    } else {
      model <- c("# fix: add a squared term (I() protects the ^ inside a formula)",
                 "model <- lm(growth ~ water + I(water^2), data = pots)")
      line <- "lines(0:60, predict(model, data.frame(water = 0:60)), col = \"#002855\", lwd = 2)"
    }
    plot <- c("par(mfrow = c(2, 3))",
              "plot(pots$water, pots$growth, col = \"steelblue\", pch = 16, xlab = \"Water (mm/week)\", ylab = \"Growth (g)\")",
              line, tail_plots)
    prnt <- common_print
  }

  if (scn == "fan") {
    data <- c(sprintf("set.seed(%d)", seed),
              "dose  <- runif(150, 0, 60)   # fertilizer, kg/ha",
              "yield <- exp(1.5 + 0.03 * dose + rnorm(150, sd = 0.25))",
              "plots <- data.frame(dose, yield)")
    if (!fix) {
      model <- "model <- lm(yield ~ dose, data = plots)"
      sc <- "plot(plots$dose, plots$yield, col = \"steelblue\", pch = 16, xlab = \"Fertilizer (kg/ha)\", ylab = \"Yield (t/ha)\")"
    } else {
      model <- c("# fix: log the response. The slope becomes a percent change.",
                 "model <- lm(log(yield) ~ dose, data = plots)")
      sc <- "plot(plots$dose, log(plots$yield), col = \"steelblue\", pch = 16, xlab = \"Fertilizer (kg/ha)\", ylab = \"log(Yield)\")"
    }
    plot <- c("par(mfrow = c(2, 3))", sc,
              "abline(model, col = \"#002855\", lwd = 2)", tail_plots)
    prnt <- common_print
    if (fix) prnt <- c(prnt, "exp(coef(model)[\"dose\"])   # x this much per kg/ha")
  }

  if (scn == "counts") {
    data <- c(sprintf("set.seed(%d)", seed),
              "size  <- runif(80, 1, 10)   # plant size",
              "seeds <- rpois(80, lambda = exp(-1 + 0.35 * size))   # a count, many zeros",
              "pods  <- data.frame(size, seeds)")
    if (!fix) {
      model <- "model <- lm(seeds ~ size, data = pods)"
      line <- "abline(model, col = \"#002855\", lwd = 2)"
    } else {
      model <- c("# fix: a model for counts, a Poisson GLM (Week 3)",
                 "model <- glm(seeds ~ size, data = pods, family = poisson)")
      line <- "lines(1:10, predict(model, data.frame(size = 1:10), type = \"response\"), col = \"#002855\", lwd = 2)"
    }
    plot <- c("par(mfrow = c(2, 3))",
              "plot(pods$size, pods$seeds, col = \"steelblue\", pch = 16, xlab = \"Plant size\", ylab = \"Seeds per plant\")",
              line, tail_plots)
    prnt <- common_print
    if (fix) prnt[2] <- "1 - model$deviance / model$null.deviance   # no R2 for a GLM: share of deviance explained"
  }

  if (scn == "typo") {
    data <- c(sprintf("wheat <- read.csv(\"%s\")", wheat_url),
              "clean <- lm(biomass ~ height, data = wheat)   # the real data")
    if (!fix) {
      data <- c(data,
                "typo <- wheat",
                sprintf("typo$biomass[7] <- %s   # the real value is 16", format(typo_val)))
      model <- "model <- lm(biomass ~ height, data = typo)"
      sc <- c("plot(typo$height, typo$biomass, col = \"steelblue\", pch = 16, xlab = \"Height (cm)\", ylab = \"Biomass (g)\")",
              "points(typo$height[7], typo$biomass[7], col = \"#FF0000\", pch = 1, cex = 2.5, lwd = 2)",
              "abline(clean, col = \"#00B050\", lwd = 2, lty = 2)")
    } else {
      model <- c("# fix: not statistical. Go back to the data sheet and restore plant 7.",
                 "model <- lm(biomass ~ height, data = wheat)")
      sc <- "plot(wheat$height, wheat$biomass, col = \"steelblue\", pch = 16, xlab = \"Height (cm)\", ylab = \"Biomass (g)\")"
    }
    plot <- c("par(mfrow = c(2, 3))", sc,
              "abline(model, col = \"#002855\", lwd = 2)", tail_plots)
    prnt <- c("coef(clean)   # the real line", common_print)
  }

  if (scn == "pots") {
    data <- c(sprintf("set.seed(%d)", seed),
              "pot_h <- rnorm(8, sd = 4)     # each pot has its own soil...",
              "pot_b <- rnorm(8, sd = 2.5)   # ...and its own luck",
              "pot   <- rep(1:8, each = 5)   # 5 tillers per pot, 8 pots",
              "height  <- 70 + pot_h[pot] + rnorm(40, sd = 3)",
              "biomass <- 2 + 0.2 * height + pot_b[pot] + rnorm(40, 0, 1)",
              "pots <- data.frame(pot, height, biomass)",
              "pot_means <- aggregate(cbind(height, biomass) ~ pot, pots, mean)")
    if (!fix) {
      model <- "model <- lm(biomass ~ height, data = pots)   # 40 tillers treated as 40 plants"
      d <- "pots"
    } else {
      model <- "model <- lm(biomass ~ height, data = pot_means)   # one value per pot, n = 8"
      d <- "pot_means"
    }
    plot <- c("par(mfrow = c(2, 3))",
              sprintf("plot(%s$height, %s$biomass, col = %s$pot, pch = 16, xlab = \"Height (cm)\", ylab = \"Biomass (g)\")", d, d, d),
              "abline(model, col = \"#002855\", lwd = 2)",
              sprintf("plot(fitted(model), resid(model), col = %s$pot, pch = 16, xlab = \"Fitted values\", ylab = \"Residuals\", main = \"Residuals vs Fitted (colour = pot)\")", d),
              "abline(h = 0, lty = 2)",
              "plot(model, which = c(2, 3, 5))",
              tail_plots[2:3])
    prnt <- c(common_print, "round(pot_means, 2)")
    if (!fix) prnt <- c(prnt, "round(coef(summary(lm(biomass ~ height, data = pot_means))), 3)   # the honest SE, n = 8")
  }

  list(data = data, model = model, plot = plot, print = prnt)
}

captions <- list(
  good   = c("Nothing to look for: a shapeless cloud, points on the Q-Q line, a flat scale-location, no point past the Cook's contours.",
             "Nothing to fix: the assumptions hold."),
  curve  = c("Look at residuals vs fitted: an arch. The line is the wrong shape (L, linearity).",
             "Fixed with a squared term: the arch is gone from residuals vs fitted."),
  fan    = c("Look at residuals vs fitted and scale-location: a fan that widens to the right (E, equal variance).",
             "Fixed by logging yield: the fan closes and scale-location goes flat."),
  counts = c("Look at the Q-Q plot and histogram: counts give lumpy, skewed residuals (N, normality).",
             "A Poisson GLM is the model for counts (Week 3). Its residual plots read differently, so judge it with Week 3 tools."),
  typo   = c("Look at residuals vs leverage: plant 7 has the largest Cook's distance and drags the line up (green dashed = the real data). Its height is near the middle, so its leverage is low; slide the typo and watch Cook's distance follow the residual.",
             "Fixed at the source: with the real value restored, no point stands out."),
  pots   = c("Look at residuals vs fitted coloured by pot: residuals clump by pot (I, independence). 40 tillers are not 40 plants.",
             "Fixed by averaging each pot: n = 8, and the slope's SE grows to its honest size.")
)

ui <- fluidPage(
  sidebarLayout(
    sidebarPanel(
      selectInput("scn", "Scenario",
                  choices = c("All good"                         = "good",
                              "1. A curve (water and growth)"    = "curve",
                              "2. A fan (fertilizer and yield)"  = "fan",
                              "3. Counts (seeds per plant)"      = "counts",
                              "4. One bad point (a typo)"        = "typo",
                              "5. Not independent (tillers in pots)" = "pots")),
      conditionalPanel("input.scn == 'typo'",
        sliderInput("typo_val", "Typo value for plant 7", 0, 200, 164, 1),
        helpText("The real value is 16 g. Slide it away from 16 and watch Cook's distance for plant 7 grow.")
      ),
      checkboxInput("fix", "Apply the fix", FALSE),
      conditionalPanel("input.scn != 'typo'",
        actionButton("resample", "New sample")
      ),
      textOutput("seed_txt"),
      hr(),
      checkboxInput("show_code", "Show R code", FALSE),
      helpText("Try: each scenario with the fix off, then on. Does the signature you found disappear?")
    ),
    mainPanel(
      plotOutput("diag_plot", height = "620px"),
      tags$p(tags$strong(textOutput("caption", inline = TRUE))),
      verbatimTextOutput("stats"),
      conditionalPanel("input.show_code == true",
        verbatimTextOutput("code_out")
      )
    )
  )
)

server <- function(input, output, session) {

  clicks <- reactiveValues(good = 0, curve = 0, fan = 0, counts = 0, pots = 0)
  observeEvent(input$resample, {
    s <- input$scn
    if (s != "typo") clicks[[s]] <- clicks[[s]] + 1
  })

  seed <- reactive({
    s <- input$scn
    if (s == "typo") NA else base_seed[[s]] + clicks[[s]]
  })

  code <- reactive({
    scenario_code(input$scn, isTRUE(input$fix), seed(), input$typo_val)
  })

  # run the same code the student sees (the wheat data are embedded here,
  # read.csv() in the student version)
  env <- reactive({
    cd <- code()
    e <- new.env(parent = environment(scenario_code))   # the app's top level, where wheat_df lives
    data_lines <- cd$data
    if (input$scn == "typo") {
      data_lines <- sub("read.csv\\(.*\\)", "wheat_df", data_lines)
    }
    eval(parse(text = c(data_lines, cd$model)), envir = e)
    e
  })

  output$seed_txt <- renderText({
    if (input$scn == "typo") "No randomness here: these are the 40 real plants."
    else paste0("Seed: ", seed())
  })

  output$diag_plot <- renderPlot({
    e <- env()
    eval(parse(text = code()$plot), envir = e)
  })

  output$caption <- renderText({
    captions[[input$scn]][if (isTRUE(input$fix)) 2 else 1]
  })

  output$stats <- renderText({
    e <- env()
    out <- character(0)
    for (ln in code()$print) {
      val <- eval(parse(text = ln), envir = e)
      out <- c(out, paste(">", ln), capture.output(print(val)))
    }
    paste(out, collapse = "\n")
  })

  output$code_out <- renderText({
    cd <- code()
    paste(c("# Residual diagnostics: paste into R to reproduce this view",
            "", cd$data, "", cd$model, "", cd$plot, "", cd$print),
          collapse = "\n")
  })
}

shinyApp(ui, server)

---

The Third Variable: Without vs With

Pick a causal structure and compare the slope of x from lm(y ~ x) (red) with the slope of x from lm(y ~ x + z) (navy), against the true direct effect (orange). At the default settings the confounder, the two paths that cancel and the collider use the same seed and n as the lecture code (Part E, sections 16 to 17b), so the numbers match. The mediator is not in the lecture code. Then move the sliders: whether adding the third variable fixes the slope or breaks it depends on which way the arrows point.

#| '!! shinylive warning !!': |
#|   shinylive does not work in self-contained HTML documents.
#|   Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 900

library(shiny)

# scenario settings: names, lecture seed, lecture n, default direct effect
scen <- list(
  confounder = list(x = "nutrients", z = "rainfall", y = "yield",
                    seed = 2, n = 200, b = 0,
                    advice = "Confounder: adjust for it. Rainfall drives both nutrients and yield; without it the nutrient slope is borrowed from rainfall."),
  cancel     = list(x = "temperature", z = "insects", y = "yield",
                    seed = 1, n = 100, b = 1,
                    advice = "Two paths that cancel: adjust for it. Heat raises yield directly but also cuts insects, which raise yield; with insects in the model we see temperature's direct effect."),
  mediator   = list(x = "fertilizer", z = "height", y = "yield",
                    seed = 7, n = 200, b = 0.3,
                    advice = "Mediator: depends on the question. Without height in the model we get the total effect of fertilizer; with height we get the direct effect (the part not passing through plant height)."),
  collider   = list(x = "x", z = "c", y = "y",
                    seed = 6, n = 300, b = 0,
                    advice = "Collider: leave it out. x and y both cause c; adjusting for c invents an x-y association.")
)

ui <- fluidPage(
  sidebarLayout(
    sidebarPanel(
      selectInput("structure", "Structure",
                  choices = c("Confounder (rainfall -> nutrients, rainfall -> yield)" = "confounder",
                              "Two paths that cancel (temperature -> insects -> yield)" = "cancel",
                              "Mediator (fertilizer -> plant height -> yield)" = "mediator",
                              "Collider (x -> c <- y)" = "collider")),
      sliderInput("b", "Direct effect of x on y (the truth)", -1, 1.5, 0, 0.1),
      sliderInput("a", "Strength of the arrows through the third variable", 0, 2, 1, 0.1),
      sliderInput("n", "n", 30, 500, 200, 10),
      actionButton("resample", "New sample"),
      hr(),
      checkboxInput("show_code", "Show R code", FALSE),
      helpText("Try: Confounder with the direct effect at 0. Without rainfall, nutrients look significant; with rainfall, the slope drops to about 0."),
      helpText("Try: Collider with the direct effect at 0. Adding c makes a strong negative slope appear from nothing.")
    ),
    mainPanel(
      plotOutput("plot", height = "430px"),
      verbatimTextOutput("stats"),
      conditionalPanel("input.show_code == true",
        verbatimTextOutput("code_out")
      )
    )
  )
)

server <- function(input, output, session) {

  # one seed per scenario, starting at the lecture's; +1 per click
  seeds <- reactiveValues(confounder = 2, cancel = 1, mediator = 7, collider = 6)

  observeEvent(input$structure, {
    s <- scen[[input$structure]]
    updateSliderInput(session, "b", value = s$b)
    updateSliderInput(session, "a", value = 1)
    updateSliderInput(session, "n", value = s$n)
  }, ignoreInit = TRUE)

  observeEvent(input$resample, {
    seeds[[input$structure]] <- seeds[[input$structure]] + 1
  })

  sim <- reactive({
    st <- input$structure
    a <- input$a; b <- input$b; n <- input$n
    seed <- seeds[[st]]
    set.seed(seed)
    if (st == "confounder") {
      z <- rnorm(n); x <- -a * z + rnorm(n); y <- b * x + a * z + rnorm(n)
    } else if (st == "cancel") {
      x <- rnorm(n); z <- rnorm(n) - a * x; y <- b * x + a * z + rnorm(n)
    } else if (st == "mediator") {
      x <- rnorm(n); z <- a * x + rnorm(n); y <- b * x + a * z + rnorm(n)
    } else {
      x <- rnorm(n); y <- b * x + rnorm(n); z <- a * (x + y) + rnorm(n, sd = 0.5)
    }
    list(x = x, y = y, z = z, a = a, b = b, n = n, seed = seed, st = st,
         m1 = lm(y ~ x), m2 = lm(y ~ x + z))
  })

  draw_dag <- function(st, a, b, s) {
    par(mar = c(1, 1, 3, 1))
    plot(NA, xlim = c(0, 1), ylim = c(0, 1), axes = FALSE, xlab = "", ylab = "",
         main = "The truth (causal arrows)")
    px <- c(0.15, 0.25); py <- c(0.85, 0.25); pz <- c(0.5, 0.8)
    edge <- function(from, to, coef, col) {
      d <- to - from; u <- d / sqrt(sum(d^2))
      p0 <- from + u * 0.13; p1 <- to - u * 0.13
      if (coef == 0) {
        arrows(p0[1], p0[2], p1[1], p1[2], col = "grey60", lty = 3, lwd = 1.5, length = 0.12)
      } else {
        arrows(p0[1], p0[2], p1[1], p1[2], col = col, lwd = 1 + 3 * abs(coef), length = 0.12)
      }
      mid <- (p0 + p1) / 2 + c(-u[2], u[1]) * 0.05
      text(mid[1], mid[2], format(round(coef, 2)), col = col, font = 2)
    }
    if (st == "confounder") {
      edge(pz, px, -a, "#7030A0"); edge(pz, py, a, "#7030A0")
    } else if (st == "cancel" || st == "mediator") {
      edge(px, pz, if (st == "cancel") -a else a, "#7030A0"); edge(pz, py, a, "#7030A0")
    } else {
      edge(px, pz, a, "#7030A0"); edge(py, pz, a, "#7030A0")
    }
    edge(px, py, b, "orange")
    box_lab <- function(p, lab) {
      rect(p[1] - 0.12, p[2] - 0.07, p[1] + 0.12, p[2] + 0.07, col = "white", border = "#002855", lwd = 2)
      text(p[1], p[2], lab, cex = 1.05)
    }
    box_lab(px, paste0(s$x, " (x)")); box_lab(py, paste0(s$y, " (y)")); box_lab(pz, paste0(s$z, " (z)"))
  }

  output$plot <- renderPlot({
    d <- sim(); s <- scen[[d$st]]
    layout(matrix(1:2, 1), widths = c(1, 1.5))
    draw_dag(d$st, d$a, d$b, s)
    par(mar = c(4.5, 4.5, 3, 1))
    pal <- colorRampPalette(c("#0070C0", "white", "#FF0000"))(100)
    cols <- pal[cut(d$z, 100)]
    plot(d$x, d$y, pch = 21, bg = cols, col = "grey40",
         xlab = s$x, ylab = s$y, main = paste0("Points coloured by ", s$z, " (blue low, red high)"))
    abline(d$m1, col = "#FF0000", lwd = 3)
    cf <- coef(d$m2)
    abline(cf[1] + cf[3] * mean(d$z), cf[2], col = "#002855", lwd = 3)
    abline(mean(d$y) - d$b * mean(d$x), d$b, col = "orange", lwd = 3, lty = 2)
    legend("topleft", bty = "n", lwd = 3, lty = c(1, 1, 2),
           col = c("#FF0000", "#002855", "orange"),
           legend = c(sprintf("without %s: slope %.2f", s$z, coef(d$m1)[2]),
                      sprintf("with %s: slope %.2f", s$z, cf[2]),
                      sprintf("truth (direct effect): %.2f", d$b)))
  })

  output$stats <- renderPrint({
    d <- sim(); s <- scen[[d$st]]
    t1 <- round(coef(summary(d$m1)), 3); rownames(t1)[2] <- s$x
    t2 <- round(coef(summary(d$m2)), 3); rownames(t2)[2:3] <- c(s$x, s$z)
    cat(sprintf("Without %s:  lm(%s ~ %s)\n", s$z, s$y, s$x)); print(t1)
    cat(sprintf("\nWith %s:  lm(%s ~ %s + %s)\n", s$z, s$y, s$x, s$z)); print(t2)
    cat(sprintf("\nTrue direct effect of %s on %s: %s\n", s$x, s$y, format(d$b)))
    if (d$st == "mediator")
      cat(sprintf("True total effect (direct + path through %s): %s\n", s$z, format(d$b + d$a^2)))
    cat(sprintf("cor(%s, %s) = %.3f\n\n", s$x, s$z, cor(d$x, d$z)))
    cat(s$advice, "\n")
  })

  output$code_out <- renderPrint({
    d <- sim(); s <- scen[[d$st]]
    X <- s$x; Y <- s$y; Z <- s$z
    sims <- switch(d$st,
      confounder = c(sprintf("%s <- rnorm(n)", Z),
                     sprintf("%s <- -a * %s + rnorm(n)", X, Z),
                     sprintf("%s <- b * %s + a * %s + rnorm(n)", Y, X, Z)),
      cancel     = c(sprintf("%s <- rnorm(n)", X),
                     sprintf("%s <- rnorm(n) - a * %s", Z, X),
                     sprintf("%s <- b * %s + a * %s + rnorm(n)", Y, X, Z)),
      mediator   = c(sprintf("%s <- rnorm(n)", X),
                     sprintf("%s <- a * %s + rnorm(n)", Z, X),
                     sprintf("%s <- b * %s + a * %s + rnorm(n)", Y, X, Z)),
      collider   = c(sprintf("%s <- rnorm(n)", X),
                     sprintf("%s <- b * %s + rnorm(n)", Y, X),
                     sprintf("%s <- a * (%s + %s) + rnorm(n, sd = 0.5)", Z, X, Y)))
    cat("# R code to reproduce this view\n\n")
    cat(sprintf("set.seed(%d)\n", d$seed))
    cat(sprintf("n <- %d\n", d$n))
    cat(sprintf("a <- %s   # strength of the arrows through %s\n", format(d$a), Z))
    cat(sprintf("b <- %s   # direct effect of %s on %s (the truth)\n\n", format(d$b), X, Y))
    cat(paste(sims, collapse = "\n"), "\n\n")
    cat(sprintf("without <- lm(%s ~ %s)\n", Y, X))
    cat(sprintf("with_z  <- lm(%s ~ %s + %s)\n", Y, X, Z))
    cat("round(coef(summary(without)), 3)\n")
    cat("round(coef(summary(with_z)), 3)\n")
    cat(sprintf("cor(%s, %s)\n\n", X, Z))
    cat(sprintf("# colour points by %s: blue low, red high\n", Z))
    cat("pal <- colorRampPalette(c(\"#0070C0\", \"white\", \"#FF0000\"))(100)\n")
    cat(sprintf("plot(%s, %s, pch = 21, col = \"grey40\", bg = pal[cut(%s, 100)])\n", X, Y, Z))
    cat(sprintf("abline(without, col = \"#FF0000\", lwd = 3)   # without %s\n", Z))
    cat(sprintf("cf <- coef(with_z)   # the %s slope, drawn at the mean of %s\n", X, Z))
    cat(sprintf("abline(cf[1] + cf[3] * mean(%s), cf[2], col = \"#002855\", lwd = 3)\n", Z))
    cat(sprintf("abline(mean(%s) - b * mean(%s), b, col = \"orange\", lwd = 3, lty = 2)   # the truth\n", Y, X))
  })
}

shinyApp(ui, server)

---

Multicollinearity and the VIF

Slide the correlation between temperature and precipitation up and watch the histogram of 200 repeat studies get wider: the model still fits, but each slope becomes unstable. Then tick the microbiome box (the lecture’s case, Part E section 18): the VIF tells us how much each standard error has been inflated, and why.

#| '!! shinylive warning !!': |
#|   shinylive does not work in self-contained HTML documents.
#|   Please set `embed-resources: false` in your metadata.
#| standalone: true
#| viewerHeight: 960

library(shiny)

ui <- fluidPage(
  sidebarLayout(
    sidebarPanel(
      sliderInput("r", "Correlation between the two predictors (r)",
                  min = 0, max = 0.99, value = 0.5, step = 0.01),
      sliderInput("n", "n", min = 20, max = 500, value = 200, step = 10),
      hr(),
      checkboxInput("micro", "Add microbiome = temperature + precipitation + small noise", FALSE),
      conditionalPanel("input.micro == true",
        sliderInput("micro_sd", "microbiome noise SD", min = 0.01, max = 1, value = 0.1, step = 0.01),
        helpText("With no noise at all, microbiome is exactly temperature + precipitation.",
                 "That is perfect collinearity: R cannot separate the three slopes and returns NA for one of them.")
      ),
      hr(),
      actionButton("resample", "New sample"),
      hr(),
      checkboxInput("show_code", "Show R code", FALSE),
      helpText("To match the lecture numbers: r = 0, n = 200, microbiome ticked with noise SD 0.1."),
      helpText("The truth: growth = 1 x temperature + 1 x precipitation + noise.",
               "Each slope's SE grows by sqrt(VIF), where VIF = 1 / (1 - R2 of that predictor on the others).")
    ),
    mainPanel(
      plotOutput("vifPlot", height = "380px"),
      verbatimTextOutput("stats"),
      conditionalPanel("input.show_code == true",
        verbatimTextOutput("code_out")
      )
    )
  )
)

server <- function(input, output, session) {

  seed <- reactive(4 + input$resample)

  sim <- reactive({
    n   <- input$n
    r   <- input$r
    mic <- isTRUE(input$micro)
    msd <- input$micro_sd
    set.seed(seed())

    make_data <- function() {
      temperature   <- rnorm(n)
      precipitation <- r * temperature + sqrt(1 - r^2) * rnorm(n)
      d <- data.frame(temperature, precipitation)
      if (mic) d$microbiome <- temperature + precipitation + rnorm(n, sd = msd)
      d$growth <- temperature + precipitation + rnorm(n)
      d
    }

    dat <- make_data()
    form <- if (mic) growth ~ temperature + precipitation + microbiome
            else     growth ~ temperature + precipitation
    model <- lm(form, data = dat)

    # 200 repeat studies at the same settings (continuing the same random stream)
    slopes <- replicate(200, {
      d <- make_data()
      X <- cbind(1, as.matrix(d[, setdiff(names(d), "growth")]))
      lm.fit(X, d$growth)$coefficients[2]
    })

    preds <- setdiff(names(dat), "growth")
    vif_tab <- t(sapply(preds, function(p) {
      others <- setdiff(preds, p)
      r2 <- summary(lm(reformulate(others, p), data = dat))$r.squared
      c(R2_j = r2, VIF = 1 / (1 - r2), sqrt_VIF = sqrt(1 / (1 - r2)))
    }))

    list(dat = dat, model = model, slopes = slopes, vif = vif_tab)
  })

  output$vifPlot <- renderPlot({
    s <- sim()
    d <- s$dat
    par(mfrow = c(1, 2), mar = c(4.5, 4.5, 3, 1))

    plot(d$temperature, d$precipitation, pch = 16,
         col = adjustcolor("steelblue", 0.6),
         xlab = "temperature", ylab = "precipitation",
         main = sprintf("r = %.2f", cor(d$temperature, d$precipitation)))

    cf  <- coef(summary(s$model))
    est <- cf["temperature", "Estimate"]
    se  <- cf["temperature", "Std. Error"]
    xl  <- range(c(s$slopes, 1, est - 2 * se, est + 2 * se))
    h <- hist(s$slopes, breaks = 25, col = adjustcolor("#002855", 0.7), border = "white",
              xlim = xl, xlab = "estimated temperature slope",
              main = "What 200 repeat studies would find")
    abline(v = 1, col = "orange", lwd = 3)
    yb <- max(h$counts) * 0.9
    segments(est - 2 * se, yb, est + 2 * se, yb, col = "#FF0000", lwd = 3)
    points(est, yb, pch = 16, col = "#FF0000", cex = 1.3)
    legend("topright", bty = "n", cex = 0.8,
           legend = c("truth (1)", "this study: estimate +/- 2 SE"),
           col = c("orange", "#FF0000"), lwd = 3)
  })

  output$stats <- renderPrint({
    s  <- sim()
    sm <- summary(s$model)
    f  <- sm$fstatistic
    cat("This sample: round(coef(summary(model)), 3)\n")
    print(round(coef(sm), 3))
    cat(sprintf("\nR2 = %.3f   F = %.1f on %d and %d df,  p = %.3g   (the model as a whole is still fine)\n\n",
                sm$r.squared, f[1], f[2], f[3], pf(f[1], f[2], f[3], lower.tail = FALSE)))
    cat("Each predictor regressed on the others:\n")
    print(round(s$vif, 3))
    cat(sprintf("\nTemperature slope: SD of the 200 repeat-study slopes = %.3f;  SE from summary() = %.3f\n",
                sd(s$slopes), coef(sm)["temperature", "Std. Error"]))
    cat("Rule of thumb: VIF > 5 look closely, > 10 drop or combine\n")
  })

  output$code_out <- renderPrint({
    n   <- input$n
    r   <- input$r
    mic <- isTRUE(input$micro)
    msd <- input$micro_sd
    preds <- c("temperature", "precipitation", if (mic) "microbiome")
    rhs   <- paste(preds, collapse = " + ")

    sim_lines <- c(
      sprintf("temperature   <- rnorm(%d)", n),
      sprintf("precipitation <- %s * temperature + sqrt(1 - %s^2) * rnorm(%d)   # cor about %s", r, r, n, r),
      if (mic) sprintf("microbiome    <- temperature + precipitation + rnorm(%d, sd = %s)   # almost their sum", n, msd),
      sprintf("growth        <- temperature + precipitation + rnorm(%d)   # true slopes are 1", n)
    )

    cat("# R code to reproduce this view\n\n")
    cat(sprintf("set.seed(%d)\n", seed()))
    cat(sim_lines, sep = "\n")
    cat(sprintf("\nmodel <- lm(growth ~ %s)\n", rhs))
    cat("summary(model)   # R2 and the F-test are still fine\n")
    cat("round(coef(summary(model)), 3)\n\n")

    cat("# VIF by hand: regress each predictor on the others\n")
    for (p in preds) {
      others <- paste(setdiff(preds, p), collapse = " + ")
      cat(sprintf("r2_%s <- summary(lm(%s ~ %s))$r.squared\n", p, p, others))
    }
    cat("c(", paste(sprintf("%s = 1 / (1 - r2_%s)", preds, preds), collapse = ",\n  "), ")\n", sep = "")
    cat("sqrt(1 / (1 - r2_temperature))   # how much the temperature SE is inflated\n")
    cat("# car::vif(model) gives the same numbers\n")
    cat("# rule of thumb: VIF > 5 look closely, > 10 drop or combine\n\n")

    cat("# 200 repeat studies at the same settings (no new set.seed: the stream continues)\n")
    cat("slopes <- replicate(200, {\n")
    cat(paste0("  ", sim_lines, collapse = "\n"), "\n", sep = "")
    cat(sprintf("  coef(lm(growth ~ %s))[2]\n", rhs))
    cat("})\n")
    cat("sd(slopes)                                    # spread across studies\n")
    cat("coef(summary(model))[\"temperature\", \"Std. Error\"]   # SE from one study\n\n")

    cat("par(mfrow = c(1, 2))\n")
    cat("plot(temperature, precipitation, pch = 16, col = \"steelblue\",\n")
    cat("     main = paste(\"r =\", round(cor(temperature, precipitation), 2)))\n")
    cat("hist(slopes, breaks = 25, col = \"#002855\", border = \"white\",\n")
    cat("     xlab = \"estimated temperature slope\", main = \"What 200 repeat studies would find\")\n")
    cat("abline(v = 1, col = \"orange\", lwd = 3)\n")
  })
}

shinyApp(ui, server)