§ INTRODUCTION

It is well known that the Poisson distribution is the most popular model to deal with count data under the framework of the generalized linear models (GLM) [Nelder and Wedderburn, 1972]. However, it is limited to equidispersed count data, i.e., when the mean of the response variable is equal to the variance. As alternatives to the Poisson model, the statistical literature for univariate count data is rich, either from overdispersed or underdispersed perspectives. For example, negative binomial (NB) type II distribution parameterized on mean and dispersion parameter known as the quadratic parametrization of the variance [Winkelmann, 2008], hurdle and zero-inflated models [Zeileis et al., 2008] and mixed Poisson regression [Winter and Bürkner, 2021] are frequent options to model overdispersed count data; Gamma Count [Zeviani et al., 2014], Conway-Maxwell-Poisson (COM-Poisson) [Shmueli et al., 2005] and the Extended Poisson Tweedie [Bonat et al., 2018] based on the Poisson Tweedie distribution [El-Shaarawi et al., 2011, Jørgensen and Kokonendji, 2016] can handle all three situations (under-, equi-, over-dispersed), even though the computation of their probability mass function (pmf) relies on numerical methods.

On the other hand, the bibliography for multivariate data is scarce. However, there is an increasing demand from researchers to analyze datasets with over one response variable. The benefit of it is to better investigate the relationship between the response variables [Bonat and Jørgensen, 2016]; at the same time that numerical methods can take advantage of it once more data is available to estimate model parameters.

Some methodologies to deal with multivariate data have been proposed and introduced by Winkelmann, 2008 with a focus on building multivariate distributions. However, it comes with the price of some practical limitations. The multivariate Poisson-lognormal regression (MPLR) and the latent Poisson-normal regression models admit only overdispersed data. The multivariate Poisson-gamma mixture model (MPGM) and multivariate NB model (MNBM) are suitable only for overdispersed data with positive correlations. Another example is the copulas framework, which allows the building of multivariate distributions. But a negative correlation between many response variables is difficult to model, especially for count data [Nikoloulopoulos and Karlis, 2009]. Inouye et al., 2017 proposed three alternatives to deal only with equi and overdispersed data constructed either via copulas or a mixture of distributions.

Distributions with no practical limitations were also developed. Famoye, 2015 proposed a multivariate generalized Poisson regression model that can deal with any kind of dispersion and correlation with an estimation based on the maximum likelihood (ML) paradigm. In a similar way, Muñoz-Pichardo et al., 2021 proposed a multivariate conditional Poisson regression model, where the relationship between response variables is measured by a coefficient in the linear predictor. In its turn, the dependence between response variables is conditional on the other response variables.

A very flexible modeling framework based on estimating functions, the Multivariate Covariance Generalized Linear Models (MCGLM) can also fit such data [Bonat, 2016]. It uses only second-order moments assumptions and they estimated the parameters based on quasi-likelihood [Wedderburn, 1974]. It can accommodate correlated data based on an approach similar to the generalized estimating equations (GEE) [Liang and Zeger, 1986], allowing both multivariate responses and correlated data. We can also cite methodologies based on Bayesian inference, such as MCMC Generalized Linear Mixed Models via  package [Hadfield, 2010] and Bayesian Regression Models using Stan -  package [BÜrkner, 2018], which comes with the price of a greater computational time. We are not going to discuss Bayesian models in this article.

Another alternative is to model the correlation between response variables in the same individual using the class of hierarchical GLM [Lee and Nelder, 1996]. This class allows to model of correlated variables or individuals via a random effect, an unobserved variable, that can follow any distribution. When the distribution of the random effect is Gaussian, we have the Generalized Linear Mixed Models (GLMM). However, GLMM is widely known and used to model the correlation between sample units, not for response variables, such a method is implemented in consolidated packages in software  [R Core Team, 2020], such as glmmTMB [Brooks et al., 2017], lme4 [Bates et al., 2015] and nlme [Pinheiro et al., 2017].

In this article, we propose to model multivariate count data under the framework of GLMMs to accommodate the correlation between response variables. Parameter estimates are obtained through the maximum likelihood method [Aldrich et al., 1997]. We use multivariate normally distributed random effects to accommodate the correlation between response variables. The estimation process is similar to the one for GLMM and was implemented in  [R Core Team, 2020] through  package [Kristensen et al., 2016].

Here, we will focus on multivariate overdispersed data, whereas da Silva et al., 2022 presented a similar approach focused on underdispersed data. This modeling approach for multivariate count data is evaluated for three different distributions of the response variables: Poisson, NB, and COM-Poisson mean parameterized [Huang, 2017]. Even though COM-Poisson does not belong to the exponential family, we refer to it as a GLMM framework, once the estimation process remains the same regardless of the distribution being used.

This article contains six sections, including this introduction. Section 2 describes the datasets used to provide illustrative applications of the model. Section 3 proposes the multivariate generalized linear mixed model (MGLMM) model. Section 4 shows the result of the simulation study to assess the estimators' properties. Section 5 presents the results of the model applied to the datasets presented in Section 2. Finally, Section 6 discusses the major contributions of this article and future work.

§ DATA SETS

§.§ Dataset I: Australian Health Survey

The AHS is the largest survey conducted in Australia concerning health issues. The data used here is available through the  package [Bonat, 2018] in the  object. The main aim of this study is to investigate whether more access to health care services and demographic covariates such as sex, age, and income are related to the number of times patients use health services.

The data has 5190 respondents and 15 variables; of which 10 were considered covariates and 5 as count response variables of a cross-sectional study with individuals over 18 years. The five response variables are Ndoc (Number of consultations with a doctor or specialist), Nndoc (Number of consultations with health professionals), Nadm (Number of admissions to a hospital, psychiatric hospital, nursing or convalescence home in the past 12 months), Nhosp (Number of nights in a hospital during the most recent admission) and Nmed (Total number of prescribed and non-prescribed medications used in the past two days). Table S1 in the supplementary material gives the description of each covariate.

<ref> shows the barplot for each response variable. The figure suggests that all variables have some overdispersion, because of the right heavy-tailed and perhaps some excess of zeros. We highlight that such interpretations are marginal in the sense it does not take into account the effects of the covariates.
<ref> presents some descriptive statistics for the 5 response variables. It is seen a positive (small in most cases) correlation between the 5 response variables, which is acceptable once they were measured in the same individual; except between Nadm and Nhosp, where the correlation is 0.996. Moreover, the mean and variance relationship between the variables shows overdispersion for all variables, which Nhosp being the most overdispersed. This is characterized by Fisher Dispersion Index (DI) greater than 1 [Fisher, 1934]. Also, the generalized dispersion index (GDI) [Kokonendji and Puig, 2018] classifies this dataset as overdispersed once it is greater than 1 and its 95% confidence interval does not contain zero. However, these results should be confirmed by model fitting.

