Total: 35 points

knitr::opts_chunk$set(error=FALSE, fig.width=6, fig.height=4)
library(tidyverse)
library(GLMsData)
library(broom)
library(statmod)
library(metR)
library(cowplot)
library(gginnards)

Data Introduction

This homework will focus on the ants dataset from the textbook “GLMs with examples in R”. Its description for Problem 10.14 on pages 417-418 is as follows:

A study of the species richness (the number of species) of ants at 22 sites in the New England region, USA, examined relationships with habitat (forest or bog), elevation (in m) and latitude. An important feature of the dataset is that there are paired observations at each site (as can be identified with a unique Site or Latitude/Elevation combination), with one observation for in the forest habitat and another in the bog habitat.

The variables are:

The chunk below will load the ants dataset into the current R session.

data(ants)

Part 1: Exploratory Data Analysis

  1. [2] Make a scatterplot with Latitude on the x-axis and Elevation on the y-axis.
ggplot(data = ants,
     aes(x = Latitude, y = Elevation)) +
geom_point()

  1. [2] Describe the scatterplot.

Answer: The scatterplot shows how elevation varies with latitude across the sites. At lower latitudes (around 42°), elevations range widely from near sea level to about 500 m. In contrast, at higher latitudes (close to 45°), sites are mostly concentrated at mid-range elevations (about 200–400 m). Although there isn’t a clear linear trend, higher elevations tend to occur more often at both the lowest and highest latitudes, while mid-latitudes (around 43°) have fewer sites and a narrower spread of elevations.

Make scatterplots with species richness on the y-axis and different colored points depending on the Habitat as follows:

  1. [2] with Latitude on the x-axis
ggplot(data = ants,
aes(x = Latitude, y = Srich, color = Habitat)) +
  geom_point()

  1. [2] with Elevation on the x-axis
ggplot(data = ants,
aes(x = Elevation, y = Srich, color = Habitat)) +
  geom_point()

  1. [3] Based on the scatterplots in Questions #2 and #3, describe the effects of Habitat, Latitude, and Elevation on species richness.

Answer: The scatterplots suggest that habitat type significantly affects species richness, with forests generally hosting more ant species than bogs at the same latitude. Richness also tends to decrease with latitude, as northern sites contain fewer species. Likewise, higher elevations are linked to lower species richness, whereas lower elevations typically support more species overall.

Fit a Poisson regression model

  1. [2] Write R code that fits the following Poisson regression model to the data:

\[Srich \sim Poisson(\mu)\] \[\log(\mu) = \beta_0 + \beta_1 Latitude + \beta_2 Elevation + \beta_3 I(Habitat=Forest)\]
where \(I(\cdot)\) is the indicator function such that \(I(A)=1\) if \(A\) is true and \(I(A)=0\) if \(A\) is false.

poisson_model <- glm(Srich ~ Latitude + Elevation + Habitat,
data = ants, family = poisson)
summary(poisson_model)
## 
## Call:
## glm(formula = Srich ~ Latitude + Elevation + Habitat, family = poisson, 
##     data = ants)
## 
## Coefficients:
##                 Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   11.9368121  2.6214970   4.553 5.28e-06 ***
## Latitude      -0.2357930  0.0616638  -3.824 0.000131 ***
## Elevation     -0.0011411  0.0003749  -3.044 0.002337 ** 
## HabitatForest  0.6354389  0.1195664   5.315 1.07e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 102.76  on 43  degrees of freedom
## Residual deviance:  40.69  on 40  degrees of freedom
## AIC: 209.04
## 
## Number of Fisher Scoring iterations: 4
b <- coef(poisson_model)
cat(sprintf(
"log(mu) = %.4f + %.4f*Latitude + %.4f*Elevation + %.4f*I(Habitat=Forest)\n",
b[1], b[2], b[3], b[4]
))
## log(mu) = 11.9368 + -0.2358*Latitude + -0.0011*Elevation + 0.6354*I(Habitat=Forest)
exp(coef(poisson_model))
##   (Intercept)      Latitude     Elevation HabitatForest 
##  1.527888e+05  7.899442e-01  9.988596e-01  1.887850e+00
  1. [2] Display a table of the estimated beta coefficients.
tidy(poisson_model)
## # A tibble: 4 × 5
##   term          estimate std.error statistic     p.value
##   <chr>            <dbl>     <dbl>     <dbl>       <dbl>
## 1 (Intercept)   11.9      2.62          4.55 0.00000528 
## 2 Latitude      -0.236    0.0617       -3.82 0.000131   
## 3 Elevation     -0.00114  0.000375     -3.04 0.00234    
## 4 HabitatForest  0.635    0.120         5.31 0.000000107
  1. [2] Write an equation that relates the estimated mean to the variables Latitude, Elevation, and Habitat.

