Multiple Regression

Bio300B Lecture 8

Richard J. Telford (Richard.Telford@uib.no)

Institutt for biovitenskap, UiB

17 August 2025

Outline

• Types of models
• Interactions
• Model selection
  • Exploratory data analysis
  • Anova for hypothesis testing
• Multi-collinearity
• Autocorrelation
• Reporting statistics

Adding more variables

# Two way anova - two categorical predictors
mod_anova2 <- lm(bill_length_mm ~ sex + species, data = penguins)

# ancova - one categorical one continuous predictor
mod_ancova <- lm(bill_length_mm ~ sex + flipper_length_mm, data = penguins)

# Multiple regression - two continuous
mod_mult <- lm(bill_length_mm ~ body_mass_g + flipper_length_mm, data = penguins)

summary(mod_anova2)

Call:
lm(formula = bill_length_mm ~ sex + species, data = penguins)

Residuals:
   Min     1Q Median     3Q    Max 
-6.087 -1.377 -0.071  1.225 11.013 

Coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
(Intercept)        36.977      0.231   160.2   <2e-16 ***
sexmale             3.694      0.255    14.5   <2e-16 ***
speciesChinstrap   10.010      0.341    29.3   <2e-16 ***
speciesGentoo       8.698      0.287    30.3   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.32 on 329 degrees of freedom
  (11 observations deleted due to missingness)
Multiple R-squared:  0.821, Adjusted R-squared:  0.819 
F-statistic:  503 on 3 and 329 DF,  p-value: <2e-16

summary(mod_ancova)

Call:
lm(formula = bill_length_mm ~ sex + flipper_length_mm, data = penguins)

Residuals:
   Min     1Q Median     3Q    Max 
-9.720 -2.812 -0.812  2.021 19.764 

Coefficients:
                  Estimate Std. Error t value Pr(>|t|)    
(Intercept)        -4.4669     3.2364   -1.38     0.17    
sexmale             2.0727     0.4568    4.54    8e-06 ***
flipper_length_mm   0.2359     0.0163   14.46   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.03 on 330 degrees of freedom
  (11 observations deleted due to missingness)
Multiple R-squared:  0.46,  Adjusted R-squared:  0.457 
F-statistic:  141 on 2 and 330 DF,  p-value: <2e-16

summary(mod_mult)

Call:
lm(formula = bill_length_mm ~ body_mass_g + flipper_length_mm, 
    data = penguins)

Residuals:
   Min     1Q Median     3Q    Max 
-8.806 -2.590 -0.705  1.991 18.829 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)       -3.436694   4.580553   -0.75     0.45    
body_mass_g        0.000662   0.000567    1.17     0.24    
flipper_length_mm  0.221865   0.032348    6.86  3.3e-11 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 4.12 on 339 degrees of freedom
  (2 observations deleted due to missingness)
Multiple R-squared:  0.433, Adjusted R-squared:  0.43 
F-statistic:  129 on 2 and 339 DF,  p-value: <2e-16

Different formula

What is an interaction

Effect of one predictor depends on value of another

Interaction between categorical predictors

mod <- lm(flipper_length_mm ~ species * sex, data = penguins)
summary(mod)

Call:
lm(formula = flipper_length_mm ~ species * sex, data = penguins)

Residuals:
    Min      1Q  Median      3Q     Max 
-15.795  -3.411   0.088   3.459  17.589 

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)    
(Intercept)               187.795      0.662  283.72  < 2e-16 ***
speciesChinstrap            3.941      1.174    3.36  0.00088 ***
speciesGentoo              24.912      0.995   25.04  < 2e-16 ***
sexmale                     4.616      0.936    4.93  1.3e-06 ***
speciesChinstrap:sexmale    3.560      1.661    2.14  0.03278 *  
speciesGentoo:sexmale       4.218      1.397    3.02  0.00274 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.66 on 327 degrees of freedom
  (11 observations deleted due to missingness)
Multiple R-squared:  0.84,  Adjusted R-squared:  0.837 
F-statistic:  342 on 5 and 327 DF,  p-value: <2e-16

Understanding coefficients

#| label: coef-explain-app
#| standalone: true
#| viewerHeight: 650

library(shiny)
library(bslib)

