# Selecting error distribution and link function for GLMMs on count data

I would like to find out how to select an appropriate error distribution and link function to model count data to determine a treatment effect using a generalized linear mixed effects model.

I have mostly followed this practical guide, which I found helpful, but I still have some questions.

Q1: In the particular case presented below, I would use the model with the negative binomial error distribution, because a) the data seems to be too overdispersed for a poisson distribution, b) the deviance is lowest using the negative binomial distribution and c) the predictions match well the observed mean values of both treatment groups. Is my choice reasonable?

Q2: What is the link function that the glmer.nb function (for negative binomial distributions) uses? I did not find any info on this in the documentation.

Q3: In the models with gamma distribution, the link function does not affect the predictions or the deviance, but the p-value. What criteria would I use to choose the correct link function?

Q4: Why do gamma models predict the same value for all blocks?

Q5: My colleague suggested that I use for count data either a poisson or a negative binomial distribution irrespective of whether a model with gamma distribution has a lower deviance because the gamma distribution is for continuous data. Is that right?

Q6: How could I graphically examine model residuals of models using a non-Gaussian distribution (similarly to the qqp plots below)?

I would be already very happy if someone would answer Q 1-2.

my_data = structure(list(treatment = structure(c(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), .Label = c("a", "b"), class = "factor"),
counts = c(118L, 121L, 133L, 28L, 1L, 1L, 1L, 1L, 15L, 43L,
33L, 14L, 55L, 2L, 1L, 7L, 29L, 187L, 228L, 208L, 20L, 1L,
4L, 1L, 47L, 73L, 50L, 62L, 0L, 4L, 135L, 1L, 126L, 2L, 25L,
39L, 2L, 143L, 2L, 1L, 2L, 53L, 345L, 12L, 1L, 2L, 133L,
1L, 95L, 85L, 139L, 91L, 10L, 1L, 1L, 29L, 133L, 358L, 181L,
27L, 8L, 141L, 92L, 6L), site = c("S01", "S01", "S01", "S01",
"S01", "S01", "S01", "S01", "S02", "S02", "S02", "S02", "S02",
"S02", "S02", "S02", "S03", "S03", "S03", "S03", "S03", "S03",
"S03", "S03", "S04", "S04", "S04", "S04", "S04", "S04", "S04",
"S04", "S05", "S05", "S05", "S05", "S05", "S05", "S05", "S05",
"S10", "S10", "S10", "S10", "S10", "S10", "S10", "S10", "S11",
"S11", "S11", "S11", "S11", "S11", "S11", "S11", "S12", "S12",
"S12", "S12", "S12", "S12", "S12", "S12"), block = c("B07",
"B07", "B07", "B07", "B17", "B17", "B17", "B17", "B02", "B02",
"B02", "B02", "B03", "B03", "B03", "B03", "B14", "B14", "B14",
"B14", "B18", "B18", "B18", "B18", "B09", "B09", "B09", "B09",
"B23", "B23", "B23", "B23", "B13", "B13", "B13", "B13", "B12",
"B12", "B12", "B12", "B16", "B16", "B16", "B16", "B04", "B04",
"B04", "B04", "B06", "B06", "B06", "B06", "B05", "B05", "B05",
"B05", "B21", "B21", "B21", "B21", "B20", "B20", "B20", "B20"
)), .Names = c("treatment", "counts", "site", "block"), row.names = c(NA,
64L), class = "data.frame")


Code

library(MASS)
library(car)
library(lme4)
library(lmerTest)
library(blmeco)
my_data_a = subset(my_data, treatment == "a")
my_data_b = subset(my_data, treatment == "b")

nd = data.frame(treatment = c("a","b"))

# Checking means and variances
mean_counts_a = mean(my_data_a$counts) mean_counts_b = mean(my_data_b$counts)

var(my_data_a$counts) var(my_data_b$counts)

my_data$counts.t = my_data$counts +1

