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:


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  

                        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!

my data

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, 

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

  • 1
    $\begingroup$ It is telling you the effect of light controlling for the effect of location and vice versa. Your data-set seems very over-dispersed relative to what you would expect from a Poisson so you might need to investigate a negative binomial model. $\endgroup$ – mdewey Jan 27 '18 at 16:25
  • $\begingroup$ I would probably treat location as a random effect, not a fixed effect, because you aren't interested in the specific locations themselves. And, as @mdewey observes, you should probably investigate a negative binomial model (glm.nb) instead of the Poisson model to take that into account, otherwise all your standard errors will be underestimated. $\endgroup$ – jbowman Jan 27 '18 at 17:12
  • $\begingroup$ Can you add your data as described here: adv-r.had.co.nz/Reproducibility.html ? That will help people to help you. $\endgroup$ – Stefan Jan 27 '18 at 17:46
  • $\begingroup$ @Stefan I have tried to follow the instructions! thank you I hope this works! $\endgroup$ – nausicaa Jan 28 '18 at 0:32
  • $\begingroup$ @mdewey thanks ill have a read of the r book and see if I can figure glm.nb out! $\endgroup$ – nausicaa Jan 28 '18 at 0:34

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.

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

enter image description here

First the two Poisson models:

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

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)
m_nb2 <- glm.nb(Passes ~ Al.N + Location + Al.N : Location, data=Data)

#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  

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

enter image description here

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:

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.

  • $\begingroup$ thank you so much for this! you have no idea how many lecturers and different methods I've been told to use without any guidance into how! This might sound stupid but what is the negative binomial model actually doing? and how do I interpret the summary tables? unfortunately I cant find access to that book and I'm rather low on funds and time to buy it $\endgroup$ – nausicaa Jan 28 '18 at 16:05
  • $\begingroup$ @nausicaa In order to make inferences based on the sample data, we want to fit this data to a distribution that best describes the data. There are many different distribution available to model your data and the negative binomial is one of those distribution that is used when there is overdispersion in the counts since it can model the variance independently of the mean (also see my answer). BTW, if I change Passes to Buzzes the story remains very similar. I assume Buzzes represents feeding behaviour? $\endgroup$ – Stefan Jan 28 '18 at 16:17
  • $\begingroup$ @nausicaa as for output it simply gives you the estimated counts on a log scale since the default link function in the Poisson model is the log. The intercept for example represents the estimate for the forest location and no light, i.e. 3.3455. Since this is on the log scale you can back transform that value with exp(3.3455) which will be 28.375. $\endgroup$ – Stefan Jan 28 '18 at 16:28
  • $\begingroup$ oh ok! so its just bc my data best fits a negative binomial distribution? yes buzzes is feeding activity! in reference to light and dark what assumptions do you make from this? as I believed passes were greater in artificial light compared to the natural light whereas buzzes were greater in the dark locations, just from looking at basic data pooling all locations, however reading the last summary table it looks to be location specific as to which was greater however I'm not sure how to interpret this entirely $\endgroup$ – nausicaa Jan 28 '18 at 16:31
  • $\begingroup$ @nausicaa Yes, it seems like feeding activity and passes depend on the location. As to why that is, is part of you being the researcher :) I have no clue about bat activity and feeding behaviour. $\endgroup$ – Stefan Jan 28 '18 at 16:38

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.