# Hacking GLMs to fit population growth models

I’m currently teaching Ecological Statistics and Data, a class I inherited from Lee Brown and Elizabeth Crone. In a lecture on population dynamics, they do some really cool things with generalized linear model—things that I don’t think are standard practice and as far as I can tell from googling, aren’t well documented. And let me tell you, I did a lot of googling to make sure I understood this stuff before teaching it. So, I thought I’d put it up on the blog for others.

# Data

I’ll be using data on the Northern Rocky Mountain grey wolf population. You can read more about the history of these wolves here.

library(tidyverse)
library(here)
wolves <- read_csv(here("static", "files", "NRMwolves.csv")) %>%
mutate(year_post = year - 1982)
head(wolves)
## # A tibble: 6 x 8
##    year num.wolves MT.wolves WY.wolves ID.wolves OR.wolves WA.wolves
##   <dbl>      <dbl>     <dbl>     <dbl>     <dbl>     <dbl>     <dbl>
## 1  1982          8         8        NA        NA        NA        NA
## 2  1983          6         6        NA        NA        NA        NA
## 3  1984          6         6        NA        NA        NA        NA
## 4  1985         13        13        NA        NA        NA        NA
## 5  1986         15        15        NA        NA        NA        NA
## 6  1987         10        10        NA        NA        NA        NA
## # … with 1 more variable: year_post <dbl>

# Exponential growth

Exponential growth describes unregulated reproduction and is described by the equation:

$N_{T} = \lambda^TN_0$

where $\lambda$ is the population growth rate, $T$ is a number of time steps (e.g. years) and $N_0$ is the population at some initial time.

## Hacking an exponential growth rate GLM

$log(N_T) = log(N_0) + log(\lambda)\times T$

Compare to a generic GLM equation with a log-link:

$log(y)= \beta_0 + \beta_1x_1$

Here’s the glm for an exponential growth model fit to the wolf data:

m_exp <- glm(num.wolves ~ year_post,
family = poisson(link = "log"), data = wolves)
exp(coef(m_exp))
## (Intercept)   year_post
##   19.615307    1.176583

The backtransformed intercept is the estimate for $N_0$, the estimated number of wolves at year_post = 0

The backtransformed coefficient for year_post = $\lambda$, the population growth rate

# Process error

The exponential growth model fit above is an observation error model. It assumes variation from predicted values are due to inaccuracies in estimating the number of wolves.

A process error model estimates a population growth that depends on current population size. This can be modeled as a rate, $N_{t+1}/N_{t}$.

## Hacking a process error GLM

Again, we can use a log-link to linearize this:

$log(N_{t+1}/N_{t}) = \beta_0$

$log(N_{t+1}) +log(N_{t}) = \beta_0$ $log(N_{t+1}) = \beta_0 + log(N_{t})$

The $log(N_t)$ term, which has no coefficient associated with it, is an offset. We can hack a glm to fit this model like so:

wolves2 <- wolves %>%
mutate(num.prev = lag(num.wolves)) %>% #create a column of lagged wolf numbers
filter(!is.na(num.prev))

m_process <- glm(num.wolves ~ 1, offset = log(num.prev),
family = poisson(link = "log"), data = wolves2)

The backtransformed intercept is the yearly rate of increase ($N_{t+1}/N_t$)

exp(coef(m_process))
## (Intercept)
##    1.108238

$N_{t+1} = e^{\beta_0} \times N_{t}$

So, if there are 13 wolves in 1985, how many would it predict in 1986?

1.108238 * 13
## [1] 14.40709

# Ricker model

Finally, the most complicated, possibly mind blowing example of hacking a GLM. This one took me quite a while to wrap my head around.

A Ricker model takes carrying capacity into account and allows growth rate to change as the population increases. It approximates logistic growth.

$N_{t+1}=N_{t}e^{r\left(1-{\frac {N_{t}}{K}}\right)}$ where $r = ln(\lambda)$ and $K$ is the carrying capacity

## Hacking a Ricker model GLM

$log(N_{t+1})=log(N_{t}) + log\left( (e^r)^{\left(1-{\frac {N_{t}}{K}}\right)}\right)$

$log(N_{t+1})=log(N_{t}) + log (e^r)\times{\left(1-{\frac {N_{t}}{K}}\right)}$

$log(N_{t+1})= r-\frac{r}{K}N_{t} + \textrm{offset}[log(N_{t})]$

We can model this with the following GLM:

m_rick <- glm(num.wolves ~ num.prev, offset = log(num.prev),
family = poisson, data = wolves2)

$r = \beta_0$

coef(m_rick)[1]
## (Intercept)
##   0.3119624

$\beta_1 = -r/K$

coef(m_rick)[2]
##      num.prev
## -0.0001683389

Which means that $K = -\beta_0/\beta_1$

-coef(m_rick)[1]/coef(m_rick)[2]
## (Intercept)
##     1853.18
emo::ji("exploding_head")
## 🤯