A Wald test for zero inflation and deflation for correlated count data from dental caries research

Abstract

Tests of homogeneity in zero-inflated models for count data have been well discussed in the literature, but the existing methodologies have relied primarily on score statistics. An often cited justification for the use of these statistics is that they do not require the model to be fitted under the alternative. But the advent of computer software with robust functions and procedures has made it easy for these alternative models to be fitted routinely in practice. In this article, we exploit this opportunity by using results generated from these analyses to develop a Wald test statistic to evaluate homogeneity in the class of zero-inflated models. We show how the proposed test can be performed with minimal programming effort for the practicing statistician. Technically, the test is based on a reparameterization of the mixing weights that naturally incorporates covariates under heterogeneity, a characteristic often ignored by existing testing procedures. A quasi-likelihood function derived from a working independence model coupled with the so-called sandwich estimator of the covariance matrix of the parameter estimates is used to accommodate correlation in the data. A simulation study is conducted to evaluate the empirical performance of the proposed Wald test and its real life applications are illustrated using correlated dental caries counts in young children.

Keywords

dental caries research Huber sandwich estimator negative binomial model Poisson model two-component mixture models working independence model zero inflation and deflation

1 Introduction

Count data with excess zeros relative to generalized linear models are common in population-based studies. A well accepted and established statistical technique to analyze such data is the use of zero-inflated models. These models view the data as being generated from a model resulting from the mixture of a distribution that places all its mass at zero and a non-degenerate distribution. Under this mixture model, the zero counts are generated from two distribution sources, resulting in a heterogeneous population. In this article, the term heterogeneity refers to this mixture in distribution, for which the extent is measured by the associated mixing weight. Likewise, the term homogeneity will be used to characterize a population generated from a non-degenerate distribution. In other words, under the homogeneous model, the mixing weight is equal to zero.

The choice of the non-degenerate distribution for the most part depends on the nature of counts under consideration. For bounded counts, a binomial distribution is typically used (Hall, 2000). In the case of unbounded counts, Poisson and negative binomial distributions are the standard (Böhning et al., 1999). Ridout et al. (1998) and references therein provide an extensive review of this literature.

In practice, it is often the case for count data with potentially extra zeros to be correlated due to the longitudinal or hierarchical structure of the study design. Extensions of zero-inflated models to accommodate correlated count data are well established. These correlated data models can be broadly classified into two categories: (1) models that explicitly model the association in the data using random effects terms, and as such allow a full likelihood approach to parameter estimation (Hall, 2000; Lee et al., 2006); and (2) models that treat association in the data as a nuisance but use the so-called sandwich estimator of the covariance matrix of the parameter estimates to correct for any misspecification of the working correlation model (Hall and Zhang, 2004; Hall and Shen, 2010).

In many applications of zero-inflated models for count data, it is often of interest to evaluate whether the inherent heterogeneity is consistent with observed data. Existing testing methodologies have relied primarily on constant alternatives under heterogeneity (Broek, 1995; Xiang et al., 2006). Although valuable, such a restrictive approach may lead to inconsistent estimates of the parameters of interest and may suffer from reduced power in detecting heterogeneity that varies with covariate profiles. Most importantly, existing testing procedures often rely on score statistics (van den Broek, 1995; Jansakul and Hinde, 2009; Todem et al., 2012). An often cited motivation for the use of score test statistics is that they only require the model to be fitted under the null hypothesis. But the advent of computer software with robust functions and procedures has made it easy for the alternative models to be fitted routinely in practice. As a result, conducting a Wald test may require significantly less programming effort than conducting a score test. Another competing testing procedure is the likelihood ratio test, which like the Wald test enjoys the ease of implementation in real life applications. But in settings where a likelihood is not assumed and a general estimating function is used for model estimation, this testing procedure is not an option. In this article we exploit estimation results generated from the computer software to develop a Wald test statistic to evaluate homogeneity in this class of models. Technically, the test is based on a reparameterization of the mixing weights that naturally incorporates covariates into the mixing weights under the alternative, a characteristic often ignored by the existing testing procedures. We extend the methodology to account for the correlation in the data using a quasi-likelihood function derived from a working independence model coupled with the so-called sandwich estimator of the covariance matrix of the parameter estimates. From a practical perspective, we show how the proposed test can be performed with minimal programming effort for the applied statistician.

The work proposed in this article is motivated by dental caries data (discussed in Section 4.2) generated from a unique three-wave longitudinal study designed to collect oral health information on low-income African–American children under the age of 6 and their main caregivers residing in the city of Detroit. Todem et al. (2012) used these data to show how a test that neglects parameter heterogeneity under the alternative could fail to reject homogeneity in the population. They conducted their analysis using the first wave data and found that these data exhibit, relatively to a negative binomial model, an inflation of zeros for younger children and a deflation of zeros for older children that are roughly of the same magnitude. Thus, by averaging over the covariate space, a testing procedure based on constant mixing weights had virtually no power in detecting heterogeneity in this population. In this article, we develop a testing procedure to detect such forms of heterogeneity across the three waves of these dental caries data.

In Section 2, we describe the marginal representation of zero-inflated models that naturally accommodates both inflation and deflation of zeros in count data, and extend the model to correlated data. A Wald test statistic for homogeneity against varying alternatives is constructed and described in Section 3. In Section 4, we conduct a numerical study to evaluate the finite sample properties of the test statistic, and apply the methodology to evaluate homogeneity of dental caries scores in children from the Detroit longitudinal study. Some concluding remarks are given in Section 5. Finally, a sample SAS programme that can be used to fit the zero-inflated model under its marginal representation and conduct the test is given in the Appendix.

2 Zero-inflated models for correlated count data

