# Matroid variant of Matiyasevich formula and its application

E.Yu. Lerner

###### Abstract

In 1977, Yu.V. Matiyasevich proposed a formula expressing the chromatic
polynomial of an arbitrary graph as a linear combination of flow polynomials
of subgraphs of the original graph. In this paper, we prove that this
representation is a particular case of one (easily verifiable) formula,
namely, the representation of the characteristic polynomial of an arbitrary
matroid as a linear combination of characteristic polynomials of dual
matroids. As applications, we consider an explicit expression for the flow
polynomial of a complete graph and a formula for the characteristic polynomial
of the matroid dual to the matroid of the projective geometry over a finite
field. We prove, in particular, that their major coefficients are defined by
the beginning of a certain row in the Pascal triangle. We study in detail the
connection with convolution formulas and other results for Tutte polynomials.

Mathematics Subject Classifications: 05B35, 05C31, 05B25.

Keywords: dual matroid, characteristic polynomial, flow polynomial,
Matiyasevich formula, $PG(n,q)$, Tutte polynomial, convolution formula.

##  1 The structure of the paper, the main theorem

The main goal of this paper is to provide a very simple derivation of the
formula that expresses the characteristic polynomial of an arbitrary matroid
as a linear combination of characteristic polynomials of dual matroids and to
consider applications of this formula. As particular cases, this
representation includes the Matiyasevich formula [12], which expresses the
flow polynomial of an arbitrary graph as a combination of chromatic
polynomials of its subgraphs, as well as its dual formula, which uses factor
graphs. We can treat both of these formulas as corollaries of one well-known
correlation between vacuum Feynman amplitudes over a residue ring [9], [8].
With the help of Feynman amplitudes one can also represent a flow polynomial
as the sum of Legendre symbols [8]. However, in this case, everything is much
simpler. For deriving a generalized variant of the Matiyasevich formula in a
more convenient matroid language, we need only elementary combinatorics.

In spite of the simplicity of derivation, the new formula appears to be useful
for applications. In particular, the standard approach to calculating the flow
polynomial of the complete graph $K_{n}$ does not allow one to express it
immediately in a simple form. At the same time, the expression for the
chromatic polynomial of the complete graph (and its factor graphs) is trivial.
Therefore, for calculating the flow polynomial of the complete graph it is
natural to use the approach that implies its expressing in terms of the
chromatic one. As a result, one needs only to calculate the sum with respect
to all partitions of a positive integer. According to the Hardy–Ramanujan
formula, the number of such partitions does no grow so fast; this, in
particular, allows one to easily calculate the flow polynomial of a complete
graph for dimensions, to which the usual approach is inapplicable.

An analogous situation occurs in calculating the characteristic polynomial of
the matroid which is dual to $PG(n,q)$, i. e., the matroid of an
$n$-dimensional projective space over a finite field $\mathbb{F}_{q}$. The
characteristic polynomial of the matroid $PG(n,q)$ has roots
$1,q,q^{2},\ldots,q^{n-1}$ and therefore the expression for it is evident, as
distinct from the characteristic polynomial of the dual matroid. In
factorization of the matroid $PG(n,q)$ in a nontrivial case there occur
matroids $PG(k,q)$, $k\leq n$; their number equals the $q$-binomial
coefficient ${n\choose k}_{q}$. As a result, we obtain an explicit expression
for the characteristic polynomial of the matroid dual to the matroid
$PG(n,q)$.

We describe these results in the next section, together with the Matiyasevich
formula (deduced as a corollary of the main theorem) and the main results
obtained in [11]. We observe that major coefficients of all characteristic
polynomials of matroids considered in Section 2 (matroids dual to matroids of
projective geometries, uniform matroids, bond matroids of complete graphs, i.
e., their flow polynomials) coincide (accurate to alternating signs) with
beginnings of rows of the Pascal triangle. In Section 3, we study the
connection of the main theorem with convolution formulas for Tutte polynomials
obtained earlier in papers [5] and [7]. In this section, we tried to
thoroughly systemize results obtained for Tutte polynomials in many papers;
their more concrete statements with certain additions have allowed us to
obtain formulas given in sections 1 and 2. The author of this paper is most
impressed by the paper [18] published by Tutte in 1967. As appeared, the
formulas for Tutte polynomials of complete graphs given in the mentioned paper
are the same as correlations for flow polynomials obtained by us in Section 2
of this paper. Moreover, Tutte considers the exponential generating function
for Tutte polynomials of complete graphs. This function allows one to obtain a
short expression for the flow polynomial of a complete graph in terms of the
differentiation operation.

Now we can easily implement two algorithms for calculating coefficients of the
flow polynomial of the complete graph with $n$ vertices that correspond to the
obtained formulas (see (13) and (27)); these algorithms are asymptotically
much faster than the standard procedure in Wolfram Mathematica. The efficiency
of the first algorithm slightly exceeds that of the second one; at the
beginning of Section 3 we discuss this property in more detail.

Let us now state the main theorem. By introducing the elementary
“$\zeta$-function” $\zeta_{q}(z)=\frac{1}{1-q^{-z}}$, $q\not\in\\{0,1\\}$,
$z\in\mathbb{C}$, we state it in the most symmetric form. The function
$\zeta_{q}$ satisfies the following trivial functional equations:

$\zeta_{q}(z)=-q^{z}\zeta_{q}(-z),\qquad\zeta_{q}(z)-1=-\zeta_{q}(-z).$ (1)

Recall [15] that we understand a matroid $M=(E,r)$ as the set $E$ (we consider
only the case of finite sets) with the function $r\equiv r_{M}$, which is
defined on all subsets $2^{E}$ of the set $E$, takes on nonnegative integer
values $\mathbb{N}_{0}$, and possesses the following properties:
(1) If $A\subseteq E$, then $0\leq r(A)\leq|A|$.
(2) If $A\subseteq B\subseteq E$, then $r(A)\leq r(B)$.
(3) If $A$ and $B$ are subsets of $E$, then

$r(A\cup B)+r(A\cap B)\leq r(A)+r(B).$

The function $r$ is the rank function of $M$, and $E$ is the ground set of
$M$.

The dual matroid $M^{*}=(E,r^{*})$ to the matroid $M=(E,r)$ is defined by the
rank function $r^{*}$ which is defined on the same set as $r$ and obeys the
correlation

$r^{*}(A)=r(E-A)+|A|-r(E),\quad\mbox{for any $A\subseteq E$.}$ (2)

Let $A\subseteq E$. The restriction of $M=(E,r)$ to $A$ is the matroid $M|A$
with the ground set $A$ and the rank function coinciding with $r$ on $2^{A}$.
Let $C=E-A$. The contraction of $M$ onto $A$ is the matroid $M.A\equiv M/C$
with the ground set $A$ and the rank function obeying the correlation

$r_{M.A}(B)=r_{M}(C\cup B)-r_{M}(C)\quad\mbox{for any $B\subseteq A$.}$

Evidently, if $B\subseteq A\subseteq E$, then $(M|A)|B=M|B$ and $(M.A).B=M.B$
. In addition, as is known [15],

$(M.A)^{*}=M^{*}|A,\quad(M|A)^{*}=M^{*}.A$ (3)

The characteristic polynomial $\chi_{M}$ of the matroid $M=(E,r)$ obeys the
formula

$\chi_{M}(x)=\sum_{A\subseteq E}(-1)^{|A|}x^{r(E)-r(A)}.$ (4)

Let us give some more definitions, which will be necessary in the following
sections of this paper. The flat of a matroid is its inclusion-maximal subset
$A$ of some fixed rank. This means that there exists no element $e\in E-A$
such that $r(A\cup\\{e\\})=r(A)$. We understand a link as an element $e\in E$
such that $r(\\{e\\})=0$. If the rank of an analogous subset in the matroid
$M^{*}$ equals zero, then $e$ is called a colink. One can easily make sure
that

if a matroid $M$ contains links, then $\chi_{M}(x)\equiv 0$. (5)

######  Theorem 1.

Let $M=(E,r)$ be an arbitrary matroid on the ground set $E$ with the rank
function $r$. The following formulas are valid:

$\displaystyle\chi_{M^{*}}(q)\,\zeta_{q}(-1)^{|E|}=\sum\limits_{A\subseteq
E}(-1)^{|E|-|A|}\,\frac{\chi_{M|A}(q)}{q^{r(A)}}\,\zeta_{q}(1)^{|A|},$ (6)
$\displaystyle\frac{\chi_{M^{*}}(q)}{q^{r^{*}(E)}}\zeta_{q}(1)^{|E|}=\sum\limits_{A\subseteq
E}\zeta_{q}(-1)^{|E|-|A|}\,\chi_{M/A}(q).$ (7)

###### Proof.

Let $A\subseteq E$. In accordance with formula (3) correlation (6) means that

$\chi_{M^{*}.A}(q)\,\zeta_{q}(-1)^{|A|}=\sum\limits_{B\subseteq
A}(-1)^{|A|-|B|}\,\frac{\chi_{M|A}(q)}{q^{r(B)}}\,\zeta_{q}(1)^{|B|}.$

Therefore, applying the Möbius inversion formula (in the case under
consideration it coincides with the inclusion-exclusion formula) to identity
(6), we can rewrite the latter as follows:

$\frac{\chi_{M}(q)}{q^{r(E)}}\zeta_{q}(1)^{|E|}=\sum_{A\subseteq
E}\zeta_{q}(-1)^{|A|}\,\chi_{M^{*}.A}(q).$ (8)

Replacing $A$ with $E-B$ and doing the matroid $M$ in formula (8) with $M^{*}$
(recall that $\left(M^{*}\right)^{*}=M$) we get formula (7). Therefore, it
suffices to prove only correlation (6). Performing substitution (2) and
putting $A=E-B$ in the left-hand side of formula (6), we turn it into the form

$\sum_{A\subseteq E}(-1)^{|E|-|A|}\,q^{|A|-r(A)}\zeta_{q}(-1)^{|E|}.$ (9)

We can write the right-hand side of formula (6) as follows:

$\displaystyle\sum_{A\subseteq E}q^{-r(A)}\zeta_{q}(1)^{|A|}\sum_{B:A\subseteq
A\subseteq E}\zeta_{q}(1)^{|B|-|A|}(-1)^{|E|-|B|}=$ (10) $\displaystyle=$
$\displaystyle\sum_{A\subseteq
E}q^{-r(A)}\zeta_{q}(1)^{|A|}(\zeta_{q}(1)-1)^{|E|-|A|}.$

Substituting correlations (1) in (10), we conclude that formula (10) is
equivalent to (9). ∎

Note that the first of correlations (1) and the equality $r^{*}(E)=|E|-r(E)$
allow us to rewrite formula (7) in a form, which is more convenient for
calculations.

######  Corollary 1.

The following formula is valid:

$\chi_{M^{*}}(x)=\frac{(-1)^{|E|}}{x^{r(E)}}\sum\limits_{A\subseteq
E}(1-x)^{|A|}\,\chi_{M/A}(x).$ (11)

Just this corollary we use in most applications.

##  2 Applications of Theorem 1 and Corollary 1

