poisson glm to observe whether effects of artificial light on the number of bat passes in each location were significant

Question

I am a university student looking for help with analysis for my dissertation, I do not have the best knowledge of R studios or statistics in general so please bear with me!

I would like to perform a Poisson Generalized Linear Model (GLM) to see if artificial light had an effect on the number of bat passes recorded in each of 5 locations. to do this I was going to use the code:

m1 <- glm(passes ~ location + source, family = poisson(link = "log"))

The summary of this model produces:

summary(m1) 

Call:
glm(formula = Passes ~ source + Location, family = poisson(link = "log"))

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-10.074   -5.342   -2.030    1.903   26.322  

Coefficients:
                        Estimate Std. Error z value Pr(>|z|)    
(Intercept)              4.26925    0.03290 129.778  < 2e-16  
sourcelight             -0.08672    0.02724  -3.183  0.00146  
LocationGrace Dieu wood -0.05823    0.04335  -1.343  0.17920    
LocationMartinshaw wood  0.03848    0.04231   0.910  0.36308    
LocationOld wood         0.21655    0.04058   5.336 9.51e-08  
LocationSwithland wood  -0.35303    0.04702  -7.508 6.01e-14
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 3597.4  on 79  degrees of freedom
Residual deviance: 3415.4  on 74  degrees of freedom
AIC: 3883.9

Number of Fisher Scoring iterations: 5

I have never performed a Poisson model with multiple explanatory variables before and so I'm not entirely sure what this is telling me.

I am also concerned this isn't comparing the effects of artificial light and location but both separately instead.

I've read about posthoc tests which may be more appropriate, however I can't find a source which explains this dumbed down enough for my level of understanding.

Is my original code correct or am I going wrong entirely?

For more information about my data, I have recorded the number of bat passes (so count data) in 5 different locations in an area with an artificial light source and in an area with a natural light source. Therefore source is a binary factor. I have 8 repeats for both light and dark in each location. Any help would be gratefully received! I apologise if this is not clear, I am trying!

Data <- structure(list(Location = structure(c(3L, 3L, 3L, 3L, 3L, 3L, 
3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 
5L, 5L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 
2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L, 
4L, 4L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 5L, 1L, 1L, 1L, 1L, 1L, 1L, 
1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L), .Label = c("Forest", 
"Grace", "Martinshaw", "Old", "Swithland"), class = "factor"), 
Al.N = structure(c(2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 
2L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 
2L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 
2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 
1L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 
1L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L), .Label = c("Dark", 
"light"), class = "factor"), Buzzes = c(2L, 4L, 3L, 2L, 88L, 
68L, 1L, 63L, 3L, 1L, 20L, 15L, 24L, 24L, 17L, 15L, 0L, 0L, 
2L, 6L, 8L, 3L, 5L, 7L, 6L, 2L, 2L, 2L, 7L, 4L, 5L, 4L, 7L, 
10L, 5L, 0L, 5L, 1L, 4L, 3L, 1L, 0L, 0L, 2L, 28L, 32L, 2L, 
21L, 2L, 2L, 6L, 3L, 17L, 22L, 31L, 29L, 0L, 9L, 3L, 3L, 
2L, 1L, 13L, 11L, 14L, 9L, 31L, 16L, 2L, 1L, 0L, 2L, 18L, 
29L, 22L, 3L, 18L, 2L, 15L, 6L), Passes = c(37L, 48L, 34L, 
28L, 279L, 216L, 7L, 198L, 29L, 17L, 154L, 120L, 68L, 134L, 
157L, 144L, 5L, 19L, 45L, 67L, 72L, 48L, 51L, 58L, 48L, 23L, 
25L, 20L, 39L, 25L, 23L, 34L, 53L, 57L, 48L, 26L, 57L, 25L, 
17L, 29L, 15L, 12L, 9L, 24L, 61L, 79L, 8L, 84L, 40L, 46L, 
55L, 46L, 50L, 98L, 104L, 99L, 24L, 93L, 74L, 54L, 15L, 39L, 
45L, 61L, 123L, 150L, 376L, 104L, 24L, 35L, 15L, 32L, 125L, 
156L, 107L, 47L, 142L, 51L, 59L, 35L), Date = structure(c(4L, 
19L, 35L, 13L, 2L, 21L, 34L, 15L, 39L, 16L, 33L, 10L, 37L, 
17L, 32L, 11L, 14L, 30L, 6L, 25L, 12L, 31L, 6L, 27L, 7L, 
24L, 38L, 18L, 5L, 26L, 36L, 20L, 9L, 28L, 3L, 22L, 8L, 29L, 
1L, 23L, 4L, 19L, 35L, 13L, 2L, 21L, 34L, 15L, 39L, 16L, 
33L, 10L, 37L, 17L, 32L, 11L, 14L, 30L, 6L, 25L, 12L, 31L, 
6L, 27L, 7L, 24L, 38L, 18L, 5L, 26L, 36L, 20L, 9L, 28L, 3L, 
22L, 8L, 29L, 1L, 23L), .Label = c("01/09/2017", "02/08/2017", 
"02/09/2017", "03/08/2017", "04/08/2017", "04/09/2017", "05/08/2017", 
"06/08/2017", "07/08/2017", "07/09/2017", "08/09/2017", "09/08/2017", 
"09/09/2017", "10/08/2017", "10/09/2017", "11/08/2017", "12/08/2017", 
"12/09/2017", "13/08/2017", "13/09/2017", "14/08/2017", "14/09/2017", 
"15/09/2017", "16/08/2017", "16/09/2017", "17/08/2017", "17/09/2017", 
"18/08/2017", "19/08/2017", "21/08/2017", "22/08/2017", "25/08/2017", 
"26/08/2017", "27/08/2017", "28/08/2017", "29/08/2017", "30/07/2017", 
"30/08/2017", "31/07/2017"), class = "factor")), .Names = c("Location", 
"Al.N", "Buzzes", "Passes", "Date"), class = "data.frame", row.names = c(NA, 
-80L))