qqp(log(my_data$counts.t), "norm") poisson <- fitdistr(my_data$counts.t, "Poisson")
qqp(my_data$counts.t, "pois", poisson$estimate)

gamma <- fitdistr(my_data$counts.t, "gamma") qqp(my_data$counts.t, "gamma", shape = gamma$estimate[[1]], rate = gamma$estimate[[2]])

nbinom <- fitdistr(my_data$counts.t, "Negative Binomial") qqp(my_data$counts.t, "nbinom", size = nbinom$estimate[[1]], mu = nbinom$estimate[[2]])

# Fitting models
lmer_counts_gaussian = lmer(log(counts+1) ~ treatment + (1|site/block), data = my_data)
summary(lmer_counts_gaussian)
plot(fitted(lmer_counts_gaussian), residuals(lmer_counts_gaussian), xlab = "Fitted Values", ylab = "Residuals")

glmer_counts_gamma = glmer(counts+1 ~ treatment + (1|site/block) , family = "Gamma"(link = "inverse"), data = my_data)
summary(glmer_counts_gamma) # deviance = 624.4
plot(fitted(glmer_counts_gamma), residuals(glmer_counts_gamma), xlab = "Fitted Values", ylab = "Residuals")
dispersion_glmer(glmer_counts_gamma) # just below the threshold of 1.4

glmer_counts_gamma_log = glmer(counts+1 ~ treatment + (1|site/block) , family = "Gamma"(link = "log"), data = my_data)
summary(glmer_counts_gamma_log) # deviance = 624.4
plot(fitted(glmer_counts_gamma_log), residuals(glmer_counts_gamma_log), xlab = "Fitted Values", ylab = "Residuals")
dispersion_glmer(glmer_counts_gamma_log) # just below the threshold of 1.4

glmer_counts_poisson = glmer(counts+1 ~ treatment + (1|site/block) , family = "poisson"(link = "log"), data = my_data)
summary(glmer_counts_poisson) # deviance = 3348.5; p-values is a bit larger
plot(fitted(glmer_counts_poisson), residuals(glmer_counts_poisson), xlab = "Fitted Values", ylab = "Residuals")
dispersion_glmer(glmer_counts_poisson) # clearly above 1.4

glmer_counts_neg_binomial = glmer.nb(counts+1 ~ treatment + (1|site/block), data = my_data)
summary(glmer_counts_neg_binomial) # deviance = 623.7
plot(fitted(glmer_counts_neg_binomial), residuals(glmer_counts_neg_binomial), xlab = "Fitted Values", ylab = "Residuals")
dispersion_glmer(glmer_counts_neg_binomial) # below 1.4

counts_means = cbind(nd, mean_real = c(mean_counts_a, mean_counts_b),
mean_gaussian = exp(predict(lmer_counts, newdata=nd, type="response", re.form=NA)),
mean_gamma = predict(glmer_counts_gamma, newdata=nd, type="response", re.form=NA),
mean_gamma_log = predict(glmer_counts_gamma_log, newdata=nd, type="response", re.form=NA),
mean_poisson = predict(glmer_counts_poisson, newdata=nd, type="response", re.form=NA),
mean_neg_binomial = predict(glmer_counts_neg_binomial, newdata=nd, type="response", re.form=NA)
)

# gamma distribution predicts mean values similar to the real means irrespective of the link function

my_data$counts_gaussian = exp(predict(lmer_counts_gaussian, newdata=my_data, type="response")) my_data$counts_gamma = predict(glmer_counts_gamma, newdata=my_data, type="response")
my_data$counts_gamma_log = predict(glmer_counts_gamma_log, newdata=my_data, type="response") my_data$counts_poisson = predict(glmer_counts_poisson, newdata=my_data, type="response")
my_data\$counts_neg_binomial = predict(glmer_counts_neg_binomial, newdata=my_data, type="response")

View(my_data)
# models with gamma distribution does not predict different values for different blocks