Let $G=(V,E)$ be an undirected graph with the vertex set $V$ and the edge set
$E$; $q\in\mathbb{N}$. The number of ways of coloring vertices of the graph in
$q$ colors so as to make adjacent vertices of the graph have various colors is
defined by the chromatic polynomial of the graph $P_{G}(q)$. The polynomial
$P_{G}(q)$ is a multiple of $q^{c(G)}$, where $c(G)$ is the number of
connectivity components of the graph $G$. Its normalized version
$P^{\prime}_{G}(q)=P_{G}(q)/q^{c(G)}$ coincides with the characteristic
polynomial of the so-called graphic matroid $M_{G}=(E,r)$. The rank function
$M_{G}$ obeys the formula $r(A)=|V(A)|-c(A)$, where $|V(A)|$ and $c(A)$ are,
correspondingly, the number of vertices and the number of connectivity
components of the subgraph induced by the subset of edges $A$, $A\subseteq E$.

The characteristic polynomial of the dual matroid $M^{*}_{G}$ is the so-called
flow polynomial $F_{G}$ [17]. Let us fix some orientation of edges in $E$ and
denote by the symbol $\mathbb{Z}_{q}$ the additive group of integers modulo
$q$. In what follows, we understand a nonzero flow as a correspondence between
each directed edge $e\in E$ and a nonzero element $k_{e}$ in $\mathbb{Z}_{q}$
such that for each vertex $v\in V$ the sum $\sum k_{e}$ for all $e$ that enter
in $v$ coincides with the analogous sum for all $e$ outgoing $v$. The number
of everywhere nonzero flows is independent of the orientation of edges; it
equals the value of the flow polynomial $F_{G}(q)$.

In the case, when $M=M_{G}$, formula (8) takes the form

$\frac{P_{G}(q)}{q^{|V|}}\zeta_{q}(1)^{|E|}=\sum_{A\subseteq
E}\zeta_{q}(-1)^{|A|}F_{H(A)}(q),$

where $H(A)$ is the subgraph of the graph $G$ induced by the subset $A$ of
edges of this graph. The latter correlation coincides with the Matiyasevich
formula [12] accurate to denotations.

Correlation (6) implies the formula that is inverse with respect to the above
one, namely,

$F_{G}(q){\zeta_{q}(-1)}^{|E|}=\sum_{A\subseteq
E}(q-1)^{|E|-|A|}\frac{P_{H(A)}(q)}{q^{|V(A)|}}\zeta_{q}(1)^{|A|},$

where $|V(A)|$ is the number of vertices of the graph $H(A)$. See this
representation of the flow polynomial as a linear combination of chromatic
polynomials at the end part of the paper [11] (see also formula (2.5) in [4]).

We can obtain a more interesting property, substituting the matroid $M=M_{G}$
in correlation (11). We can express matroids $M/A$ in terms of factor graphs
of the graph $G$ with respect to subgraphs $H(A)$. Moreover, it is evident
that if a factor graph $G/H$ contains links then $P_{G/H}(q)\equiv 0$. As a
result, we obtain a slightly modified statement of theorem 2 in the paper
[11].

######  Theorem 2.

Assume that the set of vertices $V$ is divided into subsets
$V_{1},\ldots,V_{m}$, where $m$ is an arbitrary number, but $V_{i}$ are such
that $m$ subgraphs formed by all edges in $E$, whose both endpoints belong to
one and the same subset $V_{i}$, are connected. Denote by $A$ the set of all
these edges and do by $H(A)$ the union of all subgraphs, correspondingly. Then
the following formula is valid:

$F_{G}(x)=\frac{(-1)^{|E|}}{x^{|V|}}\sum_{H(A)}(1-x)^{|A|}P_{G/H(A)}(x);$ (12)

here the sum is calculated over all subgraphs of $H$ (all such partitions of
the set $V$).

Here, similarly to formula (11) given above, and in what follows, we use the
“real” variable $x$ instead of $q$ (we will need the symbol $q$ later for the
standard denotation of the number of elements in a finite field, though
earlier it was convenient for us to use this variable as a subscript for the
“zeta-function”). Let us now apply formula (12) in the case, when the graph
$G$ is $K_{n}$, i. e., a complete graph with $n$ vertices. Let us first recall
some assertions of the partition theory [1].

A partition $\lambda$ of a positive integer $n$ is its representation as the
sum of positive integer terms (two sums, which differ in the order of addends,
are considered as equal ones). The total number of partitions $p(n)$ satisfies
the well-known formula derived by Hardy and Ramanujan (and refined by
Rademacher). The asymptotics for $p(n)$ is subexponential, in particular,
$p(50)$ is only $204\,226$. At the same time, it is impossible to calculate
(using a computer) $F_{K_{n}}(x)$ with $n=50$ with the help of standard
techniques (which are not specific for complete graphs), because the
computation time grows faster than $\exp(n)$. Formula (12) allows us to reduce
the calculation of $F_{K_{n}}(x)$ to summing $p(n)$ terms and to solve the
problem of calculating $F_{K_{n}}(x)$ for $n\sim 50$.

We can represent the partition $\lambda$ of the number $n$ in the form
$(\lambda_{1},\ldots,\lambda_{m})$, where

$\lambda_{1}\leq\lambda_{2}\leq\ldots\leq\lambda_{m},\quad\sum\nolimits_{i=1}^{m}\lambda_{i}=n,$

or in the form $(1^{n_{1}}2^{n_{2}}3^{n_{3}},\ldots)$, where $n_{i}$ is the
number of values $i$ among the collection $(\lambda_{1},\ldots,\lambda_{m})$,
correspondingly, $\sum_{i}i\,n_{i}=n$. The number of partition elements $m$ is
called the length of the partition; we denote it as $\ell(\lambda)$. Let

$\displaystyle\lambda!$ $\displaystyle=$
$\displaystyle\lambda_{1}!\lambda_{2}!\lambda_{3}!\ldots=1!^{n_{1}}2!^{n_{2}}3!^{n_{3}}\ldots,$
$\displaystyle|\lambda|!$ $\displaystyle=$ $\displaystyle
n_{1}!n_{2}!n_{3}!\ldots,\qquad f(\lambda)=\lambda!\times|\lambda|!.$

The number $f(\lambda)$ is called the Faá di Bruno coefficient of the
partition $\lambda$. We need the following assertion (see [1, theorem 13.2])
(one can prove it rather easily).

######  Lemma 1.

The number of ways for dividing a set of $n$ distinct elements into parts,
whose powers are defined by the partition $\lambda$, equals $f(\lambda)$.

######  Theorem 3.

Let $s(\lambda)=\sum_{i=1}^{\ell(\lambda)}\lambda_{i}(\lambda_{i}-1)/2$. The
following formula is valid:

$F_{K_{n}}(x)=\frac{(-1)^{\frac{n(n-1)}{2}}}{x^{n}}\sum_{\lambda}f(\lambda)(1-x)^{s(\lambda)}x(x-1)\ldots(x-\ell(\lambda)+1);$
(13)

here the sum is calculated over all partitions $\lambda$ of the number $n$.

###### Proof.

Let us apply formula (12) for $G=K_{n}$. All terms in the right-hand side for
various representations of $V$ as the union of disjoint parts
$V_{1},V_{2},\ldots,V_{m}$ are the same, provided that their powers form one
and the same partition $\lambda$ of the number $n$. By Lemma 1 the quantity of
such representations equals $f(\lambda)$. The number of edges in the subgraph
$H$ such that each their connectivity component represents a complete subgraph
with $\lambda_{i}$ vertices obeys the evident formula

$\sum\nolimits_{i=1}^{\ell(\lambda)}\frac{\lambda_{i}(\lambda_{i}-1)}{2};\quad\mbox{and}\
P_{K_{m}}(x)=x(x-1)\ldots(x-m+1).$

∎

Formula (13) allows us to easily calculate, for example, the value
$F_{K_{50}}(x)$.

######  Corollary 2.

Let $N=n(n-1)/2$ and $N^{\prime}=(n-1)(n-2)/2$. The major $n-2$ powers of the
flow polynomial $F_{K_{n}}(x)$ take the form

$x^{N^{\prime}}-Nx^{N^{\prime}-1}+{N\choose
2}x^{N^{\prime}-2}-\ldots+(-1)^{n-3}{N\choose n-2}x^{N^{\prime}-n+3}.$

###### Proof.

Consider the function $g(\lambda)=s(\lambda)+\ell(\lambda)$, where $\lambda$
is a partition of a fixed number. With $\ell(\lambda)=2$ the value of the
function $g$ decreases as so does the value $|\lambda_{1}-\lambda_{2}|$.
Moreover, with $\ell(\lambda)=1$ its value is greater than that with
$\ell(\lambda)=2$. In view of this fact we conclude that with
$\ell(\lambda)\not\in\\{1,2\\}$ we can increase the value of $g(\lambda)$ by
“merging together” two elements of the partition. As a result, the largest
value of $g(\lambda)$ is attained with $\lambda^{\prime}$ such that
$\ell(\lambda^{\prime})=1$, while the second largest one is attained with
$\lambda^{\prime\prime}=(1,n-1)$, $\ell(\lambda^{\prime\prime})=2$.

By Theorem 3 the value $(g(\lambda)-n)$ is the power of the term that
corresponds to the partition $\lambda$ in sum (13). Using calculations
performed in the previous paragraph, we conclude that coefficients of
$F_{K_{n}}(x)$ coincide with coefficients at major powers $(x-1)^{N}$,
provided that the number of major powers equals

$g(\lambda^{\prime})-g(\lambda^{\prime\prime})=(N+1)-(N^{\prime}+2)=n-2.$

Here the power of the polynomial $F_{K_{n}}(x)$ equals
$g(\lambda^{\prime})-n=N+1-n=N^{\prime}$. ∎

As an illustration, let us consider the result of the calculation of the flow
polynomial of a complete graph with 5 vertices

$F_{K_{5}}(x)=\mathbf{1}\,x^{6}-\mathbf{10}\,x^{5}+\mathbf{45}\,x^{4}-115x^{3}+175x^{2}-147x+51.$

Let us now discuss the application of the main theorem in cases when the
matroid $M$ is not, generally speaking, a graphic one. Recall that [15] the
uniform matroid $U_{m,n}$, $m\leq n$, $m,n\in{\mathbb{N}}_{0}$, is defined on
the ground set $E$, $|E|=n$, by the rank function $r$ as follows:

$\mbox{with $A\in 2^{E}$ we get the equality $r(A)=\min\\{m,|A|\\}$}.$

Note that $U^{*}_{m,n}=U_{n-m,n}$. One can easily verify that for $A\subseteq
E$, $|A|=k$,

$U_{m,n}|A=U_{\min\\{m,k\\},k}\quad\text{and}\quad
U_{m,n}.A=U_{\max\\{0,m-n+k\\},k}.$ (14)

We directly calculate the characteristic polynomial (4) of the matroid
$U_{m,n}$: $\chi_{U_{0,0}}=1$, in the rest cases (with all other values $m\leq
n$) the following formula is valid:

$\chi_{U_{m,n}}(x)=\sum_{i=0}^{m-1}(-1)^{i}{n\choose i}(x^{m-i}-1).$

In this case, sums over subsets of $E$ of power $k$ in formulas in Theorem 1
and Corollary 1 give ${n\choose k}$ equal terms. As a result, we get some
identities for binomial coefficients. In the next section, we consider a more
elegant and at the same time elementary generalization of these identities
(see Remark 1).

Let us thoroughly study the application of formula (11) in the case of finite
projective geometries, i. e., projective spaces over a finite field
${\mathbb{F}}_{q}$. Recall that elements of the ground set $E$ of the matroid
$PG(n-1,q)$ are all distinct points of the projective space of dimension $n-1$
over the field ${\mathbb{F}}_{q}$ (see details in [15, Ch. 6]). The rank of
the subset $A$ of the set $E$ is the dimension of the linear space spanned by
$A$. One can easily make sure that $|E|=\frac{q^{n}-1}{q-1}={n\choose 1}_{q}$;
here the Gaussian binomial coefficient ${n\choose k}_{q}$,
$n,k\in\mathbb{N}_{0}$, obeys the formula

