Modelling multivariate,overdispersed binomial data with additive and multiplicative random effects

Abstract

When modelling multivariate binomial data, it often occurs that it is necessary to take into consideration both clustering and overdispersion, the former arising from the dependence between data, and the latter due to the additional variability in the data not prescribed by the distribution. If interest lies in accommodating both phenomena at the same time, we can use separate sets of random effects that capture the within-cluster association and the extra variability. In particular, the random effects for overdispersion can be included in the model either additively or multiplicatively. For this purpose, we propose a series of Bayesian hierarchical models that deal simultaneously with both phenomena. The proposed models are applied to bivariate repeated prevalence data for hepatitis C virus (HCV) and human immunodeficiency virus (HIV) infection in injecting drug users in Italy from 1998 to 2007.

Keywords

Bayesian hierarchical model binomial data clustering MCMC overdispersion

1 Introduction

Clustering and overdispersion are two major issues that must be addressed when modelling data that cannot be assumed to be normally distributed, e.g., aggregated binary data and count data.

The clustering issue refers to the hierarchical structure of data, where measurements belonging to the same cluster are assumed to be associated. This issue can be accommodated using cluster-specific random effects, usually assumed to be normally distributed, which induce association between the repeated or multivariate measurements. Such models can be easily fitted within the framework of generalized linear mixed models (GLMMs, Breslow and Clayton, 1993; Molenberghs and Verbeke, 2005).

We encounter issues of overdispersion when the data present additional variability than that prescribed by the mean-variance relation of the distribution. The phenomenon of overdispersion has been widely considered in the literature, mostly in relation to the binomial and the Poisson distributions. Ignoring to account for overdispersion can lead to the underestimation of the standard errors and therefore to a wrong inference for the regression parameters. Possible solutions to this issue can be of two types (Hinde and Demétrio, 1998). The first approach consists in generalizing the variance function by including additional parameters, such as, the heterogeneity factor in overdispersed binomial data, and then estimating the regression parameters using quasi-likelihood methods (Agresti, 2002). The second approach assumes a two-stage model, where in the first stage we define for the data a distribution depending on certain parameters, whose distribution is then specified in the second stage. Examples are the beta-binomial model (Skellam, 1948) for binomial data and the negative binomial model (Breslow, 1984) for count data, but also some versions of the GLMMs can be used. A wide review of approaches that can deal with overdispersion can be found in Hinde and Demétrio (1998).

However, often interest may lie in simultaneously combining these two phenomena, clustering and overdispersion. Both marginal and random-effects models can be used to address them. If the focus is on the estimation of the fixed effects rather than on modelling the correlation structure, marginal models can be used to adjust the variance-covariance structure to accommodate clustering and overdispersion. For instance, both the GEE2 approach (Qaqish and Liang, 1992) and, especially suited for binomial data, alternating logistic regression (Carey et al., 1993) can be used to model the marginal means of the outcomes as well as the correlation between pairs of within-cluster measurements. Chen and Ahn (1997) go further in this direction and develop a marginal model for multivariate overdispersed binomial data where the mean structure depends on two multiplicative nested random effects. On the other hand, conditional models depending on random effects are a better choice if we are more interested in modelling the individual profiles rather than the population mean, as well as in estimating the correlation among the measurements. For instance, Booth et al. (2003) developed a negative binomial loglinear mixed model and they applied it to dependent and overdispersed count data. More generally, Molenberghs et al. (2007, 2010) proposed a class of extended GLMMs that accommodate clustering and overdispersion, making use of two separate sets of random effects, which then are estimated with maximum likelihood (ML) methods. These are meant to be used for modelling normal, binomial, Poisson and time-to-event data.

In this article, we focus on conditional models, and we thus extend the work presented in Molenberghs et al. (2010) by focusing on modelling multivariate, repeated and overdispersed binomial data. For this purpose, we develop a series of Bayesian hierarchical models that account for overdispersion through a set of random effects, either additively or multiplicatively included in the model, while dealing with clustering. We chose to use the Bayesian approach in order to avoid the computational difficulties encountered with GLMMs and ML estimation (Molenberghs et al., 2010). All the Bayesian hierarchical models presented in the article are fitted using Monte Carlo Markov Chain (MCMC) methods (Clayton, 1996). In this way, it is possible to specify a prior distribution for the unknown parameters, in particular for the overdispersion random effects and for their covariance matrix, and then calculate their posterior distribution through Gibbs sampling.

We apply the proposed methodology to a sample of prevalence data of hepatitis C virus (HCV) and human immunodeficiency virus (HIV) infection of injecting drug users (IDUs) in treatment from 20 Italian regions from 1998 to 2007. The article is organized as follows. In Section 2, we introduce the data. In Section 3, we describe the proposed methodology, giving details about the Bayesian hierarchical models with additive and multiplicative overdispersion parameters. Section 4 describes posterior computation, model selection and model checking. Section 5 is dedicated to the presentation of the main results, while in Section 6 we wrap up with a discussion of the proposed methodology. Model diagnostic and convergence assessment are presented in the supplementary material web appendix of the article.

2 Data

The data analyzed in the paper were reported to the European Monitoring Centre for Drugs and Drug Addiction (EMCDDA) and consist of diagnostic testing binomial data providing information about the HCV and HIV infection prevalence from samples of IDUs in treatment from 20 regions of Italy between 1998 and 2007. For each IDU, a serum specimen was taken and tested for antibodies against HCV and HIV. Further details about data collection can be found in Del Fava et al. (2011). Note that the data analyzed in this article are updated to 2007.

Figure 1 shows the observed prevalence profiles over the years for HCV and HIV infections, with a bold line representing the national prevalence profile, obtained by pooling together the regional results. We notice that the prevalence of HCV infection is much higher than the prevalence of HIV infection, reflecting the fact that HCV is reported to be about 10 times more infectious than HIV (Coutinho, 1998). In addition, the two plots reveal a pattern of large between-region and within-region variability, revealing an issue of overdispersion within regions over the years.

