Step 1 — Explore pairwise relationships

Before fitting anything, look at how each predictor relates to ozone and to each other.

pairs(
  aq[, c("Ozone","Solar.R","Wind","Temp")],
  col   = "#4a7c5e",
  pch   = 16,
  cex   = 0.65,
  upper.panel = function(x, y, ...) {
    points(x, y, ...)
    r <- round(cor(x, y, use = "complete.obs"), 2)
    legend("topleft", legend = paste0("r = ", r),
           bty = "n", cex = 0.85, text.col = "#c8832a")
  },
  main = "Airquality: pairwise relationships"
)

Pairwise scatterplot matrix. Off-diagonal plots show each variable against every other; diagonal shows the marginal distribution. Notice ozone’s right-skewed distribution and its positive correlation with temperature.

Step 2 — Fit the multiple regression model

model_aq <- lm(Ozone ~ Solar.R + Wind + Temp, data = aq)
summary(model_aq)

## 
## Call:
## lm(formula = Ozone ~ Solar.R + Wind + Temp, data = aq)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -40.485 -14.219  -3.551  10.097  95.619 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -64.34208   23.05472  -2.791  0.00623 ** 
## Solar.R       0.05982    0.02319   2.580  0.01124 *  
## Wind         -3.33359    0.65441  -5.094 1.52e-06 ***
## Temp          1.65209    0.25353   6.516 2.42e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 21.18 on 107 degrees of freedom
## Multiple R-squared:  0.6059, Adjusted R-squared:  0.5948 
## F-statistic: 54.83 on 3 and 107 DF,  p-value: < 2.2e-16

Step 3 — Interpret the coefficients

From the model output we can write:

\[\widehat{\text{Ozone}} = \hat\beta_0 + \hat\beta_1 \cdot \text{Solar.R} + \hat\beta_2 \cdot \text{Wind} + \hat\beta_3 \cdot \text{Temp}\]

Key interpretations (sign and direction matter most):

• Solar radiation: each additional Langley of solar radiation is associated with a small increase in ozone, holding wind and temperature fixed — more photons, more photochemistry.
• Wind speed: each additional mph of wind is associated with a decrease in ozone — wind dilutes and disperses the pollutant.
• Temperature: each additional °F is associated with an increase in ozone — heat accelerates the chemical reactions forming ozone.

Step 4 — Diagnostic plots

Linear regression rests on four assumptions: linearity, constant variance (homoscedasticity), normally distributed residuals, and no influential outliers. The four standard diagnostic plots test all of these at once.

par(mfrow = c(2,2), mar = c(4,4,3,1))
plot(model_aq, col = "#4a7c5e", pch = 16, cex = 0.7)

Standard regression diagnostics. Top-left: residuals vs fitted — no strong pattern is good. Top-right: Q-Q plot — points near the line mean normally distributed residuals. Bottom-left: scale-location — flat line means constant variance. Bottom-right: leverage plot — points beyond Cook’s distance contours may be unduly influential.

par(mfrow = c(1,1))

Step 5 — Visualise partial (added-variable) relationships

Partial regression plots show the unique relationship between each predictor and the response after removing the linear effects of all other predictors. They are the honest picture of what each slope actually estimates.

library(car)
avPlots(model_aq,
        col = "#4a7c5e", pch = 16, cex = 0.7,
        col.lines = "#c8832a", lwd = 2,
        main = "Partial regression plots — Ozone ~ Solar.R + Wind + Temp")

Added-variable (partial regression) plots. Each panel shows the component of ozone unexplained by the other two predictors plotted against the component of the focal predictor unexplained by the others. The slope in each panel equals the partial regression coefficient.

Step 1 — Explore the response

ggplot(birds, aes(richness, fill = disturbance)) +
  geom_histogram(binwidth = 1, colour = "white", alpha = 0.85, position = "identity") +
  scale_fill_manual(values = c("#4a7c5e","#c8832a"), name = "Disturbance") +
  labs(x = "Bird species richness", y = "Number of sites",
       title = "Distribution of riparian bird richness") +
  theme_minimal(base_size = 13) +
  theme(
    legend.position = "top",
    plot.title = element_text(family = "serif", face = "bold", colour = "#1e3a2f")
  )

