Parameters Estimation of a Bivariate Generalized Poisson Distribution with Applications to Metabolic Syndrome Data ()
1. Introduction
For many years there has been an increase in the statistical research to develop multivariate distributions for count data. There are numerous ways by which we can generalize the distributions of univariate counts to multivariate forms. The applications of these multivariate distributions are found in many fields such as engineering, bioinformatics, genetic epidemiology and biomedicine. An interested reader is referred to the recent editions of the book by Winkelmann [1] and the book by Cameron and Trivedi [2] for thorough reviews of regression models for multivariate count data and the very recent work by Tzougas et al. [3]. Needless to say, fitting these multivariate distribution to data is quite complicated and requires a great deal of computing resources. Therefor there is a need to develop more flexible bivariate models and demonstrate their applicability to real life data.
The main purpose of this article is to investigate the statistical characteristics of a class of bivariate generalized Poisson models with a correlation parameter that offers sufficient flexibility for accommodating overdispersion and accounting for the positive correlation between the variables of interest. The distribution is named “Bivariate Generalized Poisson Distribution”, BGPD, which is a special member of a larger family of distributions developed by Shenton and Consul [4]. Some of the interesting properties of this class were studies further by Jain and Sing [5], and by Shoukri [6] who investigated some of the properties of the distribution as a member of the class of bivariate modified power series distributions.
This paper has three-fold objectives. Firstly, we derive the locally powerful score test on the hypothesis that the correlation parameter is significantly different from zero, and then assess the power of this test. Secondly, construct confidence interval on the ratio of the means of the two correlated Poisson variables. Thirdly, we fit the model to metabolic syndrome data available from a random sample of spousal pairs.
2. The Bivariate Generalized Poisson Distribution (BGPD)
A bivariate random variable (X, Y) is said to have BGPD when the joint probability distribution is given by (1):
(1)
,
,
,
, and zero otherwise. The BGPD can be obtained as follows:
By expanding the bivariate probability generating function (BPGF)
, using the bivariate Lagrange expansion of an implicit function, as was shown in [4].
3. Recurrence Relation among Moments
From (1) we can write
(2)
Therefore, from (2) we can see that the distribution posses a power series expansion in terms of the parameters (
&
).
On differentiating the above equation {2} partially w.r.t.
and
respectively, dividing by
, and on summation, we get by solving for
and
.
(3)
(4)
To obtain a recurrence relation among the higher non-central product moments, we write
(5)
Define the cross product moments:
.
We can show that
(6)
By symmetry we can write:
(7)
These relationships among the moments are obtained on following Kendall [7].
In particular, from (6) and (7):
(8)
(9)
are respectively the mean and variance of X in terms of the model parameters.
And the values of
and
can be written down by symmetry. The coefficient of correlation
, between X and Y is given by
(10)
As can be seen, from the restrictions on the model parameters, the correlation (10) cannot be negative.
3.1. Maximum Likelihood Estimators
Assuming
and
to be known constants, let
be a random sample of fixed size
, taken from the BGPD family given in (1). It can be easily shown that the maximum likelihood (ML) estimators for
and
are as shown in (11) and (12):
(11)
(12)
With some considerable algebra, converting the Fisher’s information matrix we can show that the asymptotic variances and covariances of the ML estimators of
and
are given respectively in (13)-(15):
(13)
(14)
(15)
3.2. Properties of the BGPD
1: If
, then
and
are stochastically independent if and only if
.
2: If
, then
and
cannot be perfectly correlated.
Proof. The quality of
, given in (10), to unity, will imply that
. This is not true unless
, which is a contradiction to the condition enforced by the strict inequality
.
3: If
is a random sample taken from the BGPD family, the joint probability distribution of the sums
and
is also a GDPD and its probability function takes the form
(17)
4: The marginal distribution of
is the generalized Poisson distribution given by Consul and Jain [8], where their
and
are taken as
and
respectively. Thus, the marginal probability distribution of
and
are given respectively in (18) and (19):
(18)
Similarly
(19)
5: The regression equation of
on
is given by
3.3. Score Test on the Hypothesis of Independence (H0:m = 0)
In the subsequent analyses we shall assume that
. Let
be an srs. The likelihood is given by:
The log-likelihood function is:
(20)
Moreover denote the sample means by:
and
.
Differentiating the log-likelihood function (20) with respect to m, the score function
is given by:
Replacing
with their MLE under
;
, we have
(21)
Under
, from (21) we have:
Since
And under
Therefor,
, and the variance of
under the null hypothesis is:
Hence
(21)
as
the distribution of the test statistic (21) follows Chi-square distribution with one degree of freedom.
To evaluate the asymptotic power of the score test we need to find the non-null asymptotic distribution of
.
When the null hypothesis does not hold, the expected value of, denoted by
, is such that:
(22)
Moreover, denoting the variance of the test statistics
under the non-null distribution by
, then
, where,
(23)
An expression for the estimated sample size needed to verify the hypothesis of independence for given type I error rate and power 1-β is shown to be:
(24)
For example, in (24) given 5% type error rate
, and power 80%,
.
In the Appendix we give the algebraic expressions of the elements of D as functions of the population parameters.
3.4. Confidence Interval on the Ratio of BGPD Means
1) Delta Method
Let
denote the ratio of the mean of X to the mean of Y, The point estimator
. From Kendall [7], we have to the first order of approximation:
(25)
The (1 − α) 100% confidence interval on R is therefore given by:
(26)
2) Feiller’s Interval
Let
Hence, the
A point estimator of
is given by:
Therefore,
where
and
Solving the quadratic in R we get
, where
(27)
(28)
In the data analysis section, we shall compare the above intervals to the non-parametric bootstrap confidence interval.
4. Application: Modeling Metabolic Syndrome in Spousal Pairs Using GDPD
The metabolic syndrome is the co-aggregation of hypertension, impaired glucose tolerance, dyslipidemia, and abdominal obesity and is associated with an increased risk of total and cardiovascular mortality in adults [9] [10]. Genetics as well as environmental influences have been implicated in obesity and several cardiovascular risk factors [11] [12]. Family is one of the most important factors affecting metabolic risk factors in children, in that family displays an interaction between genetic and shared environmental factors [13] [14]. Recent research showed that childhood and adolescent overweight has been increasing in Asian countries due to urbanization and economic development. For example, over the past 10 years, the rates of overweight among Korean children and adolescents aged 5 - 20 years have doubled [15], which may ultimately cause an increase in adverse cardiovascular outcomes. Globally, the prevalence of the metabolic syndrome is high among obese children and adolescents, and increases with increasing obesity [16]. The widely used definition of metabolic syndrome is that of the World Health Organization [17] [18]. The components of each definition and criteria for making the diagnosis of the metabolic syndrome are summarized in Table 1.
The abbreviations of the medical official bodies given in Table 1 are:
IDF ≡ International Diabetes Federation, WHO ≡ World Health Organization, EGIR ≡ European Group for the Study of Insulin Resistance, and NCEP ≡ National Cholesterol Education Program.
Table 1. Definition of the components of metabolic syndrome [Journal of the Royal Society of Medicine Vol. 99 September 2006].
Component |
IDF |
WHO |
EGIR |
NCEP |
M |
F |
M |
F |
M |
F |
M |
F |
Central Obesity |
≥102 |
≥88 |
≥102 |
≥88 |
≥94 |
≥80 |
≥102 |
≥88 |
Raised TG |
≥1.7 |
≥1.7 |
≥1.7 |
≥1.7 |
≥2.0 |
≥2.0 |
≥1.7 |
≥1.7 |
Low HDL |
<1.03 |
<1.29 |
≤.9 |
≤1 |
≤1 |
≤1 |
≤1.03 |
≤1.29 |
Hypertension |
≥130/
85 |
≥130/
85 |
≥140/
90 |
≥140/
90 |
≥140/
90 |
≥140/
90 |
≥130/85 |
≥130/
85 |
Fasting Glucose |
≥5.6 |
≥5.6 |
≥6.1 |
≥6.1 |
≥6.1 |
≥6.1 |
≥6.1 |
≥6.1 |
The Metabolic syndrome data:
Based on a random sample of 4000 spousal pairs, the female and male counts are cross classified in Table 2 and Table 3.
The test independence had a value of U= 7.027, with a corresponding p-value= 2.11e−12.
Therefore the hypothesis of independence is not supported by the metabolic
Table 2. Cross classification of the counts of components of MS in spousal pairs.
|
Female count |
Total |
0.00 |
1.00 |
2.00 |
3.00 |
4.00 |
Male count |
0.00 |
2337 |
0 |
0 |
0 |
0 |
2337 |
1.00 |
1228 |
0 |
0 |
0 |
0 |
1228 |
2.00 |
0 |
360 |
0 |
0 |
0 |
360 |
3.00 |
0 |
0 |
44 |
12 |
0 |
56 |
4.00 |
0 |
0 |
0 |
7 |
12 |
19 |
Total |
3565 |
360 |
44 |
19 |
12 |
4000 |
Table 3. Summary statistics of counts.
|
N |
Mean |
Std. Deviation |
Variance |
Male count |
4000 |
0.5480 |
0.75421 |
0.569 |
Female count |
4000 |
0.1382 |
0.45353 |
0.206 |
syndrome data.
To fit the model to the data we shall use the method of moments to estimate the three population parameters. We shall use the notations
&
to respectively denotes the sample means of x & y, while
denotes the sample covariance between x & y.
(29)
(30)
(31)
From Equations (29)-(31), the admissible moment estimator for
is.
(32)
The moments estimators of
&
are given respectively by
&
in Equations (33) and (34).
(33)
(34)
From the metabolic syndrome data we have, for females,
, and for male,
, and
. The familiar chi-square goodness of fit, with 4 degrees of freedom had a p-value = 0.141. This shows that the model gives a reasonable fit to the data.
The following is a short R code from which we evaluate the nonparametric confidence interval and assess the large sample distribution of the ratio of correlated Poisson means.
R-CODE for bootstrapping the confidence interval on the ratio of means
x<-data$Male_count
y<-data$Female_count
df<-data.frame(x,y)
library(boot)
boot_df=boot(df,function(i,d)mean(i[d] [1])/mean(i[d] [2]),R=999)
boot_df
Bootstrap Statistics :
original bias std. error
ratio 3.964 0.0163 0.148
It is clear from Figure 1 that the quantile plot of the sample ratio of means follows a normal distribution (Table 4).
Figure 1. Histogram and normal quantile plot for 1000 bootstrap replicates of the estimated ratio of means.
Table 4. The upper and lower limits of the 95% confidence interval on the ratio of means for the three methods (delta, Feiller’s, and bootstrap).
Limits |
Delta method |
Feiller’s |
Bootstrap |
Lower limit |
3.417 |
3.467 |
3.674 |
Upper limit |
4.514 |
4.575 |
5.254 |
5. Concluding Remarks
In this article, we discussed some of the interesting applications of a bivariate generalized Poisson distribution. The model has three parameters and has provided satisfactory fit to the metabolic syndrome data collected from spousal pairs. The common parameter (m) captures overdispersion and accurately accounts for the strength of the positive correlation between the two variables.
It is worth noting that this family of models, which is a member of the Lagrange family of bivariate distributions, is suitable for many applications, and in particular in situation where queuing theory is applicable. For example, in large tertiary hospital outpatients departments, the patients usually wait for admission form queues. If we have two types of queues (e.g. males and females), then the random variables would represent the numbers of admitted patients during a busy period of the admission office (single server queue). We also note that the marginal probability distributions have been applied to a variety of data from genomics, and case fatality due to COVID [19].
Finally, it is worth noting that with the not so complicated structure of the model, it can be shown that it belongs to the class of bivariate of generalized linear models where inclusion of covariates will make this model attractive alternative in many applications.
Appendix: Higher Order Cross-Product Moments
On using the recurrence relationship in (6) we can show that: Denote
where
and
And from the recurrence relation, we have:
The
depends on the:
&
We can show that:
It should also be noted that
and
can be obtained on exchanging
by
in
and
and visa-versa.