In many bioscience studies, there is interest in simultaneously investigating two or more measurements within a single subject. In particular, it is often considered that two or more associated outcomes are affected at the same time depending on the level of certain treatments. For instance, developmental toxicity studies examine how a variety of variables are affected by pregnant mice exposed to several different levels of substances, such as rat weight, liver weight, gravid uterine weight, number of dead fetuses, presence of malformation, and fetal weight. Such variables should be jointly analyzed to take into account association between them because separate modeling can yield biased inferences (Catalano and Ryan, 1992; Regan and Catalano, 1999; Dunson, 2000; McCulloch, 2008; Hwang and Pennell, 2014; Hwang and Pennell, 2018). In the literature, several models have used latent variables to jointly analyze different types of outcomes such as continuous, binary, ordinal, unordered categorical or count variables. Catalano and Ryan (1992) proposed a joint model for binary and continuous outcomes assuming that a latent continuous variable underlying the binary outcome and continuous outcome have a bivariate normal random effects. Dunson (2000) proposed a general framework for modeling clustered mixed outcomes using a mixture of generalized linear models to describe the joint distribution of a set of underlying variables, and Dunson and Herring (2005) provided a general underlying Poisson variable model for mixed discrete outcomes. They followed a Bayesian approach to inference of the parameters and latent variables via a Markov chain Monte Carlo (MCMC) sampling algorithm. Nonparametric Bayesian methods also have been proposed for joint modeling of dicrete and continuous outcomes to describe irregular changes in distributions of the response variables (Hwang and Pennell, 2014; Fronczyk and Kottas, 2017; Hwang and Pennell, 2018).

In modeling count data in any research fields, it is common to encounter a large number of zero values. For example, developmental toxicology studies may be interested in how many birth defects the pregnant dams exposed to a toxic substance have. Then, many dams may not have any defects if the toxic level is not that high. Commonly used regression models for count data using Poisson or negative binomial distributions may not fit these count outcomes well, resulting in underestimating the zero-count probability. Such count data with excess zeros should be analyzed using discrete mixture models such as zero-inflated Poisson (ZIP) model and zero-inflated negative binomial (ZINB) model. Lambert (1992) proposed the zero-inflated Poisson (ZIP) regression model to describe discrete data with many zeros. The ZIP model then has a Poisson-Bernoulli mixture structure, which assumes two different sources for zero values: structural zeros by those subjects whose count response will always be zero, and sampling zeros that would occur by chance generated by a Poisson distribution. In addition to ZIP models, excess zeros in count data can also be accommodated by zero-inflated negative binomial models (Hall, 2000) or the hurdle models (Mullahy, 1986). Many researchers have utilized a variety of models to analyze zero-inflated count data. Hall (2000), Yau and Lee (2001), and Min and Agresti (2005) have proposed various zero-inflated models with random effects. Liu and Powers (2007) applied a growth curve model to longitudinal count data with excess zeros, and Lee et al. (2011) proposed a marginalized model approach to analyze clustered zero-inflated count data. Further, Kassahun et al. (2015) presented a joint modeling framework for hierarchical continuous and zero-inflated overdispersed count data.

There have recently been many literatures on Bayesian alternatives to fitting zero-inflated models. Rodrigues (2003) presented a Bayesian ZIP model based on the augmented data, and Ghosh et al. (2006) have proposed a Bayesian approach of zero-inflated power series (ZIPS) regression models including ZIP model for cross-sectional data using MCMC with data augmentation method. A hierarchical Bayesian model of correlated zero-inflated count data was presented by Dagne (2004), and several Bayesian models for longitudinal zero-inflated count data were also proposed by Neelon et al. (2010). Neelon et al. (2010) described Bayesian approaches to fitting repeated measures zero-inflated count data using three models: Poisson hurdle model, ZIP model, and zero-altered Poisson (ZAP) model. More recently, Neelon (2019) proposed a Bayesian ZINB regression model based on Pólya-Gamma mixtures to apply the approach to a spatiotemporal data analysis.