The general data structure of interest is similar to that of many correlated data. There are n independent clusters indexed by i (i = 1, …, n), and within each cluster i there are λ_i correlated observations indexed by j (j = 1, …, λ_i). Let Y_ij be the random variable representing the count response for subunit j in cluster i drawn from M, a mixture of an atom distribution at 0 and an unknown non-degenerate discrete distribution H governed by a finite dimensional parameter vector ,

\Pr (Y_{ij} y_{ij}) \{\begin{matrix} _{ij} (1_{ij}) h_{ij} (0;) & if y_{ij} 0 \\ (1_{ij}) h_{ij} (y_{ij};) & if y_{ij} 0 . \end{matrix}

(2.1)

Here h_ij(y_ij; θ) represents the probability mass associated with the distribution H at observed response y_ij. Parameter ϕ_ij measures the extent of heterogeneity in the population and its support depends on the representation of the mixture model. Under the hierarchical representation of the model, ϕ_ij is a probability mass requiring the distributional constraints 0 ≤ ϕ_ij ≤ 1, j = 1, …, λ_i, i = 1, …, n. But the marginal representation of the model only requires the condition 0 ≤ Pr(Y_ij = y_ij) ≤ 1 which results in constraints on mixing weights,

\frac{h_{ij} (0;)}{1 h_{ij} (0;)}_{ij} 1, j 1,, m_{i}, i 1,, n .

(2.2)

The constraint in (2.2) simply allows the model to accommodate not only inflation (positive ϕ_ij) but also deflation (negative ϕ_ij), a property that the hierarchical representation of the model, implied by the classical constraints 0 ≤ ϕ_ij ≤ 1, does not have. In this article, we will focus on the marginal representation of the model which has the ability to accommodate both inflation and deflation at zero (Dietz and Böhning, 2000). We follow Todem et al. (2012) and adopt the transformation,

_{ij} \frac{_{ij} h_{ij} (0;)}{1 h_{ij} (0;)}, j 1,, m_{i}, i 1,, n,

(2.3)

where π_ij ∈ [0, 1] represents the probability at 0 under the model in (2.1). One appealing feature of this transformation is that it naturally embeds the mixing weight constraints into the model, therefore avoiding the need for a constrained optimization.

Regarding the dependency in the data, one popular approach to accommodate the inherent within-cluster association consists of assuming random effects coupled with the so-called conditional independence assumption to tie together all observations from the same cluster (Lee et al., 2006; Xiang et al., 2006; Xiang et al., 2007). This approach is full likelihood and explicitly models the association in the data, which can be cumbersome in some practical settings. Alternatively, a quasi-likelihood that treats the within-cluster association as a nuisance but generates consistent estimates of parameters of interest is often advocated when the full likelihood is intractable. In this article, we adopt such an approach and consider a working independence model (WIM) for which the contribution of cluster i to the quasi-likelihood function is:

L_{WIM} (y_{i};)_{j 1}^{m_{i}} \Pr (Y_{ij} y_{ij}),

where y_i = (y_i₁, …, y_imi) and $R^{q}$ represents the parameters used to describe the marginal probabilities Pr(Y_ij = y_ij), j = 1, …, λ_i, i = 1, …, n. It is well known that under a correct specification of the marginal distribution Pr(Y_ij = y_ij), the WIM estimator $\overset{}{}$ is consistent for }*, the true value of }. Specifically under mild regularity conditions, $\sqrt{n} (\overset{}{}^{*})$ converges in distribution as $n$ to a centered Gaussian distribution with asymptotic variance-covariance matrix $? (^{*}) \lim_{n} n [I_{n} (^{} {)]}^{1} V_{n} (^{*})[I_{n} (^{*} {)]}^{1}$ , where $I_{n} (^{*})$ represents the expected information matrix $E [_{i 1}^{n}_{j 1}^{m_{i}} s_{ij} ()/]$ evaluated at $^{*}$ with $s_{ij} () \log (\Pr (Y_{ij} y_{ij}))/$ and $V_{n} (^{*})_{i 1}^{n} E [_{j 1}^{m_{i}} s_{ij} (^{*})]$ ${[_{j 1}^{m_{i}} s_{ij} (^{*})]}^{t}$ , where t stands for transposition (details of the derivation of the asymptotic normality, see, for example, Freedman, 2006). A consistent estimator of this asymptotic variance–covariance matrix that corrects for any misspecification of the true association structure is the so-called Huber sandwich estimator (Freedman, 2006),

\overset{}{?} (\overset{}{}) n [{\overset{}{I}}_{n} (\overset{}{} {)]}^{1} {\overset{}{V}}_{n} (\overset{}{})[{\overset{}{I}}_{n} (\overset{}{} {)]}^{1},

(2.4)

where ${\overset{}{V}}_{n} (\overset{}{})_{i 1}^{n} [_{j 1}^{m_{i}} s_{ij} (\overset{}{})] {[_{j 1}^{m_{i}} s_{ij} (\overset{}{})]}^{t}$ and ${\overset{}{I}}_{n} (\overset{}{})$ is $_{i 1}^{n}_{j 1}^{m_{i}} s_{ij} ()/$ evaluated at $\overset{}{}$ .

In general, we have $V_{n} (^{*})_{i 1}^{n} E \{[_{j 1}^{m_{i}} s_{ij} (^{})] {[_{j 1}^{m_{i}} s_{ij} (^{})]}^{t}\}$ . However, if the true correlation within-cluster is zero, we then have $E \{[_{j 1}^{m_{i}} s_{ij} (^{*})]$ ${[_{j 1}^{m_{i}} s_{ij} (^{*})]}^{t}\}_{j 1}^{m_{i}} E [s_{ij} (^{*})][s_{ij} (^{*} {)]}^{t}$ . Then $V_{n} (^{*})$ is simply equal to the expected information matrix $I_{n} (^{*})$ . As a result, the sandwich estimator reduces to $\overset{}{?} (\overset{}{}) n [{\overset{}{I}}_{n} (\overset{}{} {)]}^{1}$ which is the so-called model-based estimator.

