4
$\begingroup$

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?

$\endgroup$
5
$\begingroup$

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.