Motivated by developmental toxicity studies, we present a Bayesian joint modeling framework for continuous and zero-inflated count data. We utilize the zero-inflated Poisson (ZIP) regression model proposed by Lambert (1992) for count data with excess zeros. To account for dependency structure between continuous and count outcomes, the bivariate normal distribution will be assigned to subject-specific random effect shared across outcomes, which is expected to improve efficiency. We use the MCMC procedure with data augmentation and adaptive rejection sampling methods proposed by Ghosh et al. (2006) for parameter estimations. The proposed joint model will be applied to developmental toxicity data to estimate a benchmark dose (BMD) in a risk assessment.

The rest of the paper is organized as follows. In Section 2, we present a Bayesian joint model for continuous and zero-inflated count data, and describe the Bayesian inference procedures as well as model comparison criterion using deviance information criterion (DIC). In Section 3, we apply the proposed model to data from a developmental toxicity study of diethylhexyl phthalate (DEHP) in timed-pregnant CD-1 mice conducted by the National Toxicology Program. We conclude with a discussion of our results and future work in Section 4.

2.1. Zero-inflated Poisson regression model

Let Y = (Y₁, …, Y_n) be a vector of responses which is count data with an excess number of zeros. Lambert (1992) proposed the zero-inflated Poisson (ZIP) regression model, in which the count data was treated as a mixture of inflated zeros and Poisson counts. Let Z be an unobserved latent variable underlying the zero-inflated indicator outcome: for i = 1, 2, …, n,

Thus, the mixture of inflated zeros and Poisson counts can be written as

We can typically incorporate covariates into the model using the log link function for the expected mean count λ_i and the logit link function for p_i, which provide

where x_1i denotes a vector of explanatory variables and β₁ the vector of regression coefficients for the Poisson count stage, and x_2i denotes a vector of explanatory variables and β₂ the vector of regression coefficients for the Bernoulli zero-inflation stage. It follows that

The resulting distribution gives

We note that zero counts arise from both the Poisson count stage with parameter λ and the zero-inflation stage with parameter p. It is easy to show that

which implies that the ZIP regression model accommodates for overdispersion.

2.2. Joint model for continuous and zero-inflated count data

Let x_i be a predictor variable and Y = (Y_i1, Y_i2)^T be a 2 × 1 response vector measured on subject i (i = 1, …, n), where Y_i1 denotes the continuous response and Y_i2 denotes the count data with an excess number of zeros. For instance, in a developmental toxicity study, x_i may be the dose level of a toxic substance assigned to i^th litter, Y_i1 may be the average fetal body weight of i^th litter, and Y_i2 may be the number of malformed pups that i^th litter has. We know that the body weight and the number of malformations of the litters may be correlated and thus joint modeling can help to achieve better analysis results.

We jointly model the continuous response and the count response as follows, for i = 1, 2, …, n,

where

Note that λ_i and p_i are the parameters for Y_i2, and β_{0 j} and β_{1 j} are the fixed intercept and slope of x_i for j = 1, 2, 3. We accounted for excess zeros in Y_i2 using the logit model with the zero-inflation probability p_i. Further, we assume that ɛi~N(0,σɛ2), and U_i1 and U_i2 represent a subject-specific random effect assuming bivariate normal distribution with mean and variance-covariance matrix, independent of ε_i, U_i = (U_i1, U_i2)^T ~ N(0, ∑), where Σ=[σ12ρσ1σ2ρσ1σ2σ22]. The subject-specific random effect (U_i1, U_i2)^T models the correlation across the two outcomes, resulting in more flexible modeling and more efficient estimates (McCulloch, 2008). The complete data likelihood function of this model can be written as

where λ_i and p_i were defined in (2.1).

2.3. Bayesian inference for joint model

2.4. Bench mark dose lower limit (BMD/BMDL) estimation

We apply the proposed model into dose-response analysis in toxicology studies to estimate the benchmark dose (BMD), defined as the dose associated with a predefined small change in the response, called the benchmark response (BMR), from the background response level (Crump, 1995). Different definitions of the BMD exist depending on different types of response data, binary, continuous and count response data, and different types of changes in the response (Jensen et al., 2019; Jensen et al., 2020). Both the US Environmental Protection Agency (EPA) and the European Food Safety Authority (EFSA) recommend the definition of the BMD using relative response for count or continuous response data as follows,

where f (.) denotes the dose response model, and β is the vector of parameters in the model (US Environmental Protection Agency, 2012; Hardy et al., 2017). The US EPA and the EFSA also recommend choosing the BMR value of 5%. In practice, the BMD lower limit (BMDL), defined by the lower limit of a one-sided confidence interval of the BMD estimate, is used as reference values or limit values in regulatory applications.

