Fitting a GAM with a count response variable across multiple years

Question

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?

Answer 1

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])