Discussion from textbook

“A series of studies [22] sampled the forest biomass in Eurasia [21]. [We study] part of that data, for small-leaved lime trees (Tilia cordata).”

“A model for the foliage biomass \(y\) is sought (Foliage). The foliage mostly grows on the outer canopy, which could be crudely approximated as a spherical shape, so one possible model is that the mean foliage biomass \(\mu\) may be related to the surface area of the approximately-spherical canopy.”

[Recall: formula for surface area \(A\) of a sphere with radius \(r\) is]

\[A = 4\pi r^2\]

“In turn, the canopy diameter may be proportional to the diameter of the tree trunk (or DBH), \(d\). This suggests a model where \(\mu\) is proportional to the surface area \(4\pi(d/2)^2 = \pi d^2\); taking logs, \(\log y\) [or \(\log(\mu)\)] is related to \(\log \pi + 2 \log d\).

“In addition, the tree diameter may be related to the age (Age) of the tree. However, since diameter measures some physical quantity and is easier to measure precisely, expect the relationship between foliage biomass and DBH to be stronger than the relationship between foliage biomass and age.”

To summarize, the discussion above suggests the model

• Foliage \(\sim Gamma(\mu,\phi)\) <- random component

• \(\log(\mu) = \beta_0 + \beta_1 \log\)(DBH) <- systematic component

Plot of Foliage vs. DBH with fit line

augment(fit, type.predict = "response") |> 
  ggplot(
    mapping = aes(x = exp(log_dbh))
  ) + 
  geom_point(aes(y = Foliage)) + 
  geom_line(aes(y = .fitted)) +
  labs(x = "DBH")

Residual plot

Plot quantile residuals on y-axis and DBH on x-axis

# set.seed(123)
augment(fit, type.predict = "response") |> 
  mutate(quantile_resid = qresid(fit)) |>
  ggplot(
    mapping = aes(x = exp(log_dbh), y = quantile_resid)
  ) + 
    geom_point() + 
    geom_hline(yintercept=0) +
  labs(x = "DBH") + 
  geom_smooth()

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

Are there outliers? Yes, we have quantile residuals of almost 4 or above

Which points are these? Two points sticking up with DBH < 10

Influential points

Are they influential? (Does removing them change the fit line by a lot?) Calculate Cook’s distances: No (see below)

augment(fit) |>
  arrange(desc(.cooksd)) |>
  relocate(.cooksd)

## # A tibble: 385 × 8
##    .cooksd Foliage log_dbh .fitted .resid    .hat .sigma .std.resid
##      <dbl>   <dbl>   <dbl>   <dbl>  <dbl>   <dbl>  <dbl>      <dbl>
##  1  0.260     3.5    2.08   -0.877  3.25  0.00530  0.638      4.33 
##  2  0.239     1.6    1.79   -1.41   2.70  0.00868  0.644      3.61 
##  3  0.0636    0.36   1.39   -2.15   1.39  0.0157   0.655      1.87 
##  4  0.0427   14.0    3.24    1.26   1.79  0.00538  0.653      2.38 
##  5  0.0258   13.6    3.36    1.48   1.38  0.00666  0.655      1.84 
##  6  0.0235    0.06   0.693  -3.43   0.688 0.0340   0.658      0.931
##  7  0.0186   14.1    3.47    1.69   1.14  0.00805  0.657      1.52 
##  8  0.0171   10.6    3.27    1.32   1.26  0.00572  0.656      1.67 
##  9  0.0142    0.01   1.19   -2.51  -1.56  0.0200   0.654     -2.10 
## 10  0.0139    0.36   1.59   -1.78   0.868 0.0119   0.658      1.16 
## # ℹ 375 more rows

Fit model without them:

fit_del <- glm(
  Foliage ~ log_dbh, family = Gamma(link="log"), 
  data = lime |> 
    filter(
      Foliage != 1.6 & DBH != 6,
      Foliage != 3.5 & DBH != 8
    )
)

Compare fit lines

augment(fit, type.predict = "response") |> 
  ggplot(
    mapping = aes(x = exp(log_dbh))
  ) + 
  geom_point(aes(y = Foliage)) + 
  geom_line(aes(y = .fitted)) + 
  geom_line( # line based on deleting two observations
    data = augment(fit_del, type.predict = "response"), 
    aes(y = .fitted), 
    color = "red"
  )