3 A Wald test of homogeneity for correlated count data

We are interested in evaluating the model H against the model M in (2.1) coupled with the constraints in (2.2). We are specifically interested in the two-sided hypotheses,

H_{0} :_{ij}^{*} 0 for all i, j v s . H_{1} :_{ij}^{*} 0 for some i, j,

where $_{ij}^{*}$ is the true value of $_{ij}$ . Throughout this article, the term heterogeneity will be used to characterize the distribution M in (2.1) and that of homogeneity to characterize the non-degenerate discrete distribution H. In other words, under the homogeneous model H, the true mixing weights are all zero.

Evaluating these hypotheses is not trivial, owing to their high dimensional nature under the alternative. To address this curse of dimensionality, existing inferential procedures typically neglect parameter heterogeneity by assuming $_{ij}$ for all i, j, the so-called constancy assumption under heterogeneity (e.g., van den Broek, 1995). Although valuable, such restrictive approaches may lead to inconsistent estimators of parameters of interest and may have reduced power in detecting heterogeneity that varies across sampling units.

In this article, we consider a general formulation where the parameters ϕ_ij, j = 1, …, λ_i, i = 1, …, n are allowed to depend on covariates. Such an approach has been proposed by Jansakul and Hinde (2009) who related the mixing weight to covariates via a regression technique. A limitation however is that the identity link function was used and estimation under the alternative model required a constrained optimization to satisfy conditions in (2.2). To avoid this numerical complexity, we adopt the transformation in (2.3) which naturally incorporates covariates into the testing scheme. Under this transformation, the hypotheses under evaluation become,

H_{0} :_{ij}^{*} h_{ij} (0;^{*}) for all i, j vs . H_{1} :_{ij}^{*} h_{ij} (0;^{*}) for some i, j,

(2.5)

where $_{ij}^{*}$ and i* represent the true values of π_ij and i. The new null hypothesis is simply obtained by setting ϕ_ij = 0 for all i, j in (2.3), which is equivalent to π_ij= h_ij (0; i) for all i, j. Likewise, the new alternative hypothesis can be obtained when ϕ_ij ! 0 for some i, j in (2.3).

To allow ϕ_ij to depend on covariates, we relate π_ij and h_ij (0; i) to covariate vectors x_ij and z_ij, respectively, of dimensions r and s, observed alongside the count outcomes y_ij, j = 1, …, λ_i, i = 1, …, n. Covariate x_ij is related to the probability h_ij (0; i) through the regression model $_{ij} b^{1} (x_{ij}^{t})$ , where $_{ij}$ represents the mean of Y_ij under the distribution H, b(.) is a known link function and b is a subset of i. Likewise, parameter π_ij is related to covariate ϕ_ij via the regression $_{ij} g (z_{ij}^{t})$ , where g(.) is an increasing, differentiable and invertible function. Following Todem et al. (2012), the function g(.) is judiciously chosen to mimic h_ij (0; i) in such a way that the homogeneity hypothesis above is reduced to a hypothesis involving only a linear combination of model parameters. As an example, suppose that H is a Poisson process with b(.) = log(.), we have $_{ij} \exp {x_{i j}^{t}}$ and $h_{ij} (0;) \exp {_{ij}} \exp {\exp {x_{ij}^{t}}}$ with b = i. Thus, a natural parameterization that mimics the regression model for h_ij (0; i) is $_{ij} \exp {\exp {z_{ij}^{t}}}$ . Under these formulations, the homogeneity hypothesis simply reduces to,

H_{0} : z_{ij}^{t}^{*} x_{ij}^{t}^{*}, for all i, j,

(2.6)

where γ* and β* represent the true values of γ and β. For the Poisson model, g(.) is known, but is in general unknown. When H is a negative binomial process, g(.) is a parametric function indexed by the overdispersion parameter of the distribution H. A general form of these transformations can be found in Todem et al. (2012). In the analysis of dental caries data, we give an explicit form of this transformation when H is a negative binomial process.

Based on (2.6), the new hypothesis formulation can be further rewritten as a linear combination of elements of the parameter vector $(^{t},^{t})^{t}$ , $H_{0} : C^{*} 0,$ where C is a contrast matrix of dimension of p × q. For example, when x_ij= ϕ_ij for all i, j, we have $C^{*}^{}^{*}$ and the homogeneity hypothesis is equivalent to H₀ : b* = c*.

A Wald test statistic of homogeneity that accommodates correlated count data is,

W_{n} (\sqrt{n} C \overset{}{})^{t} {\{C \overset{}{?} (\overset{}{}) C^{t}\}}^{1} (\sqrt{n} C \overset{}{}),

where $\overset{}{?} (\overset{}{})$ is defined in (2.4). This test statistic under the null hypothesis follows a |² distribution with rank(C) degrees of freedom. The limiting distribution of the Wald test statistic as a |² distribution is guaranteed through the asymptotic normality of the maximum likelihood estimator of the parameter vector $\overset{}{}$ (see, for example, Shao, 1999, p. 433–434). The null hypothesis is then rejected in favour of the alternative for large values of the statistic W_n. Note that the sandwich estimator is the only component of the Wald test statistic that is affected when the true within-cluster correlation is zero.

4 Numerical studies

4.1 Computations

Statistical software with robust optimization functions and subroutines can be used to fit the zero-inflated model under the marginal representation described in Section 2 and to perform the proposed Wald test. For users familiar with SAS® software (SAS Institute, Cary, NC, USA), the procedure PROC NLMIXED provides a convenient module to compute the sandwich estimator of the variance-covariance matrix of the model estimates using the ‘EMPIRICAL’ option. This option, however, requires the practicing statistician to specify random effects terms in the model formulation. To circumvent this requirement, the user can elect to use one random effect term for each cluster and request that only one non-adaptive Gaussian quadrature point can be used for the numerical integration by using options ‘QPOINTS=1’ and ‘NOAD’. With these specifications and restrictions, only the point zero is used for the numerical integration of the full likelihood, resulting in the WIM likelihood. Once the model is fitted and the sandwich variance-covariance matrix estimate obtained, the Wald statistic can be readily obtained by multiplying the observed F statistic under the ‘CONTRAST’ statement by the rank of the contrast matrix. An illustration of this implementation in SAS^® is provided in the Appendix. An alternative approach is to output model estimates (parameter estimates and corresponding variance-covariance matrix estimate) to a dataset and use any matrix oriented software such as PROC IML to compute the Wald statistic and perform the test.