Barplot for the count of each response variable from Australian Health Survey (AHS) data. Ndoc (Number of consultations with a doctor or specialist), Nndoc (Number of consultations with health professionals), Nadm (Number of admissions to a hospital, psychiatric hospital, nursing or convalescence home in the past 12 months), Nhosp (Number of nights in a hospital during the most recent admission) and Nmed (Total number of prescribed and non-prescribed medications used in the past two days).

Spearmean correlation, mean, variance, and dispersion index (DI) for the Australian Health Survey (AHS) response variables. The generalized dispersion index (GDI) and standard error (SE) equal 17.944 (1.99). Ndoc (Number of consultations with a doctor or specialist), Nndoc (Number of consultations with health professionals), Nadm (Number of admissions to a hospital, psychiatric hospital, nursing or convalescence home in the past 12 months), Nhosp (Number of nights in a hospital during the most recent admission) and Nmed (Total number of prescribed and non-prescribed medications used in the past two days).
1c      4cSpearman Correlation $\rho$     1cMean     1cVariance     1cDI

(l3ptr3pt)2-5 (l3ptr3pt)6-6 (l3ptr3pt)7-7 (l3ptr3pt)8-8
Nndoc     Nadm     Nhosp     Nmed

Ndoc     0.11     0.158     0.160     0.290     0.302     0.637     2.111

Nndoc          0.107     0.113     0.165     0.215     0.932     4.341

Nadm               0.996     0.135     0.174     0.258     1.484

Nhosp                    0.139     1.334     37.455     28.083

Nmed                         1.218     2.423     1.989

§.§ Dataset II: Abundance ant species

This data was obtained from a study conducted in south-eastern Australia by Gibb et al., 2015, and it is available in the software R [R Core Team, 2020] throughout the  package [Wang et al., 2020]. The study consisted of counting the number of 41 different ant species that fell into pitfall traps during 18-day sessions in 30 different sites in south-eastern Australia between November 2007 and April 2008. The primary interest of this dataset is to investigate the relationship between the environmental variables with the occurrence of different ant species and the co-occurrence of ant species.

The response variables considered are the occurrence of each of the 41 species at the 30 sites, while the covariates comprise 5 environmental variables and we gave their full description in the supplementary material Table S2. The name of the response variables starts with an index number (1,...,41) followed by the abbreviated name of each species.

<ref> presents the barplot for each response variable from ANT data. The 41 different species present high variability. Some species were only seen one time in a single site, such as Polyrhachis and Solenopsis, while Iridomyrmex and Pheidole were seen 20 times on over one site.

Barplot for the occurrence of each ant species genus from ANT data.

<ref> presents the mean, variance, and DI for every response variable and the GDI for the dataset. The variance is higher than the mean almost for all variables, except for Polyrhachis and Solenopsis; which is described by a DI>1. This dataset is classified as overdispersed based on the GDI index = 11,543 and the lower bound of a 95% confidence interval was greater than 1.

Mean, variance, dispersion index (DI) for each genus for ANT data. Generalized dispersion index (GDI) and its standard error (SE) equals 11.54 (.92)
Response     Mean     Variance     DI     Response     Mean     Variance     DI

1.Amblyopone     0.533     1.085     2.034     21.Meranoplus.A     1.333     8.713     6.534

2.Aphaenogaster     2.033     32.999     16.229     22.Monomorium.L     4.733     32.409     6.847

3.Camponotus.Ci     0.167     0.213     1.276     23.Monomorium.Ro     3.433     20.944     6.100

4.Camponotus.Cl     1.933     7.099     3.672     24.Monomorium.Sy     1.167     6.902     5.916

5.Camponotus.Co     5.300     33.252     6.274     25.Myrmecia     0.400     0.662     1.655

6.Camponotus.Ni     0.867     3.430     3.958     26.Notoncus.Ca     0.633     2.654     4.191

7.Camponotus.Ns     0.667     1.816     2.724     27.Notoncus.Ec     2.900     12.300     4.241

8.Cardiocondyla     0.500     2.879     5.759     28.Nylanderia     3.733     25.720     6.889

9.Crematogaster     0.567     3.633     6.412     29.Ochetellus     0.433     1.082     2.496

10.Heteroponera     4.067     19.857     4.883     30.Paraparatrechina     0.700     4.976     7.108

11.Iridomyrmex.Bi     5.333     26.437     4.957     31.Pheidole.A     11.100     40.300     3.631

12.Iridomyrmex.Dr     0.600     1.834     3.057     32.Pheidole.B     0.567     2.047     3.613

13.Iridomyrmex.Mj     4.300     15.803     3.675     33.Pheidole.E     5.467     22.809     4.172

14.Iridomyrmex.Pu     2.033     20.447     10.056     34.Pheidole.J     0.200     0.372     1.862

15.Iridomyrmex.Ru     13.300     59.045     4.439     35.Polyrhachis     0.300     0.217     0.724

16.Iridomyrmex.Su     2.133     6.809     3.192     36.Rhytidoponera.A     8.300     25.528     3.076

17.Iridomyrmex.Ss     0.900     1.610     1.789     37.Rhytidoponera.B     2.400     11.834     4.931

18.Melophorus.E     0.333     0.506     1.517     38.Solenopsis     0.133     0.120     0.897

19.Melophorus.F     2.133     6.740     3.159     39.Stigmacros     0.967     3.826     3.958

20.Melophorus.H     0.900     2.783     3.092     40.Tapinoma     2.533     11.154     4.403

41.Tetramorium     1.933     7.030     3.636

<ref> explores how the occurrence of different species is related to each other by the Spearman correlation in a correlogram. The marginal correlation ranges from -.58 up to .74, which is well distributed along all potential values of the correlation parameter $\rho = \left ( -1,1 \right )$.

Correlogram of ANT species occurrence using spearman correlation

§ MULTIVARIATE GENERALIZED LINEAR MIXED MODEL (MGLMM) FOR COUNT DATA

This section follows the same description in da Silva et al., 2022 and presents the MGLMM model to completeness; even though the explanations are the same, but succinct. The fundamental difference is that in this article we will explore the model for overdispersed counts. In summary, we extended the standard GLMM [Breslow and Clayton, 1993] for multiple responses.

Let $Y_{ij}$ be a random variable for $i=1,\dots,n$ independent subjects and $r=1,\dots,k$ response variables. Also, consider that $\mathrm{x}_{irj}$ for each subject i and response variable r is the value of the j-th covariate of k known covariates. We first specify a standard GLMM joint model with only a random intercept:

\begin{equation*}
% \label{eq:dist_assumption}
\mathrm{Y}_{ir} \mid \mathrm{b}_{ir} \sim f(\mu_{ir};\phi_{ir}),
\end{equation*}

where the conditional distribution of $\mathrm{Y}_{ir}$ on the random effects follows a same distribution f allowing different mean and dispersion. Next, the linear predictor:

\begin{equation*}
\label{eq:linear_predictor}
\mathrm{g}_r(\mu_{ir}) = x_{irj}^{\top}\boldsymbol{\beta}_{r} + \mathrm{b}_{ir},
\end{equation*}