2.5. Model comparison

We also fit two other models to compare the performance of our proposed model. An independent ZIP model accounts for continuous outcome and zero-inflated count outcome independently assuming ρ = 0 in the joint model in 2.1. A joint Poisson model jointly analyzes two outcomes, but does not account for excess zeros assuming β₀₃ = β₁₃ = 0 in the joint model in 2.1. For model comparison, the deviance information criterion (DIC) would be used (Spiegelhalter et al., 2002). The DIC is defined by DIC=D(θ)¯+pD, where D(θ) = −2 log f (y|θ) is the deviance, D(θ)¯ is the posterior mean of deviance, and pD=D(θ)¯-D(θ^) is the effective number of parameters in a model. Like the Akaike information criterion (AIC), the DIC accounts for both a measure of fit (deviance) and a measure of complexity (p_D) from a Bayesian perspective, and a model with a smaller DIC fits the data better. Celeux et al. (2006) proposed a modified DIC for models with mixtures of distributions and random effects as follows,

where y is the observed data, Z are the random effects and latent variables, and f (y, Z|θ) is the complete data likelihood function. We use the modified DIC for model comparison in our application.

3.1. Data description

We describe the proposed approach using data from a developmental toxicity study of diethylhexyl phthalate (DEHP) in timed-pregnant CD-1 mice conducted by the National Toxicology Program (NTP study: TER84064). DEHP is a substance mainly used as an additive in plastics to make them more flexible, and its widespread use in everyday and medical products has raised some concerns about its safety. DEHP breaks down slowly in water, but it can remain for a long time in soil and sediment, so it accumulates in organisms living in soil and sediment. In this study, pregnant CD-1 mice were exposed to DEHP in their feed at five different dose levels(0%, 0.025%, 0.05%, 0.10%, and 0.15% in feed). All live fetuses for each litter were examined for body weights and the occurrence of any types of malformations. Table 1 summarizes the average fetal body weight and the number of litters with no malformed pups per dose group. In the table, the average fetal body weight shows a trend of decreasing, and the number of litters with malformation increases (the count of litters with no malformed fetuses decreases) as the dose level increases. It is obvious that two litter-level outcomes are related each other so that they should be jointly modeled to consider the correlation between two outcomes within each litter. Moreover, the malformed pups shows a large incidence of zero counts, which requires zero-inflated count model, zero-inflated Poisson model.

3.2. Bayesian joint model

We use the proposed joint model to analyze the average fetal weight and the number of malformation within a litter. Let x_i be the assigned dose level of i^th litter, which has five different dose values(0%, 0.025%, 0.05%, 0.10%, and 0.15% in feed), Y_i1 be the average fetal weight (g) of i^th litter, and Y_i2 be the number of malformed fetuses within i^th litter. To implement Bayesian approach, we set the noninformative priors for β, σɛ2 and ∑: β_j ~ N(0, I₂), j = 1, 2, 3, σɛ2~IG(1,1), and ∑ ~ Inverse-Wishart(Ψ = I₂, ν = 4). Results were obtained by Markov chain Monte Carlo algorithm with data augmentation (Ghosh et al., 2006) described in Section 2.3.2. A total of 10,000 MCMC samples were generated from the conditional posterior distributions of parameters, and the posterior summary was based on 5,000 samples after a burn-in of 5,000 iterations. Convergence of MCMC algorithm was confirmed via traceplots and Gelman and Rubin’s potential scale reduction factor R̂ (Gelman et al., 2014). The MCMC algorithm was implemented using R 3.6.3 (R Core Team, 2020), and particularly the adaptive rejection sampling (ARS) method for sampling of β’s was conducted via arms function of armspp package in R.

3.3. Results

Posterior summaries for parameters in our proposed model (joint ZIP model) as well as the independent ZIP model and the joint Poisson model are provided in Table 2. According to the model comparison criteria DICs, we confirmed that the joint ZIP model had the smallest DIC (−1621.4) and thus fit the DEHP data comparatively better than the other two models. Table 2 also provides parameter estimates of the three models (posterior means and 95% credible intervals). The joint ZIP model shows significant effects of dose level on fetal weight and the number of malformations, respectively: β̂₁₁ = −0.941 indicates that fetal weight decreases as the DEHP level increases, and β̂₁₂ = 4.872 indicates that the number of malformed pups increases as the DEHP level increases.

