knitr::opts_chunk$set(echo = TRUE)

Load libraries

Note: it’s a good idea to load libraries in a separate code chunk, so you don’t have to mess with the setup chunk again.

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.2.1     ✔ readr     2.1.5
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ ggplot2   4.0.3     ✔ tibble    3.3.0
## ✔ lubridate 1.9.5     ✔ tidyr     1.3.2
## ✔ purrr     1.1.0     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

Function to simulate simple linear regression data

simulate_plant_growth <- function(n = 500, 
                                  alpha = 2, 
                                  beta1 = 20, 
                                  rse = 8, 
                                  moisture_range = c(0.1, 5.0)) {
  
  soil_moisture <- runif(n, min = moisture_range[1], max = moisture_range[2])
  error <- rnorm(n, mean = 8, sd = rse)
  growth_rate <- alpha + beta1 * soil_moisture + error
  
  data.frame(
    soil_moisture = soil_moisture,
    growth_rate = growth_rate
  )
}

# Run simulation
set.seed(123)
data_single <- simulate_plant_growth()

head(data_single)

##   soil_moisture growth_rate
## 1     1.5091298    37.17777
## 2     3.9626952    84.75889
## 3     2.1039869    49.32840
## 4     4.4267853    99.25968
## 5     4.7082897   116.95386
## 6     0.3232268    15.75602

Fit linear model

fit_single <- lm(growth_rate ~ soil_moisture, data = data_single)

# Extract fitted values and residuals
data_single <- data_single %>%
  mutate(
    fitted = predict(fit_single),
    resid_model = growth_rate - fitted,         # Residuals from line (SSE)
    resid_mean  = growth_rate - mean(growth_rate) # Deviations from Y_bar (SST)
  )

Visualization and Manual \(R^2\) Calculation

# Sums of squares
SST <- sum(data_single$resid_mean^2)   # Total variance
SSE <- sum(data_single$resid_model^2)  # Residual variance

R2_manual <- 1 - (SSE / SST)
R2_lm <- summary(fit_single)$r.squared

cat(sprintf("Manual R²: %.5f\nlm() R²:     %.5f\n", R2_manual, R2_lm))

## Manual R²: 0.92365
## lm() R²:     0.92365

Visualizing Residuals from the Regression Line

ggplot(data_single, aes(x = soil_moisture, y = growth_rate)) +
  geom_point(alpha = 0.6) +
  geom_smooth(method = "lm", color = "darkblue", se = FALSE) +
  geom_segment(aes(xend = soil_moisture, yend = fitted), alpha = 0.3, color = "red") +
  theme_classic() +
  labs(
    title = "Model Residuals (Vertical Distance to Regression Line)",
    x = "Soil Moisture Content", 
    y = "Growth Rate"
  )

## `geom_smooth()` using formula = 'y ~ x'

Vectorized Two-Population Simulation

simulate_two_populations <- function(n_total = 200, moisture_range = c(0.1, 2.0)) {
  n_sub <- n_total / 2
  
  # Define group-level parameters
  params <- list(
    A = list(alpha = 2.348, beta = -8, rse = 8.019),
    B = list(alpha = 6.500, beta = 20.750, rse = 5.000)
  )
  
  # Generate categorical factor and predictor values
  pop <- rep(c("A", "B"), each = n_sub)
  moisture <- runif(n_total, min = moisture_range[1], max = moisture_range[2])
  
  # Vectorized parameter mapping
  alpha <- ifelse(pop == "A", params$A$alpha, params$B$alpha)
  beta  <- ifelse(pop == "A", params$A$beta,  params$B$beta)
  rse   <- ifelse(pop == "A", params$A$rse,   params$B$rse)
  
  growth <- alpha + beta * moisture + rnorm(n_total, mean = 0, sd = rse)
  
  data.frame(
    population = factor(pop),
    soil_moisture = moisture,
    growth_rate = growth
  )
}

set.seed(42)
data_two <- simulate_two_populations()
head(data_two)

##   population soil_moisture growth_rate
## 1          A     1.8381315   -2.726511
## 2          A     1.8804433   -4.317687
## 3          A     0.6436651  -10.846051
## 4          A     1.6778505    3.748172
## 5          A     1.3193165  -13.553388
## 6          A     1.0862823   -5.496143

m0 <- lm(growth_rate ~ 1, data = data_two)

# Model 1: Ignores population (pooled model)
m1 <- lm(growth_rate ~ soil_moisture, data = data_two)

# Model 2: Population main effect (different intercepts / additive model / ANCOVA)
m2 <- lm(growth_rate ~ soil_moisture + population, data = data_two)

# Model 3: Population interaction (different intercepts + different slopes)
m3 <- lm(growth_rate ~ soil_moisture * population, data = data_two)

# Model comparison table
data.frame(
  Model = c("Soil Moisture Only", "+ Population (Separate Intercepts)", "* Population (Interaction)"),
  R_Squared = c(summary(m1)$r.squared, summary(m2)$r.squared, summary(m3)$r.squared),
  Adj_R_Squared = c(summary(m1)$adj.r.squared, summary(m2)$adj.r.squared, summary(m3)$adj.r.squared)
)

