Multivariate generalized linear mixed models for underdispersed count data

Researchers are often interested in understanding the relationship between a set of covariates and a set of response variables. To achieve this goal, the use of regression analysis, either linear or generalized linear models, is largely applied. However, such models only allow users to model one response variable at a time. Moreover, it is not possible to directly calculate from the regression model a correlation measure between the response variables. In this article, we employed the Multivariate Generalized Linear Mixed Models framework, which allows the specification of a set of response variables and calculates the correlation between them through a random effect structure that follows a multivariate normal distribution. We used the maximum likelihood estimation framework to estimate all model parameters using Laplace approximation to integrate out the random effects. The derivatives are provided by automatic differentiation. The outer maximization was made using a general-purpose algorithm such as PORT and Broyden-Fletcher-Goldfarb-Shanno algorithm (BFGS). We delimited this problem by studying count response variables with the following distributions: Poisson, negative binomial, Conway-Maxwell-Poisson (COM-Poisson), and double Poisson. While the first distribution can model only equidispersed data, the second models equi and overdispersed, and the third and fourth models all types of dispersion (i.e. including underdispersion). The models were implemented on software R with package TMB, based on C++ templates. Besides the full specification, models with simpler structures in the covariance matrix were considered (fixed and common variance, and rho set to 0) and fixed dispersion. These models were applied to a dataset from the National Health and Nutrition Examination Survey, where two response variables are underdispersed and one can be considered equidispersed that were measured at 1281 subjects. The double Poisson full model specification overcame the other three competitors considering three goodness-of-fit measures: Akaike Information Criteria (AIC), Bayesian Information Criteria (BIC), and maximized log-likelihood. Consequently, it estimated parameters with smaller standard error and a greater number of significant correlation coefficients. Therefore, the proposed model can deal with multivariate count responses and measures the correlation between them taking into account the effects of the covariates.