ui <- page_sidebar(
  sidebar = sidebar(checkboxGroupInput("pred",
    "Predictor",
    choiceValues = c(
      "Food", "Temperature",
      "Interaction"
    ), choiceNames = c(
      "Food (A vs. B)",
      "Temperature (Low vs. High)", "Interaction"
    )
  )), card(
    textOutput("formula")
  ), card(
    card_header("Coefficients"),
    tableOutput("coef_table")
  ), card(plotOutput("plot")),
  tags$head(tags$style("#coef_table td{\n                     position:relative;\n                     };\n\n                     ")),
)
server <- function(input, output, session) {
  data <- reactive({
    b0 <- 5
    b1 <- 3
    b2 <- 2
    b3 <- 2
    food <- factor(rep(c("A", "B"), times = 24))
    temperature <- factor(rep(c("Low", "High"), each = 24),
      levels = c("Low", "High")
    )
    response <- b0 + b1 * (food == "B") + b2 * (temperature ==
      "High") + b3 * (food == "B" & temperature ==
      "High") + rnorm(48)
    data.frame(
      food = food, temperature = temperature,
      response = response
    )
  })
  form <- "response ~"
  form2 <- reactive({
    if ("Interaction" %in% input$pred) {
      paste(form, "food * temperature")
    } else if ("Food" %in% input$pred & "Temperature" %in%
      input$pred) {
      paste(form, "food + temperature")
    } else if ("Food" %in% input$pred) {
      paste(form, "food")
    } else if ("Temperature" %in% input$pred) {
      paste(form, "temperature")
    } else {
      paste(form, "1")
    }
  })
  observe({
    if ("Interaction" %in% input$pred) {
      updateCheckboxGroupInput(session, "pred", selected = c(
        "Food",
        "Temperature", "Interaction"
      ))
    }
  })
  model <- reactive({
    lm(form2(), data = data())
  })
  coefs <- reactive({
    coef(model())
  })
  coef_colours <- reactive({
    c(
      `(Intercept)` = "skyblue", foodB = "red", temperatureHigh = "green",
      `foodB:temperatureHigh` = "orange"
    )[names(coefs())]
  })
  coef_table <- reactive({
    c1 <- "<div style=\"width: 100%; height: 100%; z-index: 0; background-color: "
    c2 <- "; position:absolute; top: 0; left: 0; padding:5px;\">\n<span>"
    c3 <- "</span></div>"
    tab <- data.frame(
      Beta = paste0(
        c1, coef_colours(),
        c2, "β", seq_along(coefs()), c3
      ), Coefficent = names(coefs()),
      Estimate = coefs()
    )
    tab
  })
  output$formula <- renderText(paste("Model formula:", form2()))
  output$coef_table <- renderTable(coef_table(), sanitize.text.function = function(x) x)
  output$plot <- renderPlot({
    set.seed(1)
    f <- as.formula(form2())
    fc <- as.character(f)
    ylim <- c(0, max(data()$response))
    par(mar = c(1.5, 2.5, .5, .5), mgp = c(1.5, 0.5, 0))
    if (fc[3] == "1") {
      stripchart(data()$response,
        method = "jitter",
        jitter = 0.1, vertical = TRUE, pch = 1, ylim = ylim,
        ylab = "response"
      )
    } else {
      stripchart(f,
        data = data(), method = "jitter",
        jitter = 0.1, vertical = TRUE, pch = 1, ylim = ylim
      )
    }
    cols <- c("food", "temperature")[c(grepl(
      "food",
      fc[3]
    ), grepl("temperature", fc[3]))]
    if (length(cols) == 0) {
      cols <- "food"
    }
    pred <- predict(model(), newdata = unique(data()[,
      cols,
      drop = FALSE
    ]))
    points(seq_along(pred), pred,
      col = "#832424", pch = 16,
      cex = 4
    )
    xs <- seq_along(pred) - 0.2
    arrows(xs, rep(0, length(pred)), xs, rep(
      coefs()[1],
      length(pred)
    ),
    col = coef_colours()[1], lwd = 4,
    length = 0.1
    )
    pos <- 2 - 0.2
    if ("foodB" %in% names(coefs())) {
      arrows(pos, coefs()[1], pos, coefs()[1] + coefs()["foodB"],
        col = coef_colours()["foodB"], lwd = 4, length = 0.1
      )
      pos <- pos + 1
    }
    if ("temperatureHigh" %in% names(coefs())) {
      arrows(pos, coefs()[1], pos, coefs()[1] + coefs()["temperatureHigh"],
        col = coef_colours()["temperatureHigh"], lwd = 4,
        length = 0.1
      )
      pos <- pos + 1
    }
    if (all(c("foodB", "temperatureHigh") %in% names(coefs()))) {
      arrows(pos, coefs()[1], pos, coefs()[1] + coefs()["foodB"],
        col = coef_colours()["foodB"], lwd = 4, length = 0.1
      )
      arrows(pos, coefs()[1] + coefs()["foodB"], pos,
        coefs()[1] + coefs()["foodB"] + coefs()["temperatureHigh"],
        col = coef_colours()["temperatureHigh"], lwd = 4,
        length = 0.1
      )
      if ("foodB:temperatureHigh" %in% names(coefs())) {
        arrows(pos, coefs()[1] + coefs()["foodB"] +
          coefs()["temperatureHigh"], pos, coefs()[1] +
          coefs()["foodB"] + coefs()["temperatureHigh"] +
          coefs()["temperatureHigh"],
        col = coef_colours()["foodB:temperatureHigh"],
        lwd = 4, length = 0.1
        )
      }
    }
  })
}
shinyApp(ui, server)

