# Fitting a GAM with a count response variable across multiple years

I am trying to fit a GAM where the response variable is a count. For my particular problem, I have counts for a species, across a number of different years. So it would seem to me that the best thing to do would be to generate a mean across the three years, and use this as the response variable. However this doesn't work when the GAM is configured to use a poisson distribution; presumably as this distribution only works with the integers?

Reproducible code as follows:

set.seed(0)
dat <- data.frame(count_y1 = rpois(100,1),
count_y2 = rpois(100,1),
count_y3 = rpois(100,1),
pred1 = rnorm(100, 10, 1),
pred2 = rnorm(100, 0, 1),
pred3 = rnorm(100, 0, 1),
pred4 = rnorm(100, 0, 1))

library(dplyr)
dat <- dat %>%
rowwise() %>%
mutate(count_mn = mean(c(count_y1, count_y2, count_y3)))


If I set up the GAM as follows:

model <- gam(count_mn ~ s(pred1) + s(pred2) + s(pred3) + offset(pred4),
data = dat,
family = poisson(link="log"),
method = "REML",
select = TRUE)


(in my real world problem I have an offset that handles a variable survey area for where the species count was taken). This gives the following output:

Error in if (abs(old.score - score) > score.scale * conv.tol) { :
missing value where TRUE/FALSE needed
In addition: There were 50 or more warnings (use warnings() to see the first 50)


If I now assign quasipoisson to family, which I understand is used for zero-inflated and continuous response data, I get a model fit:

Family: quasipoisson
Link function: log

Formula:
count_mn ~ s(pred1) + s(pred2) + s(pred3) + offset(pred4)

Estimated degrees of freedom:
1.640 0.102 1.754  total = 4.5

REML score: 54.81147


So my question would be, is it possible to fit the model across the three count years , using a poisson distribution?

## 1 Answer

You could use a sum across years, not a mean. The sum of Poisson variables is Poisson; the mean isn't.

set.seed(0)
dat <- data.frame(count_y1 = rpois(100,1),
count_y2 = rpois(100,1),
count_y3 = rpois(100,1),
pred1 = rnorm(100, 10, 1),
pred2 = rnorm(100, 0, 1),
pred3 = rnorm(100, 0, 1),
pred4 = rnorm(100, 0, 1))

library(dplyr)
dat <- dat %>%
rowwise() %>%
mutate(count_sm = sum(c(count_y1, count_y2, count_y3)))


and now your gam fits fine

model <- gam(count_sm ~ s(pred1) + s(pred2) + s(pred3) + offset(pred4),
data = dat,
family = poisson(link="log"),
method = "REML",
select = TRUE)
summary(model)


You could also add another offset equal to log(3) so that the rest of the linear predictor is now a model for log(E[sum]/3) rather than log(E[sum])