I'd like to illustrate the importance of accounting for the dependency between mean and variance in inference with linear models. Is my example below a good one? Do you agree with my comments on it?

Here I simulated count data from a negative binomial distribution for 4 groups with different means and same dispersion parameter. I fitted the same model using ordinary linear regression (lm) and a generalised linear model with negative binomial error distribution (glm). One can see that:

  • The ordinary linear model returns the same error estimates for each of the four groups (and therefore the same width of confidence intervals, CI). This is because lm assumes constant variance. The glm instead adapts the error to the mean. (Note that the error estimates depend also on the group size, which is the same for all groups here).

  • Consequently, lm underestimates the error in groups with large mean (i.e. A) and overestimates in groups with small mean (C and D). In fact, the 95% CI for group A just misses the true of value of 50 and visually the CI for C and D doesn't feel right.

  • The CI from lm contain negative values which are unphysical for count data.

enter image description here

So, how does glm adapt the variance (and therefore error estimates and CI) to the mean? This is in the choice of error distribution. For count data we could choose Poisson distribution were the variance is equal to the mean. Poisson often underestimates the variance since you would assume that different observations are identical in everything and the only source of variation is due to random sampling (the urn model). If you count parasite eggs from different animals it is unlikely that Poisson fits well since animals will differ in many aspects such as exposure to infection and response to infection.

The negative binomial allows extra variation via the dispersion parameter. Like the mean, the dispersion is estimated from the data in such way that the likelihood of the data is maximised for the given mean and dispersion (in fact all glms are fitted via maximum likelihood).

R code to reproduce



a <- 4
b <- 4
c <- 4
d <- 4
count <- list(
    A=rnbinom(n=a, mu=50, size=10),
    B=rnbinom(n=b, mu=30, size=10),
    C=rnbinom(n=c, mu=5, size=10),
    D=rnbinom(n=d, mu=0.5, size=10)
dat <- data.table(
    cnd=rep(names(count), c(a, b, c, d))

norm <- lm(count ~ cnd, data=dat)
emm <- as.data.table(emmeans(norm, 'cnd'))

negbin <- glm.nb(count ~ cnd, data=dat)
emm2 <- as.data.table(emmeans(negbin, 'cnd', type='response'))

cc <- brewer.pal(3, 'Dark2')
cols <- c('lm'=cc[1], 'glm'=cc[2])
gg <- ggplot(data=dat, aes(x=cnd, y=count)) +
    geom_quasirandom(colour='grey20', width=0.1) +
    geom_hline(yintercept=0, colour='grey60', linetype='dashed') +
    geom_segment(data=emm2, aes(x=as.numeric(cnd) - 0.1, y=asymp.LCL,
    yend=asymp.UCL, colour='glm'), linewidth=0.5) +
    geom_point(data=emm2, aes(x=as.numeric(cnd) - 0.1, y=response, 
    colour='glm'), pch=4, size=2) +
    geom_segment(data=emm, aes(x=as.numeric(cnd) + 0.1, y=lower.CL, 
    yend=upper.CL, colour='lm'), linewidth=0.5) +
    geom_point(data=emm, aes(x=as.numeric(cnd) + 0.1, y=emmean, 
    colour='lm'), pch=4, size=2) +
    scale_colour_manual(name="95% CI from:", values=cols) +
    ggtitle('Data simulated from negative binomial distributions') + 
    theme_light() +
ggsave('error.pdf', width=12, height=10, units='cm')


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.