Interaction between continuous and categorical predictors

mod_interact <- lm(bill_length_mm ~ sex * flipper_length_mm, data = adelie)
summary(mod_interact)

Call:
lm(formula = bill_length_mm ~ sex * flipper_length_mm, data = adelie)

Residuals:
   Min     1Q Median     3Q    Max 
-6.395 -1.329  0.029  1.427  5.438 

Coefficients:
                          Estimate Std. Error t value Pr(>|t|)    
(Intercept)                39.7935     8.3528    4.76  4.6e-06 ***
sexmale                   -20.2163    11.0648   -1.83    0.070 .  
flipper_length_mm          -0.0135     0.0445   -0.30    0.762    
sexmale:flipper_length_mm   0.1217     0.0583    2.09    0.039 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.11 on 142 degrees of freedom
Multiple R-squared:  0.385, Adjusted R-squared:  0.372 
F-statistic: 29.6 on 3 and 142 DF,  p-value: 6.43e-15

Interaction between two continuous predictors

mod <- lm(flipper_length_mm ~ body_mass_g * bill_length_mm, data = adelie)
summary(mod)

Call:
lm(formula = flipper_length_mm ~ body_mass_g * bill_length_mm, 
    data = adelie)

Residuals:
    Min      1Q  Median      3Q     Max 
-14.653  -3.101  -0.003   3.369  16.585 

Coefficients:
                            Estimate Std. Error t value Pr(>|t|)    
(Intercept)                 2.04e+02   5.72e+01    3.58  0.00048 ***
body_mass_g                -6.63e-03   1.52e-02   -0.44  0.66330    
bill_length_mm             -9.19e-01   1.48e+00   -0.62  0.53531    
body_mass_g:bill_length_mm  3.17e-04   3.89e-04    0.82  0.41546    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 5.79 on 142 degrees of freedom
Multiple R-squared:  0.229, Adjusted R-squared:  0.212 
F-statistic:   14 on 3 and 142 DF,  p-value: 4.62e-08

Power needed for interactions

• Main effect = \(\bar{y_1}\) - \(\bar{y_2}\)
• Interaction = (\(\bar{y_1}\) - \(\bar{y_2}\)) - (\(\bar{y_3}\) - \(\bar{y_4}\))

set.seed(42)
sigma <- 10
n <- 1000
random <- tibble(
  x = sample(c("A", "B"), n, replace = TRUE),
  z = sample(factor(1:2), n, replace = TRUE),
  y = rnorm(n, mean = 0, sd = sigma)
)

lm(y ~ x * z, data = random) |> tidy()

# A tibble: 4 × 5
  term        estimate std.error statistic p.value
  <chr>          <dbl>     <dbl>     <dbl>   <dbl>
1 (Intercept)    0.281     0.615     0.457  0.647 
2 xB            -1.90      0.864    -2.20   0.0280
3 z2             0.759     0.880     0.863  0.388 
4 xB:z2          1.06      1.24      0.849  0.396 