where $\mathrm{g}_r(.)$ is a suitable link function, $\boldsymbol{\beta}_{r}$ is a $p \times 1$ vector of parameter estimates and $\mathrm{b}_{ir}$ is the random intercept value for subject i and response r. Finally, the random effects distribution is specified as:

\begin{equation*}
% \label{covmatrix}
\begin{pmatrix}
\mathrm{b}_{i1}\\
\mathrm{b}_{i2}\\
\vdots \\
\mathrm{b}_{ir}
\end{pmatrix}
\sim \mathrm{NM} \begin{pmatrix}
\begin{bmatrix}
\vdots\\
\end{bmatrix};
\underset{r\times r}{\boldsymbol{\Sigma}} = \begin{bmatrix}
\sigma^{2}_{1} & \rho_{12}\sigma_1\sigma_2 & \dots & \rho_{1r}\sigma_1\sigma_r\\
\rho_{21}\sigma_2\sigma_1 & \sigma^{2}_{2}  & \dots & \rho_{2r}\sigma_2\sigma_r \\
\vdots & \vdots & \ddots  & \vdots \\
\rho_{r1}\sigma_r\sigma_1 & \rho_{r2}\sigma_r\sigma_2 & \dots & \sigma^{2}_{r}
\end{bmatrix}
\end{pmatrix}
\end{equation*}

where each random effect has mean 0, variance $\sigma^2$, and correlation $\rho_{r{r}'} (r\neq{r}')$ between each pair of random effects. This framework is general and can be applied to any distribution f and link function $\mathrm{g}_r(.)$. Nevertheless, as we are dealing with count data, we considered only Poisson, binomial negative type II and COM-Poisson distributions, and a logarithm link function. Extending this model for different distributions for each response ($f_r$) is also possible, but we will not address it.

Maximum likelihood is the estimation procedure, and it is fully described in da Silva et al., 2022. We integrated out the random effects via numerical integration using Laplace approximation [Tierney and Kadane, 1986]. The Newton-Raphson method was efficiently implemented to perform the inner optimization. We optimized the marginal likelihood using first-order derivatives methods, such as  and  routines available in the  software. The derivatives of the joint and marginal likelihood were obtained through automatic differentiation [Baydin et al., 2017]. This was efficiently implemented by the  package, which provides  templates where an objective function must be supplied, in this case, the log-likelihood function.

The variability of the random variable in this model is being measured by two parameters, the variance of the random effect and the dispersion of the pmf. A crucial point of the model is the great flexibility to learn these two types of variances with no identifiability problems.

The code used to produce the results of this paper is available on <https://github.com/guilhermeparreira/papers/tree/master/AJS_MGLMM_Overdispersed_Count> and <http://www.leg.ufpr.br/doku.php/publications:papercompanions>.

§ SIMULATION STUDIES

In this section, we present simulation studies to assess the properties of the MLE estimators (bias and consistency). We considered a bivariate regression model for count data. We designed 12 simulation scenarios with four different sample sizes, 100, 250, 500, and 1000, and three different correlations between random effects, $\rho = -0.5, 0, 0.5$. For the regression structure, we considered only an intercept for each response, with $\beta_{01} = \text{log}(7)$ and $\beta_{02} = \text{log}(1.5)$. The variance of random effects were $\sigma^{2}_{1} = .3$ and $\sigma^{2}_{2} = .15$. The dispersion parameter for NB and COM-Poisson was equal to $\phi = 1$ and $\nu = .7$ respectively, which induces a small overdispersion.

We generated 150, 200, and 300 datasets for Poisson, COM-Poisson, and NB distribution for each design. The primary idea was to generate 100 datasets for each distribution. However, as it did not return the SE of the parameter estimates in every repetition, it was necessary to increase the number of datasets generated proportionally to the number of SE failures for each distribution to obtain at least 100 valid estimations. The following three subsection presents the results for Poisson, NB, and COM-Poisson scenarios.

§.§ Poisson scenario

<ref> presents the average bias and its confidence interval based on the average SE by sample size and simulation scenario for the bivariate Poisson regression model. It shows that, for all simulation scenarios, both the average bias and SE tend to be 0 as the sample size increases. This suggests the consistency and unbiasedness of the MLE estimator for large samples.

Average bias and its 95% confidence interval based on the average standard error (SE) by four different sample sizes, three correlation coefficients for the bivariate Poisson regression model. The true value for each parameter were with $\beta_{01} = \text{log}(7)$, $\beta_{02} = \text{log}(1.5)$, $\sigma^{2}_{1} = .3$ and $\sigma^{2}_{2} = .15$.

<ref> presents the coverage rate for each parameter by sample size and simulation scenario for the bivariate Poisson regression model. Overall, all empirical coverage rates are close to the nominal level of 95%, varying between 90% and 98% approximately. In special, the coverage rate of regression parameters $\boldsymbol{\beta}_{0r}$ are slightly greater than the nominal level. For the variance of random effects $\boldsymbol{\sigma}_{r}$, the coverage rate is slightly smaller than the nominal level. Finally, for the correlation between random effects $\rho$, the coverage rate is slightly lower when $\rho = .5$, and slightly greater when $\rho = \{0, -0.5\}$ compared to the 95% nominal level.

Coverage rate for each parameter (in the columns) by sample size and correlation coefficient for the Poisson regression model. The true value for each parameter were $\beta_{01} = \text{log}(7)$, $\beta_{02} = \text{log}(1.5)$, $\sigma^{2}_{1} = .3$ and $\sigma^{2}_{2} = .15$.

Even for the Poisson distribution, which is the simplest case to estimate because there is no dispersion parameter, there were 13 out of 1800 (3$\times$4$\times$150) simulations that did not return SE for some parameters estimates or produced extreme values due to large SEs that were not considered into the results. We used the  algorithm to estimate the model because it was more stable and faster than  in most situations. It also happened in Kristensen et al., 2016, where 9 study cases from different model settings ranging from linear regression to multivariate stochastic volatility models were considered, and  had a better performance than .

§.§ Negative Binomial scenario

<ref> suggests that for all simulation scenarios, both the average bias and SE tend to be 0 as the sample size increases. This suggests the consistency and unbiasedness of the MLE estimator for large samples. Concerning the variance of the random effect $\sigma_1$ and $\sigma_2$ were slightly underestimated.

Average bias and its 95% confidence interval based on the average standard error (SE) by four different sample sizes, three correlation coefficients for the bivariate negative binomial (NB) regression model. The true value for each parameter were $\beta_{01} = \text{log}(7)$, $\beta_{02} = \text{log}(1.5)$, $\phi = 1$, $\sigma^{2}_{1} = .3$ and $\sigma^{2}_{2} = .15$.

