Population model to model year to year dynamics

Question

My task is to assess how various environmental variables affect annual population fluctuations. For this, I would use a model like: $$ \mbox{log} ( \mu_{i,j+1} ) = \mbox{log} ( \mu_{i,j} ) + R_{j} + \sum\limits_{k} \alpha_k x_{k,j} \\ N_{i,j} \sim \mbox{Poiss} ( \mu_{i,j} ) $$ Where $N_{i,j}$ is number of observed individuals at site i in year j, $\mu_{i,j}$ is the expected number of individuals at site i in year j, $x_{k,j}$ accross all $k$ is vector of environmental variables in year $j$, $\alpha_k$ are coefficients and $R_{j}$ is model coefficient meant to handle "background" population growth (it is not measured, this is just a model coefficient). However I think I will remove the year index $j$ and use it as intercept only, so that it doesn't hide possible global effect of environmental variables.

This was the simplest version of the model - in the next stage I would like to handle overdispersion (not sure how yet) and maybe add some per-site random effect.

Questions:

1. How can I fit this model in R? I can write model in WinBUGS but prefer to have "frequentist" solution in R, because it is much faster and inference is easier (one has p-values, t-tests, F-tests...). But which function or package use to fit it? I don't think this can be implemented using GLM! I spotted that my equation can be converted to:

   $$ \log\left({\mu_{i,j+1} \over \mu_{i,j}}\right) = R_{j} + \sum\limits_{k} \alpha_k x_k$$ Which resembles logistic regression: $$\log\left({\mu_{i,j+1} \over \mu_{i,j}}\right) = \text{logit}\left({\mu_{i,j+1} \over { \mu_{i,j}} + \mu_{i,j+1}}\right) = R_{j} + \sum\limits_{k} \alpha_k x_k$$ $$N_{i,j+1} \sim \text{Binom}\left(N_{i,j} + N_{i,j+1}, p = {\mu_{i,j+1} \over { \mu_{i,j}} + \mu_{i,j+1}}\right)$$ However, I am not sure this gives equivalent result; this converts the Poisson counts to Binomial; and possibly, it would not be quite straightforward to handle overdispersion (the Poisson overdispersion for animal counts is well covered and published; it is not clear how would it work in the binomial version). For this reason, I prefer to compute the original model as it is (Poisson).

2. How to incorporate a negative density dependence (i.e. the population growth is lower where there is a lot of individuals)? Add something like $\beta * \ln(\mu_{i,j})$ to the right side? Seems little strange to me...

Answer 1

Baker 2012 (Journal of Applied Ecology) used similar model. I was asking him and he replied he uses normal glm()! He inspired me to use the following transformation (that he actually used in the linked article) - just recursively substitute the $\mbox{log} (\mu_{i,j})$, until you get this:

$$\mbox{log} ( \mu_{i,j+1} ) = \mbox{log} ( \mu_{i,1} ) + \sum\limits_{t=1}^{j} R_{t} + \sum\limits_{k} \alpha_k \sum\limits_{t=1}^{j}x_{k,t}$$

and then, $\mbox{log}(\mu_{i,1})$ can be simply taken as a site-specific intercept:

$$\begin{eqnarray} \mbox{log} ( \mu_{i,j+1} ) &=& \alpha_i + \sum\limits_{t=1}^{j} R_{t} + \sum\limits_{k} \alpha_k \sum\limits_{t=1}^{j}x_{k,t} \\ \mbox{log} ( \mu_{i,1} ) &=& \alpha_i \end{eqnarray}$$

so this can be easily solved by classic GLM. It is trivial to see that the transformed equations are equivalent to the original model. I do not trivially see that the whole fit proccess including poisson errors will also be equivalent, but this is probably more limitation of my brain than an actual problem :).

The transformed model is of course very easily fitted using glm()! Including overdispersion using the quasipoisson family.

Answer 2

Population growth models often use Poisson modeling framework. In R, fitting a Poisson GLM is easy. See ?glm. An example is:

f <- glm(N ~ x + R, family=poisson)

To estimate: $\log(\mu | x) = \beta_0 + \beta_1 x + \beta_2 R$. $\beta_1$ is interpreted as a relative rate comparing the rate (or Poisson lambda) for $N$ differing by 1 unit in $x$. This may be desirable when $x$ is an experimental condition with controlled values, or completely unconditional. I would agree it doesn't immediately suggest to me a good model for $x$ being time, because inferring the previous year's population from the current is much different than inferring the next year's population from the last. Let's hold on to the question of whether a fixed effect is adequate for modeling baseline time trends.

Intuitively, we have some notion that there is correlation in these data since population growth is exponential. In horseshoe crabs, for instance, higher population means more mingling and higher fertility rates (see Agresti Categorical Data Analysis 2nd ed). This motivates the use of a quasilikelihood, or quasipoisson model to account for overdispersion if our interest lies specifically in measuring the relative rates for some exposure accounting for time.

However, neither of these approaches really answer the question we're truly interested in: "What is the relative rate for $N$ comparing values of $x$ differing by 1 unit for a standardized population?". That is, for every thousand subjects, say, if we were to observe $x$ taking a different value, how would that affect seasonal trends? The way to account for this is the use of an offset.

Offsets are constrained parameters in the model that can account for lagged effects in time and standardize a denominator of growth. For instance, consider the following linear model:

$\log \left( N_{i} / N_{i-1} \right) = \beta_0 + \beta_1 x$

$N_{i-1}$ (the previous year's rates) are considered fixed and known. However, if $x$ is correlated with $N_{i-1}$, adjusting for the lagged effect in the linear model will lead to biased estimates of $\beta_1$. So we rewrite the model as:

$\log N_{i} = \beta_0 + \beta_1 x + 1 \cdot \log N_{i-1}$.

This allows us to account for the previous years' volume while estimating the association between $x$ and the current years' incremental growth.