Formula for interactions

y ~ x + z + x:z

y ~ x * z

y ~ (x + z)^2

Use y ~ x + I(x^2) to get a quadratic. Or y ~ poly(x, 2)

mod1 <- lm(flipper_length_mm ~ species + sex + species:sex, data = penguins)
mod2 <- lm(flipper_length_mm ~ species * sex, data = penguins)
mod3 <- lm(flipper_length_mm ~ (species + sex)^2, data = penguins)

coef(mod1)

             (Intercept)         speciesChinstrap            speciesGentoo 
                 187.795                    3.941                   24.912 
                 sexmale speciesChinstrap:sexmale    speciesGentoo:sexmale 
                   4.616                    3.560                    4.218 

coef(mod2)

             (Intercept)         speciesChinstrap            speciesGentoo 
                 187.795                    3.941                   24.912 
                 sexmale speciesChinstrap:sexmale    speciesGentoo:sexmale 
                   4.616                    3.560                    4.218 

coef(mod3)

             (Intercept)         speciesChinstrap            speciesGentoo 
                 187.795                    3.941                   24.912 
                 sexmale speciesChinstrap:sexmale    speciesGentoo:sexmale 
                   4.616                    3.560                    4.218 

Model Selection

You want to the best model!

The best model for what?

• Exploratory data analysis
• Inference (hypothesis testing)
• Predictions

Tredennick et al 2021 A practical guide to selecting models for exploration, inference, and prediction in ecology Ecology. https://doi.org/10.1002/ecy.3336

Exploratory analysis

• Consider all plausible models
• P-values not meaningful
• High type I error rate
• Suggest hypotheses for hypothesis testing with independent data

The more biology you include in the model

Automatic model selection

Last resort - many problems and biases

• Forward selection

• Backwards selection

• All possible models

library(MuMIn)
library(conflicted)
conflict_prefer("filter", "dplyr")
conflict_prefer("select", "dplyr")

penguins2 <- penguins |>
  drop_na() |>
  select(-year)
full <- lm(body_mass_g ~ ., data = penguins2, na.action = "na.fail")

AIC

AIC Akaike information criterion \(2k - 2 \times log(likelihood)\)
(k the number of parameters)

Measure of how well the model fit the data

Penalised for model complexity

AICc correction for small sample sizes

AIC weights - probability model is best of those tested

dredge(full) |> gt::gt()