We can examine zero-inflation from the Bernoulli model component in the joint model in 2.1. The significant positive coefficient of the intercept (β̂₀₃ = 1.415) indicated strong evidence of zero-inflation although the coefficient of the slope was not significant (β̂₁₃ = −0.244, CI=(−2.203, 1.718)). Considering zero-inflation aspect in the model greatly reduced DIC value from the joint Poisson model, meaning that the model fit was greatly improved. Parameter estimates for the random effects can explain the performance of jointly modeling the data. Obviously, fetal weight and number of malformed pups showed a negative association as evidenced by significant negative correlation of the random effects ( ρ̂ = −0.127). The joint ZIP model had better model fit (smaller DIC value) than the independent ZIP model separately modeling two outcomes, although they have very similar parameter estimates as a whole.

Figure 1 represents the estimated relative response curves along with 95% credible bands for fetal weight and number of malformed pups, respectlvely for the joint ZIP model as defined in Section 2.4. As dose level increases from 0% to 0.15%, the relative response curves linearly decreases for fetal weight and increases for number of malformations in a quadratic way. When we used the BMR value of 5% according to the recommendations of the US EPA and the EFSA, we obtained the estimated BMD and BMDL, 0.051% and 0.033%, respectively for fetal weight, and 0.011% and 0.008%, respectively for the number of malformations in the joint ZIP model. Table 3 gives the estimated BMDs and BMDLs for each model. The estimated BMD and BMDL of the joint ZIP model were comparable to or slightly smaller than the other two models. We know that the joint ZIP model can be used to determine the BMD for regulatory purposes, and suggests regulatory limits which are more protective (smaller BMDL) in a risk assessment.

We have proposed a Bayesian joint zero-inflated Poisson model for continuous and count data with excess zeros, in which the association between the two outcomes was taken into account via the subject-specific random effect. When we applied the proposed model to DEHP dataset from developmental toxicity study, the model showed better model fit than independent zero-inflated Poisson model and joint Poisson model based on DIC measures. We were able to confirm the negative effect of DEHP level on fetal weight and the positive effect of DEHP level on the number of malformations, respectively. It also turned out that the DEHP data should be analyzed using models dealing with excess zeros and jointly considering two outcomes. Further, our proposed model can be applied to a developmental toxicity study to estimate BMD and BMDL in dose-response analysis. We have confirmed the joint ZIP model suggests more protective regulatory limits in a risk assessment.

In terms of estimation, we presented a Bayesian approach to fitting the joint zero-inflated Poisson model based on the MCMC procedure with data augmentation method. The Bayesian approach has several advantages including use of prior information, easy interpretation, accounting for uncertainty in the parameter estimates and suitability for small sample data compared to maximum likelihood procedures. The data augmentation method used for count data with excess zeros greatly helped to improve efficiency of estimation compared to commonly used Metropolis-Hastings algorithm (Chip and Greenberg, 1995).

In the proposed model, ZIP model was used to deal with excess zeros in count outcome. However, the ZIP model can often be unstable if there is zero-deflated counts at a particular level of a categorical predictor, even if there are many zeros overall (Min and Agresti, 2005). The zero-inflated negative binomial (ZINB) model or the hurdle model may be considered to be used to describe zero-inflated count data instead of ZIP model to address unstableness of the model. In the proposed model (equation (2.1)), it was assumed that the fixed effects had a linear trend with dose level. The fixed effects could be generalized using a nonparametric method, e.g., regression splines. The random effects could also be extended in several ways. Nonparametric Bayesian methods such as Dirichlet process priors can be used for the random effect distribution instead of bivariate normal distribution, which will result in more flexible modeling (Hwang and Pennell, 2014). Another extension may include allowing random slopes in addition to random intercepts. As seen in Table 1, litter size is noticeably small at a high level of DEHP(0.1% and 0.15%) in Section 3. To see how litter size affects modeling results, we can take into account joint modeling framework of informative litter size and outcomes (Dunson et al., 2003; Hwang and Pennell, 2018).