# Population model to model year to year dynamics

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...

• "I tried to use binomial family - logistic regression" -- why? You stated above this that it was conditionally Poisson. Wouldn't you use a Poisson model if you think it's Poisson? (I'm not saying that would be the solution to all your problems there ... it just seems a rather odd choice not to use the distribution you have in your model) – Glen_b Jun 17 '13 at 9:27
• @Glen_b you are absolutely right. It was just a thought - because the equation is equivalent - however I know it is not enough for the models to be equivalent. Anyways, I don't know how to properly evaluate the poisson model in R. It is not ordinary poisson regression, as you have the term also on the right side of the equation. – Curious Jun 17 '13 at 9:29
• Can you clarify your notation? Are you measuring population annually or monthly or...? To what degree of granularity is $R_t$ measured? If it's annual, is it subscripted by $j$? Or is it a lagged effect like $R_{j-1}$? – AdamO Jun 18 '13 at 16:00
• @AdamO, it is annually (see expl. above: "site i, year j"). $R_t$ is not measured, only $N_{i,j}$ is measured. It is just model coefficient. You are right, it should be $R_j$, thanks! – Curious Jun 18 '13 at 16:10
• Your $N_{i,j+1} \sim \text{Binom}(N_{i,j+1} + N_{i,j}, p)$ includes $N_{i,j+1}$ on both its RHS and its LHS. Is that what you meant to write? – Josh O'Brien Jun 18 '13 at 18:50

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.

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.

• Adam, I don't see your model - that's what I was asking for. Your model f <- glm(N ~ x + R, family=poisson) is not correct! – Curious Jun 18 '13 at 16:17
• You've only given us abstract definitions of your variables, I can't give you exact syntax without a description of the data. As I suggested, did you read ?glm? – AdamO Jun 18 '13 at 16:27
• Or else you can point out exactly where you're finding the difficulty. – AdamO Jun 18 '13 at 16:36
• Adam, I tried to give precise specification of my variables and the model. If something is not clear, please ask me. Anyway, I am not quite sure why you are pointing me to manual. I wonder what makes you think I am a beginner who needs help with the basics? Maybe you didn't understand why this problem is not a simple one? – Curious Jun 18 '13 at 16:43
• As for the difficulty - watch the $\mbox{log}(\mu_{i,j})$ in the predictor... – Curious Jun 18 '13 at 16:48