(Intercept)	bill_depth_mm	bill_length_mm	flipper_length_mm	island	sex	species	df	logLik	AICc	delta	weight
-1461.0	67.22	18.204	15.95	NA	+	+	8	-2354	4724	0.000	7.682e-01
-1500.0	67.58	18.189	16.24	+	+	+	10	-2354	4728	3.543	1.307e-01
-1211.5	74.38	NA	17.54	NA	+	+	7	-2357	4729	4.537	7.948e-02
-1245.4	74.51	NA	17.85	+	+	+	9	-2357	4732	8.027	1.388e-02
-759.1	NA	21.633	17.85	NA	+	+	7	-2360	4734	9.536	6.526e-03
-785.6	NA	21.528	18.14	+	+	+	9	-2359	4738	13.153	1.070e-03
-365.8	NA	NA	20.02	NA	+	+	6	-2364	4741	16.649	1.863e-04
-390.0	NA	NA	20.32	+	+	+	8	-2364	4745	20.118	3.287e-05
844.0	87.93	26.537	NA	NA	+	+	7	-2369	4752	27.263	9.232e-07
829.8	88.33	26.664	NA	+	+	+	9	-2369	4756	31.392	1.172e-07
1577.0	102.04	NA	NA	NA	+	+	6	-2375	4763	38.571	3.236e-09
1576.2	102.15	NA	NA	+	+	+	8	-2375	4767	42.751	4.000e-10
2169.3	NA	32.537	NA	NA	+	+	6	-2378	4768	43.901	2.252e-10
2167.9	NA	32.535	NA	+	+	+	8	-2378	4772	48.086	2.778e-11
-4282.1	141.77	39.718	20.23	NA	NA	+	7	-2385	4784	59.631	8.643e-14
3372.4	NA	NA	NA	NA	+	+	5	-2388	4786	61.394	3.580e-14
-4324.0	141.82	39.599	20.58	+	NA	+	9	-2384	4786	62.049	2.580e-14
3375.8	NA	NA	NA	+	+	+	7	-2388	4790	65.442	4.731e-15
-2219.3	-50.03	NA	35.69	+	+	NA	7	-2397	4807	82.977	7.364e-19
-2115.6	-50.34	6.052	33.94	+	+	NA	8	-2396	4808	83.566	5.488e-19
-4484.7	181.95	NA	25.53	NA	NA	+	6	-2400	4811	86.962	1.004e-19
-4524.1	181.54	NA	25.95	+	NA	+	8	-2399	4814	89.247	3.204e-20
-3720.1	NA	NA	39.30	+	+	NA	6	-2401	4815	90.409	1.792e-20
-3629.6	NA	5.808	37.64	+	+	NA	7	-2401	4816	91.143	1.242e-20
-1709.5	180.78	54.061	NA	NA	NA	+	6	-2405	4822	97.444	5.318e-22
-1735.5	181.30	54.252	NA	+	NA	+	8	-2404	4825	100.777	1.005e-22
-3864.1	NA	60.117	27.54	NA	NA	+	6	-2411	4833	109.078	1.583e-24
-2246.8	-86.95	NA	38.19	NA	+	NA	5	-2412	4834	109.703	1.158e-24
-2288.5	-86.09	-2.329	38.83	NA	+	NA	6	-2412	4836	111.527	4.653e-25
-3902.9	NA	59.868	27.96	+	NA	+	8	-2410	4836	111.692	4.284e-25
-5410.3	NA	NA	46.98	NA	+	NA	4	-2427	4863	138.223	7.424e-31
-5433.5	NA	-5.201	48.21	NA	+	NA	5	-2427	4864	139.128	4.723e-31
-1000.8	256.55	NA	NA	NA	NA	+	5	-2431	4871	146.983	9.299e-33
-1014.6	257.16	NA	NA	+	NA	+	7	-2430	4875	150.703	1.448e-33
-5707.1	51.24	13.748	42.90	+	NA	NA	7	-2438	4891	166.751	4.741e-37
-6114.7	56.67	NA	47.40	+	NA	NA	6	-2441	4895	170.817	6.209e-38
-4013.2	NA	NA	40.61	NA	NA	+	5	-2443	4895	171.067	5.478e-38
-4047.8	NA	NA	41.09	+	NA	+	7	-2442	4897	173.113	1.970e-38
200.5	NA	90.298	NA	NA	NA	+	5	-2444	4899	174.468	1.001e-38
-4355.0	NA	16.883	39.66	+	NA	NA	6	-2445	4902	177.784	1.906e-39
182.0	NA	90.510	NA	+	NA	+	7	-2444	4903	178.295	1.477e-39
-4687.9	NA	NA	44.89	+	NA	NA	5	-2450	4909	184.809	5.685e-41
-6537.6	19.88	NA	51.77	NA	NA	NA	4	-2460	4928	203.674	4.552e-45
-5872.1	NA	NA	50.15	NA	NA	NA	3	-2461	4928	203.836	4.198e-45
-5836.3	NA	4.959	48.89	NA	NA	NA	4	-2461	4929	204.973	2.377e-45
-6445.5	17.84	3.293	50.76	NA	NA	NA	5	-2460	4930	205.355	1.965e-45
4784.5	-155.12	44.769	NA	+	+	NA	7	-2470	4955	230.236	7.767e-51
1747.4	NA	60.616	NA	+	+	NA	6	-2505	5022	297.540	1.886e-65
6551.4	-257.31	36.519	NA	NA	+	NA	5	-2507	5023	299.035	8.929e-66
7580.2	-211.52	NA	NA	+	+	NA	6	-2506	5024	299.241	8.057e-66
3706.2	NA	NA	NA	NA	NA	+	4	-2513	5035	310.285	3.221e-68
3709.7	NA	NA	NA	+	NA	+	6	-2513	5039	314.411	4.093e-69
8779.1	-299.34	NA	NA	NA	+	NA	4	-2529	5066	341.864	4.474e-75
1697.5	-26.23	76.000	NA	+	NA	NA	6	-2537	5086	361.234	2.782e-79
1239.3	NA	76.905	NA	+	NA	NA	5	-2538	5086	361.346	2.632e-79
4375.8	NA	NA	NA	+	+	NA	5	-2562	5133	409.061	1.145e-89
3413.5	-145.51	74.813	NA	NA	NA	NA	4	-2595	5199	474.205	8.187e-104
746.1	NA	74.025	NA	NA	+	NA	4	-2614	5236	512.010	5.056e-112
5590.3	-54.76	NA	NA	+	NA	NA	5	-2615	5240	516.032	6.767e-113
4719.2	NA	NA	NA	+	NA	NA	4	-2618	5244	519.963	9.482e-114
388.8	NA	86.792	NA	NA	NA	NA	3	-2629	5264	539.833	4.594e-118
7520.0	-193.01	NA	NA	NA	NA	NA	3	-2658	5322	598.047	1.050e-130
3862.3	NA	NA	NA	NA	+	NA	3	-2667	5340	615.648	1.582e-134
4207.1	NA	NA	NA	NA	NA	NA	2	-2700	5404	679.945	1.727e-148