${n\choose
k}_{q}=\frac{(q^{n}-1)(q^{n-1}-1)\cdots(q^{n-k+1}-1)}{(q^{k}-1)(q^{k-1}-1)\cdots(q-1)}.$

The following assertion is a generalization of the remark about the number of
elements in the ground set ([15, Proposition 6.1.4]):

The number of rank-$k$ flats in $PG(n-1,q)$ equals ${n\choose k}_{q}$. (15)

In addition, it is well known that (see [21, Proposition 7.5.3] and references
to fundamental works by Richard Stanley mentioned there) the characteristic
polynomial of the matroid $PG(n-1,q)$ takes the form

$\chi_{PG(n-1,q)}(x)=(x-1)(x-q)(x-q^{2})\cdots(x-q^{n-1}).$ (16)

The factor space of a linear space of dimension $n$ with respect to its
subspace of dimension $k$, $k<n$, is its subspace of dimension $n-k$
orthogonal to the subspace, which we use for the factorization. In terms of
matroids, this means that if $A,A\subseteq E$, is the rank-$k$ flat of
$PG(n-1,q)$, then the matroid $PG(n-1,q)/A$ is isomorphic to $PG(n-k-1,q)$.
But if $A$ is not a flat of $PG(n-1,q)$, then the matroid $PG(n-1,q)/A$
contains links. Therefore, by using formulas (5), (15), and (16), we conclude
that Corollary 1 in the case under consideration can be formulated as follows:

######  Theorem 4.

The characteristic polynomial of the matroid dual to $PG(n-1,q)$ obeys the
formula

$\chi_{PG(n-1,q)^{*}}(x)=\frac{(-1)^{\frac{q^{n}-1}{q-1}}}{x^{n}}\sum_{k=0}^{n}{n\choose
k}_{q}\ (1-x)^{\frac{q^{k}-1}{q-1}}\,\prod_{i=0}^{n-k-1}(x-q^{i}).$

Here the empty product equals one. Correspondingly, the latter term in the sum
equals $(1-x)^{\frac{q^{n}-1}{q-1}}$, while the second largest one is the
polynomial of $x$ of the power $\frac{q^{n-1}-1}{q-1}+1$. Hence, we get the
following analog of Corollary 2:

######  Corollary 3.

Let $N=\frac{q^{n}-1}{q-1}$, $N^{\prime}=N-n$, $k=q^{n-1}-1$. The $k$ major
powers of the polynomial $\chi_{PG(n-1,q)^{*}}(x)$ take the form

$x^{N^{\prime}}-Nx^{N^{\prime}-1}+{N\choose
2}x^{N^{\prime}-2}-\ldots+(-1)^{k-1}{N\choose k}x^{N^{\prime}-k+1}.$

For example, the matroid dual to the Fano matroid $F_{7}=PG(2,2)$ satisfies
the equality

$\chi_{F_{7}^{*}}(x)=\mathbf{1}\,x^{4}-\mathbf{7}\,x^{3}+\mathbf{21}\,x^{2}-28x+13.$

##  3 The Tutte polynomial and the matroid variant of the Matiyasevich formula

In this section, we prove that one can obtain formulas given in sections 1 and
2 by refining results obtained for the Tutte polynomial. The author of this
paper is most impressed by paper [18], where Tutte generalizes Theorem 3.
Moreover, in the mentioned paper, Tutte also obtains the exponential
generating function for Tutte polynomials of complete graphs. It implies that
one can get the flow polynomial of the complete graph by multiple
differentiations (see Corollary 4 below). The mathematical formula in this
case is more compact than formula (13) in Theorem 3 (one can also obtain the
latter formula as a corollary of the Tutte decomposition of the Tutte
polynomial of the complete graph $K_{n}$ in all partitions $\lambda$ of the
number $n$). Note that the direct calculation of the flow polynomial of a
complete graph by means of differentiation (using a computer) requires
slightly more time than its calculation by decomposition in all partitions (in
essence, both techniques are equivalent, they differ only in computation
aspects and intermediate results storing methods; the latter are not specified
in the shorter procedure which implies the use of the computation technique
mentioned in Corollary 4).

See fundamentals of the theory of Tutte polynomials, in particular, in [20],
[3], [21], and in the book [19] by Tutte (in the case of graphic matroids).
Specific properties of characteristic polynomials that are of interest for us
are described in classical papers [17] and [21]. The review [13] has appeared
to be very useful for us as a reference book of formulas for various Tutte
polynomials.

Recall that the Tutte polynomial of the matroid $M=(E,r)$ obeys the formula

$T_{M}(x,y)=\sum_{A\subseteq E}(x-1)^{r(E)-r(A)}(y-1)^{|E|-r(A)}.$

By definition, $T_{M^{*}}(x,y)=T_{M}(y,x)$. In addition,

$\chi_{M}(z)=(-1)^{r(E)}T_{M}(1-z,0),\qquad\chi_{M^{*}}(z)=(-1)^{|E|-r(E)}T_{M}(0,1-z).$
(17)

In certain (exceptional) cases one can easily calculate the Tutte polynomial.
For example, in the case of a uniform matroid $U_{m,n}$, it is evident that

$T_{U_{m,n}}(x,y)=\sum_{i=0}^{m-1}{n\choose i}(x-1)^{m-i}+{n\choose
m}+\sum_{i=m+1}^{n}{n\choose i}(y-1)^{i-m}.$

In a general case, for any matroid with $(x,y)\in H_{1}$, where $H_{1}$ is the
so-called critical hyperbola $(x-1)(y-1)=1$, the following formula is valid:

$\displaystyle T_{M}(x,y(x))$ $\displaystyle=$ $\displaystyle\sum_{A\subseteq
E}(x-1)^{r(E)-r(A)}(x-1)^{r(A)-|A|}=\sum_{A\subseteq
E}\frac{(x-1)^{r(E)}}{(x-1)^{|A|}}=$ (18) $\displaystyle=$
$\displaystyle(x-1)^{r(E)}\left(1+\frac{1}{(x-1)^{|E|}}\right)=x^{|E|}(x-1)^{r(E)-|E|}.$

If the ground set $E$, $|E|=n$, consists of $i$ links and $j$ colinks (the
base case), then $T_{M}(x,y)=x^{j}y^{i}$, $i+j=n$. In particular, in this
case, the characteristic polynomial satisfies the correlation

$\chi_{M}(z)=z^{n}\quad\mbox{if
$j=n$;}\qquad\chi_{M}(z)=0\quad\mbox{otherwise}.$ (19)

Assume that an edge $e$ is neither a link nor a colink, and let $A=E-e$
(hereinafter we identify the set consisting of one edge and the edge itself).
Denote the matroid $M|A$ as $M\setminus e$. Any matroid and any $e\in E$
different from a colink satisfy the formula $T_{M}(x,y)=T_{M\setminus
e}(x,y)+T_{M/e}(x,y)$ (the Tutte polynomial is the Tutte–Grothendieck
invariant).

The latter correlation, together with the base case, uniquely defines
$T_{M}(x,y)$ for any $M$. In the case of the characteristic polynomial, the
recurrent correlation for it takes the form

$\chi_{M}(z)=\chi_{M\setminus e}(z)-\chi_{M/e}(z).$ (20)

Correlations (19) and (20) allow us to prove formula (6) by induction with
respect to the number of edges different from colinks. Below we discuss a more
natural proof of Theorem 1 which is based on the application of the
convolution formula for Tutte polynomials. This important correlation for
Tutte polynomials established by W. Kook, V. Reiner, and D. Atanton allows
several generalizations.

######  Proposition 1 ([5]).

The following formula is valid:

$T_{M}(x,y)=\sum_{A\subseteq E}T_{M|A}(0,y)T_{M/A}(x,0).$ (21)

For example, for uniform matroids, the following remark is valid.

######  Remark 1.

Let $n\in\mathbb{N}$, $m\in\mathbb{N}_{0}$, $m\leq n$. For the matroid
$M=U_{m,n}$ identity (21) takes the form

$T_{M}(x,y)=T_{M}(x,0)+T_{M}(0,y).$

###### Proof.

Let us use formula (14) and the fact that the sign alternating sum of binomial
coefficients equals zero. With $k\in\mathbb{N}$ we get equalities

$T_{U_{k,k}}(0,y)=0;\qquad T_{U_{k,0}}(x,0)=0.$

Therefore, with $0<|A|\leq m$ we get the equality $T_{U_{m,n}|A}(0,y)=0$, and
with $m\leq|A|<n$ we conclude that $T_{U_{m,n}/A}(x,0)=0$. Therefore, among
all terms in the sum in the right-hand side of the equality, only those ones
that correspond to cases $A=\emptyset$ and $A=E$ differ from zero. ∎

Certainly, one can easily immediately verify the formula in Remark 1 without
referring to Proposition 1.

We need one more notion which is closely connected with the Tutte polynomial.
Recall that the Whitney rank generating function $R_{M}(u,v)$ of the matroid
$M=(E,r)$ obeys the formula

$R_{M}(u,v)=\sum_{A\subseteq E}u^{r(E)-r(A)}v^{|A|-r(A)},$

i. e., $R_{M}(u-1,v-1)=T_{M}(u,v)$. Evidently,

$\chi_{M}(x)=(-1)^{r(E)}R_{M}(-x,-1),\quad\chi_{M^{*}}(x)=(-1)^{|E|-r(E)}R_{M}(-1,-x).$
(22)

We can rewrite formula (18) as follows:

$R_{M}(x,1/x)=(x+1)^{|E|}x^{r(E)-|E|}.$ (23)

In the paper [16], V. Reiner generalizes equality (21) to a certain two-
parameter family of equalities, where formula (21) is a limit correlation. In
[6], J. P. A. Kung establishes some two-parameter identities for
characteristic polynomials of matroids. Finally, in [7] J. P. S. Kung obtains
several identities which generalize identities established by him earlier and
the convolution formula (21). Identity 3 in this paper takes the form

$R_{M}(\lambda\xi,xy)=\sum_{A\subseteq
E}\lambda^{r(E)-r(A)}(-y)^{|A|-r(A)}R_{M|A}(-\lambda,-x)R_{M/A}(\xi,y)$ (24)

(we use denotations introduced by J. P. S. Kung). Now we represent formulas
(6) and (7) as particular cases of identity (24).

###### Proof of Theorem 1 as a corollary of equality (24).

Putting $xy=-q$ and $\lambda\xi=-1$, in accordance with (22) in the left-hand
side of equality (24) we obtain the desired expression
$(-1)^{|E|-r(E)}\chi_{M^{*}}(q)$. For establishing (6) it suffices to
concretize $x=1$, $\lambda=q$ and, correspondingly, $y=-q$, $\xi=-1/q$; then
to perform substitutions (22) and (23), and sum degrees of $q$ and $(-1)$.
Formula (7) can be obtained analogously, provided that $y=-1$, $\xi=-q$ and,
correspondingly, $x=q$ and $\lambda=1/q$. ∎

Note that Theorem 1 without referring to Matiyasevich (the author of this
paper was not aware of his works) was first stated by us in preprint [10] (see
theorem 2 therein). The proof described above is adduced in the mentioned
preprint.