4.2 Application to the Detroit Dental Caries data

4.2.1 Study description

We apply the proposed methodology to data generated from a longitudinal survey conducted by the Detroit Center for Research on Oral Health Disparities, hosted at the University of Michigan, Ann Arbor. This is a unique three-wave study designed to collect oral health information on low-income African–American children (0–5 years) and their main caregivers (14+ years), living in the city of Detroit. The overall goal of the study is to promote oral health and reduce its disparities within the community of low-income African Americans through the understanding of the social, familial, biological, and neighbourhood determinants of dental caries and periodontal disease. Data on oral examinations involving assessment of dental caries, periodontal disease, oral cancer and oral hygiene were collected during each of the three waves of examinations and interviews completed between 2002–2007. For our illustration, we focus on dental caries scores which represent the cumulative severity of tooth decay for each surveyed child. Specifically, we focus on the number of decayed surfaces (DS) showing signs of clinically detectable enamel lesions comprising both noncavitated and cavitated lesions. These dental caries scores have well documented shortcomings but they continue to be instrumental in evaluating and comparing the risks of dental caries across population groups (Lewsey and Thomson, 2004). The DS score represents the generic response variable Y_ij introduced in earlier sections. Considered covariates include AgeB (child’s age at baseline in years), SC (child’s daily sugar consumption, a wave dependent variable), Wave2 (second wave indicator), and Wave3 (third wave indicator). A more detailed description of the study can be found in the papers of Sohn et al. (2007) and Ismail et al. (2011).

4.2.2 Analysis of the Detroit Caries Indices in early childhood

Our numerical computations are based on data collected on 1,021 surveyed children. Dental caries scores collected across children are assumed independent but those collected on the same child across the three waves are correlated. It is a common practice to use zero-inflated negative binomial models to describe dental caries scores in children (Böhning et al., 1999). In this article, we are interested in evaluating whether or not the inclusion of mixing weights is consistent with observed data. We first fit a zero-inflated negative binomial model with constant mixing weight $_{ij}$ coupled with the sandwich estimate for the variance-covariance matrix of the parameter estimates. The associated negative binomial model has mean $_{ij} \exp {x_{ij}^{t}}$ , where x_ij = (1, AgeB_i, SC_ij, Wave2_ij, Wave3_ij, AgeB_i*Wave2_ij, AgeB_i*Wave3_ij, AgeB_i* SC_ij)^t and $(_{0},_{1},,_{7})^{t}$ , and overdispersion parameter l. The constant mixing weight test fails to reject the homogeneity hypothesis $^{*} 0$ , with $^{}$ being the true value of , for the outcome of interest with an observed statistic of 2.37 and a p-value of 0.12 (the limiting distribution is a |² distribution with one degree of freedom). A thorough evaluation of this model reveals that the constant mixing weight does not provide an appropriate representation of zeros in this population. In Figure 1, we compare the observed proportions to the predicted probabilities from a homogeneous negative binomial model for each age group throughout the three study waves. This informal model assessment reveals that for younger children there is an inflation of zeros and a deflation for the older age group throughout the three study waves. One may argue that a zero-inflated negative binomial model with a heterogeneous overdispersion parameter may well accommodate both inflation and deflation of zeros in the data without allowing ϕ_ij to be negative. We evaluated such a claim and found that by allowing the overdispersion parameter to depend on available covariates using the following log-linear model, $_{ij} \exp {x_{ij}^{t}}$ , where a = (a₀, a₁, …, a₇)^t. A Wald test with 7 degrees of freedom revealed that these covariates were not important in predicting the overdispersion parameter (Wald test statistic=11.2, p-value = 0.131).

Figure 1

Source: Authors’ own.

The analyses above then draw into question the efficiency of models that ignore this varying heterogeneity under alternatives. We fit the proposed zero-inflated negative binomial model with a covariate-dependent mixing weight,

_{ij} \frac{\{1 \exp \{x_{ij}^{t}\}}^{1 /} \{1 \exp \{x_{ij}^{t}\}}^{1 /}}{1 \{1 \exp \{x_{ij}^{t}\}}^{1 /}},

where γ = (γ₀, γ₁, …, γ₇)^t. Note here that ϕ_ij = x_ij and the regression model $_{ij} {1 \exp {x_{ij}^{t} {}}}^{1 /}$ mimics the form of $h_{ij} (0;) {1 \exp {x_{ij}^{t} {}}}^{1 /}$ , and the function g(.) depends on l. Under the null hypothesis $x_{ij}^{t}^{} x_{ij}^{t}^{*} for all i, j$ , where c* and b* represent the true values of γ and b. Compared to a $^{2}$ distribution with 8 degrees of freedom, the proposed Wald test statistic W_n rejects the homogeneity hypothesis at 5% significance level with an observed value of 627 (p-value < 0.0001), providing a strong evidence of heterogeneity in this population. This is a good example where a constant mixing weight model is unable to capture varying heterogeneity in the population, but its counterpart that accommodates inflation and deflation at zero greatly improves efficiency in detecting the underlying heterogeneity. Model parameter estimates, corresponding standard errors, and model summary statistics are displayed in Table 1. The QIC statistic (Pan, 2001) provides additional evidence that the proposed model coupled with the transformation (2.3) fits the data reasonably well compared to the null model and the constant mixing weight model. These analyses however by no means imply that the proposed two-component model provides the best representation for these data.