Distribution of bird species richness across 80 riparian survey sites. The right-skewed, non-negative integer distribution is characteristic of count data — unsuitable for ordinary linear regression.

Step 2 — Fit the Poisson GLM

poisson_mod <- glm(richness ~ canopy_pct + disturbance,
                   family = poisson,
                   data   = birds)
summary(poisson_mod)

## 
## Call:
## glm(formula = richness ~ canopy_pct + disturbance, family = poisson, 
##     data = birds)
## 
## Coefficients:
##                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)    0.480071   0.131667   3.646 0.000266 ***
## canopy_pct     0.022169   0.001729  12.824  < 2e-16 ***
## disturbanceLow 0.677422   0.076459   8.860  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 393.792  on 79  degrees of freedom
## Residual deviance:  72.458  on 77  degrees of freedom
## AIC: 395.56
## 
## Number of Fisher Scoring iterations: 4

Step 3 — Interpret coefficients

GLM coefficients are on the log scale (because the link is log). Exponentiating converts them to rate ratios — multiplicative effects on expected counts.

exp(coef(poisson_mod))

##    (Intercept)     canopy_pct disturbanceLow 
##       1.616189       1.022417       1.968796

exp(confint(poisson_mod))

##                   2.5 %   97.5 %
## (Intercept)    1.243314 2.083474
## canopy_pct     1.018992 1.025923
## disturbanceLow 1.697198 2.290641

Step 4 — Check for overdispersion

A Poisson model assumes the variance equals the mean. Environmental count data often have variance greater than the mean (overdispersion), which inflates Type I error if uncorrected.

# Quick check: residual deviance / residual df should be near 1
with(summary(poisson_mod), deviance / df.residual)

## [1] 0.9410136

If this ratio is substantially above 2, switch to family = quasipoisson or a negative binomial model (MASS::glm.nb()).

Step 5 — Visualise predictions

pred_grid <- expand.grid(
  canopy_pct  = seq(20, 95, length.out = 120),
  disturbance = c("Low","High")
)
pred_grid$fit <- predict(poisson_mod, newdata = pred_grid, type = "response")

# SE for ribbon
pred_se <- predict(poisson_mod, newdata = pred_grid, type = "link", se.fit = TRUE)
pred_grid$lwr <- exp(pred_se$fit - 1.96 * pred_se$se.fit)
pred_grid$upr <- exp(pred_se$fit + 1.96 * pred_se$se.fit)

ggplot() +
  geom_jitter(data = birds,
              aes(canopy_pct, richness, colour = disturbance),
              width = 0, height = 0.25, alpha = 0.45, size = 1.8) +
  geom_ribbon(data = pred_grid,
              aes(canopy_pct, ymin = lwr, ymax = upr, fill = disturbance),
              alpha = 0.18) +
  geom_line(data = pred_grid,
            aes(canopy_pct, fit, colour = disturbance), linewidth = 1) +
  scale_colour_manual(values = c("#4a7c5e","#c8832a"), name = "Disturbance") +
  scale_fill_manual(values  = c("#4a7c5e","#c8832a"), name = "Disturbance") +
  labs(x = "Canopy cover (%)", y = "Expected bird species richness",
       title = "Poisson GLM: riparian bird richness") +
  theme_minimal(base_size = 13) +
  theme(
    legend.position = "top",
    plot.title = element_text(family = "serif", face = "bold", colour = "#1e3a2f")
  )

Predicted bird richness as a function of canopy cover for low- and high-disturbance sites (Poisson GLM). Shaded bands show 95% confidence envelopes on the predicted mean. Points are observed counts, jittered slightly to reduce overplotting.

Step 1 — Visualise raw detection rates

Before modelling, it helps to bin the continuous predictor and plot the observed proportion of detections in each bin.