Inference

Test small number of a priori hypothesis.

Use anova() to compare nested models

H₀ there is no interaction between sex and species for predicting body mass

mod1 <- lm(body_mass_g ~ species + sex, data = penguins)
mod2 <- lm(body_mass_g ~ species * sex, data = penguins)
anova(mod1, mod2)

Analysis of Variance Table

Model 1: body_mass_g ~ species + sex
Model 2: body_mass_g ~ species * sex
  Res.Df      RSS Df Sum of Sq    F Pr(>F)    
1    329 32979185                             
2    327 31302628  2   1676557 8.76  2e-04 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

car::Anova vs anova

mod2 <- lm(body_mass_g ~ species + sex, data = penguins)
mod3 <- lm(body_mass_g ~ sex + species, data = penguins)

• anova - tests terms sequentially
• car::Anova - marginal test

anova(mod2)

Analysis of Variance Table

Response: body_mass_g
           Df   Sum Sq  Mean Sq F value Pr(>F)    
species     2 1.45e+08 72595110     724 <2e-16 ***
sex         1 3.71e+07 37090262     370 <2e-16 ***
Residuals 329 3.30e+07   100241                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

anova(mod3)

Analysis of Variance Table

Response: body_mass_g
           Df   Sum Sq  Mean Sq F value Pr(>F)    
sex         1 3.89e+07 38878897     388 <2e-16 ***
species     2 1.43e+08 71700792     715 <2e-16 ***
Residuals 329 3.30e+07   100241                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

car::Anova(mod2)

Anova Table (Type II tests)

Response: body_mass_g
            Sum Sq  Df F value Pr(>F)    
species   1.43e+08   2     715 <2e-16 ***
sex       3.71e+07   1     370 <2e-16 ***
Residuals 3.30e+07 329                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

car::Anova(mod3)

Anova Table (Type II tests)

Response: body_mass_g
            Sum Sq  Df F value Pr(>F)    
sex       3.71e+07   1     370 <2e-16 ***
species   1.43e+08   2     715 <2e-16 ***
Residuals 3.30e+07 329                   
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Multicollinearity

• Two or more predictor variables in a multiple regression model are highly correlated. Example: pH and calcium
• Coefficient estimates are unstable
  • erratic change in response to small changes in the model or the data.
• Solve by having lots of data

col <- performance::check_collinearity(full)
col

# Check for Multicollinearity

Low Correlation

   Term  VIF     VIF 95% CI adj. VIF Tolerance Tolerance 95% CI
 island 3.73 [ 3.14,  4.49]     1.39      0.27     [0.22, 0.32]
    sex 2.33 [ 2.00,  2.77]     1.53      0.43     [0.36, 0.50]

Moderate Correlation

              Term  VIF     VIF 95% CI adj. VIF Tolerance Tolerance 95% CI
    bill_length_mm 6.10 [ 5.06,  7.40]     2.47      0.16     [0.14, 0.20]
     bill_depth_mm 6.10 [ 5.06,  7.41]     2.47      0.16     [0.14, 0.20]
 flipper_length_mm 6.80 [ 5.63,  8.26]     2.61      0.15     [0.12, 0.18]

