Understanding when to use a negative binomial GLMM

Question

I have 16 birds (191978,191984, 191977, 191980, 191986, 201446, 191983, 201447, 211598, 211590, 211595, 191981, 211591, 201441, 201445, 211592). There are 6 males and 10 females. The dataset is called Gbirds_sex. I have the number of times they revisit their release site (visitIdx). I also have (timeInside) column to show the residence time within the release site. I want to do statistical analysis in R to see if sex affects number of revisits (visitIdx) and residence time (timeInside). I am confused as to which test works best and the R code for it. LLM, GLM. Should I use a random intercept?

I have tried these:

#1
#visitIdx is a count variable so
# Fit a GLMM with Poisson distribution
library(lme4)
model_glmm <- glmer(visitIdx ~ sex + (1 | id), family =
poisson(), data = Gbirds_sex)
summary(model_glmm)
#2
#residenceTime is continuous - use the Gaussian family.
# Fit the GLMM with Gaussian response distribution and random      
# intercept for each individual
model_timeInside_Gaussian <- lm(timeInside ~ sex, data =       
Gbirds_sex)
summary(model_timeInside_Gaussian)

However, dispersion was high (2.8)

Should I now use:

library(glmmTMB)
model_glmm_nb <- glmmTMB(visitIdx ~ sex + (1 | id), 
   family = nbinom2(), data = Gbirds_sex)
summary(model_glmm_nb)

In relation to timeInside: Firstly, I ran a GLM:

Fit GLM for 'timeInside' with 'sex' as a predictor using Gaussian family

glm_timeInside <- glm(timeInside ~ sex, data = Gbirds_sex)
summary(glm_timeInside)

Secondly, I checked for overdispersion:

residual_deviance <- deviance(glm_timeInside)
df_residual <- df.residual(glm_timeInside)

Compute overdispersion statistic

overdispersion_statistic <- residual_deviance / df_residual
print(overdispersion_statistic)

Thirdly, as dispersion was high(7), I ran a GLMM with timeInside as the response variable, sex as the fixed effect, and id as a random effect

glmm_timeInside <- lmer(timeInside ~ sex + (1|id), data = Gbirds_sex)
summary(glmm_timeInside)

This didn't fit

boundary (singular) fit: see help('isSingular')

Is the next step a Quasi-Gaussian GLM? Is this okay that it doesn't give an AIC value ("NA")?

Answer 1

First off, the code for the second model's lm() function only fits a standard ordinary least squares (OLS) regression in R, which won't account for the clustered data in your model. If you wanted to do that, it would be better to use lmer() instead, which is simply the Gaussian equivalent to glmer().

model_timeInside_Gaussian <- lmer(timeInside ~ sex + (1|id), data = Gbirds_sex)

To your primary question about whether or not one should fit to a negative binomial GLMM using glmmTMB, I would say that this should be the default for count models, as Poisson count models almost never meet the assumption of equi-dispersion (which clearly your model seems to violate anyway).

As a side note, I noticed that you only have 16 birds in this example. Hopefully you have a lot of observations per cluster (ID), or else the parameters from this model will be somewhat untrustworthy. Mixed models are a bit data hungry, so it is important to have enough to work with in your model.

Edit

Its difficult to know how accurate the parameters will be given what you've noted in the comments about how many observations there are. Perhaps the effect isn't very noisy or is strong enough to not matter. Just be aware that you will likely get questions about this, and in the future one should perform an a priori power analysis to determine if this sort of model will be reliable. Some discussion about how many clusters or observations per cluster are necessary can be found in this answer.

In any case, there is no reason to care about dispersion with continuous data (the confusing wording of the new edit doesn't make it clear now which model you are referring to here). This only matters for discrete data, where dispersion is more important for properly modeling a Poisson-like response (either by itself or via a mixture distribution).

As for the singular matrix issue, you need to do some digging to see what is happening, most importantly to determine if there is enough random effect variance to estimate this model. If the output of your mixed model shows that this variance is near zero (visible at the beginning of the model output summary), then this typically causes mixed models to crash. A very simple check is to plot the response against the clusters or look at the caterpillar plot of the random intercepts to see if they vary enough to estimate (see a coded example in the answer here).