Table 1

The fitted^‡ homogeneous, constant mixing weight and covariate-dependent mixing weight models to longitudinal dental caries data

	Null model		Constant mixing weight model		Covariate-dependent mixing weight model
Parameters	Estimates	S.E.	Estimates	S.E.	Estimates	S.E.
β ₀	–1.03*	0.18	–0.92*	0.20	1.60*	0.17
β ₁	4.36*	0.29	4.26*	0.31	0.62*	0.28
β ₂	2.23*	0.58	2.05*	0.61	0.17	0.51
β ₃	2.14*	0.17	2.16*	0.17	0.44*	0.17
β ₄	2.72*	0.19	2.67*	0.19	0.51*	0.18
β ₅	–2.99*	0.27	–3.05*	0.28	–0.74*	0.27
β ₆	–4.14*	0.30	–4.08*	0.31	–1.00*	0.29
β ₇	–2.95*	0.94	–2.68*	0.98	0.38	0.84
κ	1.58*	0.07	1.39*	0.13	0.86*	0.05
	–	–	0.04	0.03	–	–
γ ₀	–	–	–	–	–3.24*	0.19
γ ₁	–	–	–	–	6.37*	0.44
γ ₂	–	–	–	–	5.70*	1.22
γ ₃	–	–	–	–	2.28*	0.22
γ ₄	–	–	–	–	4.07*	0.30
γ ₅	–	–	–	–	–1.95*	0.45
γ ₆	–	–	–	–	–5.49*	0.55
γ ₇	–	–	–	–	–8.04*	2.00
QIC	12852.70		12853.19		12311.71

Source: Authors’ own.

Notes: ‡: AgeB and SC are rescaled to range 0–1

*: p-value < 0.05, ϕ-test

T: constant mixing weight

QIC $2_{i 1}^{n} \log {L_{WIM} (y_{i}; \overset{}{})} 2 trace {({\overset{}{I}}_{n} (\overset{}{} {))}^{1} n^{1} \overset{}{?} (\overset{}{})}$ .

4.3 Simulations

We conducted a numerical study to evaluate the finite sample performances of the proposed Wald test statistic W_n for homogeneity. These performances are compared to those of the Wald test statistic $W_{n}^{c}$ that assumes $_{ij}$ , the constancy assumption under heterogeneity. The two methods are further compared using the integrated mean squared error (IMSE) of the estimated mixing weight relative to its true counterpart throughout the covariate profile. Throughout the simulations, two true mixing weight configurations are considered, the constant and the covariate-dependent mixing weights.

To keep the simulation simple, unlike in the real data analysis, the homogeneous model is a Poisson process. To simplify the data generating mechanism for the correlation model, the systematic components of the data such as the true mean of the homogeneous Poisson model and the true mixing weights may vary with each cluster i but not with the subunit j. In order words, the covariates are allowed to be cluster level but not individual level. For example, the true mean of the Poisson process takes the form $_{ij}^{*} \exp {_{0}^{*} 0.5 x_{i}}$ , where the intercept $_{0}^{*}$ takes values in the set {0.25, 0.5, 0.75} and x_i, assumed to be constant across all j = 1, …, m, is a covariate generated from a uniform distribution on (0, 1). In the sequel, we will just write $_{i}^{*}$ and $_{i}^{*}$ to emphasize that these quantities are constant across j. Independent counts from the zero-inflated Poisson model can be generated from the following model,