All theorems in Section 2 are more elementary corollaries of results on Tutte
polynomials which were obtained earlier. We deduce Theorem 4 by substituting
values $u=0$ and $v=1-x$ in the formula for the Tutte polynomial in the
projective geometry

$T_{PG(n-1,q)}(u,v)=\frac{1}{(u-1)^{n}}\sum_{k=0}^{n}{n\choose k}_{q}\
v^{\frac{q^{k}-1}{q-1}}\,\prod_{i=0}^{n-k-1}((u-1)(v-1)-q^{i})$ (25)

and by multiplying it by $(-1)^{|E|-n}$ in accordance with correlation (17).
Expression (25) is established in paper [14].

Another way to derive the same expression (25) together with the exponential
generating function for Tutte polynomials of projective geometries over the
field $\mathbb{F}_{q}$ is proposed in the paper [2]. Evidently, the value of
the $n$th derivative of the exponential generating function at zero coincides
with the $n$th term of the series. Below we apply such an approach for
deriving an explicit formula for the characteristic polynomial, calculating
the flow polynomial of the complete graph in a different way. The technique
used by us is based on results obtained by Tutte in [18].

In the mentioned paper, Tutte does not immediately consider the flow
polynomial, but he elegantly treats the polynomial of two variables
$Q_{G}(u,v)$ which is closely connected with the chromatic polynomial
$P_{G}(x)$. One can express this dichromatic Tutte polynomial $Q_{G}(u,v)$ in
terms of the Whitney rank function as follows:

$Q_{G}(u,v)=u^{c(G)}R_{M_{G}}(u,v);$

here $c(G)$ is the number of connectivity components of the graph $G$. For the
dichromatic polynomial of a complete graph, Tutte has obtained an equation
([18, equation (16)]), which in terms used in Section 2 takes the form

$Q_{K_{n}}(u,v)=\frac{1}{v^{n}}\sum_{\lambda}f(\lambda)(1+v)^{s(\lambda)}uv(uv-1)\ldots(uv-\ell(\lambda)+1).$

By using the second formula in (22) we obtain the statement of Theorem 3 as a
particular case of the above correlation.

The previous correlation was strengthened by Tutte as the following formula
for the exponential generating function for $Q_{K_{n}}(u,v)v^{n}$ ([18,
equation (17)]):

$\sum_{n=0}^{\infty}\frac{z^{n}v^{n}}{n!}Q_{K_{n}}(u,v)=\Bigl{\\{}\sum_{n=0}^{\infty}\frac{z^{n}}{n!}(1+v)^{\frac{1}{2}n(n-1)}\Bigr{\\}}^{uv}.$
(26)

######  Corollary 4 (Corollary of the Tutte formula (26)).

The flow polynomial of a complete graph obeys the formula

$F_{K_{n}}(x)=\frac{(-1)^{\frac{1}{2}n(n-1)}}{x^{n}}\left.\left(\frac{\partial^{n}}{\partial
z^{n}}\Bigl{\\{}\sum_{i=0}^{n}\frac{z^{i}}{i!}(1-x)^{\frac{1}{2}i(i-1)}\Bigr{\\}}^{x}\,\right)\right|_{z=0}.$
(27)

###### Proof.

Substituting (22) in formula (26), we get the correlation

$\sum_{i=0}^{\infty}\frac{z^{i}x^{i}(-1)^{\frac{1}{2}i(i-1)}}{i!}F_{K_{i}}(x)=\Bigl{\\{}\sum_{i=0}^{\infty}\frac{z^{i}}{i!}(1-x)^{\frac{1}{2}i(i-1)}\Bigr{\\}}^{x}.$

By applying the operator $\left.\left(\frac{\partial^{n}}{\partial
z^{n}}\,\cdot\,\right)\right|_{z=0}$ to both sides of the latter equality, we
obtain formula (27). ∎

As we have mentioned, the major coefficients of characteristic polynomials of
all matroids $M^{*}_{K_{n}}$, $U_{m,n}$, and $PG(n,q)^{*}$ considered in
Section 2 coincide (accurate to alternating signs) with the beginning of the
row of the Pascal triangle, whose number equals the quantity of elements in
the corresponding matroid. We have established this property as a corollary of
the matroid form of the Matiyasevich formula. The classical assertion about
combinatorial sense of coefficients $\chi_{M}$ stated by G. C. Rota [17]
reveals a deeper reason for the result mentioned above. Namely, the reason is
that the complete graph $K_{n}$, the projective geometry $PG(n,q)$, and
especially the uniform matroid have “rather thick cuts” (cocycles). Let us
discuss this property.

Recall that the set $A$: $A\subseteq E$ of the matroid $M=(E,r)$ is said to be
independent, if $r(A)=|A|$. If, in addition, $r(A)=r(E)$, then the independent
set is called a base. The minimal set, which is not independent, is called a
cycle of the matroid $M$. A cycle of the matroid $M^{*}$ is called a cocycle.
Let us associate all elements $e$, $e\in E$, with distinct positive integers
$l_{e}$. A cycle $C$ without an element $e^{\prime}$ such that
$l_{e^{\prime}}=\max_{e\in C}l_{e}$ is called a broken circuit. In [17], G. C.
Rota has obtained the following result ([17, P. 359, Corollary 1]):

######  Proposition 2.

Assume that the characteristic polynomial of the matroid $M=(E,r)$ takes the
form

$\chi_{M}(x)=x^{r(E)}+m_{1}x^{r(E)-1}+m_{2}x^{r(E)-2}+\cdots+m_{r(E)}.$

Then $(-1)^{k}m_{k}$ is the number of independent subsets $E$ consisting of
$k$ elements without broken circuits.

As a result, we obtain the following assertion. If any cocycle of the matroid
$M$ contains at least $t$ elements, then for any $k\leq t-2$ the absolute
value $(-1)^{k}m^{\prime}_{k}$ of the $k$th coefficient of the characteristic
polynomial $\chi_{M^{*}}$

$\chi_{M^{*}}(x)=x^{r^{*}(E)}+m^{\prime}_{1}x^{r^{*}(E)-1}+\cdots+m^{\prime}_{r^{*}(E)}\quad\mbox{(here
$r^{*}(E)=|E|-r(E)$)}$

coincides with ${|E|\choose k}$. Really, in this case, any subset consisting
of $k$ elements of the set $|E|$ is necessarily independent in the matroid
$M^{*}$ and contains no broken circuit.

Since any cut in the complete graph $K_{n}$ consists of at least ${n-1}$
edges, the mentioned property implies Corollary 2.

We can use the same property for finding coefficients of the characteristic
polynomial of a uniform matroid or the dual one. Recall that the latter also
is a uniform matroid.

Really, any cycle of a uniform matroid $U_{m,n}$, by definition, contains
$m+1$ elements, and the power of its characteristic polynomial equals $m$.
Consequently, all coefficients of the characteristic polynomial, except,
possibly, the free term, coincide (accurate to alternating signs) with the
beginning of the $n$th row in the Pascal triangle.

Finally, for describing cocycles in $PG(n-1,q)$ we can make use of the fact
that a cocycle is the minimal set of elements which has at least one common
element with any base. Hence, one can easily understand cocycles of the
matroid $PG(n-1,q)$ geometrically. Namely, they represent sets of points in
the projective space which belong to some hyperplane ${\mathbb{F}}_{q}^{n}$
(which does not contain the origin of coordinates). The number of points in
any hyperplane ${\mathbb{F}}_{q}^{n}$ equals $q^{n-1}$. Hence we get Corollary
3.

##  4 Conclusion

In this paper, we propose a matroid form of one variant of the Matiyasevich
formula and describe its application for deriving characteristic polynomials
of bond matroids and matroids dual to matroids of finite projective
geometries. One can deduce most results from correlations for the
corresponding Tutte polynomials. However, the use of one variable in the case,
when two of them are not necessary, seems to give an easier way to obtain such
formulas. On the other hand, this way allows us to better understand (in a
sense) the interconnection of two variables in the Tutte polynomial. It also
allows us to clarify the structure of major powers of considered polynomials.

In what follows, we are going to consider “conditional” flow and chromatic
polynomials and their generalizations. In the simplest case, we are interested
in the number of regular colorings of the graph, when colors of some vertices
are fixed a priori. For calculating this number, similarly to the conditional
probability formula, we can write down the ratio of two chromatic polynomials,
which explains the terminology used by us. It makes sense to consider the
matroid variant of such correlations.

### Acknowledgements

The research was supported by RSF (project No. 24-21-00158).

## References

  * [1] G. E. Andrews, The Theory of Partitions, Cambridge Univ. Press, Cambridge, 1976.
  * [2] M. Barany and V. Reiner, The Tutte polynomial of a finite projective space, (2005) preprint.
  * [3] T. Brylawski and J. G. Oxley, The Tutte polynomial and its applications, in Matroid Applications (ed. N. White), Encyclopedia of Mathematics and Its Applications, 40, Cambridge Univ. Press, 1992.
  * [4] B. Bychkov, A. Kazakov, and D. Talalaev, Functional Relations on Anisotropic Potts Models: from Biggs Formula to the Tetrahedron Equation, SIGMA 17, (2021) 035. DOI: 10.3842/SIGMA.2021.035
  * [5] W. Kook, V. Reiner and D. Stanton, A convolution formula for the Tutte polynomial, J. Combin. Theory Ser. B 76 (1999), no. 2, 297–300. MR 1699230
  * [6] J. P. S. Kung, A multiplication identity for characteristic polynomials of matroids, Adv. in Appl. Math. 32 (2004), no. 1-2, 319–326. MR 2037633
  * [7] J. P. S. Kung, Convolution-multiplication identities for Tutte polynomials of graphs and matroids, J. Combin. Theory Ser. B 100 (2010), no. 6, 617–624. MR 2718681
  * [8] A. P. Kuptsov, E. Yu. Lerner, S. A. Mukhamedjanova, Flow polynomials as Feynman amplitudes and their $\alpha$-representation, Electron. J. Combin., 24 (2017), no. 1, paper 11, 19 pp. DOI: 10.37236/6396
  * [9] E. Yu. Lerner, Feynman integrals of $p$-adic argument in the momentum space. II. Explicit expressions, Teoret. Mat. Fiz. 104 (1995), no. 3, 371–392; translation in Theoret. and Math. Phys. 104 (1995), no. 3, 1061–1077 (1996). MR 1607005
  * [10] E. Yu. Lerner, The $\alpha$-representation for the characteristic function of a matroid, arXiv:1611.02746, DOI: 10.48550/arXiv.1611.02746 (2016)
  * [11] E. Yu. Lerner and S. A. Mukhamedjanova, Matiyasevich Formula for Chromatic and Flow Polynomials and Feynman Amplitudes, Lobachevskii Journal of Mathematics 43 (2022), no. 12, 3552–3561 MR 4567102
  * [12] Yu. V. Matiyasevich, On a certain representation of the chromatic polynomial, Diskretnyi Analiz 91, issue 31, 61–70 (1977). MR 543806, arXiv:0903.1213
  * [13] C. Merino, M. Ramirez-Ibañez, G. Rodriguez-Sànchez The Tutte polynomial of some matroids, 2012 (2012) Article ID 430859 DOI: 10.1155/2012/430859
  * [14] E. G. Mphako. Tutte Polynomials of Perfect Matroid Designs, Combinatorics, Probability and Computing bf 9 (2000) 363–367 DOI: 10.1017/S0963548300004338
  * [15] J. G. Oxley, Matroid Theory, Oxford University Press, Oxford, New York, 1992. MR 1207587
  * [16] V. Reiner, An interpretation for the Tutte polynomial, European J. Combin. 20 (1999), no. 2, 149–161 MR 1676189
  * [17] G. C. Rota On the foundations of combinatorial theory I. Theory of Möbius Functions. Z. Wahrscheinlichkeitstheorie verw Gebiete 2 (1964) 340–-368. DOI: 10.1007/BF00531932
  * [18] W. T. Tutte, On Dichromatic Polynomials, Journal of Combinatorial Theory, Series A, 2 (1967), 301–320. DOI: 10.1016/S0021-9800(67)80032-2
  * [19] W. T. Tutte, Graph theory, reprint of the 1984 original, Encyclopedia of Mathematics and its Applications, 21, Cambridge Univ. Press, Cambridge, 2001. MR 1813436
  * [20] D. Welsh, The Tutte polynomial, Random Structures Algorithms 15 (1999), no. 3-4, 210–228. MR 1716762
  * [21] T. Zaslavsky, The Möbius Function and the Characteristic Polynomial, Combinatorial geometries, Encyclopedia Math. Appl. 29 (1987) 114-–138 MR 921071