High Correlation

    Term   VIF     VIF 95% CI adj. VIF Tolerance Tolerance 95% CI
 species 63.52 [51.73, 78.06]     2.82      0.02     [0.01, 0.02]

plot(col)

Error: Package `see` required to plot collinearity-check.
  Please install it by running `install.packages("see")`.

Autocorrelation

Linear models assume residuals are independent

If data are spatially or temporally structured, residuals may be correlated.

Positive autocorrelation

• confidence intervals are too narrow
• increases risk of false positive (Type I error)

Luteinising Hormone concentration

library(ggfortify)

Error in library(package, pos = pos, lib.loc = lib.loc, character.only = TRUE, : there is no package called 'ggfortify'

lh <- fortify(lh) |>
  mutate(time = Index * 10) |>
  rename(concentration = Data)

Error in `fortify()`:
! `data` must be a <data.frame>, or an object coercible by `fortify()`,
  or a valid <data.frame>-like object coercible by `as.data.frame()`.
Caused by error in `.prevalidate_data_frame_like_object()`:
! `dim(data)` must return an <integer> of length 2.

mod_lh <- lm(concentration ~ time, data = lh)

Error in eval(predvars, data, env): object 'concentration' not found

augment(mod_lh, interval = "confidence") |>
  ggplot(aes(x = time, y = concentration)) +
  geom_line() +
  geom_point() +
  geom_ribbon(aes(ymin = .lower, ymax = .upper), alpha = 0.2, fill = "#ee5050") +
  geom_line(aes(y = .fitted), colour = "#ee5050") +
  labs(x = "Time, minutes", y = "Hormone concentration")

Error: object 'mod_lh' not found

Detecting autocorrelation

For time series with equally-spaced observations, use autocorrelation function (ACF)

autoplot(acf(mod_lh$residuals, plot = FALSE))

Error: object 'mod_lh' not found

Detecting autocorrelation

Tests for equally-spaced observations

Durbin-Watson test

performance::check_autocorrelation(mod_lh)

Error: object 'mod_lh' not found

Other methods for spatial data and non-equally spaced data

Solution - use a model that accounts for autocorrelation

• generalised least squares (nlme::gls())

Reporting regression results - tables

Make a table - estimates, confidence intervals, p-values

Can calculate confidence intervals with broom::tidy() and format output

library(gt)
tidy(mod_mult, conf.int = TRUE) |>
  select(term, estimate, conf.low, conf.high, p.value) |>
  mutate(p.value = format.pval(p.value, eps = 0.001)) |>
  gt() |>
  fmt_number(
    columns = c(estimate, conf.low, conf.high),
    decimals = 4
  )

term	estimate	conf.low	conf.high	p.value
(Intercept)	−3.4367	−12.4466	5.5732	0.45
body_mass_g	0.0007	−0.0005	0.0018	0.24
flipper_length_mm	0.2219	0.1582	0.2855	<0.001

gtsummary

Or make table with gtsummary

library(gtsummary)

Error in library(package, pos = pos, lib.loc = lib.loc, character.only = TRUE, : there is no package called 'gtsummary'

tbl <- tbl_regression(mod_mult,
  label = list(
    body_mass_g = "Body mass g",
    flipper_length_mm = "Flipper length mm"
  ),
  estimate_fun = label_style_sigfig(digits = 4)
)

Error in tbl_regression(mod_mult, label = list(body_mass_g = "Body mass g", : could not find function "tbl_regression"

tbl

function (src, ...) 
{
    UseMethod("tbl")
}
<bytecode: 0x6161da6799a8>
<environment: namespace:dplyr>

Reporting regression results - inline

Can extract values directly from regression .

• estimate
• confidence intervals
• p-value, degrees of freedom, statistic

Longer flipper length is associated with longer bill length (estimate = `r round(coef(mod_mult)["flipper_length_mm"], 2)` mm/mm, 95% CI = `r round(confint(mod_mult)["flipper_length_mm", ], 2)` p = `r format.pval(tidy(mod_mult) |> filter(term == "flipper_length_mm") |> select("p.value"), eps = 0.001)`).

Longer flipper length is associated with longer bill length (estimate = 0.22 mm/mm, 95% CI = 0.16, 0.29 p = <0.001).

Also function inline_text() from gtsummary

See APA guide for reporting statistics