I've tried to follow the instructions, apologies if this is incorrect!

Answer 1

After looking at your raw data, you probably also want a model that includes the interaction of Al.N : Location. As you can see below, Passes are different among Location.

require(ggplot2)
ggplot(Data, aes(Al.N, Passes)) + geom_point() + facet_wrap(~ Location)

First the two Poisson models:

m_p1 <- glm(Passes ~ Al.N + Location, family = poisson, data = Data)
summary(m_p1)
m_p2 <- glm(Passes ~ Al.N + Location + Al.N : Location, family = poisson, data=Data)
summary(m_p2)

By looking at the AIC, it can also be seen that m_p2 is better (the smaller the better):

AIC(m_p1, m_p2)
#     df      AIC
#m_p1  6 3883.873
#m_p2 10 2871.577

One property of the Poisson distribution is that the mean $\lambda$ equals the variance, which means that there is no additional parameter available to model the variance independently from the mean (such as $\sigma$ in the normal distribution). So if the the variance in the model is larger then $\lambda$, we have overdispersion, and the model will give you biased standard errors. To check how the dispersion of the fitted models compares to 1, I adapted the following code from here (please read the caveats of this approach in the link): http://bbolker.github.io/mixedmodels-misc/glmmFAQ.html#overdispersion

df_resid <- nrow(model.frame(m_p1)) - length(coef(m_p1))
pearson_resid <- residuals(m_p1, type = "pearson")
pearson_sq <- sum(pearson_resid^2)
pearson_sq / df_resid 
#[1] 58.84301 

df_resid <- nrow(model.frame(m_p2)) - length(coef(m_p2))
pearson_resid <- residuals(m_p2, type = "pearson")
pearson_sq <- sum(pearson_resid^2)
pearson_sq / df_resid 
#[1] 35.96936

Although adding the interaction to the model reduces overdispersion, it is still way too high compared to the expected 1.

As @mdewey and @jbowman already pointed out, a negative binomial model may be more appropriate here. It has one parameter more than the Poisson distribution and hence allows to adjust the extra variance independently from the mean. To do so, you can use the glm.nb() function from the MASS package.

require(MASS) # needed for glm.nb() function
# same two models
m_nb1 <- glm.nb(Passes ~ Al.N + Location, data=Data)
summary(m_nb1)
m_nb2 <- glm.nb(Passes ~ Al.N + Location + Al.N : Location, data=Data)
summary(m_nb2)