# Multivariate generalized linear mixed models for underdispersed count data

Guilherme Parreira da Silvaa, Henrique Aparecido Laureanoa, Ricardo Rasmussen
Petterleb, Paulo Justiniano Ribeiro Júniora, Wagner Hugo Bonata CONTACT
Guilherme Parreira da Silva. Email<EMAIL_ADDRESS>aLaboratory of Statistics and Geoinformation, Department of Statistics, Paraná
Federal University, Curitiba, Brazil; bDepartment of Integrative Medicine,
Paraná Federal University, Curitiba, Brazil

###### Abstract

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 BFGS. We delimited this problem by
studying only count response variables with the following distributions:
Poisson, negative binomial (NB) and COM-Poisson. While the first distribution
can model only equidispersed data, the second models equi and overdispersed,
and the third models all types of dispersion. 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, 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: AIC, BIC and
likelihood. As a consequence, 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.

###### keywords:

Regression models. Automatic differentiation. Multivariate models. Template
Model Builder. Optimization. Laplace Approximation.

††articletype: ARTICLE

##  1 Introduction

Researchers are often interested to understand the relationship between a set
of response variables and a set of covariates in their different areas of
study. For example, doctors are interested to know whether polypharmacy is
related to complications after surgery; veterinarians may be interested to
know whether animal welfare is related to meat quality; administrators may be
interested to know whether the usage of the new policy has improved social
indicators. When we have a set of covariates and one specific objective
(response variable), these situations can be addressed in the statistical
literature by regression models. Certainly, one of the most known and widely
applied is the linear regression (LM) model [1].

LM is widely used due to its simplicity and the general ordinary least squares
estimation procedure, which is covered in different textbooks in different
areas such as business, numerical optimization, agronomy, among others. To
correct apply it, it is necessary to verify whether the residuals are
independent, not autocorrelated and with homogeneous variance. These
assumptions can be too restrictive depending on the context.

As an alternative, the generalized linear model (GLM) is a more flexible
approach [2]. GLM is a class of models that generalizes LM by supporting
response variables that belong to the exponential family. The distributions
that belong to the exponential family are binomial, gamma, inverse Gaussian,
normal, and Poisson. Moreover, it is built upon a link function that connects
the linear predictor to the expectation of the response variable, and its
variance can be related to its mean.

A count data represents the number of times that an event occurs in a fixed
interval, such as time, space, distance, area, among others. Therefore, it is
finite and non-negative. One example of such data is the number of times ear,
body posture and head orientation changes in ewes after brushing (a treatment
proposed to increase welfare in animals) [3]. The Poisson distribution is
widely used for this purpose but relies on the fact that the variance of the
data is equal to its mean, which is known as equidispersion. However, this
assumption is too restrictive, and a different mean-variance relationship can
be found, such as overdispersion and underdispersion. Overdispersion occurs
when the variance of the data is greater than the mean, and it is often found
in practice. It usually happens due to excess of zeros, heavy-tailed
distribution, or absence of a covariate to model the data [4]. On the other
hand, underdispersion occurs when the variance is smaller than the mean and it
is less usual than the overdispersion case [5]. An underdispersed random
variable is characterized by a smaller range of observed counts compared to an
overdispersed count data.

Different distributions have been proposed to model count data. When the
variance is equal to the mean, after considering the effects of all
covariates, Poisson is the most obvious choice. When the data are
overdispersed, negative binomial (NB) type II under the same framework of GLM
is a good choice. The Extended Poisson Tweedie [6] based on the Poisson
Tweedie distribution [7, 8], Conway-Maxwell-Poisson (COM-Poisson) [9] and
Gamma Count [5] distributions can be used to model either under-, equi- or
overdispersed data. The drawback of them is that the probability mass function
(pmf) does not have a closed-form expression, making the process of inference
time consuming for procedures that rely on the pmf, such as likelihood.

Beyond different distributions, hurdle and the zero-inflated models [10] can
be used to model count data, especially when the data is zero-inflated [11].
Two widely-known alternatives are the hurdle and the zero-inflated models. The
disadvantage of them is that interpretability becomes cumbersome especially
for hurdle model.

All the strategies pointed out here consider that we have available only one
response variable. However, it is not difficult to find in the literature
datasets where the researchers possess more than 1 response variable for the
same study. Usually, the analysis is made by each response individually due to
the lack of alternatives in statistical software. Nevertheless, there is an
increasing interest in the literature to develop models or distributions that
can handle multivariate responses, that is, when there is more than one
response variable [12].

One approach to model more than one response variable simultaneously is to
construct multivariate distributions for count data. [13] present three
alternatives to model multivariate count data. The first assumes that the
marginal distribution is Poisson, and a multivariate distribution is build
under copulas or multivariate distribution theory [14]. The second uses a
mixture of independent Poisson. The third method generalizes the first one,
where the conditional distributions are also Poisson. However, none of them
deals with under-dispersed data. [15] proposed a multivariate generalized
Poisson regression model based on the multivariate generalized Poisson
distribution (MGPD) that can deal with equi-, under- or overdispersed data,
the correlation estimates can either be positive or negative, and the
estimation is made under the maximum likelihood (ML) framework.

[16] provides an overview of different distributions/models for count data. It
presents the multivariate NB model (MNBM) and the multivariate Poisson-gamma
mixture model (MPGM), which allow only for overdispersion and nonnegative
correlation. It also presents the multivariate Poisson-lognormal regression
(MPLR) model and the latent Poisson-normal regression model, which allows
positive and negative correlation. However, it is suitable only for
overdispersed data.

[17] proposed the Multivariate Covariance Generalized Linear Models (MCGLM),
which is a class of models based on quasi-likelihood [18] and generalized
estimating equations (GEE) [19], that allows fitting multivariate models using
only second-order moments assumptions with correlated data. Bayesian
Regression Models using Stan - brms package [20] and MCMC Generalized Linear
Mixed Models - MCMCglmm package [21] provide a framework to model multivariate
models via Bayesian inference [22].

Another alternative is to model the correlation between response variables for
the same individual using the class of hierarchical GLM [23]. This class
allows to model 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 correlation between
sample units, not for response variables, such methodology is implemented in
consolidated packages in software R [24], such as, glmmTMB [25], lme4 [26] and
nlme [27].

In this article, we propose to model multivariate underdispersed count data
under the framework of GLMMs to accommodate correlation between response
variables.

This article contains six sections including this introduction. Section 2
describes the dataset used as a model example of application. Section 3
presents a literature review of the GLMM model. Section 4 proposes the MGLMM
model along with the estimation procedure. Section 5 presents the results of
the model applied to the data from Section 2. Finally, Section 6 discusses the
main contributions of this article and future work is also pointed.

##  2 DATASET: NATIONAL HEALTH AND NUTRITION EXAMINATION SURVEY

The National Health and Nutrition Examination Survey (NHANES) is a program
that studies the health and nutrition status of adults and children in the
United States. This survey is being conducted every year since the early
1960s. Nowadays, it examines a nationally representative sample of about 5,000
subjects each year [28]. Among different types of data collected, the main
objective of this analysis was to investigate whether demographical variables
influence sexual behaviour.

We used the sex behaviour and demographic datasets. From the first, it was
selected three response variables Nmsp (Number of male sex partners in the
past year), Nmosp (Number of male oral sex partners in the past year) and
Nspfy (Number of sex partners who are five years older in the past year). From
the second dataset, it was used the covariates race (0 = Others, 1 = White),
education level (range from 1 = Less Than $9^{\text{th}}$ Grade to 5 = College
Graduate or above) and marital status (0 = Marital, 1 = Others). After
deleting those respondents who had missing data, the sample consisted of 1281
women, with ages ranging from 18 to 80, 57% were married, 43% white and 31%
had some college or associates (AA) degree.

One way to analyse this dataset is via regression models. Once all three
response variables are count data, we need to choose a suitable probability
mass function for this data. In this paper, we compared the Poisson, NB and
COM-Poisson probability distributions. As this is a cross-sectional study, we
do not have responses correlated over time. However, we may have correlated
responses variables, once they were measured in the same individual. The model
proposed in section 4 can accommodate this situation.

Table 1 presents the mean, variance, Fisher Dispersion Index (DI) [29] for
every response variable and the generalized dispersion index (GDI) [30] for
the dataset. The main reason for choosing those response variables is the fact
that Nmsp and Nmosp can be considered marginally underdispersed, once the
sample variance is smaller than the sample mean. Moreover, Nspfy is only a
little overdispersed. It is easily seen from the DI, which is calculated by
dividing the variance of the variable by its mean. A variable with
$\text{DI}>1$ can be said as overdispersed, DI = 1 as equidispersed and
$\text{DI}<1$ underdispersed.

Moreover, the GDI is a recently proposed multivariate dispersion index. It is
calculated based on the expectation of each variable and the covariance
between them. When the number of variables is equal to 1, it is just the
classical Fisher DI. A $\text{GDI}>1$ classifies the multivariate responses as
overdispersed, GDI = 1 as equidispersed and $\text{GDI}<1$ underdispersed. The
standard error (SE) for GDI is calculated using the asymptotic behaviour of
the estimator. According to this index, a 95% confidence interval based on SE
suggests that this dataset is equidispersed as 1 is included on it.

Table 1: DESCRIPTIVE MEASUREMENTS FOR NHANES RESPONSE VARIABLES | Spearman Correlation $\rho$ | Mean | Variance | DI | GDI(SE)
---|---|---|---|---|---
| Nmsp | Nmosp | Nspfy |  |  |  |
Nmsp |  | 0.033 | 0.222 | 1.313 | 1.084 | 0.826 |
Nmosp |  |  | 0.038 | 0.372 | 0.350 | 0.939 | 1.092(0.22)
Nspfy |  |  |  | 0.368 | 0.411 | 1.117 |

The three response variables show no or small (Nmsp and Nspfy) correlation
between them. Figure 1 shows the barplot of each response variable. We can see
that there is a higher frequency for non-occurrence of events, rather than the
occurrence.

Figure 1: BARPLOT FOR EACH RESPONSE VARIABLE FROM NHANES DATA

##  3 GENERALIZED LINEAR MIXED MODELS (GLMM)

The GLMM class can be used to address dependent observations, overdispersion,
among others [31]. In addition to the fixed effect from the GLM, it
accommodates extra variability through an unobservable variable, the random
effect, and is calculated for each set of dependent observations. Once the
random effect is an unobserved random variable, we can write the GLMM in a
hierarchical structure. The following notation is valid when the dependent
observations come for repeated measurements and longitudinal studies:

$\displaystyle\mathrm{Y}_{il}\mid\mathrm{\mathbf{b}}_{i\cdot}\sim
f(\mu_{il},\phi)$
$\displaystyle\mathrm{g}(\mu_{il})=x_{ij}^{\top}\bm{\beta}+z_{ij}^{\top}\mathrm{\mathbf{b}}_{i\cdot}$
(1)
$\displaystyle\mathrm{\mathbf{b}}_{i\cdot}\sim\mathcal{N}\left(\mathbf{0},\bm{\Sigma}\right),$

where $\mathrm{Y}_{il}$ is the i-th unit sample measured in the l-th times (in
a longitudinal study) or in the l-th group (in a repeated measurements study),
$\mathrm{\mathbf{b}}_{i\cdot}$ is a vector of random effects for each sample
unit that follows a multivariate normal distribution with mean $\mathbf{0}$
and covariance matrix $\bm{\Sigma}$ to accommodate the variance of each random
effect and covariance among them, $f$ is the conditional distribution of
$\mathrm{Y}_{il}$ on $\mathrm{\mathbf{b}}_{i\cdot}$, $\mu_{il}$ is the mean
for every sample unit $i$ and measurement $l$, $\mathrm{g}(.)$ is a suitable
link function, and $\mathrm{z}_{ij}^{\top}$ is the design matrix associated to
the random effect.

##  4 MULTIVARIATE GENERALIZED LINEAR MIXED MODEL (MGLMM) FOR COUNT DATA

Let $\mathrm{Y}_{ir}$ be the multivariate outcome for subject _i_ , where
$i=1,\dots,n$ and response variable _r_ , where $r=1,\dots,k$. Suppose that a
set of _p_ known covariates is available for each response _r_ , therefore,
$\mathrm{x}_{irj}$ is the value of the _j-th_ covariate for individual _i_ and
response _r_. The purpose of this article is to provide a joint model for a
set of response variables. We start from the standard GLMM model specification
with only a random intercept. The first component we need to assume is the
conditional distribution of the response variables:

$\mathrm{Y}_{ir}\mid\mathrm{b}_{ir}\sim f(\mu_{ir};\phi_{r}),$

where all response variables share the same distribution. However, with
different mean and dispersion parameters. The second component we need to
specify is the linear predictor:

$\mathrm{g}_{r}(\mu_{ir})=x_{irj}^{\top}\bm{\beta}_{r}+\mathrm{b}_{ir},$

where $\mathrm{g}_{r}(.)$ is a suitable link function, $\bm{\beta}_{r}$ is a
$p\times 1$ vector of parameter estimates and $\mathrm{b}_{ir}$ is the random
intercept value for each individual and response variable. Lastly, the
distribution of the random effects is specified by:

$\begin{pmatrix}\mathrm{b}_{i1}\\\ \mathrm{b}_{i2}\\\ \vdots\\\
\mathrm{b}_{ir}\end{pmatrix}\sim\mathrm{NM}\begin{pmatrix}\begin{bmatrix}0\\\
0\\\ \vdots\\\ 0\end{bmatrix};\underset{r\times
r}{\bm{\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},$

where each random effect has mean 0, variance $\sigma^{2}$, and
$\rho_{r{r}^{\prime}}(r\neq{r}^{\prime})$ measures the correlation between
each pair of random effects. Despite this is a general modelling framework,
this article addresses only the same distribution (either Poisson, binomial
negative or COM-Poisson) and link function (logarithm) for all responses
variables.

An important question for this type of model is whether it is possible to
estimate simultaneously the dispersion parameters of the pmf and the variance
parameters of the random effects, once they measure variability related to the
same random variable, and can cause identifiability problems.

###  4.1 INFERENCE AND ESTIMATION

In this section, we present the estimation procedure used to obtain the
parameter estimates based on the likelihood function. In order to estimate the
model presented in section 4 under the maximum likelihood (ML) paradigm, it is
necessary to construct the joint distribution for both random variables
($\mathrm{Y}_{ir}\text{ and }\mathrm{b}_{ir}$). Once we have a hierarchical
structure, the joint distribution of $\mathrm{Y}_{ir}\text{ and
}\mathrm{b}_{ir}$ can be factored into
$f(\mathrm{Y}_{ir},\mathrm{b}_{ir})=f(\mathrm{Y}_{ir}\mid\mathrm{b}_{ir})f(\mathrm{b}_{ir})$.
As only $\mathrm{Y}$ is observed, we are interested on the marginal
distribution of $\mathrm{Y}$ that can be obtained integrating out
$\mathrm{b}_{ir}$, that is, $f(\mathrm{Y}_{ir})=\int
f(\mathrm{Y}_{ir}\mid\mathrm{b}_{ir})f(\mathrm{b}_{ir})\mathrm{db}_{ir}$,
which in practical terms means that the joint distribution is averaged over
the $\mathrm{b}_{ir}$ terms.

Now, we particularize the ML estimation method to the model described in
section 4. The objective is to estimate the parameter vector
$\bm{\theta}=\\{\bm{\beta},\bm{\phi},\bm{\sigma^{2}},\bm{\rho}\\}$, where
$\bm{\beta}$ is the regression parameter vector, $\bm{\phi}$ the extra
parameters of each distribution (dispersion in most of cases), and
$\bm{\sigma^{2}}\text{ and }\bm{\rho}$ are the parameters that compose the
variance covariance matrix $\bm{\Sigma}$. The marginal likelihood function for
the model described in section 4 for each sample unit is

$\mathrm{L}_{i}(\bm{\beta},\bm{\Sigma},\bm{\phi}\mid\mathrm{\mathbf{y}})=\int\prod_{r=1}^{k}f(\mathrm{y}_{r}\mid\mathrm{\mathbf{b}},\bm{\beta},\bm{\phi})f(\mathrm{\mathbf{b}}\mid\bm{\Sigma})\mathrm{d}\mathrm{\mathbf{b}},$
(2)

where $\mathrm{\mathbf{y}}$ is the $k$-response vector, $\mathrm{y_{r}}$ is
each response and $\mathrm{\mathbf{b}}$ is the random effects vector. The full
likelihood for $\bm{\theta}$ is given by

$\mathrm{L}(\bm{\beta},\bm{\Sigma},\bm{\phi})=\prod_{i=1}^{N}\mathrm{L}_{i}(\bm{\beta},\bm{\Sigma},\bm{\phi}\mid\mathrm{\mathbf{y}}),$

where _N_ is the total number of sample units. Under independence between
sample units, we have:

$\mathrm{\mathbf{L}}(\bm{\beta},\bm{\Sigma},\bm{\phi})=\prod_{i=1}^{N}\int\prod_{r=1}^{k}f(\mathrm{y}_{r}\mid\mathrm{\mathbf{b}},\bm{\beta},\bm{\phi})f(\mathrm{\mathbf{b}}\mid\bm{\Sigma})\mathrm{d}\mathrm{\mathbf{b}}.$
(3)

The Equation 3 is composed of the product of two probability distributions.
The first one is the probability distribution of the sample units, while the
second is the probability distribution of the random effects. The distribution
of the random effects is assumed to be a multivariate normal distribution. As
the probability distribution of the response is not normal, the integral does
not present a closed form solution expression, thus numerical methods are
required to solve such integral.

####  4.1.1 Numerical integration via Laplace approximation

Once we are dealing with non-normal data, one possibility is to use numerical
integration methods to solve the integral in Equation 2 for each set of sample
units. Other possibilities include Penalized Quasi Likelihood (PQL) and Monte
Carlo techniques that will not be covered here. It is important to note that
the integral is of dimension $k$, and $k$ can vary from two to $r$ (the number
of response variables). Methods based on numerical quadratures such as
trapezium, $1/3$ Simpson, Gauss-Hermite or adaptative Gauss-Hermite use too
many points to evaluate the integrand and its complexity is proportional to
the dimension of the integral [32]. To solve the integral in Equation 2 we
used the LA, which is a special case of the adaptive Gauss-Hermite quadrature
(AGHQ) when it is used only one integration point [33]. Although LA is faster
than AGHQ, it is less accurate [34].

The idea of LA is to replace an integral with a tractable closed-form
expression. This resulting expression is maximized concerning the variable
that we wanted to integrate and the integral is solved. Applying the LA to the
integral in Equation 2 results in the marginal likelihood function. The
approximation is obtained by

$\int_{\mathbb{R}^{k}}\text{exp}\\{\mathrm{\mathbf{Q(b)}}\\}\text{d}\mathrm{\mathbf{b}}\approx(2\pi)^{n_{\mathrm{\mathbf{b}}}/2}\left|-\mathbf{Q^{\prime\prime}(\hat{b})}\right|^{-1/2}\text{exp}\\{\mathbf{Q(\hat{b})}\\},$
(4)

where $\mathrm{\mathbf{Q(b)}}$ is a uni-modal and bounded function of the
variable $\mathrm{\mathbf{b}}$ and $n_{\mathrm{\mathbf{b}}}$ is the dimension
of the integral, that it is $k$ in this case. Careful accounting of the
approximation error shows it to generally be $O(n^{-1})$ where $n$ is the
sample size (assuming a fixed length for $\mathrm{\mathbf{b}}$) [35]. In this
case, $\mathrm{\mathbf{Q(b)}}$ is the log of the integrand in Equation 2,
$\mathrm{\mathbf{b}}$ the random effect, $\mathbf{\hat{b}}$ the maximized
random effect values, while all other $\bm{\theta}$ parameters remain
constant. Therefore, $\mathbf{Q(\hat{b})}$ is the maximum value of the
function with respect to $\mathrm{\mathbf{b}}$ and
$\mathbf{Q^{\prime\prime}(\hat{b})}$ is the curvature of the function in the
maximum, or in other words, it corresponds to the value of the second
derivative (Hessian) of the function in the maximum.

It is important to say that the integral of order $k$ in Equation 2 has to be
solved repeatedly in every step of the ML method for each one of the $N$
sample units. Therefore, we have two optimizations procedures: one that is
called the outer maximization, which is the maximization of the marginal
likelihood (result from the integral), and another one that is called the
inner maximization, which is the maximization required inside the LA.

For each optimization procedure, we have different situations. In the inner
process is feasible to calculate the Hessian matrix efficiently (once we know
the integrand) and to choose a good initial guess for $\mathrm{\mathbf{b}}$.
Therefore, the inner process maximization can be efficiently implemented by
the Newton-Raphson (NR) method (that requires second-order derivatives). On
the other hand, in the outer maximization, it is not feasible to calculate the
Hessian matrix efficiently, neither we have a good initial guess after the
integral is solved. Therefore, it will be maximized with algorithms that
require only first-order derivatives and do not depend strongly on the initial
guess as NR does. The derivatives were obtained through automatic
differentiation [36].

##  5 SOFTWARE IMPLEMENTATION

In this section, we present the software implementation employed to fit the
model presented in section 4. The software used was R version 4.0.2 [24] along
the Template Model Builder - TMB [37] package. TMB is an R package that offers
a collection of tools to build complex random effects statistical models
through C++ templates. TMB is based on state-of-the art software: CppAD [38],
Eigen C++ [39], BLAS [40], among others libraries written in C++, which are
responsible for obtaining the derivatives through AD, linear algebra
computations and parallelization, respectively.

Once the user supplies an objective function in a C++ template file (which,
usually is the negative log-likelihood function), it obtains the marginal
likelihood via LA integrating out the $\mathrm{\mathbf{b}}$ random effects.
The internal optimization of the LA is made via Newton’s method, with first
and second derivatives provided via AD from TMB. After that, the marginal can
be optimized to obtain the MLE parameters with any general-purpose algorithm
(using first derivative information provided from TMB), such as PORT or BFGS
implemented in nlminb and optim routine in R. Moreover, the standard deviation
of the parameters (or a function of it) can be obtained via the Delta method
[41]; profiling is available too.

A implementation of the Poisson MGLMM is available on
http://www.leg.ufpr.br/doku.php/publications:papercompanions.

###  5.1 Reparametrization

We had to use reparametrization techniques for those parameters that vary in a
subset of $\mathbb{R}$. Dispersion parameters for both NB and COM-Poisson
distributions were log reparametrized (ln to be more precise) into the
dispersion parameter as input, allowing it to will vary on all $\mathbb{R}$.
For $\bm{\sigma}$ and $\rho$ we used an efficient reparametrization from TMB
under the normal multivariate specification.

Instead of estimating $\bm{\rho}$ directly in $\bm{\Sigma}$, the TMB
reparametrization consists to estimate a different unrestricted parameter
$\bm{\varrho}$ in a lower triangular matrix with unit diagonal $L$. An
unstructured symmetric positive definite correlation matrix $\bm{\Omega}$ can
be obtained via
$\bm{\Omega}=\mathbf{D}^{-\frac{1}{2}}\mathbf{L}\mathbf{L}^{\top}\mathbf{D}^{-\frac{1}{2}}$,
where $\mathbf{D}=\text{diag}(\mathbf{L}\mathbf{L}^{\top})$.

The $\bm{\Omega}$ matrix is only one part of TMB’s $\bm{\Sigma}$
decomposition. Consider the case we have 4 response variables, so, we want to
estimate 6 correlation parameters. Thus, the lower triangular matrix
$\mathbf{L}$ has order 4 and is filled row-wise, such that:

$\mathbf{L}=\begin{pmatrix}1&&&\\\ \varrho_{0}&1&&\\\
\varrho_{1}&\varrho_{2}&1&\\\
\varrho_{3}&\varrho_{4}&\varrho_{5}&1\end{pmatrix}.$

According to the the restrictions imposed into $\bm{\Omega}$, and the way the
normal multivariate was coded into TMB, the variance parameters of the random
effects need to be supplied as standard deviations. To obtain the full
variance-covariance matrix $\bm{\Sigma}$ we can use
$\bm{\Sigma}=\mathbf{W}\bm{\Omega}\mathbf{W}$ where $\mathbf{W}$ is a diagonal
matrix with entries being equal to the random effects’s standard deviation.

##  6 DATA ANALYSIS

In this section we apply the proposed model in section 4 to analyze the NHANES
data set presented in section 2. Initial values had to be chosen carefully for
the applied models. Firstly, it was fitted a MCGLM model [17] with log link
function, power parameter fixed, and variance Tweedie, which corresponds to
fit a Poisson Model via quasi-likelihood, in order to obtain initial parameter
estimates for the regression and variance parameters (based on the variance of
the residuals). On the other hand, the correlation parameter was set to 0. We
did not use MCGLM to obtain initial estimates for the correlation parameter
due to the difference in methodologies. This parameter estimates values were
considered as initial values to the Poisson model.

Then, for every distribution, the estimation was made in two steps. In the
first, the model was fitted based on a random sample (SRS) of size 350. In the
second step, the initial values from the sample model were used to fit the
model with the whole dataset. After that, the first estimation was carried out
with the PORT routine, followed by BFGS and PORT. For the dataset analyzed,
PORT was the fastest and produced a more consistent optimization over BFGS in
most scenarios. In every step, initial values were used from the last
optimization, and the choice of the intermediate models and the final results
reported were based on the maximum value of the logLik.

The initial values for the NB model were the final estimates from the Poisson
model (with initial value to the dispersion parameter $\phi=1$, i.e., small
overdispersion). In the same way, the initial values for the COM-Poisson model
were the final estimates from the NB model (with initial value to the
dispersion parameter $\nu=1$, i.e., equidispersion). We used the mean-
parametrization of the COM-Poisson model, where $\nu>1$ is considered
underdispersion, and $\nu<1$, overdispersion [42].

This workflow is described in Figure 2. The reported results are only from the
full data, as the samples were used to obtain initial estimates.

Figure 2: Model estimation process for each dataset MCGLM Poisson sample
Poisson full data NB sample 1ºPORT 2ºBFGS 3ºPORT NB full data COM-Poisson
sample 1ºPORT 2ºBFGS 3ºPORT COM-Poisson full data 1ºPORT 2ºBFGS 3ºPORT

In addition to the models presented in section 4, simpler versions of them
were also fitted. We used some different scenarios concerning the variance of
the random effects, one for the dispersion parameter and another for the
correlation parameter. The first case considered a model with a common
variance specification, where the variance of the random effect was equal to
all response variables. In the second case, we used a fixed variance
specification, where the variance was fixed at $1$. The third scenario was
when the dispersion parameter was fixed to $\phi=1$ for the NB and $\nu=1.5$
for the COM-Poisson model indicating small underdispersion for this
distribution. The final and fourth scenario sets the correlation parameter to
zero, which aims to reproduce the case where every response distribution would
be fitted separately. As the Poisson distribution does not have a dispersion
parameter, only the full specification of the model was used and when the
correlation parameter was set to zero. This workflow and the linear predictor
selection that is explained in the next subsection is presented in Figure 3.

Figure 3: Model structure and linear predictor selection Best model Poisson
Best Poisson1 Poisson5 NB Best NB1 NB2 NB3 NB4 NB5 COM-Poisson Best COM-
Poisson1 COM-Poisson2 COM-Poisson3 COM-Poisson4 COM-Poisson5 Choosing model
structureLinear Predictor Selection1 - Full model 2 - Fixed variance 3 - Comum
variance 4 - Fixed dispersion 5 - Correlation zero

We reported the results of the models presented in Figure 2 and their
variations with log-likelihood value (logLik \- larger the best), Akaike
Information Criterion (AIC - lower the best), Bayesian Information Criterion
(BIC - lower the best), number of parameters estimates (np) and whether the SE
was returned or not for the parameter estimates.

###  6.1 NHANES data

Table 2 presents the goodness-of-fit measures for NHANES data from different
distributions and specifications.

Table 2: Goodenss-of-fit measures for NHANES data from different distributions and specifications Model | np | AIC | BIC | Loglik | SE
---|---|---|---|---|---
Poisson | 21 | 7145.1 | 7253.4 | -3551.6 | ✓
Rho Zero Poisson | 18 | 7169.8 | 7262.6 | -3566.9 | ✗
NB | 24 | 7150.2 | 7273.9 | -3551.1 | ✓
Fixed Dispersion NB | 21 | 8092.4 | 8200.6 | -4025.2 | ✗
Comum Variance NB | 22 | 7203.4 | 7316.8 | -3579.7 | ✗
Fixed Variance NB | 21 | 7924.2 | 8032.5 | -3941.1 | ✓
Rho Zero NB | 21 | 7176.8 | 7285.1 | -3567.4 | ✗
COM-Poisson | 24 | 4615.9 | 4739.6 | -2284.0 | ✓
Fixed Dispersion COM-Poisson | 21 | 6991.7 | 7099.9 | -3474.8 | ✗
Comum Variance COM-Poisson | 22 | 5608.2 | 5721.6 | -2782.1 | ✓
Fixed Variance COM-Poisson | 21 | 6649.8 | 6758.0 | -3303.9 | ✓
Rho Zero COM-Poisson | 21 | 5526.3 | 5634.5 | -2742.1 | ✓

Overall, the model which best fitted the data with respect to the three fit
measures was the COM-Poisson with full specification. We can also note a large
difference in absolute numbers against the competitors. The second best
specification was with rho zero for the COM-Poisson model. Moreover, we can
note a very close logLik between the NB and the Poisson: it happened due to a
large value for the dispersion parameter of the NB, approaching for the
Poisson model. Among the smaller specifications, the fixed dispersion was the
worst scenario for both NB and COM-Poisson models.

These results agree with the underdispersion characteristic of the data
presented in Table 1. It is important to note that beforehand only COM-Poisson
is able to deal with underdispersion data due to the $\nu$ parameter. For
Poisson and NB models, the random effects structure only added some extra
variance for each response variable. For Poisson, the extra variance parameter
gives the flexibility to model only overdispersion; while for NB, it allows
modelling even greater variability than the traditional NB model. Therefore,
the analysis of this dataset confirmed the expected results by the specified
model and distributions.

The checkmark ✓ in the last column specifies that the SE of the parameter
estimate was returned at least for one parameter of the model. The ✗ indicates
that it was not possible to return the SE for any parameter estimate of the
model. In order to deal with this problem, we modified the linear predictor
excluding those covariates that their respective parameter estimates did not
have the SEs returned for each response variable based on the best model, and
it was refitted. The same linear predictor was used to fit the best model from
each distribution, in this case, Poisson, and NB. The results are presented in
Table 3.

Table 3: Goodness-of-fit measures for NHANES data from the best specification for each distribution Model | np | AIC | BIC | logLik | SE
---|---|---|---|---|---
Poisson | 16 | 7147.8 | 7230.3 | -3557.9 | ✓
NB | 19 | 7153.1 | 7251.0 | -3557.5 | ✓
COM-Poisson | 19 | 4448.6 | 4546.6 | -2205.3 | ✓

We can see that reducing the linear predictor improved the COM-Poisson model
fit, while it worsened the Poisson and NB in terms of logLik and AIC, but
improved the BIC because of a greater penalty due to the number of parameters.
In order to better understand the behaviour of the parameter estimates, Figure
4 compares the regression estimates, Table 4 compares the dispersion and
Equation 5 presents a matrix where out of the diagonal are presented
correlation estimates and the SE of each random effect, along with the
standard deviation and SE in the diagonal.

Figure 4: Regression parameter estimates and 95% confidence intervals by
outcome and final model

From Figure 4 we can see that for most of the estimates, the COM-Poisson
confidence intervals were smaller than NB and Poisson models. As these
distributions can not handle underdispersed data, they tend to overestimate
the SE of the parameter estimates. The confidence interval for Nmosp intercept
based on COM-Poisson was the only larger than NB and Poisson models. Moreover,
we see that the intercept point estimate was not close between COM-Poisson and
the other two models in every scenario.

Regarding the relationship between the covariates and the response variables,
the models considered, mostly agree in the direction of the relationship to
each response variable. Race (0 = Others, 1 = White) covariate had a null
effect for Nmsp and Nmosp; for Nspfy had a negative effect for Poisson and NB
models and null for COM-Poisson model. Marital status (0 = Married, 1 =
Others) had a positive effect for all three response variables and Poisson and
NB models; while COM-Poisson model had only a positive effect for Nmsp.
Finally, education level (range from 1 = Less Than $9^{th}$ Grade to 5 =
College Graduate or above) had a negative effect on Nspfy for NB and Poisson
models and null effect for COM-Poisson. In those models that can not handle
underdispersion, more information was captured by the regression parameters
than in the COM-Poisson model.

Table 4: Dispersion parameter estimates and SEs for each model and outcome of NHANES data | NB ($\phi$) | COM-Poisson ($\hat{\nu}$)
---|---|---
Outcome | Estimate | SE | Estimate | SE
Nmsp | 4958.17 | 21947.1 | 46.513 | 0.059
Nmosp | 753.63 | 3878.4 | 20.964 | 5.829
Nspfy | 1996.65 | 12993.9 | 19.675 | 1.286

From Table 4 we can see that $\phi$ was large enough to approximate the NB to
a Poisson model followed by an even larger SE; it justifies the similar model
measures presented in Table 2 and Table 3. For COM-Poisson, large $\hat{\nu}$
values indicate underdispersion and a small SE indicates that $\hat{\nu}$ is
not one, so, it cannot be considered equidispersed.

$\displaystyle\bm{{\Sigma}^{\prime}_{\text{Poisson}}}$
$\displaystyle=\begin{bmatrix}0.18(0.03)^{*}&0.94(0.1)^{*}&0.97(0.04)^{*}\\\
&0.23(0.08)^{*}&0.92(0.13)^{*}\\\ &&0.63(0.07)^{*}\end{bmatrix}$ (5)
$\displaystyle\bm{{\Sigma}^{\prime}_{\text{NB}}}$
$\displaystyle=\begin{bmatrix}0.18(0.03)^{*}&0.97(0.08)^{*}&0.99(0.02)^{*}\\\
&0.23(0.08)^{*}&0.97(0.09)^{*}\\\ &&0.63(0.07)^{*}\end{bmatrix}$
$\displaystyle\bm{{\Sigma}^{\prime}_{\text{COM-Poisson}}}$
$\displaystyle=\begin{bmatrix}0.48(<0.01)^{*}&<-0.01(<0.01)^{*}&<0.01(<0.01)\\\
&1.21(0.15)^{*}&<-0.01(0.01)\\\ &&1.25(0.06)^{*}\end{bmatrix}$

Equation 5 presents the correlation between random effects in the upper
diagonal and the standard deviation in the diagonal. Stars represent
statistical significance at 5% level. It was necessary to make a distinction
between $\bm{\Sigma}$ and $\bm{{\Sigma}^{\prime}}$ because the last one uses
standard deviation in the diagonal, and the first use variance.

Equation 5 shows that the correlation estimates for Poisson and NB model are
close to one and are significant at 5% level. In contrast, the correlation
estimates for the COM-Poisson model are close to zero and only one is
significant (between Nmsp and Nmosp). It seems that the underdispersion
information not modelled by Poisson and NB models are somewhat present in the
correlation estimates. Moreover, the standard deviation of the random effect
was larger for COM-Poisson than was the Poisson and NB models. It may occur
because at the same time that the COM-Poisson has greater variability due to
the random effect standard deviation, it balances out with smaller dispersion
of $\hat{\nu}$ parameter.

##  7 Discussion

The main focus of this article was to propose MGLMM for underdispersed count
data. With this new class of statistical models, we can model more than one
response variable in the same framework accommodating a random effect that
follows a multivariate normal distribution (MN). The parameters of the
covariance matrix from the MN allow us to account for the correlation between
the random effects and the variance of them. This extra feature is not
available in a GLMM model, where it is necessary to model one response
variable at a time.

In particular, we addressed only count data problems, considering Poisson, NB
and COM-Poisson distributions. This model was implemented using the TMB
package in R. It was estimated under the ML paradigm using numerical
integration via LA with inner and outer optimization based on Newton’s method
and general-purpose algorithm respectively, such as BFGS and PORT routines,
which derivatives are provided through AD.

A dataset from the NHANES survey was analysed by each model and variations of
them: rho set to zero, fixed dispersion, fixed variance and common variance
for all random effects. It comprised of three response variables, being one
equidispersed and two underdispersed. The COM-Poisson was the best model
compared to their counterparts according to logLik, AIC and BIC, especially
because of its ability to model underdispersion. The SEs were smaller compared
to the Poisson and NB models. The dispersion parameter $\nu$ indicated
underdispersion, which is expected due to the nature of the data and the
correlation parameter $\rho$ was almost zero.

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

Future work could include the exploration of the residuals for this model and
propose a correlation coefficient between the response variables, not the
random effects, that it is based on the correlation of the random effect.

## Disclosure statement

No potential competing interest was reported by the authors.

## Funding

The corresponding author was funded with a scholarship provided by the
Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES).