3 Methodology

3.1 A joint model for HCV and HIV infection

The data consist of multivariate repeated binomial measurements collected in a period of 10 years, from 1998 to 2007. Let Y _θ = ( Y _θ₁, Y _θ₂) be the response vector representing the number of reported cases of HCV and HIV infection, respectively, in region θ. In turn, let Y _ik = (Y_i_1k,…, Y_iJk) be the response vector containing the repeated measurements per infection k within the θth region in year j. Del Fava et al. (2011) discussed a joint Bayesian hierarchical model for HCV and HIV infection prevalence that merely took into account the regional clustering of data. In the first stage of the hierarchical model, they assumed that the distribution of Y_ijk is binomial, with sample size equal to n_ijk: Y_ijk(3.1)Binπ_ijkn_ijkjk

Figure 1

Source: Authors’ own.

The primary interest is in the estimation of π_ij₁ and π_ij₂, which are the prevalence of HCV and HIV infections in the θth region in year j, respectively, and in the association between the two infections at the population level, ρ(π_ij₁, π_ij₂).

In the second stage of the hierarchical model, for infection k, they specified a logistic model for the prevalence π_ij. It contains unstructured fixed effects adjusting for time, α and β_j, and region-specific random intercepts, γ _θ = (γ_i₁, γ_i₂):

