Day 4: Linear and Multiple Regression

0. Overview

In this practical example, we will use two factors (i) number of factors (Lv) and (ii) factor dissimilarity (dissim) to predict soil decomposition rate (Decom), with the dataset MultiFactor_data.csv.

Goal

  • Explain what linear regressions estimate
  • Fit simple linear regression and extend to multiple regression
  • Compare models (ANOVA, AIC/BIC, ΔR²) and report results clearly

What we use from the data

  • Decom: Soil decomposition rate (response variable).
  • Lv: factor level: number of factors acting simultaneously (numerical predictor).
  • dissim: factor dissimilarity, an index measures how dissimilar the co-acting factors are (continuous predictor)
##Import dataset
# Read the combined dataset (ensure the CSV file is in your working directory)
dat <- read.csv("MultiFactor_data.csv")

Install and load packages

library(dplyr)
library(ggplot2)
library(tidyverse)
library(broom)
library(car)
library(lmtest)

1. Subsetting data for regression analysis

# Keep only variables used in regression analysis and drop missing values
dat2 <- dat %>% 
  dplyr::select(Decom, Lv, dissim) %>% 
  drop_na()

nrow(dat2)
## [1] 343
summary(dat2)
##      Decom              Lv            dissim      
##  Min.   :0.4242   Min.   :0.000   Min.   :0.0000  
##  1st Qu.:0.6342   1st Qu.:1.000   1st Qu.:0.0000  
##  Median :0.6570   Median :1.000   Median :0.0000  
##  Mean   :0.6507   Mean   :2.603   Mean   :0.2165  
##  3rd Qu.:0.6769   3rd Qu.:5.000   3rd Qu.:0.4515  
##  Max.   :0.7735   Max.   :8.000   Max.   :1.0000

2. Understand data patterns

# Distribution of Decom
ggplot(dat2, aes(Decom)) +
  geom_histogram(bins = 30) +
  labs(title = "Distribution of Decomposition Rate (Decom)")

# Scatter with linear smoother: Decom ~ n_factors
ggplot(dat2, aes(Lv, Decom)) +
  geom_point(alpha = 0.6, position = position_jitter(width = 0.05, height = 0)) +
  geom_smooth(method = "lm", se = TRUE) +
  labs(x = "Number of factors", y = "Decom",
       title = "Decomposition vs. Number of Factors")
## `geom_smooth()` using formula = 'y ~ x'

# Optional: pairs plot to see relationships & correlations
GGally::ggpairs(dat2, columns = c("dissim","Lv","Decom"))
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2

3. Build simple linear model

Model formular:

Decom ~ Lv Predict soil decomposition rate from factor number lv.

m1 <- lm(Decom ~ Lv, data = dat2)

summary(m1)     # Coefs, p-values, R²
## 
## Call:
## lm(formula = Decom ~ Lv, data = dat2)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.192035 -0.021474 -0.000894  0.021596  0.138152 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.6673085  0.0030802 216.644  < 2e-16 ***
## Lv          -0.0063842  0.0008416  -7.586 3.17e-13 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.0401 on 341 degrees of freedom
## Multiple R-squared:  0.1444, Adjusted R-squared:  0.1419 
## F-statistic: 57.55 on 1 and 341 DF,  p-value: 3.168e-13
broom::glance(m1)[c("r.squared","adj.r.squared","AIC","BIC")]
## # A tibble: 1 × 4
##   r.squared adj.r.squared    AIC    BIC
##       <dbl>         <dbl>  <dbl>  <dbl>
## 1     0.144         0.142 -1229. -1218.

Interpret model results

  • estimate: For every 1-unit increase in Lv, how much Decom increases on average. The sign of estimate indicates positive or nagative correlation between independent and dependent variables.

  • Model fit: r.squared tells you the proportion of variance in Decom explained by Lv. adj.r.squared is slightly penalized version of R² — more useful for comparing models with different numbers of predictors.

4. Multiple regression (add factor dissimilarity index)

Model formular:

Decom ~ Lv + dissim Predict soil decomposition rate with two predictors lv and dissim.

m2 <- lm(Decom ~ Lv + dissim, data = dat2) # build multiple regression model from two predictors

summary(m2)
## 
## Call:
## lm(formula = Decom ~ Lv + dissim, data = dat2)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.191693 -0.022074  0.000972  0.019407  0.145637 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.665849   0.003018 220.598  < 2e-16 ***
## Lv          -0.001760   0.001331  -1.322    0.187    
## dissim      -0.048874   0.011086  -4.408  1.4e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03905 on 340 degrees of freedom
## Multiple R-squared:  0.1907, Adjusted R-squared:  0.1859 
## F-statistic: 40.05 on 2 and 340 DF,  p-value: 2.412e-16
broom::glance(m2)[c("r.squared","adj.r.squared","AIC","BIC")]
## # A tibble: 1 × 4
##   r.squared adj.r.squared    AIC    BIC
##       <dbl>         <dbl>  <dbl>  <dbl>
## 1     0.191         0.186 -1246. -1231.

Compare m1 vs m2

Does adding another predictor dissim improve explanatory power beyond Lv?

# Anova test comparing two models
anova(m1, m2)   # nested F-test: significant p => m2 improves fit
## Analysis of Variance Table
## 
## Model 1: Decom ~ Lv
## Model 2: Decom ~ Lv + dissim
##   Res.Df     RSS Df Sum of Sq      F    Pr(>F)    
## 1    341 0.54823                                  
## 2    340 0.51859  1  0.029643 19.434 1.398e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Comparing model AIC (Akaike Information Criterion)
delta_AIC <- broom::glance(m2)["AIC"] - broom::glance(m1)["AIC"]

# Increased r2 of model_2 compared to model_1
delta_r2 <- glance(m2)$r.squared - glance(m1)$r.squared
delta_r2
## [1] 0.04626221

Interpret model comparison

  • ANOVA testing: If additional predictors in the more complex model explain a significant amount of extra variance in the outcome variable.
  • Comparing model AIC: AIC measures the trade-off between model fit and model complexity. The lower AIC, the better the model (better balance of fit and simplicity). The absolute AIC value doesn’t matter — only differences between models do. To compare AIC values, compute them for each model and pick the one with the lower AIC. If the difference (ΔAIC) is > 2, the lower-AIC model is meaningfully better.
  • delta_r2: How much r2 has been improved when adding another predictor.