Solved – $R^2$ from a generalized linear mixed-effects models (GLMM) using a negative binomial distribution

generalized linear modelglmmmixed modelnegative-binomial-distributionr

I try to compute the marginal and conditional $R^2$ for a GLMM using a negative binomial distribution by following the procedure recommended by Nakagawa & Schielzeth (2013) . Unfortunately, the supplementary material of their article does not include an example of a negative binomial distribution (see the online version of the article stated below, I also added their code below).
I fitted my model using the glmmPQL function from the MASS package.

full_model  <- glmmPQL ( Y~ a + b + c,  random = ~ 1 +  A | location  
, family = negative.binomial (1.4 ) ,data= mydata

In particular, I do have the following problems:

First, I need to extract the fixed effect design matrix of my model. However, full_model @X or model.matrix(full_model) does not work. I also tried to set the argument x=TRUE before extracting the matrix. Well, this should not be too tricky, but the following problems are.
Second, I need to specify the distribution-speciﬁc variance of my model. Examples in the article (see table 2 & and the supplementary R code of the online article) specify this for a binomial and a Poisson distribution. With some deeper statistical knowledge, it should not be difficult to specify this for a negative binomial distribution.
Finally, I would need to know if glmmPQL uses additive dispersion or to multiplicative dispersion. In the paper, they state: "we only consider additive dispersion implementation of GLMMs although the formulae that we present below can be easily modiﬁed for the use with GLMMs that apply to multiplicative dispersion. " Thus, in case glmmPQL uses multiplicative dispersion, I would need further help to adjust the formula.

Can anybody help?

Thanks, best
Philipp

P.S. R-code is welcome.

Nakagawa & Schielzeth (2013) A general and simple method for obtaining R 2 from generalized linear mixed-effects models. Methods in Ecology and Evolution 2013, 4, 133–142. doi: 10.1111/j.2041-210x.2012.00261.x

Their R script:

  #A general and simple method for obtaining R2 from generalized linear mixed-effects models
  #Shinichi Nakagawa1,2 and Holger Schielzeth3
  #1 National Centre of Growth and Development, Department of Zoology, University of    Otago, Dunedin, New Zealand
  #2 Department of Behavioral Ecology and Evolutionary Genetics, Max Planck Institute for Ornithology, Seewiesen, Germany
  #3 Department of Evolutionary Biology, Bielefeld University, Bielefeld, Germany
  #Running head: Variance explained by GLMMs
  #Correspondence:
  #S. Nakagawa; Department of Zoology, University of Otago, 340 Great King Street,    Dunedin, 9054, New Zealand
  #Tel:  +64 (0)3 479 5046
  #Fax: +64 (0)3 479 7584
  #e-mail: shinichi.nakagawa@otago.ac.nz 


  ####################################################
  # A. Preparation
  ####################################################
  # Note that data generation appears below the analysis section.
  # You can use the simulated data table from the supplementary files to reproduce exactly the same results as presented in the paper.

  # Set the work directy that is used for rading/saving data tables
  # setwd("/Users/R2")

  # load R required packages
  # If this is done for the first time, it might need to first download and install the package
  # install.package("arm")
  library(arm)
  # install.package("lme4")
  library(lme4)


  ####################################################
  # B. Analysis
  ####################################################

  # 1. Analysis of body size (Gaussian mixed models)
  #---------------------------------------------------

  # Clear memory
  rm(list = ls())

  # Read body length data (Gaussian, available for both sexes)
  Data <- read.csv("BeetlesBody.csv")

  # Fit null model without fixed effects (but including all random effects)
  m0 <- lmer(BodyL ~ 1 + (1 | Population) + (1 | Container), data = Data)

  # Fit alternative model including fixed and all random effects
  mF <- lmer(BodyL ~ Sex + Treatment + Condition + (1 | Population) + (1 | Container), data = Data)

  # View model fits for both models
  summary(m0)
  summary(mF)

  # Extraction of fitted value for the alternative model
  # fixef() extracts coefficents for fixed effects
  # mF@X returns fixed effect design matrix
  Fixed <- fixef(mF)[2] * mF@X[, 2] + fixef(mF)[3] * mF@X[, 3] + fixef(mF)[4] * mF@X[, 4]

  # Calculation of the variance in fitted values
  VarF <- var(Fixed)

  # An alternative way for getting the same result
  VarF <- var(as.vector(fixef(mF) %*% t(mF@X)))

  # R2GLMM(m) - marginal R2GLMM
  # Equ. 26, 29 and 30
  # VarCorr() extracts variance components
  # attr(VarCorr(lmer.model),'sc')^2 extracts the residual variance
  VarF/(VarF + VarCorr(mF)$Container[1] + VarCorr(mF)$Population[1] + attr(VarCorr(mF), "sc")^2)

  # R2GLMM(c) - conditional R2GLMM for full model
  # Equ. XXX, XXX
  (VarF + VarCorr(mF)$Container[1] + VarCorr(mF)$Population[1])/(VarF +    VarCorr(mF)$Container[1] + VarCorr(mF)$Population[1] + (attr(VarCorr(mF), "sc")^2))

  # AIC and BIC needs to be calcualted with ML not REML in body size models
  m0ML <- lmer(BodyL ~ 1 + (1 | Population) + (1 | Container), data = Data, REML = FALSE)
  mFML <- lmer(BodyL ~ Sex + Treatment + Condition + (1 | Population) + (1 | Container), data = Data, REML = FALSE)

  # View model fits for both models fitted by ML
  summary(m0ML)
  summary(mFML)


  # 2. Analysis of colour morphs (Binomial mixed models)
  #---------------------------------------------------

  # Clear memory
  rm(list = ls())
  # Read colour morph data (Binary, available for males only)
  Data <- read.csv("BeetlesMale.csv")

  # Fit null model without fixed effects (but including all random effects)
  m0 <- lmer(Colour ~ 1 + (1 | Population) + (1 | Container), family = "binomial", data = Data)

  # Fit alternative model including fixed and all random effects
  mF <- lmer(Colour ~ Treatment + Condition + (1 | Population) + (1 | Container), family = "binomial", data = Data)

  # View model fits for both models
  summary(m0)
  summary(mF)

  # Extraction of fitted value for the alternative model 
  # fixef() extracts coefficents for fixed effects 
  # mF@X returns fixed effect design matrix
  Fixed <- fixef(mF)[2] * mF@X[, 2] + fixef(mF)[3] * mF@X[, 3]

  # Calculation of the variance in fitted values
  VarF <- var(Fixed)

  # An alternative way for getting the same result
  VarF <- var(as.vector(fixef(mF) %*% t(mF@X)))

  # R2GLMM(m) - marginal R2GLMM
  # see Equ. 29 and 30 and Table 2
  VarF/(VarF + VarCorr(mF)$Container[1] + VarCorr(mF)$Population[1] + pi^2/3)

  # R2GLMM(c) - conditional R2GLMM for full model
  # Equ. XXX, XXX
  (VarF + VarCorr(mF)$Container[1] + VarCorr(mF)$Population[1])/(VarF +     VarCorr(mF)$Container[1] + VarCorr(mF)$Population[1] + pi^2/3)


  # 3. Analysis of fecundity (Poisson mixed models)
  #---------------------------------------------------

  # Clear memory
  rm(list = ls())

  # Read fecundity data (Poisson, available for females only)
  Data <- read.csv("BeetlesFemale.csv")

  # Creating a dummy variable that allows estimating additive dispersion in lmer 
  # This triggers a warning message when fitting the model
  Unit <- factor(1:length(Data$Egg))

  # Fit null model without fixed effects (but including all random effects)
  m0 <- lmer(Egg ~ 1 + (1 | Population) + (1 | Container) + (1 | Unit), family = "poisson", data = Data)

  # Fit alternative model including fixed and all random effects
  mF <- lmer(Egg ~ Treatment + Condition + (1 | Population) + (1 | Container) + (1 | Unit), family = "poisson", data = Data)

  # View model fits for both models
  summary(m0)
  summary(mF)

  # Extraction of fitted value for the alternative model 
  # fixef() extracts coefficents for fixed effects 
  # mF@X returns fixed effect design matrix
  Fixed <- fixef(mF)[2] * mF@X[, 2] + fixef(mF)[3] * mF@X[, 3]

  # Calculation of the variance in fitted values
  VarF <- var(Fixed)

  # An alternative way for getting the same result
  VarF <- var(as.vector(fixef(mF) %*% t(mF@X)))

  # R2GLMM(m) - marginal R2GLMM 
  # see Equ. 29 and 30 and Table 2 
  # fixef(m0) returns the estimate for the intercept of null model
  VarF/(VarF + VarCorr(mF)$Container[1] + VarCorr(mF)$Population[1] + VarCorr(mF)$Unit[1] + log(1 + 1/exp(as.numeric(fixef(m0)))))

  # R2GLMM(c) - conditional R2GLMM for full model
  # Equ. XXX, XXX
  (VarF + VarCorr(mF)$Container[1] + VarCorr(mF)$Population[1])/(VarF +    VarCorr(mF)$Container[1] + VarCorr(mF)$Population[1] + VarCorr(mF)$Unit[1] + log(1 + 
                                                                       1/exp(as.numeric(fixef(m0)))))


  ####################################################
  # C. Data generation
  ####################################################

  # 1. Design matrices 
  #---------------------------------------------------

  # Clear memory
  rm(list = ls())

  # 12 different populations n = 960
  Population <- gl(12, 80, 960)

  # 120 containers (8 individuals in each container)
  Container <- gl(120, 8, 960)

  # Sex of the individuals. Uni-sex within each container (individuals are sorted at the pupa stage)
  Sex <- factor(rep(rep(c("Female", "Male"), each = 8), 60))

  # Condition at the collection site: dry or wet soil (four indiviudal from each condition in each container)
  Condition <- factor(rep(rep(c("dry", "wet"), each = 4), 120))

  # Food treatment at the larval stage: special food ('Exp') or standard food ('Cont')
  Treatment <- factor(rep(c("Cont", "Exp"), 480))

  # Data combined in a dataframe
  Data <- data.frame(Population = Population, Container = Container, Sex = Sex, Condition = Condition, Treatment = Treatment)


  # 2. Gaussian response: body length (both sexes)
  #---------------------------------------------------

  # simulation of the underlying random effects (Population and Container with variance of 1.3 and 0.3, respectively)
  PopulationE <- rnorm(12, 0, sqrt(1.3))
  ContainerE <- rnorm(120, 0, sqrt(0.3))

  # data generation based on fixed effects, random effects and random residuals errors
  Data$BodyL <- 15 - 3 * (as.numeric(Sex) - 1) + 0.4 * (as.numeric(Treatment) - 1) + 0.15 * (as.numeric(Condition) - 1) + PopulationE[Population] + ContainerE[Container] + 
    rnorm(960, 0, sqrt(1.2))

  # save data (to current work directory)
  write.csv(Data, file = "BeetlesBody.csv", row.names = F)


  # 3. Binomial response: colour morph (males only)
  #---------------------------------------------------

  # Subset the design matrix (only males express colour morphs)
  DataM <- subset(Data, Sex == "Male")

  # simulation of the underlying random effects (Population and Container with variance of 1.2 and 0.2, respectively)
  PopulationE <- rnorm(12, 0, sqrt(1.2))
  ContainerE <- rnorm(120, 0, sqrt(0.2))

  # generation of response values on link scale (!) based on fixed effects and random effects
  ColourLink <- with(DataM, 0.8 * (-1) + 0.8 * (as.numeric(Treatment) - 1) + 0.5 *    (as.numeric(Condition) - 1) + PopulationE[Population] + ContainerE[Container])

  # data generation (on data scale!) based on negative binomial distribution
  DataM$Colour <- rbinom(length(ColourLink), 1, invlogit(ColourLink))

  # save data (to current work directory)
  write.csv(DataM, file = "BeetlesMale.csv", row.names = F)


  # 4. Poisson response: fecundity (females only)
  #---------------------------------------------------

  # Subset the design matrix (only females express colour morphs)
  DataF <- Data[Data$Sex == "Female", ]

  # random effects
  PopulationE <- rnorm(12, 0, sqrt(0.4))
  ContainerE <- rnorm(120, 0, sqrt(0.05))

  # generation of response values on link scale (!) based on fixed effects, random effects and residual errors
  EggLink <- with(DataF, 1.1 + 0.5 * (as.numeric(Treatment) - 1) + 0.1 *   (as.numeric(Condition) - 1) + PopulationE[Population] + ContainerE[Container] +   rnorm(480, 
                                                                                                                                                         0, sqrt(0.1)))

  # data generation (on data scale!) based on Poisson distribution
  DataF$Egg <- rpois(length(EggLink), exp(EggLink))

  # save data (to current work directory)
  write.csv(DataF, file = "BeetlesFemale.csv", row.names = F)

Best Answer

For Q1 check out the following paper http://onlinelibrary.wiley.com/doi/10.1111/2041-210X.12225/full It says their code is now included in the package "MuMIn". HTH.

Related Solutions

Solved – Help interpreting count data GLMM using lme4 glmer and glmer.nb – Negative binomial versus Poisson

I believe there are some important problems to be addressed with your estimation.

From what I gathered by examining your data, your units are not geographically grouped, i.e. census tracts within counties. Thus, using tracts as a grouping factor is not appropriate to capture spatial heterogeneity as this means that you have the same number of individuals as groups (or to put another way, all your groups have only one observation each). Using a multilevel modelling strategy allows us to estimate individual-level variance, while controlling for between-group variance. Since your groups have only one individual each, your between-group variance is the same as your individual-level variance, thus defeating the purpose of the multilevel approach.

On the other hand, the grouping factor can represent repeated measurements over time. For example, in the case of a longitudinal study, an individual's "maths" scores may be recored yearly, thus we would have a yearly value for each student for n years (in this case, the grouping factor is the student as we have n number of observations "nested" within students). In your case, you have repeated measures of each census tract by decade. Thus, you could use your TRTID10 variable as a grouping factor to capture "between decade variance". This leads to 3142 observations nested in 635 tracts, with approximately 4 and 5 observations per census tract.

As mentioned in a comment, using decade as a grouping factor is not very appropriate, as you have only around 5 decades for each census tract, and their effect can be better captured introducing decade as a covariate.

Second, to determine whether your data ought to be modelled using a poisson or negative binomial model (or a zero inflated approach). Consider the amount of overdispersion in your data. The fundamental characteristic of a Poisson distribution is equidispersion, meaning that the mean is equal to the variance of the distribution. Looking at your data, it is pretty clear that there is much overdispersion. Variances are much much greater than the means.

library(dplyr)    
 dispersionstats <- scaled.mydata %>%
 + group_by(decade) %>%
 + summarise(
 + means = mean(R_VAC),
 + variances = var(R_VAC),
 + ratio = variances/means)

##   dispersionstats
##   # A tibble: 5 x 5
##   decade     means variances     ratio 
##    <int>     <dbl>     <dbl>     <dbl> 
## 1   1970  45.43513   4110.89  90.47822 
## 2   1980 103.52365  17323.34 167.33707 
## 3   1990 177.68038  62129.65 349.67087 
## 4   2000 190.23150  91059.60 478.67784 
## 5   2010 247.68246 126265.60 509.78821

Nonetheless, to determine if the negative binomial is more appropriate statistically, a standard method is to do a likelihood ratio test between a Poisson and a negative binomial model, which suggests that the negbin is a better fit.

library(MASS)
library(lmtest)

modelformula <- formula(R_VAC ~ factor(decade) + P_NONWHT * a_hinc + offset(HU_ln))

poismodel <- glm(modelformula, data = scaled.mydata, family = "poisson")   
nbmodel <- glm.nb(modelformula, data = scaled.mydata)

lrtest(poismodel, nbmodel)

## Likelihood ratio test

##  Model 1: R_VAC ~ factor(decade) + P_NONWHT * a_hinc + offset(HU_ln)  
## Model 2: R_VAC ~ factor(decade) + P_NONWHT * a_hinc + offset(HU_ln)
##   #Df  LogLik Df  Chisq Pr(>Chisq)
## 1   8 -154269
## 2   9  -17452  1 273634  < 2.2e-16 ***
##  ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

After establishing this, a next test could consider whether the multilevel (mixed model) approach is warranted using similar approach, which suggests that the multilevel version provides a better fit. (A similar test could be used to compare a glmer fit assuming a poisson distribution to the glmer.nb object, as long as the models are otherwise the same.)

library(lme4)

glmmformula <- update(modelformula, . ~ . + (1|TRTID10))

nbglmm <- glmer.nb(glmmformula, data = scaled.mydata)

lrtest(nbmodel, nbglmm)

## Model 1: R_VAC ~ factor(decade) + P_NONWHT * a_hinc + offset(HU_ln)
## Model 2: R_VAC ~ factor(decade) + P_NONWHT + a_hinc + (1 | TRTID10) +
##     P_NONWHT:a_hinc + offset(HU_ln)
##   #Df LogLik Df Chisq Pr(>Chisq)
## 1   9 -17452
## 2  10 -17332  1 239.3  < 2.2e-16 ***
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Regarding the estimates of the poisson and nb models they are actually expected to be very similar to each other, with the main distinction being the standard errors, ie if overdispersion is present, the poisson model tends to provide biased standard errors. Taking your data as an example:

poissonglmm <- glmer(glmmformula, data = scaled.mydata)
summary(poissonglmm)

## Random effects:
##  Groups  Name        Variance Std.Dev.
## TRTID10 (Intercept) 0.2001   0.4473
## Number of obs: 3142, groups:  TRTID10, 635

## Fixed effects:
##                     Estimate Std. Error z value Pr(>|z|)
## (Intercept)        -2.876013   0.020602 -139.60   <2e-16 ***
## factor(decade)1980  0.092597   0.007602   12.18   <2e-16 ***
## factor(decade)1990  0.903543   0.007045  128.26   <2e-16 ***
## factor(decade)2000  0.854821   0.006913  123.65   <2e-16 ***
## factor(decade)2010  0.986126   0.006723  146.67   <2e-16 ***
## P_NONWHT           -0.125500   0.014007   -8.96   <2e-16 ***
## a_hinc             -0.107335   0.001480  -72.52   <2e-16 ***
## P_NONWHT:a_hinc     0.160937   0.003117   51.64   <2e-16 ***
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

summary(nbglmm)
## Random effects:
##  Groups  Name        Variance Std.Dev.
##  TRTID10 (Intercept) 0.09073  0.3012
## Number of obs: 3142, groups:  TRTID10, 635

## Fixed effects:
##                     Estimate Std. Error z value Pr(>|z|)
## (Intercept)        -2.797861   0.056214  -49.77  < 2e-16 ***
## factor(decade)1980  0.118588   0.039589    3.00  0.00274 **
## factor(decade)1990  0.903440   0.038255   23.62  < 2e-16 ***
## factor(decade)2000  0.843949   0.038172   22.11  < 2e-16 ***
## factor(decade)2010  1.068025   0.037376   28.58  < 2e-16 ***
## P_NONWHT            0.020012   0.089224    0.22  0.82253
## a_hinc             -0.129094   0.008109  -15.92  < 2e-16 ***
## P_NONWHT:a_hinc     0.149223   0.018967    7.87 3.61e-15 ***
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Notice how the coefficient estimates are all very similar, the main difference being only the significance of one of your covariates, as well as the difference in the random effects variance, which suggests that the unit-level variance captured by the overdispersion parameter in the nb model (the theta value in the glmer.nb object) captures some of the between tract variance captured by the random effects.

Regarding exponentiated coefficients (and associated confidence intervals), you can use the following:

fixed <- fixef(nbglmm)
confnitfixed <- confint(nbglmm, parm = "beta_", method = "Wald") # Beware: The Wald method is less accurate but much, much faster.

# The exponentiated coefficients are also known as Incidence Rate Ratios (IRR)
IRR <- exp(cbind(fixed, confintfixed)
IRR
##                         fixed      2.5 %     97.5 %
## (Intercept)        0.06094028 0.05458271 0.06803835
## factor(decade)1980 1.12590641 1.04184825 1.21674652
## factor(decade)1990 2.46807856 2.28979339 2.66024515
## factor(decade)2000 2.32553168 2.15789585 2.50619029
## factor(decade)2010 2.90962703 2.70410073 3.13077444
## P_NONWHT           1.02021383 0.85653208 1.21517487
## a_hinc             0.87889172 0.86503341 0.89297205
## P_NONWHT:a_hinc    1.16093170 1.11856742 1.20490048

Final thoughts, regarding zero inflation. There is no multilevel implementation (at least that I am aware of) of a zero inflated poisson or negbin model that allows you to specify an equation for zero inflated component of the mixture. the glmmADMB model lets you estimate a constant zero inflation parameter. An alternative approach would be to use the zeroinfl function in the pscl package, though this does not support multilevel models. Thus, you could compare the fit of a single level negative binomial, to the single level zero inflated negative binomial. Chances are that if zero inflation is not significant for single level models, it is likely that it would not be significant for the multilevel specification.

Addendum

If you are concerned about spatial autocorrelation, you could control for this using some form of geographical weighted regression (though I believe this uses point data, not areas). Alternatively, you could group your census tracts by an additional grouping factor (states, counties) and include this as a random effect. Lastly, and I am not sure if this is entirely feasible, it may be possible to incorporate spatial dependence using, for example, the average count of R_VAC in first order neighbours as a covariate. In any case, prior to such approaches, it would be sensible to determine if spatial autocorrelation is indeed present (using Global Moran's I, LISA tests, and similar approaches).

Best Answer

Related Solutions

Solved – Help interpreting count data GLMM using lme4 glmer and glmer.nb – Negative binomial versus Poisson

Addendum

Related Question