Answer: The mean species richness is estimated using the equation \[ \hat{\mu} = \exp\left( 11.94 - 0.236 \times \text{Latitude} - 0.00114 \times \text{Elevation} + 0.635 \times I(\text{Habitat} = \text{Forest}) \right) \] In this model, latitude is in degrees and elevation is in meters. The term 𝐼(Habitat = Forest) is an indicator that equals 1 for forest sites and 0 for bogs. According to the coefficients, species richness tends to decrease with higher latitude and elevation, but forest sites are expected to support more species than bogs.

Interpret coefficients of the estimated model

  1. [3] Interpret the coefficient for Latitude.

Answer: The estimated coefficient for Latitude is –0.2358. This tells us that moving one degree farther north is associated with a 21% decline in expected species richness, after controlling for elevation and habitat. \[ e^{-0.2358} = 0.79 \] So, for each one-degree increase in latitude, the expected number of species is about 79% of the previous value, which is roughly a 21% decline.

exp(coef(poisson_model)["Latitude"])
##  Latitude 
## 0.7899442
  1. [3] Interpret the coefficient of Elevation in terms of changes of 100 m in elevation.

Answer: The estimated coefficient for Elevation is –0.0011. This implies that for every 1-meter rise in elevation, the log of expected species richness decreases slightly, holding latitude and habitat constant. For a 100-meter increase: \[ e^{-0.0011 \times 100} = e^{-0.11} \approx 0.90 \] This means species richness drops by about 10% for every 100-meter gain in elevation.

exp(coef(poisson_model)["Elevation"] * 100)
## Elevation 
## 0.8921617
  1. [3] Interpret the coefficient of \(I(Habitat=Forest)\).

Answer: The mean number of species in forest habitats is 1.89 times the mean number of species in bog habitats, when Latitude and Elevation are held constant.

exp(coef(poisson_model)["Elevation"] * 100)
## Elevation 
## 0.8921617

Visualizing the fitted model with a contour plot

As we did previously in LM 2, we can Visualize the fitted mean values over two explanatory variables with a contour plot.

The contours on the plots below (as shown by the color scale and the black numeric labels) show the fitted mean counts of species richness for different values Latitude, Elevation, and Habitat.

The white numbers show the observed counts.

fit <- glm(Srich ~ Latitude + Elevation + Habitat,
data = ants,
family = poisson)

For ease of reading the contour plots, let’s remove the observed counts.

Now, the white points on the plots below show the location of Latitude = 43 degrees and Elevation = 200 m. (Note: this is not a site in the dataset.)

  1. [3] Use the augment() function to find the fitted mean count of ants species at Latitude = 43 degrees and Elevation = 200 m for (a) Habitat = Forest and (b) Habitat = Bog. (Hint: You will need to supply these x values in the newdata argument.)
# for Habitat = Forest
augment(
poisson_model,
type.predict = "response",
newdata = tibble(Latitude = 43, Elevation = 200, Habitat = "Forest")
)
## # A tibble: 1 × 4
##   Latitude Elevation Habitat .fitted
##      <dbl>     <dbl> <chr>     <dbl>
## 1       43       200 Forest     9.07
# for Habitat = Bog
augment(
poisson_model,
type.predict = "response",
newdata = tibble(Latitude = 43, Elevation = 200, Habitat = "Forest")
)
## # A tibble: 1 × 4
##   Latitude Elevation Habitat .fitted
##      <dbl>     <dbl> <chr>     <dbl>
## 1       43       200 Forest     9.07
  1. [3] Explain how the fitted mean counts calculated in the previous problem is connected with coefficient estimate for \(I(Habitat=Forest)\).

Answer: The expected mean counts for Forest and Bog habitats are related through the coefficient of 𝐼(Habitat = Forest). According to the model: \(\hat{\mu}_{Forest} = e^{\beta_3} \hat{\mu}_{Bog}\)

# Compare fitted values for Forest vs. Bog habitats
fit_forest <- augment(
poisson_model,
newdata = tibble(Latitude = 43, Elevation = 200, Habitat = "Forest"),
type.predict = "response"
)
fit_bog <- augment(
poisson_model,
newdata = tibble(Latitude = 43, Elevation = 200, Habitat = "Bog"),
type.predict = "response"
)
# Ratio of predicted means
fit_forest$.fitted / fit_bog$.fitted
## [1] 1.88785
# Exponentiated coefficient for habitat effect
exp(coef(poisson_model)[["HabitatForest"]])
## [1] 1.88785

Plugging in the estimated coefficient:

\(\hat{\mu}_{Forest} = 1.88785 \times \hat{\mu}_{Bog}\)

This means that, for the same latitude and elevation, the expected species richness in forests is about 1.89 times that in bog habitats. 14. [3] Recall that the contours on the plots in Activity 2.3 were evenly spaced – i.e., each change in x or y produced the same change in the contour value. Explain why that is not the case for the contours shown here.

Answer:In Activity 2.3, the Multiple Linear Regression model produced a linear relationship, so equal changes in the predictors resulted in equal changes in the fitted values, with evenly spaced contour lines. In this homework, the Poisson regression with a log link created a nonlinear relationship. Changes in the predictors affected the mean multiplicatively rather than additively, leading to unevenly spaced contour lines.