#Call:
#glm.nb(formula = Passes ~ Al.N + Location + Al.N:Location, data = Data, 
#    init.theta = 2.453412501, link = log)

#Deviance Residuals: 
#    Min       1Q   Median       3Q      Max  
#-2.9174  -0.8198  -0.1160   0.4445   2.4269  

#Coefficients:
#                             Estimate Std. Error z value Pr(>|z|)    
#(Intercept)                    3.3455     0.2353  14.220  < 2e-16 ***
#Al.Nlight                      1.3424     0.3278   4.095 4.22e-05 ***
#LocationGrace                  0.6033     0.3297   1.830  0.06727 .  
#LocationMartinshaw             1.4124     0.3277   4.310 1.63e-05 ***
#LocationOld                    1.3250     0.3278   4.042 5.31e-05 ***
#LocationSwithland              0.5386     0.3300   1.632  0.10260    
#Al.Nlight:LocationGrace       -0.9426     0.4619  -2.040  0.04130 *  
#Al.Nlight:LocationMartinshaw  -2.8470     0.4640  -6.136 8.44e-10 ***
#Al.Nlight:LocationOld         -1.8638     0.4610  -4.043 5.27e-05 ***
#Al.Nlight:LocationSwithland   -1.3632     0.4632  -2.943  0.00325 ** 
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

#(Dispersion parameter for Negative Binomial(2.4534) family taken to be 1)

#    Null deviance: 130.287  on 79  degrees of freedom
#Residual deviance:  84.657  on 70  degrees of freedom
#AIC: 810.93

#Number of Fisher Scoring iterations: 1


#              Theta:  2.453 
#          Std. Err.:  0.384 

# 2 x log-likelihood:  -788.930

Checking AIC:

AIC(m_nb1, m_nb2)
#      df      AIC
#m_nb1  7 835.5290
#m_nb2 11 810.9303

Now dispersion:

df_resid <- nrow(model.frame(m_nb1)) - length(coef(m_nb1)+1)
pearson_resid <- residuals(m_nb1, type = "pearson")
pearson_sq <- sum(pearson_resid^2)
pearson_sq / df_resid 
#[1] 1.490963

df_resid <- nrow(model.frame(m_nb2)) - length(coef(m_nb2)+1)
pearson_resid <- residuals(m_nb2, type = "pearson")
pearson_sq <- sum(pearson_resid^2)
pearson_sq / df_resid 
#[1] 1.055097

Check the residuals for non-random patterns:

plot(pearson_resid ~ fitted(m_nb2))

That looks good, so m_nb2 looks like a decent model for your data.

To get mean estimates and 95% confidence intervals, you could do this via the emmeans() function, in the emmeans package:

require(emmeans)
emmeans(m_nb2, ~ Al.N * Location, type = "response")
# Al.N  Location   response        SE  df asymp.LCL asymp.UCL
# Dark  Forest       28.375  6.675952 Inf  17.89240  44.99902
# light Forest      108.625 24.794163 Inf  69.44469 169.91061
# Dark  Grace        51.875 11.982910 Inf  32.98634  81.57971
# light Grace        77.375 17.739803 Inf  49.36797 121.27075
# Dark  Martinshaw  116.500 26.571810 Inf  74.50412 182.16777
# light Martinshaw   25.875  6.111123 Inf  16.28714  41.10702
# Dark  Old         106.750 24.370911 Inf  68.24007 166.99224
# light Old          63.375 14.579257 Inf  40.37395  99.47976
# Dark  Swithland    48.625 11.249112 Inf  30.89862  76.52091
# light Swithland    47.625 11.023322 Inf  30.25626  74.96435

# Confidence level used: 0.95 
# Intervals are back-transformed from the log scale

So it definitely looks like that the effect of the light source varies among the different locations with respect to passes of bats.

If you want all pairwise comparisoins between your counts, then run:

emmeans(m_nb2, pairwise ~ Al.N * Location, type = "response")

You should also have a look at the Zuur et al. books in the "A Beginner's Guide to ..." series. They have lots of biological datasets and cover pretty much all possible analyses.