tortoise %>%
  mutate(shrub_bin = cut(shrub_density, breaks = 6)) %>%
  group_by(shrub_bin) %>%
  summarise(
    prop     = mean(detected),
    se       = sqrt(prop * (1 - prop) / n()),
    midpoint = mean(shrub_density),
    .groups  = "drop"
  ) %>%
  ggplot(aes(midpoint, prop)) +
  geom_errorbar(aes(ymin = pmax(prop - se, 0), ymax = pmin(prop + se, 1)),
                width = 2, colour = "#8aad96") +
  geom_point(size = 3.5, colour = "#4a7c5e") +
  geom_smooth(method = "glm",
              method.args = list(family = binomial),
              se = FALSE, colour = "#c8832a", linewidth = 1) +
  scale_y_continuous(labels = scales::percent_format(),
                     limits = c(0, 1)) +
  labs(x = "Shrub density (plants/0.1 ha)", y = "Detection probability",
       title = "Desert tortoise detection vs. shrub density") +
  theme_minimal(base_size = 13) +
  theme(plot.title = element_text(family = "serif", face = "bold", colour = "#1e3a2f"))

Observed detection rate across binned shrub-density classes. Detection frequency increases with shrub density — tortoises are more likely detected in denser shrublands — consistent with the known habitat preference. Error bars show ±1 SE.

Step 2 — Fit the logistic regression

logit_mod <- glm(detected ~ shrub_density + slope_pct,
                 family = binomial,
                 data   = tortoise)
summary(logit_mod)

## 
## Call:
## glm(formula = detected ~ shrub_density + slope_pct, family = binomial, 
##     data = tortoise)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -1.71683    0.58647  -2.927 0.003418 ** 
## shrub_density  0.04045    0.01043   3.878 0.000105 ***
## slope_pct     -0.03255    0.02203  -1.477 0.139560    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 157.72  on 119  degrees of freedom
## Residual deviance: 138.96  on 117  degrees of freedom
## AIC: 144.96
## 
## Number of Fisher Scoring iterations: 3

Step 3 — Interpret coefficients as odds ratios

exp(coef(logit_mod))

##   (Intercept) shrub_density     slope_pct 
##     0.1796353     1.0412780     0.9679758

exp(confint(logit_mod))

##                    2.5 %    97.5 %
## (Intercept)   0.05349237 0.5420687
## shrub_density 1.02113020 1.0640447
## slope_pct     0.92579727 1.0099260

Step 4 — Predicted probability surface

pred2 <- expand.grid(
  shrub_density = seq(5, 80, length.out = 80),
  slope_pct     = seq(0, 35, length.out = 60)
)
pred2$prob <- predict(logit_mod, newdata = pred2, type = "response")

ggplot(pred2, aes(shrub_density, slope_pct, fill = prob)) +
  geom_raster(interpolate = TRUE) +
  geom_contour(aes(z = prob), colour = "white", linewidth = 0.6,
               breaks = c(0.25, 0.50, 0.75)) +
  geom_point(data = tortoise,
             aes(shrub_density, slope_pct, shape = factor(detected)),
             colour = "white", size = 1.5, inherit.aes = FALSE, alpha = 0.6) +
  scale_fill_gradient2(low = "#1e3a2f", mid = "#8aad96", high = "#c8832a",
                       midpoint = 0.5, name = "P(detect)") +
  scale_shape_manual(values = c(1, 16), name = "Detected",
                     labels = c("No","Yes")) +
  labs(x = "Shrub density (plants/0.1 ha)", y = "Slope (%)",
       title = "Desert tortoise: predicted detection probability") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(family = "serif", face = "bold", colour = "#1e3a2f"))

Predicted detection probability surface across shrub density and slope values. Warm colours (orange) indicate higher predicted probability. White contour lines mark the 0.25, 0.50, and 0.75 probability thresholds.

Step 5 — Model fit for GLMs

For GLMs, R² has no direct equivalent, but McFadden’s pseudo-R² and AIC are standard:

# McFadden pseudo-R²
logit_null <- glm(detected ~ 1, family = binomial, data = tortoise)

mcfadden_r2 <- 1 - (logDeviance <- deviance(logit_mod) / deviance(logit_null))
cat("McFadden R² =", round(1 - deviance(logit_mod)/deviance(logit_null), 3), "\n")

## McFadden R² = 0.119

# AIC comparison
AIC(logit_null, logit_mod)

##            df      AIC
## logit_null  1 159.7179
## logit_mod   3 144.9648