Overall, <ref> shows that the empirical coverage rates are close to the nominal level of 95%. In particular, the coverage rate of regression parameters $\boldsymbol{\beta}_{0r}$, the variance of random effects $\boldsymbol{\sigma}_{r}$ and the correlation between random effects $\rho$ are slightly greater than the nominal level. On the opposite, there was a coverage rate close to 80% for $\rho$ estimated when the sample size was equal to 100 and $\rho = 0$.

Coverage rate for each parameter (in the columns) by sample size and correlation coefficient for bivariate negative binomial (NB) regression model. The true value for each parameter were $\beta_{01} = \text{log}(7)$, $\beta_{02} = \text{log}(1.5)$, $\phi = 1$, $\sigma^{2}_{1} = .3$ and $\sigma^{2}_{2} = .15$.

Estimation problems were more severe for the bivariate NB regression model compared to the Poisson and occurred in 740 out of 3600 (3$\times$4$\times$300) simulations (20,6%). It was also necessary to remove those iterations when the dispersion parameter $\phi$ was greater than 5 from the results (we simulated at 1).

§.§ COM-Poisson scenario

<ref> shows that, for all simulation scenarios, the SE tends to be 0 as the sample size increases, but the bias does not. This suggests the non-consistency of the MLE estimator. Moreover, almost every estimate is biased. While $\nu_2$ is overestimated, $\nu_1$ is underestimated (and the bias does not decrease with an increase in sample size). The correlation parameter $\rho$ has a typical behavior: when the data was simulated with $\rho=0$ there was almost no bias, for $\rho=.5$ a negative bias and for $\rho=-.5$ a positive bias; the model forces the correlation parameter towards zero. The standard deviation of random effect $\sigma_2$ is overestimated, while $\sigma_1$ is slightly underestimated. Regarding the regression parameters, $\beta_{02}$ is slightly underestimated and $\beta_{01}$ suggests a negligible bias.

An interesting possible relationship among parameters is that when $\nu$ was overestimated (making the model more underdispersed: $\nu>1$), $\sigma$ was also overestimated (increasing the variance of the model). In contrast, when $\nu$ was underestimated (making the model even more overdispersed: $\nu<1$), $\sigma$ was also underestimated (decreasing the variance of the model). It seems that when $\nu$ makes the dispersion of the model greater, $\sigma$ makes the variance lower, and vice versa. It may be related because of a possible non-orthogonality between these parameters.

Average bias and its 95% confidence interval based on the average standard error (SE) by four different sample sizes, three correlation coefficients for the bivariate COM-Poisson regression model. The true value for each parameter were $\beta_{01} = \text{log}(7)$, $\beta_{02} = \text{log}(1.5)$, $\nu = .7$, $\sigma^{2}_{1} = .3$ and $\sigma^{2}_{2} = .15$.

Overall, <ref> shows that empirical coverage rates are not close to the nominal level of 95%, which agrees with the results presented in <ref>, where the bias did not decrease even for a higher sample size. In particular, the coverage rate of $\sigma_{2}$ had the worst results among all parameters (the coverage rate decreases as sample sizes increase), followed by $\nu_2$. Not surprisingly, these parameters had the two largest biases in <ref>. In contrast, results for $\beta_{01}$, $\sigma_{1}$ and $\rho$ (when the data was generated with $\rho=0$) had coverage rates close to 95% level.

Coverage rate for each parameter (in the columns) by sample size and correlation coefficient for bivariate COM-Poisson regression model. The true value for each parameter were $\beta_{01} = \text{log}(7)$, $\beta_{02} = \text{log}(1.5)$, $\nu = .7$, $\sigma^{2}_{1} = .3$ and $\sigma^{2}_{2} = .15$.

Estimation problems were more severe for the COM-Poisson regression model compared to the Poisson and less severe if compared to NB, accounting for 245 (10,2%) out of 2400 (3$\times$4$\times$200). Besides the rules used for Poisson to classify extreme values, it was necessary to remove those models when the dispersion parameter $\nu$ was greater than 4 and the SE of $\rho$ estimate was greater than 2.

§ DATA ANALYSES

This section presents the data analyses of the two datasets presented from <ref>. Initial values for the model in <ref> were chosen carefully in the following way. Firstly, and for every dataset, we fitted a Poisson model via quasi-likelihood in  package [Bonat, 2016] to obtain initial parameter estimates for the regression and variance parameters (based on the variance of the residuals). Because of the difference in methodologies, we set the correlation parameter to 0. These were initial values for the Poisson model; the final estimates from the Poisson became the initial values for the NB model, and from the NB, the initial estimates of COM-Poisson. The initial values to the dispersion parameter were $\phi = 1$ for NB, i.e., small overdispersion, and $\nu = 1$, i.e., equidispersion to the COM-Poisson model.

The estimation process started with a simple random sample (SRS) of size 300 (seed 2390) for AHS data. The parameters obtained from the sample were used as initial values to fit the model to the whole dataset. For ANT data, we did not use this step because of the small size of the dataset (30). We then optimize the same model more than once iterating over optimizers. In the first round we used  routines and its final estimates were used as initial estimates to optimize the same model with ; in its turn, the final estimates from  were used as initial estimates to ; and another final round to . We used this rationale because these models are complex to estimate and this procedure produced reasonable results.

Besides the full model presented in <ref>, we also tried to estimate four simpler versions to check whether the model was over-specified. We made the simplifications either by fixing values from the $\boldsymbol{\Sigma}$ or the dispersion parameter. The first one aimed to reproduce independent response variables setting the correlation parameter to 0. In the second version, we fixed $\phi = 1$ and $\nu=1.5$ for NB and COM-Poisson models inducing small overdispersion and underdispersion, respectively. In the third simplified model, we set the variance to $1$. The last strategy used only a single variance parameter to estimate all response variables. The last three simpler versions described aimed to answer whether we can estimate both dispersion and variance parameters in the same model. We only tested them for NB and COM-Poisson models since they have dispersion parameters.

The reported results are only from the full data and the best model version for each distribution after removing covariates that did not have the SEs estimated. We then presented the maximized log-likelihood value, Akaike and Bayesian Information Criterion (AIC and BIC), and the number of parameters estimated (np).

§.§ AHS data

<ref> presents the goodness-of-fit measures for AHS data from different distributions and specifications.

Goodness-of-fit measures for AHS data from different distributions and specifications. Number of parameters estimated (np), Akaike and Bayesian Information Criterion (AIC and BIC), and log-likelihood value (Loglik).
Model     np     AIC     BIC     Loglik

Poisson     70     33814     34272     -16837

NB     75     33826     34318     -16838

COM-Poisson     75     33106     33598     -16478

The COM-Poisson full model was the best among all distributions and the simpler versions of each distribution.

<ref> presents the regression parameter estimates and 95% confidence intervals by outcome and final model. We can see that the confidence interval amplitude is almost the same for all models, but slightly smaller for the COM-Poisson models over the competitors for almost all estimates for the Ndoc response variable; for the other response variables, we saw little differences. The point estimates were very close between Poisson and NB models, with a difference for the COM-Poisson model in some estimates, for example, freepoor, age, sex, illness, and chcond limited for Ndoc. We may relate this to the bias found in <ref>.