{\begin{array}{l} logit (π_{i j 1}) = α_{1} + β_{j 1} + γ_{i 1}, \\ logit (π_{i j 2}) = α_{2} + β_{j 2} + γ_{i 2} . \end{array}

(3.2)

The time-specific coefficients, α_j and β_j, are assumed to be independently normally distributed with mean zero and variance depending on the hyperparameter $σ_{β}^{2}$ :

α_{j}, β_{j} \sim N (0, σ_{β}^{2}) .

(3.3)

We then specify a uniform distribution for the variance hyperparameter σ_b: σ_b(3.4)U

Gelman (2006) argued that this choice of prior on σ is preferable over other possible priors, e.g., the inverse-Gamma(ε, ε) or the half-Cauchy, because it is, in practice, less ‘informative’. The same prior will be therefore used for all the univariate variance hyperparameters in the other models.

The random intercepts γ_ik, which are assumed to follow a bivariate normal distribution, account for the association between the repeated measurements within region, independently of time. Depending on the structure of their covariance matrix D _γ , the random intercepts might even account for the association between the two infections, ρ(π_ij₁, π_ij₂). Del Fava et al. (2011) fitted a series of joint models with different variance structures and concluded that the model with the best goodness-of-fit is the so-called correlated random-intercept model. For this model, the region-specific random intercepts are assumed to follow a bivariate normal distribution with a mean vector of zeros and unstructured covariance matrix D _γ :

D_{γ} = (\begin{matrix} σ_{γ_{1}}^{2} & ρ_{γ_{1} γ_{2}} σ_{γ_{1}} σ_{γ_{2}} \\ ρ_{γ_{1} γ_{2}} σ_{γ_{1}} σ_{γ_{2}} & σ_{γ_{2}}^{2} \end{matrix}) .

(3.5)

The infection-specific variances of the random intercepts, $σ_{γ_{1}}^{2}$ and $σ_{γ_{2}}^{2}$ , account for the within-region association between the repeated measurements, whereas the parameter ρ_γ_1c2 is the correlation coefficient between the two infections at the level of the linear predictor, being therefore a measure of the association between HCV and HIV infection, ρ(π_ij₁, π_ij₂). The variance–covariance matrix D _γ is given a standard inverse-Wishart prior distribution, a distribution that generates random positive definite matrices, with degrees of freedom υ = 2 and scale matrix ψ equal to a two-by-two identity matrix, as it is usually done in the literature, i.e.,

D_{γ} \sim W^{- 1} (Ψ, 2), with Ψ = (\begin{matrix} 1 & 0 \\ 0 & 1 \end{matrix}) .

(3.6)

Prior distributions for the other unknown parameters can be found in Del Fava et al. (2011), where this basic model is introduced.

3.2 Joint models with additive overdispersion parameters

In this section, we propose extensions to the random-intercept joint model (3.2) to deal simultaneously with clustering and overdispersion. This is achieved by including separate sets of random effects for overdispersion and for clustering. We opt for a set of random effects θ_ijk, which, for convenience, are assumed to be independent of the random intercepts γ_ik. In this section, we focus on additive overdispersion parameters θ_ijk (McLachlan, 1997), which are introduced on the same scale of the linear predictor. We consider four possible situations of interest for the overdispersion random effects: (i) they are shared by the two infections; (ii) they are differentiated by infection and independent; (iii) they are differentiated by infection and allowed to be dependent; (iv) they are differentiated by infection and by year, leading thus to time-specific covariance matrices, D _θj .

3.2.1 Shared overdispersion parameters

The first model we consider is the shared overdispersion model, where the overdispersion parameters are shared between the infections, i.e., θ_ij₁ = θ_ij₂ = θ_ij:

{\begin{array}{l} logit (π_{i j 1}) = α_{1} + β_{j 1} + γ_{i 1} + θ_{i j}, \\ logit (π_{i j 2}) = α_{2} + β_{j 2} + γ_{i 2} + δ θ_{i j} . \end{array}

(3.7)

We assume that θ_ij has a normal prior distribution, $θ_{i j} \sim N (0, σ_{θ}^{2})$ , and the hyperparameter σ_i is distributed according to U(0,100). The underlying assumption behind the shared overdispersion model (3.7) is that the correlation between the infections described by the overdispersion parameters is equal to one. We use the parameter δ to relax the assumption of common variance between the overdis-persion parameters of HCV and HIV infections, since $σ_{θ_{H I V}}^{2} = δ^{2} σ_{θ_{H C V}}^{2}$ . Note that the case with $σ_{θ}^{2} = 0$ implies the absence of region-specific evolution patterns during the years.

3.2.2 Independent overdispersion parameters

Another possible model is the independent overdispersion model, which, unlike model (3.7), includes infection-specific random effects, θ_ijk:

{\begin{array}{l} logit (π_{i j 1}) = α_{1} + β_{j 1} + γ_{i 1} + θ_{i j 1}, \\ logit (π_{i j 2}) = α_{2} + β_{j 2} + γ_{i 2} + θ_{i j 2} . \end{array}

(3.8)

For θ_ijk we now assume a bivariate normal prior distribution, with covariance between θ_ij₁ and θ_ij₂ equal to zero:

(\begin{matrix} θ_{i j 1} \\ θ_{i j 2} \end{matrix}) \sim M V N [(\begin{matrix} 0 \\ 0 \end{matrix}), D_{θ} = (\begin{matrix} σ_{θ_{1}}^{2} & 0 \\ 0 & σ_{θ_{2}}^{2} \end{matrix})] .

(3.9)

Regarding the variance components for the overdispersion parameters θ_ijk, we assume that $σ_{θ_{1}}$ and $σ_{θ_{2}}$ are independently distributed according to a uniform prior distribution within the range [0,100].

This model assumes that, although there is overdispersion in the time evolution of prevalence among the regions, all correlation between HCV and HIV infections at the regional level is captured fully by the random intercepts, not by the overdispersion parameters.

3.2.3 Correlated overdispersion parameters

As a further extension, model (3.8) can be expressed as a correlated overdispersion model. For this model, the joint distribution of θ_ij₁ and θ_ij₂ can be rewritten as

(\begin{matrix} θ_{i j 1} \\ θ_{i j 2} \end{matrix}) \sim M V N [(\begin{matrix} 0 \\ 0 \end{matrix}), D_{θ} = (\begin{matrix} σ_{θ_{1}}^{2} & ρ_{θ_{1} θ_{2}} σ_{θ_{1}} σ_{θ_{2}} \\ ρ_{θ_{1} θ_{2}} σ_{θ_{1}} σ_{θ_{2}} & σ_{θ_{2}}^{2} \end{matrix})] .

(3.10)

For the covariance matrix D _θ, similarly to what we did for D _γ , we specify an inverse-Wishart prior distribution, $D_{θ} \sim W^{- 1} (Ψ, 2)$ , where $Ψ$ is a two-by-two identity matrix.

Assuming an unstructured covariance matrix for D _θ, it is possible to estimate an additional correlation between the infections, ρ_θ_1θ2. This implies that, after having accounted for the correlation between the infections using the region-specific random intercepts γ_ik, there is still correlation between HCV and HIV infections in the time evolution of prevalence that is captured by the overdispersion parameters.

3.2.4 Correlated overdispersion parameters with time-specific correlation

The last additive model that we consider extends model (3.10) relaxing the assumption of a constant correlation between HCV and HIV infections captured by the overdispersion parameters. We now let the covariance matrix D _θ change across years:

(\begin{matrix} θ_{i j 1} \\ θ_{i j 2} \end{matrix}) \sim M V N [(\begin{matrix} 0 \\ 0 \end{matrix}), D_{θ j} = (\begin{matrix} σ_{θ_{1} | j}^{2} & ρ_{θ_{1} θ_{2} | j} σ_{θ_{1} | j} σ_{θ_{2} | j} \\ ρ_{θ_{1} θ_{2} | j} σ_{θ_{1} | j} σ_{θ_{2} | j} & σ_{θ_{2} | j}^{2} \end{matrix})] .

(3.11)

For the covariance matrix D _θj , we specify a time-specific inverse-Wishart prior distribution, $D_{θ j} \sim W_{j}^{- 1} (Ψ, 2)$ , where $Ψ$ is a two-by-two identity matrix and j = 1, …, 10.

3.3 Joint models with multiplicative overdispersion parameters

We consider a setting in which we account for overdispersion using multiplicative effects (McLachlan, 1997; Molenberghs et al., 2010). While the random intercepts γ_ik induce association between the clustered measurements, the parameters θ_ijk account for additional overdispersion. Hence, in this section, we assume that π_ijk = θ_ijkκ_ijk and

{\begin{array}{l} Y_{i j k} \sim B i n (θ_{i j k} κ_{i j k}, n_{i j k}), \\ logit (κ_{i j 1}) = α_{1} + β_{j 1} + γ_{i 1}, \\ logit (κ_{i j 2}) = α_{2} + β_{j 2} + γ_{i 2} . \end{array}

(3.12)

Note that 0 ≤ θ_ijk ≤ 1 must hold to ensure that 0 ≤ θ_ijkκ_ijk ≤ 1.

We specify a beta prior distribution for θ_ijk. As a special case, we use a beta distribution with parameters equal to 1, equivalent to a uniform distribution over the range (0,1), which implies a flat prior for the overdispersion parameters:

\begin{array}{l} θ_{i j 1} \sim B e (1, 1), \\ θ_{i j 2} \sim B e (1, 1) . \end{array}

(3.13)

However, in general, we can assume that the parameters of the beta prior distributions are hyperparameters to estimate: θ_ijk(3.14)Beabk

Then, the infection-specific variance of the overdispersion random effects can be calculated as

σ^{2} = \frac{a b}{{(a + b)}^{2} (a + b + 1)} .

(3.15)

For the prior distribution of the hyperparameters a and b in (3.15), we choose independent flat uniform distributions within the range [0,100], i.e., a, b ~ U(0, 100).

4 Posterior computation, model selection and model checking

4.1 Posterior computation

The hierarchical Bayesian models presented in the previous section are fitted to data using MCMC methods, specifically Gibbs sampling implemented through the JAGS software (Plummer, 2003). For each model, we used three chains of 250 000 iterations each, burn-in of 125000 and thinning of 125.

Initial values for each of the three chains were generated randomly according to prespecified distribution. For instance, initial values for the coefficients α_j and β_j, as well as for the random effects γ_ik and θ_ijk, were drawn from independent standardized normal distributions; initial values for the variance hyperparameter σ_b were drawn from a uniform distribution within the range [0,1]; finally, initial values for the covariance matrices D _γ and D _θ were drawn from inverse-Wishart distributions with degrees of freedom υ = 2 and a two-by-two identity scale matrix.

For convergence assessment, we used diagnostic tools contained in the R package CODA (Plummer et al., 2006), namely, the potential scale reduction factor $\hat{R}$ (Gelman and Rubin, 1992) and Geweke’s Z-diagnostic (Geweke, 1992). We monitored the main parameters of each model, which are the time-specific coefficients α_j and β_j, the random effects γ_ik and θ_ijk, the covariance matrices D _γ and D _θ_{, the correlation coefficients ρ_γ and ρ_θ, and further scaling parameters, such as δ for model (3.7). Traceplots of the main parameters were employed in order to check the mixing properties of the three chains. Posterior means with credible intervals as well as convergence statistics (the upper bound of the CI for the potential scale reduction factor and Geweke’s Z-score) for all the variance components are presented in Section 1 of the supplementary material web appendix.}

4.2 Model selection

A hierarchical mixed-effects model may be seen as a missing-data problem, where the random effects are regarded as the missing information. When dealing with missing data problems, Celeux et al. (2006) showed that the deviance information criterion (DIC, Spiegelhalter et al., 2002), which is the typical selection criterion used for Bayesian models, does not work properly with distributions outside the exponential family. Moreover, Plummer (2008) suggested that the approximation used routinely to compute the DIC is valid only when the effective number of parameters (the penalty pD) is much smaller than the number of independent observations μ.

In their work, Celeux et al. (2006) suggest some alternative versions of the original DIC, which vary among them depending on the focus of the analysis. The most straightforward among these alternative is a ‘conditional’ DIC, denoted as DIC₇ by the authors, which considers the random effects, γ and θ , as additional parameters of the model, so we do not average them out (Neelon et al., 2010). Practically, to compute DIC₇, we use the following formula (Celeux et al., 2006):

D I C_{7} = 2 \bar{D (ζ)} - D (\bar{ζ}),

(4.1)

where ζ = ( γ , α , β , θ ). The first term is easily obtained from the model, since the deviance can be monitored similarly to other parameters; the second term can be calculated by plugging the posterior means of the parameters in the formula for the deviance.

Moreover, keeping in mind the specific concerns raised against the DIC, we decided to use two further criteria to select the best model: the penalized expected deviance (PED, Plummer, 2008), and the difference in posterior deviance (Aitkin et al., 2009; Aitkin, 2010).

The PED can be considered as a loss function when predicting the data Y using the same data Y. The issue of using the data twice (for estimation as well as for prediction) makes the expected deviance $\bar{D (θ)}$ too optimistic (Plummer, 2008). Thus, the PED penalizes it with a measure for model complexity, p_opt:

P E D = \bar{D (θ)} + p_{o p t} .

(4.2)

Even though outside the exponential family it is hard to calculate the optimism parameter p_opt; it can be estimated for general models using MCMC methods. According to Plummer (2011), the software JAGS uses importance sampling to estimate the parameter p_opt. However, the author warns that the estimates may turn out to be numerically unstable when the effective number of parameters is high, as it typically occurs with random-effects models. Similar to the DIC, the smaller the PED, the better.

The difference in posterior deviances (Aitkin et al., 2009; Aitkin, 2010) permits to compare pairs of models to select the best one. This approach is based on the observation that models with growing numbers of parameters are automatically penalized by the increasing diffuseness of the posterior distributions in their parameters. Thus, we can base the selection between two models on the difference between the whole posterior distributions of their deviances,

{D_{1, 2}^{(m)} = D_{1}^{(m)} - D_{2}^{(m)} : m = 1, . . ., M},

(4.3)

where M is the length of the MCMC chain. We can derive the posterior probability that Model 1 is better than Model 2,

P (D_{1, 2}^{(m)} < 0) = \frac{1}{M} \sum_{m = 1}^{M} I (D_{1, 2}^{(m)} < K),

(4.4)

where the constant K is assumed to be equal to -2log(9) = 4.39, a value that is calibrated to correspond to a likelihood ratio test favouring Model 1 with a posterior probability of 0.9 (Aitkin, 2010). It is also possible to derive 95% credible intervals (CI) for the difference in deviances: a totally negative 95% CI implies that we favour Model 1 over Model 2.

4.3 Model checking

Once we have selected the best models, it is necessary to assess the goodness-of-fit of these chosen models. Indeed, it might be that the best fitting model still fits the data poorly. Here we follow the approach of conducting posterior predictive checks, by comparing the observed data with the replicated data obtained from the posterior predictive distribution (Gelman et al., 1996; Neelon et al., 2010). The idea behind this approach is that the distance between the observed data Y and the replicated data Y *, accounted by some chosen discrepancy measures, should be as small as possible. This discrepancy measure, T = T( Y , θ ), should take an extreme value if the model is in contrast with the observed data. From these discrepancy measures, it is possible to compute Bayesian p-values (P_B) that denote the probability that the discrepancy measure based on the replicated data, T* = T( Y *, θ ), is more extreme than the measure based on the observed data (Gelman et al., 1996). A P_B close to 0.5 as well as, more generally, lying between 0.2 and 0.8 implies adequate fit, whereas values close to 0 or 1 indicate poor fit (Neelon et al., 2010). These Bayesian p-values can be computed using MCMC by evaluating the proportion of draws for which it stands that T* > T. For our models, we specifically consider two statistics for T, namely, the sample skewness, c₁, and the coefficient of variation, γ_y = σ/μ.

In addition to the Bayesian p-value, we also compute the mean squared error (MSE) for the in-sample and the out-sample predictions. In particular, for the out-sample predictions, we remove the observations for the last two years (2006 and 2007) and then we predict them using the model based on the remaining observations. The MSEs for both types of predictions are given by the following formulae:

M S E_{i n} = \frac{\sum_{i = 1}^{20} \sum_{j = 1}^{10} \sum_{k = 1}^{2} {(Y_{i j k}^{*} - Y_{i j k})}^{2}}{\sum_{i = 1}^{20} \sum_{j = 1}^{10} \sum_{k = 1}^{2} n_{i j k}}

(4.5)

and

M S E_{o u t} = \frac{\sum_{i = 1}^{20} \sum_{j = 9}^{10} \sum_{k = 1}^{2} {(Y_{i j k}^{*} - Y_{i j k})}^{2}}{\sum_{i = 1}^{20} \sum_{j = 9}^{10} \sum_{k = 1}^{2} n_{i j k}} .

(4.6)

Of course, the smaller the MSE, the closer the replicated data to the observed data and thus the better the goodness-of-fit of the model.

5 Application to the data

Table 1 presents the information criteria for model selection, discussed in Section 4.2, while MSEs and Bayesian p-values, discussed in Section 4.3, are presented in Table 2. We refer to the supplementary material web appendix for the posterior means and respective 95% CI of the variance components of each of the models.

Table 1 shows a comparison of all the fitted models in terms of DIC₇, PED (the latter based on further 5000 iterations of the MCMC at convergence), as well as according to the difference between couples of models in their posterior deviances. We notice that the PED and the difference in posterior deviances lead us to favour different models. As expected, the worst model is the basic joint model (model (3.2)). Among the models with additive overdispersion parameters, the PED indicates that the shared random-effects model is the best one (model (3.7)), while the DIC₇ and the difference in posterior deviances favour the model with correlated overdispersion parameters and time-specific correlation coefficients between HCV and HIV infection (model (3.11)).

Table 1
Information criteria for model selection

Type
Model
DIC₇
PED
Diff.

Basic (3.2) 17745 10350 –

Additive (3.7) 6859 6033 $D_{3.7, 3.2} = \bar{D_{3.7}} - \bar{D_{3.2}}$
–5442(–5492, –5390)
P(D_3.7,3.2 < 4.39) = 1

(3.8) 4251 7025 $D_{3.8, 3.7} = \bar{D_{3.8}} - \bar{D_{3.7}}$
–1305(–1381, –1226)
P(D_3.8,3.7 < 4.39) = 1

(3.10) 4250 7033 $D_{3.10, 3.8} = \bar{D_{3.10}} - \bar{D_{3.8}}$
–0.58(–88, 85)
P(D_3.10,3.8 < 4.39) = 0.46

(3.11)
3923
7151
$D_{3.11, 3.8} = \bar{D_{3.11}} - \bar{D_{3.10}}$
–72(–155, 12)
P(D_3.11,3.8 < 4.39) = 0.95

Multiplicative (3.13) 4092 7726 $D_{3.13, 3.11} = \bar{D_{3.13}} - \bar{D_{3.11}}$
–7(–89, 73)
P(D_3.13,3.11 < 4.39) = 0.53

(3.14) 4321 6941 $D_{3.13, 3.14} = \bar{D_{3.13}} - \bar{D_{3.14}}$
–115(–208, 26)
P(D_3.13,3.14 < 4.39) = 0.99

Type	Model	DIC₇	PED	Diff.
Additive	(3.7)	6859	6033	$D_{3.7, 3.2} = \bar{D_{3.7}} - \bar{D_{3.2}}$ –5442(–5492, –5390) P(D_3.7,3.2 < 4.39) = 1
(3.8)	4251	7025	$D_{3.8, 3.7} = \bar{D_{3.8}} - \bar{D_{3.7}}$ –1305(–1381, –1226) P(D_3.8,3.7 < 4.39) = 1
(3.10)	4250	7033	$D_{3.10, 3.8} = \bar{D_{3.10}} - \bar{D_{3.8}}$ –0.58(–88, 85) P(D_3.10,3.8 < 4.39) = 0.46
(3.11)	3923	7151	$D_{3.11, 3.8} = \bar{D_{3.11}} - \bar{D_{3.10}}$ –72(–155, 12) P(D_3.11,3.8 < 4.39) = 0.95
Multiplicative	(3.13)	4092	7726	$D_{3.13, 3.11} = \bar{D_{3.13}} - \bar{D_{3.11}}$ –7(–89, 73) P(D_3.13,3.11 < 4.39) = 0.53
(3.14)	4321	6941	$D_{3.13, 3.14} = \bar{D_{3.13}} - \bar{D_{3.14}}$ –115(–208, 26) P(D_3.13,3.14 < 4.39) = 0.99

Source: Authors’ own.

Table 2

MSE for in-sample and out-sample predictions and Bayesian p-values for model checking

Type	Model	MSE in-sample	MSE out-sample	P_Bγ₁	P_BC_v
Basic	(3.2)	12426	1177	1	0.06
Additive	(3.7)	2646	1168	0.42	0.47
	(3.8)	1227	1170	0.49	0.53
	(3.10)	1224	1184	0.50	0.53
	(3.11)	1229	890	0.50	0.53
Multiplicative	(3.13)	1227	1002	0.49	0.57
Multiplicative	(3.14)	1257	1206	0.49	0.50

Source: Authors’ own

Therefore, we discard both the models with independent overdispersion parameters (model (3.8)), which rules out any correlation between the overdispersion parameters between the two infections, as well as the correlated overdispersion mode, which assumes a constant correlation of 0.36 between the overdispersion parameters for HCV and HIV infection. As a matter of fact, from Figure 2 (panel d), which shows the posterior means of the time-dependent correlation coefficients with their respective 95% CI, we observe that the correlation is never significantly different from zero, except for 2004, when it is borderline positive, and for 2006, when it becomes instead significantly positive.

For the multiplicative models, according to the difference in posterior deviances and DIC₇, model (3.13) outperforms model (3.14), while the PED ranks the models the other way around.

Finally, when comparing the best additive model and the best multiplicative model, the PED favours the additive model with shared random effects, while, according to the difference in posterior deviance, we do not have enough confidence to choose between model (3.11) and model (3.13). The DIC₇ favours the additive time-specific correlated model (3.11).

If we consider the MSE estimates shown in Table 2, those for the in-sample predictions allow to label the basic model the additive shared model and, at a lesser extent, the multiplicative Be(a,b) model as poorly fitting models. More helpful, the MSEs for the out-sample predictions highlight the additive time-specific correlated model (3.11) as good predicting model. In this sense, the latter MSEs agree with the conclusions drawn by the DIC₇ and the difference in posterior deviances. Concerning the Bayesian p-values, the chosen discrepancy measures do not allow us to discriminate between the various additive and multiplicative models, since all the P_B are close to 0.50.

Figure 2

Source: Authors’ own.

Table 3 gives us more details about the estimated parameters, since it shows, for the basic model and the three selected overdispersion models, the so-called ‘time evolution odds ratios’ of year j, j = 2, …, 10, with reference year 1998. These odds ratios are computed by exponentiating the unstructured coefficients α_j (for HCV infection) and β_j (for HIV infection), that is to say, OR_HCV_|j,1998 = exp(α_j) and OR_HIV_|j,1998 = exp(β_j). The posterior means of the odds ratio are accompanied by the posterior means of the respective standard errors. These time evolution odds ratios compare the prevalence odds of each year compared to 1998. We notice that, generally, the prevalence odds are smaller than one, implying that the prevalence of HCV and HIV infection decreased over the years, mostly after 2004. The estimated odds ratios are more or less similar among all the models, but not their standard errors, which are larger for the additive overdispersed model and, to a greater extent, for the multiplicative one. This occurs because these models take into account the additional variability in the data, not considered by the basic model.

Table 3

Time evolution odds ratios for HCV and HIV infection, comparing prevalence odds of year j, j = 2, …, 10, with reference year 1998

Type	Model	Year	OR_HCV\|j,1998 (post. se)	OR_HCV\|j,1998 (post. se)
		1999	0.99 (0.011)	0.95 (0.015)
		2000	1.07 (0.012)	1.02 (0.015)
		2001	0.96 (0.011)	0.98 (0.014)
		2002	0.89 (0.0097)	0.98 (0.015)
Basic	3.2	2003	0.90 (0.0098)	0.94 (0.015)
		2004	0.86 (0.0093)	0.94 (0.015)
		2005	0.77 (0.0085)	0.91 (0.015)
		2006	0.80 (0.0088)	0.82 (0.013)
		2007	0.76 (0.0083)	0.81 (0.013)
		1999	0.98 (0.089)	0.92 (0.055)
		2000	1.02 (0.091)	0.95 (0.056)
		2001	0.88 (0.080)	0.89 (0.053)
		2002	0.94 (0.086)	0.97 (0.058)
Additive	3.7	2003	0.95 (0.084)	0.93 (0.054)
		2004	0.92 (0.083)	0.92 (0.054)
		2005	0.77 (0.071)	0.86 (0.053)
		2006	0.80 (0.073)	0.81 (0.048)
		2007	0.77 (0.070)	0.79 (0.048)
		1999	1.07 (0.099)	0.97 (0.089)
		2000	1.06 (0.12)	1.02 (0.087)
		2001	0.94 (0.090)	1.02 (0.087)
		2002	1.02 (0.11)	0.98 (0.085)
	3.11	2003	1.02 (0.099)	0.96 (0.083)
		2004	1.00 (0.12)	0.97 (0.091)
		2005	0.82 (0.084)	0.94 (0.089)
		2006	0.89 (0.090)	0.88 (0.085)
		2007	0.83 (0.084)	0.96 (0.13)
		1999	1.05 (0.37)	1.00 (0.14)
		2000	1.28 (0.45)	0.92 (0.11)
		2001	0.95 (0.31)	0.91 (0.12)
		2002	1.00 (0.38)	0.87 (0.12)
Multiplicative	3.13	2003	0.93 (0.33)	0.83 (0.10)
		2004	1.83 (0.96)	0.86 (0.11)
		2005	0.81 (0.29)	0.86 (0.12)
		2006	0.80 (0.28)	0.82 (0.10)
		2007	0.85 (0.32)	0.89 (0.16)

Figure 2 (panels a–c), for each fitted model, presents the posterior means of the variance components (with 95% CI) for HCV and HIV infection, respectively, and the correlation for the random effects γ_ik. For the models with additive random effects only (models (3.2)–(3.11)), we notice that misspecifying the overdis-persion random effects or even ignoring them affects neither the estimates of the variances components for the clustering random effects, nor the length of their credible intervals (panels a and b). This may be in relation to the fact the γ_ik and θ_ik are assumed to be independent and are both introduced on the same scale of the linear predictor; therefore, the former accommodate the within-region association all over the years, while the latter capture the unexplained additional variability. However, the same argument does not hold for the models with multiplicative random effects. For instance, for model (3.13) and, to a lesser extent, for model (3.14), we notice that ${\hat{σ}}_{γ_{1}}^{2}$ (but not ${\hat{σ}}_{γ_{2}}^{2}$ ) is larger with a wider CI. This is reflected in smaller values of the correlation $ρ_{γ_{1} γ_{2}}$ and longer CI, as can be observed in Figure 2 (panel c).

We refer to Figures 3–4 for a graphical representation of the results. We notice that model (3.2) in Figure 3 (panel a) shows parallel regional profiles, because only the cluster-specific random effects are specified; thus, the model fails to account for the exceeding variability in the data. This is instead accomplished by the best additive and multiplicative models (according to PED and the difference in posterior deviances) plotted in panels b–d. What is most striking from the graphical representation is that the fitted regional profiles from the pair of additive models look very similar, as also happens with the pair of multiplicative models. Finally, Figure 4 displays the marginal prevalence of HCV and HIV infection with the respective 95% CI for the basic model as well as the overdispersion models with additive and multiplicative random effects. For comparison, we plotted also the national prevalence per year, obtained by pooling together all the regional results. We notice that the five estimated prevalence profiles are fairly close to the national observed prevalence. What really distinguishes the marginal prevalence estimates from each other is their credible intervals, which account for all the variability shown by the regional profiles.

Figure 3

Source: Authors’ own.

6 Discussion

In this article, we discussed a set of hierarchical models that account simultaneously with clustering and overdispersion for binomial data. This objective has been achieved using two separate sets of random effects, one to induce association among the within-region measurements and between the infections, the other to account for the extra binomial variability in the data. For the latter set of random effects, two different settings have been considered: (1) the overdispersion parameters are included additively into the model, on the same scale of the linear predictor; (2) the overdispersion parameters are included multiplicatively in the model, correcting the prevalence in order to encompass the additional variability. All fitted combined models outperform the so-called basic joint model (Del Fava et al., 2011), which accommodates only the clustered nature of the data. This result was expected because the random-intercept model (3.2) cannot capture the complete association structure in the data.

Figure 4

Source: Authors’ own.

Under the additive approach, the overdispersion parameters can be seen as random slopes adjusting for clustering, as it happens with normal outcomes (Molenberghs et al., 2010), where overdispersion and cluster-specific random effects coincide. Indeed, the region-specific random intercepts γ_ik induce association averaging out the yearly measurements, while the overdispersion random parameters θ_ijk further adjust for the variation within each year. Under the multiplicative approach, the random effects deal differently with the clustering and the overdispersion. While the random intercepts γ_ik induce the within-region association and are responsible for shifting the regional profiles, the overdispersion random effects θ_ijk act directly on the prevalence π_ijk and adjust the standard mean-variance relation of the binomial distribution by inflating the variance of the estimated prevalence according to the time-specific variations.

Several models have been proposed for each setting according to the dependence structure among the overdispersion random effects and therefore to their covariance matrix D _θ. In both settings, we parameterize simultaneously the correlation between HCV and HIV infection over time (via the joint distribution of the random intercepts γ_ik) and within a specific time point (via the joint distribution of the overdispersion parameters θ_ijk).

A very critical point with the presented models is that model selection can be quite a difficult task to accomplish, considering that the small differences among the models, all related to the parameterization of the overdispersion parameters. The original DIC (Spiegelhalter et al., 2002) cannot be considered a valid option when the condition p_D << n does not hold (Plummer, 2008), because it tends to under-penalize more complex models. The other selection criteria here used, i.e., the DIC₇, the PED and the difference in posterior deviances, may lead to different conclusions. To our experience, it seems that PED, whose penalization is almost twice the penalty for DIC (Plummer, 2008), tends to favour less complex models, while the DIC₇ and the difference in posterior deviances select more complex models. In our case, the posterior means estimated by the different best models look similar, which is apparent from the graphical representation of prevalence. For a more detailed discussion of model selection in the Bayesian context, we refer to Lesaffre and Lawson (2012).

Furthermore, another critical issue of these Bayesian hierarchical models is clearly the modelling of the variance components. In the article, we used for the covariance matrices a prior distribution based on a standard inverse-Wishart distribution with degrees of freedom equal to k. However, recent literature highlights that the standard inverse-Wishart distribution is not exactly ‘non-informative’, since it shows a strong dependence between the variance and correlation: the higher the variance, the higher the correlation. Besides being dependent on each other, the variances and the correlation cannot be estimated independently. A possible solution is using a ‘scaled’ version of the inverse-Wishart distribution (O’Malley and Zaslavsky, 2008), where the parameterization keeps separate the correlation and the variances, so that they can be estimated independently from data. Future research is needed to incorporate such a prior in our modelling framework.

Finally, we note that, apart from temporal association and overdispersion, the design of our study, consisting of 20 Italian regions, suggests the possible presence of spatial correlation. Nonetheless, given the typical way of transmission of these infections (by injecting), we expect their transmission to occur more locally, in small injecting networks, and to be less influenced by transmission in neighbouring regions (in contrast to airborne infectious diseases, for which the spatial patterns matter). From a more general point of view, while such a source of association could well be weaker than the ones incorporated in the model, it should not be excluded a priori, even though it will require substantial additional work. It is a topic of future research.

Even though in this article we presented several extensions to the basic model of Del Fava et al. (2011), some enhancements are still possible. A possible extension for the model with unstructured time-dependent overdispersion parameters could be a parametric model for the correlation between HCV and HIV infections at the level of the overdispersion parameters over time. Furthermore, we could model the correlation between the overdispersion parameters for HCV and HIV using information collected at regional level, concerning risk factors related to injecting drug use (percentage of sharing syringes or other paraphernalia, etc.) or results of interventions (percentage of supplies of clean drug injection equipment, opioid substitution and other forms of drug dependence treatments, antiviral treatments for HIV, health promotion, etc.). An investigation of these type of models is a subject for an ongoing research and will be published in the near future.

Footnotes

Appendix

This appendix provides supplementary material for the analyses presented in Del Fava et al. (2013), giving more details about the output of the presented models. Section A1 gives details about the convergence assessment of the variance components for the models discussed in Section 3 in Del Fava et al. (2013). Section A2 presents individual prevalence and marginal prevalence proﬁles, estimated by the models discarded by the selection criteria. In Section A3 we discuss the implementation of the models in the JAGS software and in Section A4 the raw data are presented. An Excel ﬁle with the data is available upon request from the ﬁrst author at the following e-mail address: emanuele.delfava@unibocconi.it.

References

Agresti

(2002) Categorical data analysis. 2nd edition. New York: John Wiley & Sons.

Aitkin

Liu

Chadwick

(2009) Bayesian model comparison and model averaging for small-area estimation. Annals of Applied Statistics, 3, 199–221.

Aitkin

(2010) Statistical inference. An integrated Bayesian/likelihood approach. Boca Raton, FL: Chapman & Hall.

Booth

Casella

Friedl

Hobert

(2003) Negative binomial loglinear mixed models. Statistical Modelling, 3, 179–91.

Breslow

(1984) Extra-Poisson variation in log-linear models. Applied Statistics, 33, 38–44.

Breslow

Clayton

(1993) Approximate inference in generalized linear mixed models. Journal of the American Statistical Association, 88, 9–25.

Carey

Zeger

Diggle

(1993) Modelling multivariate binary data with alternating logistic regression. Biometrika, 80, 517–26.

Celeux

Forbes

Robert

Titterington

(2006) Deviance information criteria for missing data models. Bayesian Analysis, 1, 651–74.

Chen

Ahn

(1997) Marginal models with multiplicative variance components for overdispersed binomial data. Journal of Agricultural, Biological, and Environmental Statistics, 2, 440–50.

10.

Clayton

(1996) Generalized linear mixed models. In Gilks

Richardson

Spiegelhalter

(eds), Markov Chain Monte Carlo in practice. London: Chapman & Hall.

11.

Coutinho

(1998). HIV and hepatitis C among injecting drug users. BMJ (Clinical research ed.), 317(7156), 424–25.

12.

Del Fava

Kasim

Usman

Shkedy

Hens

Aerts

Bollaerts

Scalia Tomba

Vickerman

Sutton

Wiessing

Kretzschmar

(2011) Joint modeling of HCV and HIV infections among injecting drug users in italy using repeated cross-sectional prevalence data. Statistical Communications in Infectious Diseases, 3, 1–24.

13.

Gelman

Rubin

(1992) Inference from iterative simulation using multiple sequences. Statistical Science, 7, 457–511.

14.

Gelman

Meng

Stern

(1996) Posterior predictive assessment of model fitness via realized discrepancies. Statistica Sinica, 6, 733–807.

15.

Gelman

(2006) Prior distributions for variance parameters in hierarchical models. Bayesian Analysis, 1, 515–33.

16.

Geweke

(1992) Evaluating the accuracy of sampling-based approaches to calculating posterior moments. In Bernado

Berger

Dawid

Smith

AFM

(eds), Bayesian Statistics 4. Oxford, UK: Clarendon Press.

17.

Hinde

Demétrio

CGB

(1998) Overdis-persion: models and estimation. Computational Statistics and Data Analysis, 27, 151–70.

18.

Lesaffre

Lawson

(2012) Bayesian biostatistics. Hoboken, NJ: John Wiley & Sons.

19.

McLachlan

(1997) On the EM algorithm for overdispersed count data. Statistical Methods in Medical Research, 6, 76–98.

20.

Molenberghs

Verbeke

(2005) Models for discrete longitudinal data. Berlin: Springer-Verlag.

21.

Molenberghs

Verbeke

Demétrio

CGB

(2007) An extended random-effects approach to modeling repeated, overdispersed count data. Lifetime Data Analysis, 13, 513–31.

22.

Molenberghs

Verbeke

Demétrio

CGB

Vieira

AMC

(2010) A family of generalized linear models for repeated measures with normal and conjugate random effects. Statistical Science, 25, 325–47.

23.

Neelon

O’Malley

Normand

(2010) A Bayesian model for repeated measures zero-inflated count data with application to outpatient psychiatric service use. Statistical Modelling, 10, 421–39.

24.

O’Malley

Zaslavsky

(2008) Domain-level covariance analysis for multilevel survey data with structured nonresponse. Journal of the American Statistical Association, 103, 1405–18.

25.

Plummer

(2003) JAGS: a program for analysis of Bayesian graphical models using Gibbs sampling. Proceedings of the 3rd International Workshop on Distributed Statistical Computing (DSC 2003), March 20–22, Vienna, Austria. ISSN 1609–395X.

26.

Plummer

(2011) JAGS Version 3.3. 0 user manual. International Agency for Research on Cancer.

27.

Plummer

Best

Cowles

Vines

(2006) CODA: convergence diagnosis and output analysis for MCMC, R News, 6, 7–11.

28.

Plummer

(2008) Penalized loss functions for Bayesian model comparison. Biostatistics, 9, 523–39.

29.

Smith

B.J.

(2007) Boa: An R package for MCMC output convergence assessment and posterior inference. Journal of Statistical Software, 21, 1–37.

30.

Qaqish

Liang

K-Y

(1992) Marginal models for correlated binary responses with multiple classes and multiple levels of nesting. Biometrics, 48, 939–50.

31.

Skellam

(1948) A probability distribution derived from the binomial distribution by regarding the probability of success as variable between the sets of trials. Journal of the Royal Statistical Society, Series B, 10, 257–61.

32.

Spiegelhalter

Best

Carlin

Van der Linde

(2002) Bayesian measures of model complexity and fit. Journal of the Royal Statistical Society, series B, 64, 583–640.