## References

  * [1] Galton F. Regression towards mediocrity in hereditary stature. The Journal of the Anthropological Institute of Great Britain and Ireland. 1886;15:246–263.
  * [2] Nelder JA, Wedderburn RW. Generalized linear models. Journal of the Royal Statistical Society: Series A (General). 1972;135(3):370–384.
  * [3] Tamioso PR, Maiolino Molento CF, Boivin X, et al. Inducing positive emotions: Behavioural and cardiac responses to human and brushing in ewes selected for high vs low social reactivity. Applied Animal Behaviour Science. 2018;208:56 – 65. Available from: http://www.sciencedirect.com/science/article/pii/S0168159118304490.
  * [4] Grunwald GK, Bruce SL, Jiang L, et al. A statistical model for under-or overdispersed clustered and longitudinal count data. Biometrical Journal. 2011;53(4):578–594.
  * [5] Zeviani WM, Ribeiro Jr PJ, Bonat WH, et al. The gamma-count distribution in the analysis of experimental underdispersed data. Journal of Applied Statistics. 2014;41(12):2616–2626.
  * [6] Bonat WH, Jørgensen B, Kokonendji CC, et al. Extended Poisson–Tweedie: Properties and regression models for count data. Statistical Modelling. 2018;18(1):24–49.
  * [7] El-Shaarawi AH, Zhu R, Joe H. Modelling species abundance using the Poisson–Tweedie family. Environmetrics. 2011;22(2):152–164.
  * [8] Jørgensen B, Kokonendji CC. Discrete dispersion models and their Tweedie asymptotics. AStA Advances in Statistical Analysis. 2016;100(1):43–78.
  * [9] Shmueli G, Minka TP, Kadane JB, et al. A useful distribution for fitting discrete data: revival of the Conway–Maxwell–Poisson distribution. Journal of the Royal Statistical Society: Series C (Applied Statistics). 2005;54(1):127–142.
  * [10] Zeileis A, Kleiber C, Jackman S. Regression models for count data in R. Journal of statistical software. 2008;27(8):1–25.
  * [11] Ridout M, Demétrio CG, Hinde J. Models for count data with many zeros. In: Proceedings of the XIXth international biometric conference; Vol. 19; International Biometric Society Invited Papers Cape Town, South Africa; 1998. p. 179–192.
  * [12] Bonat WH, Jørgensen B. Multivariate covariance generalized linear models. Journal of the Royal Statistical Society Series C (Applied Statistics). 2016;:649–675.
  * [13] Inouye DI, Yang E, Allen GI, et al. A review of multivariate distributions for count data derived from the Poisson distribution. Wiley Interdisciplinary Reviews: Computational Statistics. 2017;9(3):e1398.
  * [14] Campbell JT. The Poisson correlation function. Proceedings of the Edinburgh Mathematical Society. 1934;4(1):18–26.
  * [15] Famoye F. A multivariate generalized Poisson regression model. Comm Statist Theory Methods. 2015 Feb;44(3):497–511.
  * [16] Winkelmann R. Econometric analysis of count data. Springer Science & Business Media; 2008.
  * [17] Bonat WH. mcglm: Multivariate covariance generalized linear models; 2016. R package version 0.3.0; Available from: https://CRAN.R-project.org/package=mcglm.
  * [18] Wedderburn RW. Quasi-likelihood functions, generalized linear models, and the gauss—newton method. Biometrika. 1974;61(3):439–447.
  * [19] Liang KY, Zeger SL. Longitudinal data analysis using generalized linear models. Biometrika. 1986;73(1):13–22.
  * [20] BÜrkner PC. Advanced Bayesian multilevel modeling with the R package brms. The R Journal. 2018;10(1):395–411.
  * [21] Hadfield JD. Mcmc methods for multi-response generalized linear mixed models: The MCMCglmm R package. Journal of Statistical Software. 2010;33(2):1–22. Available from: http://www.jstatsoft.org/v33/i02/.
  * [22] Dempster AP. A generalization of bayesian inference. Journal of the Royal Statistical Society: Series B (Methodological). 1968;30(2):205–232.
  * [23] Lee Y, Nelder JA. Hierarchical generalized linear models. Journal of the Royal Statistical Society Series B (Methodological). 1996;58(4):619–678. Available from: http://www.jstor.org/stable/2346105.
  * [24] R Core Team. R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing; 2020. Available from: https://www.R-project.org/.
  * [25] Brooks ME, Kristensen K, van Benthem KJ, et al. glmmTMB balances speed and flexibility among packages for zero-inflated generalized linear mixed modeling. The R Journal. 2017;9(2):378–400. Available from: https://journal.r-project.org/archive/2017/RJ-2017-066/index.html.
  * [26] Bates D, Mächler M, Bolker B, et al. Fitting linear mixed-effects models using lme4. Journal of Statistical Software. 2015;67(1):1–48.
  * [27] Pinheiro J, Bates D, DebRoy S, et al. nlme: Linear and nonlinear mixed effects models; 2017. R package version 3.1-131; Available from: https://CRAN.R-project.org/package=nlme.
  * [28] Nhanes [https://www.cdc.gov/nchs/nhanes/about_nhanes.htm]; 2007. Accessed: 2020-08-14.
  * [29] Fisher RA. The effect of methods of ascertainment upon the estimation of frequencies. Annals of eugenics. 1934;6(1):13–25.
  * [30] Kokonendji CC, Puig P. Fisher dispersion index for multivariate count distributions: A review and a new proposal. Journal of Multivariate Analysis. 2018;165:180–193.
  * [31] Breslow NE, Clayton DG. Approximate inference in generalized linear mixed models. Journal of the American statistical Association. 1993;88(421):9–25.
  * [32] Bonat WH, Ribeiro-Jr PJ. Practical likelihood analysis for spatial generalized linear mixed models. Environmetrics. 2016;27(1):83–89.
  * [33] Tierney L, Kadane JB. Accurate approximations for posterior moments and marginal densities. Journal of the American Statistical Association. 1986;81(393):82–86. Available from: https://www.tandfonline.com/doi/abs/10.1080/01621459.1986.10478240.
  * [34] Signorelli M, Spitali P, Tsonaka R. Poisson–Tweedie mixed-effects model: A flexible approach for the analysis of longitudinal rna-seq data. Statistical Modelling. 2020;:1471082X20936017.
  * [35] Wood SN. Core statistics. Cambridge University Press; 2015. 6.
  * [36] Baydin AG, Pearlmutter BA, Radul AA, et al. Automatic differentiation in machine learning: a survey. The Journal of Machine Learning Research. 2017;18(1):5595–5637.
  * [37] Kristensen K, Nielsen A, Berg CW, et al. TMB: Automatic differentiation and Laplace approximation. Journal of Statistical Software. 2016;70(5):1–21.
  * [38] Bell BM. Cppad: a package for c++ algorithmic differentiation; 2005. Available from: http://www.coin-or.org/CppAD.
  * [39] Guennebaud G, Jacob B, et al. Eigen v3 [Http://eigen.tuxfamily.org]; 2010.
  * [40] Blackford LS, Petitet A, Pozo R, et al. An updated set of basic linear algebra subprograms (blas). ACM Transactions on Mathematical Software. 2002;28(2):135–151.
  * [41] Thygesen UH, Albertsen CM, Berg CW, et al. Validation of ecological state space models using the Laplace approximation. Environmental and Ecological Statistics. 2017;24(2):317–339.
  * [42] Huang A. Mean-parametrized Conway–Maxwell–Poisson regression models for dispersed counts. Statistical Modelling. 2017;17(6):359–380.