\Pr (Y_{ij}^{IND} k) \{\begin{matrix} _{i}^{*} (1_{i}^{*}) \exp {_{i}^{*}} & if k 0, \\ (1_{i}^{*}) \frac{\exp \{_{i}^{*}\} (_{i}^{*})^{k}}{k!} & if k > 0 . \end{matrix}

To induce correlation among the data Y_i₁, Y_i₂, …, Y_im from the same unit i, we use independent Bernoulli variables coupled with independent zero-inflated Poisson processes. Let U_i be an independent random variable with the same distribution as Y_ij^IND and a_ij be an independent Bernoulli variable with success probability t ∈ (0, 1). Define $Y_{ij} (1 a_{ij}) Y_{ij}^{IND} a_{ij} U_{i}$ , where a_ij and U_i, j = 1, …, m, i = 1, …, n are all independent across the variables and indexes. Note that observations from the same cluster i, share the same U_i, which therefore induces correlation in the data. It can be readily shown that $E (Y_{ij}) (1_{i}^{*})_{i}^{*}$ , $Var (Y_{ij}) (1_{i}^{*})(1_{i}^{*}_{i}^{*})_{i}^{*}$ , where $_{i}^{*}$ and $_{i}^{*}$ are constant across j and $Corr (Y_{ij}, Y_{i j^{}})^{2}$ , $j j^{}$ . In our simulations, the correlation is set to t² = 0.49.

A popular alternative for introducing correlation in the data is the use of random effects coupled with the so-called conditional independence assumption (see, for example, Hedeker et al., 1994; Lee et al., 2006; Xiang et al., 2006). This approach is not chosen here because of its potential effects on the mean structure when random effects are integrated out. The mechanism adopted here introduces correlation in the data without affecting the marginal mean structure. This is only achieved in our scheme when the mean structure is constant across the index j.

In our simulations, we perform the proposed Wald test assuming the working mixing weight $_{i} {_{i} \exp {_{i}}}{1 \exp {_{i} {}}}^{1}$ under the alternative, where $_{i} \exp {_{0}_{1} x_{i}}$ and $_{i} \exp {\exp {_{0}_{1} x_{i}}}$ . The constant mixing weight test is performed assuming $_{i}$ under the alternative. Under these parameterizations, the null hypotheses are given by $H_{0} :_{0}^{*}_{0}^{*},_{1}^{*}_{1}^{*}$ for our model and $H_{0} :^{*} 0$ for the constant mixing weight model. Note that the estimation under the constant mixing weight model is subject to constraints in (2.2). Since the covariate x_i is between 0 and 1, these constraints reduce to $\max {(1 \exp {\exp {_{0} {}})}^{1},(1 \exp {\exp {_{0}_{1} {}})}^{1}} 1$ .

Finally, all simulations are replicated 1,000 times for sample sizes n = 100, 200 and 400 with m = 3 replicates within each cluster.

To study the empirical size of the proposed tests, we generate data under the null hypothesis by setting $_{i}^{*}^{*} 0$ . All considered Wald tests have well controlled type I error rates at 0.05% significance level (Table 2). We further study the empirical power of the tests to detect various forms of heterogeneity in the population. First, we consider the situation where the true mixing weights are constant. As expected, the Wald test based on constant mixing weights outperforms its covariate-dependent counterpart. The IMSEs also appear to be smaller for the constant mixing weight estimate relatively to the covariate-dependent mixing weight estimate. Second, we consider the situation where the true mixing weights vary with the covariate profile under various transformations. Results of these analyses are presented in Table 3. The covariate-dependent Wald test outperforms the covariate-independent Wald test in detecting the considered varying alternatives. Furthermore, the IMSEs of the constant mixing weight estimates are higher than those of the covariate-dependent mixing weight estimates.

Table 2
Empirical size and power of Wald test statistics to detect constant mixing weight $_{i}^{}^{}$ at 5% significance level, and integrated mean squared errors (IMSEs) for estimating $^{}$ using the constant mixing weight model and the proposed covariate-dependent mixing weight model

$_{0}^{}$ n = 100
n = 200
n = 400

0.25 0.5 0.75 0.25 0.5 0.75 0.25 0.5 0.75

$_{i}^{}^{} 0$

Size $W_{n}^{c}$ 0.046 0.043 0.056 0.044 0.049 0.047 0.049 0.052 0.054

$W_{n}$ 0.056 0.047 0.058 0.050 0.051 0.056 0.056 0.056 0.045

IMSE $\overset{}{}$ 0.009 0.005 0.004 0.004 0.002 0.001 0.002 0.001 0.001

${\overset{}{}}_{i}$ 0.023 0.011 0.011 0.011 0.005 0.003 0.005 0.003 0.001

$_{i}^{}^{} 0.1$

Power $W_{n}^{c}$ 0.182 0.293 0.432 0.338 0.484 0.743 0.580 0.778 0.952

$W_{n}$ 0.132 0.216 0.341 0.239 0.361 0.656 0.447 0.674 0.908

IIMSE $\overset{}{}$ 0.010 0.005 0.003 0.004 0.003 0.001 0.002 0.001 0.001

${\overset{}{}}_{i}$ 0.023 0.012 0.007 0.010 0.006 0.003 0.005 0.003 0.002

$_{i}^{}^{} 0.2$

Power $W_{n}^{c}$ 0.551 0.737 0.927 0.796 0.958 0.997 0.982 1.000 1.000

$W_{n}$ 0.414 0.627 0.851 0.683 0.905 0.992 0.943 0.996 1.000

IMSE $\overset{}{}$ 0.009 0.005 0.003 0.004 0.002 0.002 0.002 0.001 0.001

${\overset{}{}}_{i}$ 0.022 0.012 0.007 0.010 0.006 0.003 0.004 0.003 0.002

	$_{0}^{*}$	n = 100	n = 200	n = 400
		$_{i}^{}^{} 0$
Size	$W_{n}^{c}$	0.046	0.043	0.056	0.044	0.049	0.047	0.049	0.052	0.054
	$W_{n}$	0.056	0.047	0.058	0.050	0.051	0.056	0.056	0.056	0.045
IMSE	$\overset{}{}$	0.009	0.005	0.004	0.004	0.002	0.001	0.002	0.001	0.001
	${\overset{}{}}_{i}$	0.023	0.011	0.011	0.011	0.005	0.003	0.005	0.003	0.001
		$_{i}^{}^{} 0.1$
Power	$W_{n}^{c}$	0.182	0.293	0.432	0.338	0.484	0.743	0.580	0.778	0.952
	$W_{n}$	0.132	0.216	0.341	0.239	0.361	0.656	0.447	0.674	0.908
IIMSE	$\overset{}{}$	0.010	0.005	0.003	0.004	0.003	0.001	0.002	0.001	0.001
	${\overset{}{}}_{i}$	0.023	0.012	0.007	0.010	0.006	0.003	0.005	0.003	0.002
		$_{i}^{}^{} 0.2$
Power	$W_{n}^{c}$	0.551	0.737	0.927	0.796	0.958	0.997	0.982	1.000	1.000
	$W_{n}$	0.414	0.627	0.851	0.683	0.905	0.992	0.943	0.996	1.000
IMSE	$\overset{}{}$	0.009	0.005	0.003	0.004	0.002	0.002	0.002	0.001	0.001
	${\overset{}{}}_{i}$	0.022	0.012	0.007	0.010	0.006	0.003	0.004	0.003	0.002

Source: Authors’ own.

Notes: $W_{n}^{c}$ : covariate-independent test, with df=1

$W_{n}$ : proposed covariate-dependent test, with df=2

$\overset{}{}$ : constant mixing weight estimate

${\overset{}{}}_{i} \frac{\exp \{\exp \{{\overset{}{}}_{0} {\overset{}{}}_{1} x_{i}\}\} \exp \{\exp \{{\overset{}{}}_{0} {\overset{}{}}_{1} x_{i}\}\}}{1 \exp \{\exp \{{\overset{}{}}_{0} {\overset{}{}}_{1} x_{i}\}\}}$ : covariate-dependent mixing weight estimate.

Table 3

Empirical size and power of Wald test statistics to detect various specifications of $_{i}^{*}$ at 5% significance level, and integrated mean squared errors (IMSEs) for estimating the true mixing weight $_{i}^{*}$ using the constant mixing weight model and the proposed covariate-dependent mixing weight model

	$_{0}^{*}$	n = 100			n = 200			n = 400
	$_{0}^{*}$	0.25	0.5	0.75	0.25	0.5	0.75	0.25	0.5	0.75
		$_{i}^{*} 0.2 0.4 x_{i}$
Power	$W_{n}^{c}$	0.071	0.064	0.054	0.057	0.075	0.068	0.056	0.081	0.078
	$W_{n}$	0.198	0.313	0.530	0.320	0.530	0.819	0.559	0.841	0.985
IMSE	$\overset{}{}$	0.024	0.019	0.016	0.018	0.016	0.015	0.016	0.015	0.014
	${\overset{}{}}_{i}$	0.026	0.013	0.006	0.011	0.006	0.003	0.005	0.003	0.002
		$_{i}^{*} \exp {\exp {0.1 4 x_{i}}}$
Power	$W_{n}^{c}$	0.143	0.220	0.274	0.202	0.328	0.528	0.367	0.557	0.837
	$W_{n}$	0.188	0.264	0.445	0.308	0.510	0.768	0.550	0.821	0.982
IMSE	$\overset{}{}$	0.021	0.017	0.015	0.016	0.014	0.013	0.014	0.013	0.013
	${\overset{}{}}_{i}$	0.031	0.018	0.012	0.017	0.012	0.008	0.010	0.008	0.007
		$_{i}^{*} \exp {4 5 x_{i}}{1 \exp {4 5 x_{i} {}}}^{1}$
Power	$W_{n}^{c}$	0.261	0.318	0.464	0.426	0.589	0.787	0.697	0.880	0.972
	$W_{n}$	0.598	0.778	0.951	0.855	0.978	0.999	0.991	1.000	1.000
IMSE	$\overset{}{}$	0.079	0.074	0.070	0.072	0.069	0.066	0.068	0.068	0.065
	${\overset{}{}}_{i}$	0.023	0.014	0.009	0.014	0.009	0.006	0.009	0.006	0.004
		$_{i}^{*} (4 5 x_{i})$
Power	$W_{n}^{c}$	0.122	0.145	0.170	0.217	0.286	0.303	0.358	0.473	0.593
	$W_{n}$	0.525	0.748	0.925	0.826	0.966	0.999	0.985	1.000	1.000
IMSE	$\overset{}{}$	0.105	0.099	0.099	0.097	0.096	0.095	0.093	0.094	0.094
	${\overset{}{}}_{i}$	0.035	0.023	0.016	0.023	0.018	0.014	0.019	0.016	0.012

Source: Authors’ own.

Notes: $W_{n}^{c}$ : covariate-independent test (df=1)

$W_{n}$ : proposed covariate-dependent test (df=2)

$\overset{}{}$ : constant mixing weight estimate

It is worth noting that the power of the tests also improves when the two components are well separated (as b* increases) and with increasing sample sizes.

Finally, to evaluate the impact of the degree of correlation on the Wald tests, we perform the tests using data generated under different within-cluster correlations (t²), taking value in the set {0.25, 0.49, 0.81}. The results are given in Table 4. The empirical sizes appear not to be affected by the degree of correlation, however as we expected, the empirical powers of all considered tests are affected by t². The empirical power increases as the true correlation decreases, regardless of the specification of the true mixing weight. Similar findings have been reported in the literature (see, for example, Fitzmaurice, 1995).

Table 4

Empirical size and power of Wald test statistics to detect various specifications of $_{i}^{*}$ at 5% significance level under different degree of correlation t²

Correlation	ρ²	ρ²=0.25			ρ²=0.49			ρ²=0.81
Correlation	ρ²	n = 100	n = 200	n = 400	n = 100	n = 200	n = 400	n = 100	n = 200	n = 400
		$_{i}^{}^{} 0$
Size	$W_{n}^{c}$	0.043	0.046	0.045	0.043	0.049	0.052	0.038	0.050	0.043
	$W_{n}$	0.057	0.060	0.044	0.047	0.051	0.056	0.050	0.050	0.051
		$_{i}^{}^{} 0.1$
Power	$W_{n}^{c}$	0.357	0.627	0.881	0.293	0.484	0.778	0.219	0.396	0.661
	$W_{n}$	0.248	0.497	0.804	0.216	0.361	0.674	0.147	0.280	0.540
		$_{i}^{*} \exp {4 5 x_{i}}{1 \exp {4 5 x_{i} {}}}^{1}$
Power	$W_{n}^{c}$	0.443	0.721	0.948	0.318	0.589	0.880	0.263	0.498	0.772
	$W_{n}$	0.887	0.996	1.000	0.778	0.978	1.000	0.662	0.941	0.997
		$_{i}^{*} (4 5 x_{i})$
Power	$W_{n}^{c}$	0.192	0.343	0.611	0.145	0.286	0.473	0.146	0.193	0.365
	$W_{n}$	0.865	0.992	1.000	0.748	0.966	1.000	0.623	0.888	0.999

Source: Authors’ own.

Notes: $W_{n}^{c}$ : covariate-independent test, with df=1

$W_{n}$ : proposed covariate-dependent test, with df=2

$\overset{}{}$ : constant mixing weight estimate

is the CDF of a standard normal distribution.

5 Discussion

The issue of evaluating the homogeneity hypothesis in zero-inflated models has attracted considerable attention in the literature, but most of this debate has focused primarily on restrictive assumptions on the form of the alternative. One common limitation is that constant mixing weights are often assumed under heterogeneity. This article extends the literature by developing a Wald test statistic to detect heterogeneity that varies with the covariate profile in correlated count data. Most importantly, it shows how the constraints on the mixing weight under the proposed zero-inflated model that accommodates both inflation and deflation of zeros in the population can be naturally embedded into the testing procedure. We note however that there may be a slight deterioration of power when the true mixing weight does not depend on covariates. Given that the true data generating mechanism is usually unknown to the analyst, a general approach that captures various forms of heterogeneity in the population is the most conservative data analytic strategy.

It is worth noting that the zero-inflated model coupled with the transformation in (2.3) simply reduces to a two-part model where one part is a binary outcome model and another part is a truncated count model (see Cameron and Trivedi, 1998, pp. 123–25). However, unlike in the classical hurdle model, which is also a two-part model, the two components of the proposed model potentially share some parameters. For example, in the analysis of the Detroit data, the ‘hurdle’ probability π_ij depends on the overdispersion parameter of the negative binomial model, which violates the assumption of typical hurdle models. But as discussed in the article, such a restriction has the elegant property that both inflation and deflation can be evaluated against the homogeneous model.

Finally, one fundamental assumption of our testing methodology is that the mixing parameters are linearly related to covariates through the proposed transformation. Such an assumption may not be tenable in settings where the mixing weights are related to covariates through an unspecified function. For such a complication, smoothing techniques such as generalized additive models and spline models can be used to judiciously estimate the underlying relationship between the mixing weight and covariates (see, for example, Engle et al., 1986; Barry and Welsh, 2002; Xue et al., 2004; Lam et al., 2006). This extension and other generalizations of the methodology are outside of the scope of this article and are subject of further research.

Footnotes

Acknowledgements

The authors thank the two referees for their constructive comments and suggestions. We are also grateful to Amid Ismail for his permission to use the dental caries data. This work was supported by Dr Todem’s NCI/NIH K-award, 1K01 CA131259 and its supplement from the 2009 ARRA funding mechanism.

SAS Code Sample

References

Barry

Welsh

(2002) Generalized additive modelling and zero inflated count data. Ecological Modelling, 157(2), 179–88.

Böhning

Dietz

Schlattmann

Mendonca

Kirchner

(1999) The zero-inflated poisson model and the decayed, missing and filled teeth index in dental epidemiology. Journal of the Royal Statistical Society: Series A (Statistics in Society), 162(2), 195–209.

Cameron

Trivedi

(1998) Regression analysis of count data. New York: Cambridge University Press.

Dietz

Böhning

(2000) On estimation of the Poisson parameter in zero-modified Poisson models. Computational Statistics and Data Analysis, 34(4), 441–59.

Engle

Granger

CWJ

Rice

Weiss

(1986) Semiparametric estimates of the relation between weather and electricity sales. Journal of the American Statistical Association, 81, 310–20.

Fitzmaurice

(1995) A caveat concerning independence estimating equations with multivariate binary data. Biometrics, 51, 309–17.

Freedman

(2006) On the so-called ‘Huber sandwich estimator’ and ‘robust standard errors’. The American Statistician, 60(4), 299–302.

Hall

(2000) Zero-inflated Poisson and binomial regression with random effects: A case study. Biometrics, 56(4), 1030–39.

Hall

Shen

(2010) Robust estimation for zero-inflated poisson regression. Scandinavian Journal of Statistics, 37(2), 237–52.

10.

Hall

Zhang

(2004) Marginal models for zero inflated clustered data. Statistical Modelling, 4(3), 161–80.

11.

Hedeker

Gibbons

Flay

(1994) Random-effects regression models for clustered data with an example from smoking prevention research. Journal of Consulting and Clinical Psychology, 62(4), 757.

12.

Ismail

Lim

Sohn

(2011) A transition scoring system of caries increment with adjustment of reversals in longitudinal study: evaluation using primary tooth surface data. Community Dentistry and Oral Epidemiology, 39(1), 61–68.

13.

Jansakul

Hinde

(2009) Score tests for extra-zero models in zeroinflated negative binomial models. Communications in Statistics-Simulation and Computation, 38(1), 92–108.

14.

Lam

Xue

Bun

Cheung Y

(2006) Semiparametric analysis of zero-inflated count data. Biometrics, 62(4), 996–1003.

15.

Lee

Wang

Scott

Yau

McLachlan

(2006) Multi-level zero-inflated poisson regression modelling of correlated count data with excess zeros. Statistical Methods in Medical Research, 15(1), 47–61.

16.

Lewsey

Thomson

(2004) The utility of the zero-inflated poisson and zero-inflated negative binomial models: a case study of cross-sectional and longitudinal dmf data examining the effect of socio-economic status. Community Dentistry and Oral Epidemiology, 32, 183–89.

17.

Pan

(2001) Akaike’s information criterion in generalized estimating equations. Biometrics, 57(1), 120–25.

18.

Ridout

Demetrio

CGB

Hinde

(1998) Models for count data with many zeros. In Proceedings of international biometric conference (pp. 179–92). Cape Town, South Africa.

19.

Shao

(1999) Mathematical statistics. New York: Springer.

20.

Sohn

Ismail

Amaya

Lepkowski

(2007) Determinants of dental care visits among low-income African-American children. The Journal of the American Dental Association, 138(3), 309–18.

21.

Todem

Hsu

W-W

Kim

(2012) On the efficiency of score tests for homogeneity in two-component parametric models for discrete data. Biometrics, 68(3), 975–82.

22.

van den Broek

(1995) A score test for zero inflation in a Poisson distribution. Biometrics, 51, 738–43.

23.

Xiang

Lee

Yau

McLachlan

(2007) A score test for overdispersion in zero-inflated poisson mixed regression model. Statistics in Medicine, 26(7), 1608–22.

24.

Xiang

Lee

Yau

KKW

McLachlan

(2006) A score test for zero-inflation in correlated count data. Statistics in Medicine, 25(10), 1660–71.

25.

Xue

Lam

(2004) Sieve maximum likelihood estimator for semiparametric regression models with current status data. Journal of the American Statistical Association, 99(466), 346–56.