Regression parameter estimates and 95% confidence intervals by outcome and final model for the Australian Health Survey (AHS) data.

<ref> presents the dispersion estimates for each model and outcome. Even though this data can be considered as marginally overdispersed according to the $\text{GDI} = 17.94$ presented in <ref>, we see that $\phi$ approaches the infinity (suggest an equidispersed model by NB distribution), and $\nu$ is greater than 1 for all response variables and indicates underdispersion. The expected result would be $\phi$ and $\nu$ approaching zero, which results in an overdispersion of both distributions. However, as explored in <ref>, this behavior may occur because of the variance of the random effect. For the COM-Poisson model, the only $\nu$ estimate which indicates overdispersion was for Nmed, which had a sample DI equal to 1.99 (small overdispersion).

Dispersion of parameter estimates and standard errors (SEs) for each model and outcome of the Australian Health Survey (AHS) data
1c      2cNB($\phi$)     2cCOM-Poisson($\nu$)

(l3ptr3pt)2-3 (l3ptr3pt)4-5
Outcome     Estimate     SE     Estimate     SE

Ndoc     1.7e+04     1.8e+05     9.160     0.432

Nndoc     5.9e+03     6.0e+04     2.837     0.336

Nmed     3.6e+03     2.4e+04     0.674     0.032

Nhosp     8.5e+03     1.2e+05     5.110     0.615

Nadm     8.5e+03     8.1e+04     6.665     0.431

We presented the correlation coefficient estimates and their SEs in parentheses in (1) for the COM-Poisson model. Among the 10 correlation coefficients returned, COM-Poisson had 6 significant correlation coefficients. The stars in the matrix represent significant coefficients at the 5% level. Even though we do not have a primary interest in interpreting each correlation coefficient estimate, it is important to know which of them is significant, because it may be related to a smaller SE.

We can see an almost perfect significant correlation between Nadm and Nhosp, and a strong correlation between Nmed and Ndoc. We had also seen this pattern in <ref>.

\begin{multline}
\begin{aligned}
\label{eq:corrmatrixahs}
\text{COM-Poisson}_\rho&=
\begin{blockarray}{*{5}{c} l}
\begin{block}{*{5}{>{$\footnotesize}c<{$}} l}
Ndoc & Nndoc & Nmed & Nhosp & Nadm & \\
\end{block}
\begin{block}{[*{5}{c}]>{$\footnotesize}l<{$}}
1 & -0.07(0.03)^{*} & 0.78(0.09)^{*} & -0.02(0.02) & -0.02(0.02) \bigstrut[t]& Ndoc\\
& 1 & 0.34(0.14)^{*} & -0.04(0.03) & -0.05(0.03) & Nndoc\\
&  & 1 & 0.25(0.09)^{*} & 0.25(0.09)^{*} & Nmed\\
&  &  & 1 & 0.99(<.01)^{*} & Nhosp\\
&  &  &  & 1 & Nadm \\
\end{block}
\end{blockarray}
\end{aligned}
\end{multline}

\subsection{ANT data}
\label{cap:resultANT}

We presented the best model for each distribution in \autoref{tab:antfit2}. Here, we included two models for the COM-Poisson distribution because they were equivalent in the intermediate models based on the adopted criterion.

\begin{table}[H]

\caption{\label{tab:antfit2}Goodness-of-fit measures for ANT data from the best parametrization for each distribution. Number of parameters estimated (np), Akaike and Bayesian Information Criterion (AIC and BIC), and log-likelihood value (Loglik).}
\centering
\begin{tabular}[t]{lccccc}
\toprule
Model & np & AIC & BIC & logLik\\
\midrule
Poisson & 1057 & 5085.08 & 6566.15 & -1485.54\\
NB & 1098 & 5166.98 & 6705.50 & -1485.49\\
COM-Poisson & 1098 & 4965.87 & 6504.38 & -1384.93\\
Fixed dispersion COM-Poisson & 1057 & 5016.70 & 6497.77 & -1451.35\\
% Poisson & 1057 & 5085.08 & 6566.15 & -1485.54 & \cmark\\
% NB & 1098 & 5166.98 & 6705.50 & -1485.49 & \cmark\\
% COM-Poisson & 1098 & 4965.87 & 6504.38 & -1384.93 & \cmark\\
% Fixed dispersion COM-Poisson & 1057 & 5016.70 & 6497.77 & -1451.35 & \cmark\\
\bottomrule
\end{tabular}
\end{table}

From \autoref{tab:antfit2} we can see that the best models were from COM-Poisson distribution; Poisson and NB models had similar \texttt{logLik} results. The COM-Poisson models had higher \texttt{logLik} and \texttt{BIC}, and smaller \texttt{AIC} than the fixed dispersion model. A logLik ratio test (LRT) between these two COM-Poisson models resulted in $p<.00001$ ($\chi_{41}=132.84$), which gives evidence in favor of the fully specified model. Therefore, we will present the estimates of the full models for each distribution.

\autoref{fig:antbeta1} presents the regression parameter estimates and 95\% confidence intervals by outcome and final model for the first 12 response variables. We presented the same graphic for the remaining response variables in the supplementary material, Figures S1, S2, and S3, once we found similar patterns for all response variables. First, we can see that not all response variables share the same linear predictor (covariate feral mammal dung is not presented in the second, third and ninth response variables). Second, there were still some regression estimates' SEs that were not returned by the model and are presented by a hollow circle. It was necessary to make the distinction when the SE was returned and when it was tiny (filled circle, such as bare ground for COM-Poisson and response variables 1-8 and 10-12).

Overall, the confidence intervals for the COM-Poisson model were smaller than its counterparts and the point estimates were close among all models. The feral mammal dung covariate had the largest confidence intervals compared to the other covariates. There are still regression coefficients that are very close to zero, suggesting that a better variable selection method may improve the model fit.

\begin{figure}[H]
\centering
\includegraphics[width=\textwidth]{ant_beta_final_model_1_12_v_artigo.pdf}
\caption{Regression parameter estimates and 95\% confidence intervals by outcome and final model for ANT data.}\label{fig:antbeta1}
\end{figure}

The correlation coefficient for the COM-Poisson model is presented in \autoref{fig:corantcmp}. Among the 820 correlation coefficients returned, COM-Poisson had 333 significant correlation coefficients, NB 278, and Poisson 83. It shows that the correlation coefficient from COM-Poisson had a smaller SE than their counterparts. The stars in the graphic represent the significant coefficients at the 5\% level. The correlation patterns presented in this figure differ somewhat from the ones found in the marginal correlation. It shows the importance of calculating the correlation coefficient in a model that accounts for the effects of the linear predictor. For example, the sample correlation was nearly zero between responses 1 and 2, 2 and 3; while the COM-Poisson correlation between the random effect of these two response variables was negative and significant at 5\% level, for NB only one was significant, and for Poisson neither was significant possibly because of a high SE, once the $\hat\rho$ was equal to -.7 between the random effects of Y1 and Y2, and -.54 between the random effects of Y2 and Y3.