##                                Model  R_Squared Adj_R_Squared
## 1                 Soil Moisture Only 0.01862786    0.01367144
## 2 + Population (Separate Intercepts) 0.76192562    0.75950862
## 3         * Population (Interaction) 0.91197041    0.91062302

#summary(m0)
#summary(m1)
#summary(m2)
#summary(m3)



anova(m0, m1, m2, m3)

## Analysis of Variance Table
## 
## Model 1: growth_rate ~ 1
## Model 2: growth_rate ~ soil_moisture
## Model 3: growth_rate ~ soil_moisture + population
## Model 4: growth_rate ~ soil_moisture * population
##   Res.Df   RSS Df Sum of Sq        F    Pr(>F)    
## 1    199 88119                                    
## 2    198 86477  1      1641   41.475 9.012e-10 ***
## 3    197 20979  1     65499 1654.970 < 2.2e-16 ***
## 4    196  7757  1     13222  334.078 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

What could tell us Residual Sum of Squares (RSS or SSE) and Sum of Squares (Sum of Sq)?

\(R^2 = SSR/SST\) or \(R^2 = 1 - SSE/SST\)

SSE_m1 <- sum(residuals(m1)^2) # or deviance(m1)
SSE_m2 <- sum(residuals(m2)^2) # or deviance(m2)
SSE_m3 <- sum(residuals(m3)^2) # or deviance(m3)

# Calculate SSR for each model (SSR = SST - SSE)
SSR_m1 <- SST - SSE_m1
SSR_m2 <- SST - SSE_m2
SSR_m3 <- SST - SSE_m3

# Verify R^2 manually (matches summary()$r.squared exactly)
R2_m1 <- SSR_m1 / SST
R2_m2 <- SSR_m2 / SST
R2_m3 <- SSR_m3 / SST


R2_m1

## [1] 0.7954379

# Visualizing regression lines per population
data_two <- data_two %>%
  mutate(fitted = predict(m3))

# 3. Plot points, fitted parallel lines, and residual segments
ggplot(data_two, aes(x = soil_moisture, y = growth_rate, color = population)) +
  geom_point(alpha = 0.6) +
  # Draw fitted lines for model m2 (parallel slopes)
  geom_line(aes(y = fitted), linewidth = 1) + 
  # Draw residual line segments connecting each point to its fitted line
  geom_segment(aes(xend = soil_moisture, yend = fitted), alpha = 0.3) +
  theme_minimal() +
  labs(
    title = "",
    x = "Soil Moisture Content", 
    y = "Growth Rate",
    color = "Population"
  )

Key Discussion Points & ConceptsWhat does the slope (\(\beta_1\)) represent?

It represents the expected average change in plant growth rate per unit increase in soil moisture content.

Why does \(R^2\) increase when adding population?\(R^2\) increases (from \(\approx 0.708\) to \(\approx 0.733\)) because adding population accounts for baseline group differences (intercepts), reducing the total unexplained variance (\(\text{SSE}\)).

Why does the interaction term (soil.moisture:populationB) not improve \(R^2\)?

Understanding the Simulation Parameters

The simulated data follow a simple linear regression model:

\[ Y_i = \alpha + \beta X_i + \varepsilon_i \]

where:

\(Y_i\) = Plant growth rate
\(X_i\) = Soil moisture
\(\alpha\) = Intercept
\(\beta_#\) = Slope
\(\varepsilon_i\) = Random error, assumed to follow \(N(0, \sigma^2)\)

Simulation Parameters

Parameter	Meaning	What it controls	Effect on the plot
`n_total`	Sample size	Total number of observations (\(N\)).	Number of points displayed.
`alpha` (a)	Intercept	Expected growth when moisture = 0.	Moves the regression line up or down.
`beta` (b)	Slope	Change in expected growth for each unit increase in moisture.	Changes the slope and direction of the line.
`rse` (s)	Error standard deviation (\(\sigma\))	Amount of random variation around the regression line.	Controls how spread out the points are.

Variables in the Dataset

Variable	Type	Description	Effect on the plot
`pop`	Categorical	Population group (A or B).	Determines colour or grouping.
`moisture`	Numeric (X)	Soil moisture content (predictor).	Determines horizontal position.
`growth`	Numeric (Y)	Plant growth rate (response).	Determines vertical position.

Key Concepts to Remember

Intercept (\(\alpha\)): Where the regression line crosses the Y-axis.
Slope (\(\beta\)): How much the expected response changes when X increases by one unit.
Error SD (\(\sigma\)): How much individual observations vary around the expected regression line.
Sample size (\(N\)): How many observations are available for analysis.
Population (pop): Identifies groups; it only affects the fitted model if explicitly included as a predictor.

Beckerman Chapter 5

Morbey

10/28/2025