  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 \texttt{PORT} and \texttt{BFGS}. We delimited
this problem by studying only count response variables with the following
distributions: Poisson, negative binomial (NB) and COM-Poisson. The models were
implemented on software \texttt{R} with package \texttt{TMB}. Besides the full
specification, models with simpler structures in the covariance matrix were
considered (fixed and common variance, fixed dispersion, $\rho$ set to 0).
These models were applied to a dataset from the National Health and Nutrition
Examination Survey, where three underdispersed response variables were measured
at 1281 subjects. The COM-Poisson model full specified overcome the other two
competitors considering three goodness-of-fit indexes. Therefore, the proposed
model can deal with multivariate count responses and measures the correlation
between them taking into account the effects of the covariates.