\begin{figure}[H]
\vspace{0.35cm}
\setlength{\abovecaptionskip}{.0001pt}
\caption{Correlogram of ant species occurrence from COM-Poisson Model. Stars represent significant correlation at 5\% level}\label{fig:corantcmp}
\includegraphics[width=0.95\textwidth]{cor_ant_cmp.pdf}
\centering
\end{figure}

\section{Discussion}

The focus of this article was to propose MGLMM for count data. The great advantage of this model is the ability to deal with multiple outcomes. We specified the framework based on a GLMM with a random intercept, where we structured a joint model whose random effect follows a multivariate normal distribution. As a result, we can measure the variance and correlation of the random effects motivated by the multivariate set of responses. The distributions used for variables of counts were Poisson, NB, and COM-Poisson. The estimation process is the same as a GLMM model. We used the \texttt{TMB} package to implement the model once it provides state-of-art \texttt{C++} libraries that handle automatic differentiation, \texttt{CppAD} [Bell BM}, 2005], linear algebra, \texttt{Eigen C++} [Guennebaud \emph{et~al.}, 2010] and parallel computation, \texttt{BLAS} [Blackford \emph{et~al.}, 2002], among others.

In order to evaluate the properties of the ML estimators, we conducted simulation studies for each distribution, considering three different values for the correlation parameter and four different sample sizes. They were all evaluated through average bias and confidence interval based on the average SE and coverage rate with a nominal level of 95\%. For Poisson distribution, we achieved unbiased and consistent estimators for large samples with intervals for bias ranging at most between (-.2;.2), except for sample size equal to 100 and values for the parameter $\rho$. The coverage rate was close to 95\% in all scenarios considered, varying between 90\% and 98\%. We evidenced a greater variability for NB compared to the Poisson distribution, especially because of the dispersion parameter $\phi$ in NB, necessary to model overdispersion. For NB distribution, we also achieved unbiased and consistent estimators for large samples with bias intervals ranging in most cases between (-.5;.5). In particular, we noticed a greater confidence interval width for the correlation and $\phi_2$ parameter estimates; and small bias for $\sigma_2$ and $\nu_2$ parameter estimates. The coverage rate in most cases was equal to or greater than 95\%; with a coverage rate below 80\% for the correlation parameter when $\rho = 0$ and the sample size equal to 100.

For the COM-Poisson model, the parameter estimators were neither consistent nor unbiased. The regression parameter was underestimated for $\beta_{02}$ and showed no bias for $\beta_{01}$. The correlation parameter $\rho$ was always estimated towards zero: when $\rho=-.5$ it was overestimated, $\rho=+.5$ was underestimated, and for $\rho=0$ it showed no bias. The standard deviation of random effect $\sigma_2$ was overestimated while $\sigma_1$ was underestimated. This behavior may be correlated to the dispersion parameter $\nu$: when more variance was captured from $\sigma_2$ less dispersion was captured from $\nu_2$; on the other hand, when less variance was captured from $\sigma_1$, more dispersion was captured from $\nu_1$.

% The dispersion parameter $\nu_2$ was overestimated whereas $\nu_1$ was underestimated in all scenarios. This bias may be due to the reparametrization used by Huang, 2017, wherein 1 out of 3 examples the dispersion parameter was overestimated by 1.6\% (1.754 compared to 1.727 original value); unfortunately, no simulation study was provided. On the other hand, there is a recently proposed mean parametrization of the COM-Poisson distribution by \citet{cmprepeduardo} that shows a simulation study where there is a small bias for the dispersion parameter in the overdispersion case for a sample size of 1000. A performance comparison between these two reparametrizations shows that \citet{cmprepeduardo} is more efficient than Huang, 2017, because the new parametrization only requires simple algebra, without adding any extra computational complexity to the original parametrization. Huang, 2017 requires an extra solution of a non-linear system that contains an infinite sum on it.

%Já tem uma soma infinita da COM-Poisson em si, para calcular o Z
% Agora, o que o Huang fez foi propor a resolução de um sistema para lambda com mais uma soma infinita para que a reparametrização na média dele funcione.

As most parameters were biased, the coverage rate was not close to 95\% in all scenarios. For $\sigma_2$, $\beta_{02}$ and $\nu_r$ the coverage rate were between 70\% and 90\%; for $\rho$ the coverage rate were close to 80\% in 2 out of 3 scenarios ($\rho=\{-0.5,0.5\}$). The other parameters, $\beta_{01}$, $\sigma_{1}$ and $\rho$ when $\rho=0$ had a coverage rate close to 95\%.

% Along with the simulation studies, some estimates did not have the SEs returned or produced extreme values, and therefore, they were discarded. In the Poisson model, it occurred in .72\% cases, for NB was in 20.56\% and COM-Poisson in 10.21\%. These percentages were calculated by the number of problems divided by the total number of simulations, which will change according to the number of replications for each distribution. For example, for COM-Poisson, the total number of simulations is 200 times 12 (the number of scenarios). It may be related because of the complexity of the model proposed model, which involves two parameters that measure variance (one for the observed variable, and another for the random effect) and the usage of LA because of a gain speed instead of a more precise numerical integration technique.

After that, the two datasets were analyzed by each model and variation of them: no correlation, fixed dispersion, fixed variance, and common variance for all random effects. The first dataset was from the AHS survey with five random variables and 5190 participants. According to the fit measures used, the COM-Poisson model provided the best fit. The SE of $\beta$ and $\rho$ estimates were similar among the COM-Poisson, NB, and Poisson models. The $\nu$ parameter captured 4 out of 5 underdispersed response variables and one overdispersed response variable, followed by a small $\sigma$. The COM-Poisson model produced more significant correlation values than its counterparts.

The second dataset was the ANT, which contains 41 response variables that count for the number of ANT species at 30 sites in Australia. The multivariate response can be considered as overdispersed by the GDI. The COM-Poisson model was also the best model regarding \texttt{AIC} and \texttt{logLik}; the model with the best \texttt{BIC} was the COM-Poisson with fixed dispersion. In almost all comparisons, the SEs of the COM-Poisson model was smaller than the other models. The $\nu$ parameter was greater than 1 for all response variables, indicating underdispersion.
% This may be because of a greater variance in the random effects, balancing these two parameters corresponding to a negative empirical correlation between them.

Therefore, we suggest using the MGLMM model framework for count data. In particular, the best results were obtained with the COM-Poisson model in two real datasets. The major advantage of it is the possibility of modeling all response variables at the same time and measuring the correlation between the random effects.

The estimation process of this model was computationally intensive, being a computational challenge to implement the model. For example, the estimation of the COM-Poisson model for AHS data was cumbersome. It took 5 days to fit in a Debian GNU/LINUX 8 (jessie) 92 GB RAM server with an AMD Opteron 6136 processor using two threads. Improving the computational implementation or even using other frameworks are examples of future work. Improving the estimation procedure for the COM-Poisson model in this context is in high demand as well as understanding whether the COM-Poisson distribution can replace the NB distribution for overdispersed data sets.

% Future work could. Choosing the best number of cores in advance is the key to the estimation being successful, as a high number of cores can cause overflow.

% \section*{Disclosure statement}
% %
% No potential competing interest was reported by the authors.
% %
% %
% \section*{Funding}
% %
% The corresponding author was funded with a scholarship provided by the Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES).

\begin{thebibliography}{39}
\newcommand{\enquote}[1]{``#1''}
\providecommand{\natexlab}[1]{#1}
\providecommand{\url}[1]{\texttt{#1}}
\providecommand{\urlprefix}{URL }
\expandafter\ifx\csname urlstyle\endcsname\relax
\providecommand{\doi}[1]{doi:\discretionary{}{}{}#1}\else
\providecommand{\doi}{doi:\discretionary{}{}{}\begingroup
\urlstyle{rm}\Url}\fi
\providecommand{\eprint}[2][]{\url{#2}}

[Aldrich \emph{et~al.}, 1997]
Aldrich J, \emph{et~al.} (1997).
\newblock \enquote{{RA} {F}isher and the making of maximum likelihood
\newblock \emph{Statistical science}, \textbf{12}(3), 162--176.

[Bates \emph{et~al.}, 2015]
Bates D, M{\"a}chler M, Bolker B, Walker S (2015).
\newblock \enquote{Fitting Linear Mixed-Effects Models Using {lme4}.}
\newblock \emph{Journal of Statistical Software}, \textbf{67}(1), 1--48.
\newblock \doi{10.18637/jss.v067.i01}.

[Baydin \emph{et~al.}, 2017]
Baydin AG, Pearlmutter BA, Radul AA, Siskind JM (2017).
\newblock \enquote{Automatic differentiation in machine learning: a survey.}
\newblock \emph{The Journal of Machine Learning Research}, \textbf{18}(1),

[Bell BM}, 2005]
{Bell BM} (2005).
\newblock \emph{CppAD: a package for C++ algorithmic differentiation}.
\newblock \urlprefix\url{http://www.coin-or.org/CppAD}.

[Blackford \emph{et~al.}, 2002]
Blackford LS, Petitet A, Pozo R, Remington K, Whaley RC, Demmel J, Dongarra J,
Duff I, Hammarling S, Henry G, \emph{et~al.} (2002).
\newblock \enquote{An updated set of basic linear algebra subprograms (BLAS).}
\newblock \emph{ACM Transactions on Mathematical Software}, \textbf{28}(2),

[Bonat, 2016]
Bonat WH (2016).
\newblock \emph{mcglm: Multivariate Covariance Generalized Linear Models}.
\newblock R package version 0.3.0,
\urlprefix\url{https://CRAN.R-project.org/package=mcglm}.

[Bonat, 2018]
Bonat WH (2018).
\newblock \enquote{Multiple Response Variables Regression Models in {R}: The
{mcglm} Package.}
\newblock \emph{Journal of Statistical Software}, \textbf{84}(4), 1--30.
\newblock \doi{10.18637/jss.v084.i04}.

[Bonat and J{\o}rgensen, 2016]
Bonat WH, J{\o}rgensen B (2016).
\newblock \enquote{Multivariate covariance generalized linear models.}
\newblock \emph{Journal of the Royal Statistical Society. Series C (Applied
Statistics)}, pp. 649--675.

[Bonat \emph{et~al.}, 2018]
Bonat WH, J{\o}rgensen B, Kokonendji CC, Hinde J, Dem{\'e}trio CG (2018).
\newblock \enquote{Extended {P}oisson--{T}weedie: Properties and regression
models for count data.}
\newblock \emph{Statistical Modelling}, \textbf{18}(1), 24--49.

[Breslow and Clayton, 1993]
Breslow NE, Clayton DG (1993).
\newblock \enquote{Approximate inference in generalized linear mixed models.}
\newblock \emph{Journal of the American statistical Association},
\textbf{88}(421), 9--25.

[Brooks \emph{et~al.}, 2017]
Brooks ME, Kristensen K, {van Benthem} KJ, Magnusson A, Berg CW, Nielsen A,
Skaug HJ, Maechler M, Bolker BM (2017).
\newblock \enquote{{glmmTMB} Balances Speed and Flexibility Among Packages for
Zero-inflated Generalized Linear Mixed Modeling.}
\newblock \emph{The R Journal}, \textbf{9}(2), 378--400.
\newblock
\urlprefix\url{https://journal.r-project.org/archive/2017/RJ-2017-066/index.html}.

[B\"{U}rkner, 2018]
B\"{U}rkner PC (2018).
\newblock \enquote{Advanced {Bayesian} Multilevel Modeling with the {R} Package
\newblock \emph{The R Journal}, \textbf{10}(1), 395--411.
\newblock \doi{10.32614/RJ-2018-017}.

[da~Silva \emph{et~al.}, 2022]
da~Silva GP, Laureano HA, Petterle RR, Júnior PJR, Bonat WH (2022).
\newblock \enquote{Multivariate generalized linear mixed models for
underdispersed count data.}
\newblock \doi{10.48550/ARXIV.2205.10486}.
\newblock \urlprefix\url{https://arxiv.org/abs/2205.10486}.

[El-Shaarawi \emph{et~al.}, 2011]
El-Shaarawi AH, Zhu R, Joe H (2011).
\newblock \enquote{Modelling species abundance using the {P}oisson--{T}weedie
\newblock \emph{Environmetrics}, \textbf{22}(2), 152--164.

[Famoye, 2015]
Famoye F (2015).
\newblock \enquote{A Multivariate Generalized {P}oisson Regression Model.}
\newblock \emph{Comm. Statist. Theory Methods}, \textbf{44}(3), 497--511.
\newblock ISSN 0361-0926.
\newblock \doi{10.1080/03610926.2012.743565}.

[Fisher, 1934]
Fisher RA (1934).
\newblock \enquote{The effect of methods of ascertainment upon the estimation
of frequencies.}
\newblock \emph{Annals of eugenics}, \textbf{6}(1), 13--25.

[Gibb \emph{et~al.}, 2015]
Gibb H, Stoklosa J, Warton DI, Brown A, Andrew NR, Cunningham S (2015).
\newblock \enquote{Does morphology predict trophic position and habitat use of
ant species and assemblages?}
\newblock \emph{Oecologia}, \textbf{177}(2), 519--531.

[Guennebaud \emph{et~al.}, 2010]
Guennebaud G, Jacob B, \emph{et~al.} (2010).
\newblock \enquote{Eigen v3.}
\newblock http://eigen.tuxfamily.org.

[Hadfield, 2010]
Hadfield JD (2010).
\newblock \enquote{MCMC Methods for Multi-Response Generalized Linear Mixed
Models: The {MCMCglmm} {R} Package.}
\newblock \emph{Journal of Statistical Software}, \textbf{33}(2), 1--22.
\newblock \urlprefix\url{http://www.jstatsoft.org/v33/i02/}.

[Huang, 2017]
Huang A (2017).
\newblock \enquote{Mean-parametrized {C}onway--{M}axwell--{P}oisson regression
models for dispersed counts.}
\newblock \emph{Statistical Modelling}, \textbf{17}(6), 359--380.

[Inouye \emph{et~al.}, 2017]
Inouye DI, Yang E, Allen GI, Ravikumar P (2017).
\newblock \enquote{A review of multivariate distributions for count data
derived from the {P}oisson distribution.}
\newblock \emph{Wiley Interdisciplinary Reviews: Computational Statistics},
\textbf{9}(3), e1398.

[J{\o}rgensen and Kokonendji, 2016]
J{\o}rgensen B, Kokonendji CC (2016).
\newblock \enquote{Discrete dispersion models and their {T}weedie asymptotics.}
\newblock \emph{AStA Advances in Statistical Analysis}, \textbf{100}(1),

[Kokonendji and Puig, 2018]
Kokonendji CC, Puig P (2018).
\newblock \enquote{Fisher dispersion index for multivariate count
distributions: A review and a new proposal.}
\newblock \emph{Journal of Multivariate Analysis}, \textbf{165}, 180--193.

[Kristensen \emph{et~al.}, 2016]
Kristensen K, Nielsen A, Berg CW, Skaug H, Bell BM (2016).
\newblock \enquote{{TMB}: Automatic Differentiation and {L}aplace
\newblock \emph{Journal of Statistical Software}, \textbf{70}(5), 1--21.
\newblock \doi{10.18637/jss.v070.i05}.

[Lee and Nelder, 1996]
Lee Y, Nelder JA (1996).
\newblock \enquote{Hierarchical Generalized Linear Models.}
\newblock \emph{Journal of the Royal Statistical Society. Series B
(Methodological)}, \textbf{58}(4), 619--678.
\newblock ISSN 00359246.
\newblock \urlprefix\url{http://www.jstor.org/stable/2346105}.

[Liang and Zeger, 1986]
Liang KY, Zeger SL (1986).
\newblock \enquote{Longitudinal data analysis using generalized linear models.}
\newblock \emph{Biometrika}, \textbf{73}(1), 13--22.

[Muñoz-Pichardo \emph{et~al.}, 2021]
Muñoz-Pichardo JM, Pino-Mejías R, García-Heras J, Ruiz-Muñoz F,
González-Regalado ML (2021).
\newblock \enquote{A multivariate {P}oisson regression model for count data.}
\newblock \emph{Journal of Applied Statistics}, \textbf{0}(0), 1--17.
\newblock \doi{10.1080/02664763.2021.1877637}.
\newblock \eprint{https://doi.org/10.1080/02664763.2021.1877637},
\urlprefix\url{https://doi.org/10.1080/02664763.2021.1877637}.

[Nelder and Wedderburn, 1972]
Nelder JA, Wedderburn RW (1972).
\newblock \enquote{Generalized linear models.}
\newblock \emph{Journal of the Royal Statistical Society: Series A (General)},
\textbf{135}(3), 370--384.

[Nikoloulopoulos and Karlis, 2009]
Nikoloulopoulos AK, Karlis D (2009).
\newblock \enquote{Modeling multivariate count data using copulas.}
\newblock \emph{Communications in Statistics-Simulation and Computation},
\textbf{39}(1), 172--187.

[Pinheiro \emph{et~al.}, 2017]
Pinheiro J, Bates D, DebRoy S, Sarkar D, {R Core Team} (2017).
\newblock \emph{{nlme}: Linear and Nonlinear Mixed Effects Models}.
\newblock R package version 3.1-131,
\urlprefix\url{https://CRAN.R-project.org/package=nlme}.

[R Core Team}, 2020]
{R Core Team} (2020).
\newblock \emph{R: A Language and Environment for Statistical Computing}.
\newblock R Foundation for Statistical Computing, Vienna, Austria.
\newblock \urlprefix\url{https://www.R-project.org/}.

[Shmueli \emph{et~al.}, 2005]
Shmueli G, Minka TP, Kadane JB, Borle S, Boatwright P (2005).
\newblock \enquote{A useful distribution for fitting discrete data: revival of
the {C}onway--{M}axwell--{P}oisson distribution.}
\newblock \emph{Journal of the Royal Statistical Society: Series C (Applied
Statistics)}, \textbf{54}(1), 127--142.

[Tierney and Kadane, 1986]
Tierney L, Kadane JB (1986).
\newblock \enquote{Accurate Approximations for Posterior Moments and Marginal
\newblock \emph{Journal of the American Statistical Association},
\textbf{81}(393), 82--86.
\newblock \doi{10.1080/01621459.1986.10478240}.
\newblock
\eprint{https://www.tandfonline.com/doi/pdf/10.1080/01621459.1986.10478240},
\urlprefix\url{https://www.tandfonline.com/doi/abs/10.1080/01621459.1986.10478240}.

[Wang \emph{et~al.}, 2020]
Wang Y, Naumann U, Eddelbuettel D, Wilshire J, Warton D (2020).
\newblock \emph{mvabund: Statistical Methods for Analysing Multivariate
Abundance Data}.
\newblock R package version 4.1.3,
\urlprefix\url{https://CRAN.R-project.org/package=mvabund}.

[Wedderburn, 1974]
Wedderburn RW (1974).
\newblock \enquote{Quasi-likelihood functions, generalized linear models, and
the Gauss—Newton method.}
\newblock \emph{Biometrika}, \textbf{61}(3), 439--447.

[Winkelmann, 2008]
Winkelmann R (2008).
\newblock \emph{Econometric analysis of count data}.
\newblock Springer Science \& Business Media.

[Winter and B{\"u}rkner, 2021]
Winter B, B{\"u}rkner PC (2021).
\newblock \enquote{Poisson regression for linguists: A tutorial introduction to
modelling count data with brms.}
\newblock \emph{Language and Linguistics Compass}, \textbf{15}(11), e12439.

[Zeileis \emph{et~al.}, 2008]
Zeileis A, Kleiber C, Jackman S (2008).
\newblock \enquote{Regression models for count data in {R}.}
\newblock \emph{Journal of statistical software}, \textbf{27}(8), 1--25.

[Zeviani \emph{et~al.}, 2014]
Zeviani WM, Ribeiro~Jr PJ, Bonat WH, Shimakura SE, Muniz JA (2014).
\newblock \enquote{The Gamma-count distribution in the analysis of experimental
underdispersed data.}
\newblock \emph{Journal of Applied Statistics}, \textbf{41}(12), 2616--2626.

\end{thebibliography}

\end{document}