A Bayesian approach of analysing semi-continuous longitudinal data with monotone missingness

Abstract

There is a rich literature on the analysis of longitudinal data with missing values. However, the analysis becomes complex for semi-continuous (zero-inflated) longitudinal response with missingness. In this article, we propose a partially varying coefficients regression model for analysing such data. We use a two-part model, where in the first part we propose a latent dynamic model for accounting a ‘zero’ or a ‘non-zero’ response, and in the second part we use another dynamic model for estimating the mean trajectories of non-zero responses. Two dynamic models are linked through subject-specific random effects. The missing covariates are imputed repeatedly based on their respective posterior predictive distributions and the missing responses are imputed using the working model under different identifying restrictions. We analyse data from the Health and Retirement Study (HRS) for aged individuals and develop a dynamic model for predicting out-of-pocket medical expenditures (OOPME) containing excess zeros. The operating characteristics of the proposed model are investigated through extensive simulation studies.

Keywords

Health insurance Legendre orthogonal polynomials Probit model semi-continuous longitudinal response varying coefficients model

1 Introduction

Analysis of longitudinal data is typically challenging due to inherent dependence among measurements from the same subject over different time points. Despite the difficulty, various methods have been proposed for modelling mean trajectory and covariance structure for univariate and multivariate longitudinal data. Random effects models proposed by Laird and Ware (1982) have been successfully used for analysing longitudinal response. Pourahmadi (1999, 2000), Pan and Mackenzie (2003) proposed a joint mean-covariance model for univariate longitudinal data. Das et al. (2013) generalized Pourahmadi's approach to bivariate longitudinal outcomes. Flexible semi-parametric Bayesian models have also been proposed for grouped longitudinal measurements by Das and Daniels (2014) and Das et al. (2015).

In many real applications, we have longitudinal measurements with excess zeros. Statistical methods based on mixture models have been used for analysing longitudinal measurements with ‘excess zeros’. Lambert (1992) proposed a zero-inflated Poisson (ZIP) model which can handle such excess zeros in count data. Hall (2000) proposed a zero-inflated binomial (ZIB) model for upper-bounded count data containing a significant proportion of zeros. Min and Agresti (2002) provided a nice tutorial on different methods for handling data with excess zeros. Most popularly used methods for analysing zero-inflated data include two-part models (Duan et al., 1983), Tobit models (Amemiya, 1984) and finite mixture models (Aitkin and Rubin 1985). Ghosh and Tu (2008), Boone et al. (2012), Buu et al. (2012), Alfo and Maruotti (2010) and Yang et al. (2016) developed models for zero-inflated longitudinal count data. Similar models have also been proposed for continuous longitudinal measurements with excess zeros (Olsen and Schafer, 2001; Mukherji et al., 2016, Farewell et al., 2017). In this article, we propose a Bayesian approach of analysing zero-inflated continuous longitudinal response with incomplete covariates and response. Specifically, we consider monotone incompleteness which is quite common in clinical trials and/or biomedical experiments. In longitudinal studies, incompleteness typically occurs due to dropouts. In the presence of dropouts, marginal distribution of the response variable is not identified in general (Little, 1993). Selection models (Diggle and Kenward, 1994) and pattern mixture models (PMM; Little, 1993, 1995) are commonly used for modelling incomplete longitudinal data. In a selection model all subjects contribute to the model at all time points, but in a PMM subjects only contribute to the model prior to the respective dropout times. Thus PMMs are usually under-identified (Little, 1993). Little (1995) proposed several identifying restrictions for overcoming the under-identification problem which are discussed with examples in Daniels and Hogan (2008). A PMM for handling zero-inflation and incompleteness in a longitudinal setting is proposed by Maruotti (2001). In our work, however, we do not model the missingness but impute missing covariates and missing response using a working model under different identifying restrictions (Molenberghs et al., 1998). We analyse data from the Health and Retirement Study (HRS) conducted by University of Michigan for an ageing population. Subjects were either retired or about to retire (with age 50 years or above), and they were followed for 16 years with a follow-up frequency of 2 years which we refer as waves. The goal of this study was to understand the challenges of ageing in terms of health condition and economy. According to a joint report by the US Department of State and the National Institute on Aging, nearly 500 million people in the world were 65 years or above in 2006 (Dobriansky et al., 2007). By 2030, this number is expected to increase to 1 billion. In the United States life expectancy has increased from 49 years to 78 years in the last century (Arias, 2006), and in China it has increased from 40 years to 76 years in the last 50 years. Thus, it is expected that the growth of the older population will have a global impact on the health care system, economy and insurance.

In the HRS data, variables related to physical health, financial health and demographies are measured from the subjects longitudinally. This dataset has been analysed by several authors from various perspective. Brown et al. (2015) proposed a Bayesian two-part model for household debt based on the HRS data. Bhuyan et al. (2018) proposed a two-stage regression model for handling endogeneity and incompleteness due to censoring for the HRS data. They used a Bayesian Tobit model for the response and the endogenous covariate. Das et al. (2017) proposed a dynamic hurdle model for predicting the rate of hospital admission of aged individuals. Our current work is different from the previous works (based on the HRS data) since we consider the problem of monotone missingness and zero-inflation in a longitudinal setting. Our variable of interest is out-of-pocket medical expenditure (OOPME). Note that at each wave, OOPME is the total cost of all medical services excluding the cost that is reimbursed or paid through health insurance. The OOPME is an accurate measure of the actual financial burden of medical expenditure incurred by individuals (Mukherji et al., 2016), and contains excess zeros. The objective of our study is to model and predict OOPME in terms of the covariates measuring physical and financial health. We use a two-part model for our analysis. In the first part, we consider a latent variable for the binary indicator of zero/non-zero response, and use a dynamic probit model similar to the model proposed in Albert and Chib (1993). In the second part, we propose a partially varying coefficients dynamic model for non-zero response. Missing covariates are imputed from respective posterior predictive distributions, and missing responses are imputed based on the working model under different identifying restrictions. Model parameters are estimated by Markov Chain Monte Carlo (MCMC) in a multiple imputation framework.

The rest of our article is organized as the following. In Section 2, we propose two dynamic models for predicting zero-inflated longitudinal response under monotone missingness. Results from the analysis of the HRS data are shown in Section 3. In Section 4, we show the results from the simulation study. Some concluding remarks are given in Section 5.

2 Proposed model and methods

We consider a non-negative zero-inflated longitudinal response measured from N subjects at T different time points. The set of covariates contains $J$ covariates with time-varying effects and additional $J^{'}$ covariates with fixed effects on the response. Note that all the N subjects are not observed at all the T time points, and the missingness is monotone in the sense that if a subject is missing at time T, then it is missing at time $t^{'}$ , for $t^{'} > t$ . We also note that when a subject is missing, we observe neither the response nor the time-varying covariates for that subject. Let $Y_{it}$ be the response from the ith subject at time T. We consider a two-part model where in the first part, we model the probability of zero response and in the second part, we consider a dynamic model for non-zero response.

2.1 Dynamic model for the probability of a zero response

Define binary random variables $E_{it}$ as the following:

E_{it} = \{\begin{matrix} 1, & for Y_{it} > 0; \\ 0, & for Y_{it} = 0 . \end{matrix}

We model $E_{it}$ as a Bernoulli random variable with success probability $= p_{it}$ . Thus, assume $E_{it} | p_{it} \sim$ Bernoulli $(p_{it})$ , where $p_{it} = P (E_{it} = 1)$ . We introduce a latent (continuous) random variable $R_{it}$ as the following:

E_{it} = \{\begin{matrix} 1, & for R_{it} > 0; \\ 0, & for R_{it} \leq 0, \end{matrix}

and model $R_{it}$ using the following partially varying coefficients linear mixed model:

R_{it} = \sum_{j = 1}^{J} α_{j} (t) x_{ij} (t) + \sum_{j^{'} = 1}^{J^{'}} β_{j^{'}} z_{{ij}^{'}} (t) + b_{i} + e_{it},

(2.1)

where $α_{j} (t)$ denotes time-varying effect of the covariate $x_{j}$ , $β_{j^{'}}$ denotes fixed effect of the covariate $z_{j^{'}}$ , and the residual errors $e_{it}$ are assumed to be identically and independently distributed as $N (0, 1)$ . The subject-specific random effects $b_{i}$ capture the correlation among the measurements from the same subject at different time points. We also assume that $b_{i}$ s and $e_{it}$ s are independently distributed. In (2.1) we essentially use the following probit model for $p_{it}$ :

Φ^{- 1} (p_{it}) = \sum_{j = 1}^{J} α_{j} (t) x_{ij} (t) + \sum_{j^{'} = 1}^{J^{'}} β_{j^{'}} z_{{ij}^{'}} (t) + b_{i},

(2.2)

where Φ denotes the cumulative distribution function of a standard normal random variable. One can model the time-varying coefficients $α_{j} (t)$ in Equation (2.1) either parametrically or non-parametrically. A non-parametric modelling is more appealing since the functional forms of $α_{j} (t)$ are not known in advance. We consider a non-parametric modelling based on Legendre orthogonal polynomials (LOP). These polynomials are well-behaved basis functions (Hušková and Sen, 1985; Meyer, 2000) and have been used in non-parametric regression (Cui et al., 2006; Bhuyan et al., 2018).

A Legendre polynomial of order K as a function of time T is expressed as the following:

P_{k} (t) = \sum_{l = 0}^{L} (- 1)^{l} \frac{(2 k - 2 l)!}{2^{k} l! (k - l)! (k - 2 l)!} t^{k - 2 l},

(2.3)

where $L = \frac{k}{2}$ or $\frac{k - 1}{2}$ , whichever is an integer. LOPs are defined over [−1, 1] and are orthogonal in the sense that $\int_{- 1}^{1} P_{r} (t) P_{s} (t) dt = 0$ , for $r \neq s$ . First few LOPs are given as the following:

$P_{0} (t) = 1, P_{1} (t) = t$ and $P_{2} (t) = \frac{1}{2} (3 t^{2} - 1) .$ The original time points are transformed to the LOP range [−1, 1] by defining $t^{'} = 1 + 2 (\frac{t - t_{\min}}{t_{\max} - t_{\min}})$ , where $t_{\min}$ and $t_{\max}$ are the smallest and the highest time points, respectively. Let $P^{(k)} (t) = [P_{0} (t), P_{1} (t), \dots, P_{k} (t)]^{T}$ denote the family of the first $k + 1$ basis functions, and we express $α_{j} (t)$ as a linear combination of the basis functions: $α_{j} (t) = a_{j}^{T} P^{(k)} (t)$ , where $a_{j} = (a_{j 0}, a_{j 1}, \dots, a_{jk})^{T}$ is the coefficient vector. The optimal order K is obtained from the information criteria, for example, AIC/BIC.

2.2 Dynamic model for non-zero longitudinal response

Next we model the non-zero responses for each subject at different time points. We note that in the original data different subjects are non-zero at different time points. Let the ith subject give non-zero responses at $T_{i}$ different time points where $T_{i} = \sum_{t = 1}^{T} E_{it}$ . Let $Y_{i} (t_{il})$ denote the non-zero response from the ith subject at the lth time point, where $l = 1, 2, \dots, T_{i}$ , and $i = 1, 2, \dots, N$ . We use the following partially varying coefficients model for $Y_{i} (t_{il})$ :

Y_{i} (t_{il}) = \sum_{j = 1}^{J} γ_{j} (t_{il}) x_{ij} (t_{il}) + \sum_{j^{'} = 1}^{J^{'}} δ_{j^{'}} z_{{ij}^{'}} (t_{il}) + c_{i} + ε_{i} (t_{il}),

(2.4)

where $γ_{j} (t_{il})$ denotes time-varying effect of the covariate $x_{j}$ , $δ_{j^{'}}$ denotes fixed effect of the covariate $z_{j^{'}}$ and the residual errors $ε_{i} (t_{il})$ are assumed to be identically and independently distributed as $N (0, σ_{ε}^{2})$ . The functions $γ_{j} (t_{il})$ are modelled non-parametrically using LOPs as discussed in Section 2.1.

In our application, a positive response at one time point can affect the probability of a zero response at other time points. Such dependence can be modelled by assuming a joint distribution for ( $b_{i}$ , $c_{i}$ ). Here, we assume that $(b_{i}, c_{i}$ ) follows a bivariate normal distribution with mean vector $= (0, 0)$ and unknown covariance matrix $= Σ$ . However, we assume that the latent random variables $R_{it}$ and the non-zero measurements $Y_{i} (t_{il})$ are conditionally independent given the random effects $b_{i}$ and $c_{i}$ . Two dynamic models are linked through the joint distribution of the random effects.

2.3 Likelihood and parameter estimation

Let $Θ_{1}$ and $Θ_{2}$ denote the set of all model parameters from the first dynamic model (given in Equation 2.1) and the second dynamic model (given in Equation 2.4), respectively. Also, let $f (b_{i}, c_{i})$ denote the bivariate normal density of $(b_{i}, c_{i})$ . Then, the joint marginal likelihood (after integrating out the random effects) can be expressed as the following:

\begin{matrix} L (Θ_{1}, Θ_{2} | R, Y) = \int \prod_{i = 1}^{N} [L_{1 i} (Θ_{1} | R_{i}, b_{i}) \times L_{2 i} (Θ_{2} | Y_{i}, c_{i})] \times f (b_{i}, c_{i}) d (b_{i}, c_{i}), \end{matrix}

(2.5)

where $L_{1 i} (Θ_{1} | R_{i}, b_{i}) = [\prod_{t = 1}^{T} \{1 (R_{it} > 0) 1 (E_{it} = 1) + 1 (R_{it} \leq 0) 1 (E_{it} = 0)\} \times ϕ (R_{it} | b_{i})]$ and $L_{2 i} (Θ_{2} | Y_{i}, c_{i}) = [\prod_{l = 1}^{T_{i}} ϕ (Y_{i} (t_{il}) | c_{i})]$ . Here $I (A)$ denotes an indicator variable which takes value 1 if A occurs, and 0 otherwise. Also $ϕ (R_{it} | b_{i})$ and $ϕ (Y_{i} (t_{il}) | c_{i})$ denote the normal density of $R_{it}$ conditional on $b_{i}$ (from Equation 2.1) and the normal density of $Y_{i} (t_{il})$ conditional on $c_{i}$ (from Equation 2.4), respectively.

Assuming that the prior distributions for $Θ_{1}$ and $Θ_{2}$ are independent, joint posterior density of $Θ = (Θ_{1}, Θ_{2})$ can be expressed as:

π (Θ | R, Y) \propto L (Θ_{1}, Θ_{2} | R, Y) \times π (Θ_{1}) \times π (Θ_{2}) .

(2.6)

We compute full conditional distributions of the model parameters from the aforementioned joint posterior density (details are given in the Appendix). Model parameters are estimated by their respective sample posterior means. Convergence of the Markov chains are monitored graphically and also by computing the multivariate potential scale reduction factors as suggested in Brooks and Gelman (1998).

2.4 Imputation of the missing values

Xie and Paik (1997), Schafer and Yucel (2002), Searman et al. (2012), Chen and Fu (2011) and Maruotti (2011) proposed methods for modelling or imputing the missing covariates. For a set of correlated variables, one can use a joint model and update the missing values based on the respective full conditional distributions. However, in a regression setting with a set of uncorrelated predictors, joint modelling is not meaningful. Also joint modelling becomes cumbersome for larger number of variables under consideration. Here we propose an alternative multiple imputation method where univariate models are used for imputing the predictors, and the missing responses are imputed based on the working model under a set of identifying restrictions. In real applications (Section 3), we first verify the absence of multicollinearity in our data by computing variance inflation factor (VIF) for each covariate at each time point. Additionally, we assume that the covariates and the response are ‘missing at random’ (Rubin, 1976).

In our proposed dynamic models, we have $J + J^{'}$ covariates. Consider a particular covariate, say $X_{1}$ , for which we assume a probability distribution with pmf/pdf $P (X_{1} | θ)$ , where $θ$ denotes the set of parameters for the probability distribution of $X_{1}$ . The missing values for the covariate $X_{1}$ are imputed based on the following posterior predictive distribution:

P (X_{1}^{(mis)} | X_{1}^{(obs)}) = \int P (X_{1}^{(mis)} | θ) P (θ | X_{1}^{(obs)}) d θ .

(2.7)

In other words, we first consider a suitable prior distribution for θ, and then compute the posterior distribution of θ for the given observed data ( $X_{1}^{(obs)}$ ). Then we sample from this posterior density, and for each sampled θ we impute the missing data ( $X_{1}^{(mis)}$ ) by directly sampling from $P (X_{1} | θ)$ .

Once the missing covariates are imputed, we use the dynamic models given in Equations (2.2) and (2.4) for imputing the missing response values. We start with some initial parameter values for these two dynamic models. Typically we fit two models on the subjects which are fully observed (no missing values) and compute the maximum likelihood estimates (mle) of the model parameters. These estimates are taken as the initial choice of the model parameters. For one set of imputed covariates, we impute the missing $E_{it}$ values based on the dynamic model given in Equation (2.2). Thus, we replace a missing $Y_{it}$ by 0 with probability $1 - {\hat{p}}_{it}$ . Next, the non-zero missing values (with probabilities ${\hat{p}}_{it}$ ) are imputed based on the dynamic model given in Equation (2.4) under certain identifying restrictions discussed in the following.

Considering the model in Equation (2.4) as our working model, let $f_{t}^{*} (Y_{s} | Y_{1}, \dots, Y_{s - 1})$ ; for $s = t + 1, \dots, T$ denotes the conditional distribution of $Y_{s}$ for the given measurements upto time $s - 1$ for a subject who is observed till time T. Here, we ignore the subscript i for the ease of presentation. For longitudinal measurements with monotone missingness, these conditional distributions are used for imputing the missing outcomes. Note that in longitudinal study, there is an inherent dependence among the outcomes from the same subject, and hence the imputed missing outcomes must be correlated to the observed outcomes from the same subject. By imputing the missing outcomes from these conditional distributions, we ensure such dependence. However, these conditional distributions are non-identified (Daniels and Hogan, 2008), and under this setting marginal distribution of response is also not identified (Little, 1993). We can identify these non-identified distributions by imposing some restrictions (Little, 1994; Molenberghs et al., 1998) as the following.

Under complete case missing value (CCMV) restriction, we have

f_{t}^{*} (Y_{s} | Y_{1}, \dots, Y_{s - 1}) = f_{T}^{*} (Y_{s} | Y_{1}, \dots, Y_{s - 1}),

(2.8)

for $s = t + 1, \dots, T$ . Hence under CCMV restriction, the conditional density of an unobserved data is equated to the corresponding conditional density for the subjects who completed the study (no missing value). We fit a T-variate distribution based on the subjects with no missing value and use the fitted distribution for imputing the missing outcomes at different time points. In practice, we may have only a few subjects with no missing value, and then imputation based on CCMV will be inefficient.

For available case missing value (ACMV) restriction, we have

f_{t}^{*} (Y_{s} | Y_{1}, \dots, Y_{s - 1}) = f_{\geq s}^{*} (Y_{s} | Y_{1}, \dots, Y_{s - 1}),

(2.9)

for $s = t + 1, \dots, T$ . Here the conditional density of an unobserved data (at time s) is equated to the corresponding conditional density of all the subjects who are observed at least till time s. Thus we have different distributions for imputing the missing outcomes at different time points.

The ACMV restriction in (2.9) can further be expressed as the following:

\begin{matrix} f_{\geq s}^{*} (Y_{s} | Y_{1}, \dots, Y_{s - 1}) & = & \sum_{k = s}^{T} \frac{ψ_{k} f_{k}^{*} (Y_{1}, \dots, Y_{s - 1})}{\sum_{k = s}^{T} ψ_{k} f_{k}^{*} (Y_{1}, \dots, Y_{s - 1})} f_{k}^{*} (Y_{s} | Y_{1}, \dots, Y_{s - 1}) \\ = & \sum_{k = s}^{T} ω_{sk} f_{k}^{*} (Y_{s} | Y_{1}, \dots, Y_{s - 1}), \end{matrix}

where $ψ_{k}$ denotes the proportion of individuals observed at time K and $ω_{sk} = \frac{ψ_{k} f_{k}^{*} (Y_{1}, \dots, Y_{s - 1})}{\sum_{k = s}^{T} ψ_{k} f_{k}^{*} (Y_{1}, \dots, Y_{s - 1})}$ . Thus the missing responses can be imputed from a mixture normal density under ACMV restriction.

Alternatively, under neighbouring case missing value (NCMV) restriction we have

f_{t}^{*} (Y_{s} | Y_{1}, \dots, Y_{s - 1}) = f_{s}^{*} (Y_{s} | Y_{1}, \dots, Y_{s - 1}),

(2.10)

for $s = t + 1, \dots, T$ . Here the nearest identified pattern is considered for imputing an unobserved response at a particular time point. This method is computationally easy compared to CCMV and ACMV.

Once the missing non-zero responses are imputed, we obtain one set of complete data $Z = (Y, E)$ corresponding to one set of imputed predictors. We fit the dynamic models in Equations (2.2) and (2.4) on the complete dataset and jointly estimate the model parameters (based on the joint likelihood in Equation [2.5]) using Gibbs sampler as discussed in Section 2.3. This process is repeated 1 000 times, that is, we consider 1 000 complete datasets (by imputing missing covariates and response 1 000 times) and obtain 1 000 estimates of the model parameters. The final estimates are obtained by averaging 1 000 estimated parameter values, and the corresponding Monte Carlo Standard Errors (MCSE) are also computed based on 1 000 replicates of the complete dataset.

3 HRS data analysis

We consider data from the HRS for our analysis. This study started in 1992 with 1 929 subjects of age 50 years or above, and ended in 2008 with a follow-up frequency of two years which we refer as waves. The goal of the study was to understand the challenges of ageing related to physical health condition and additional cost due to hospital admission. We note that for 1 227 subjects, we have the complete data (no missing observation) while for the others we have missingness at one or more waves but the missingness is monotone. At any particular wave, the medical expenditure of a subject is covered by exactly one health insurance.

For our analysis, we focus on modelling and predicting OOPME. Note that OOPME includes costs related to medical services (e.g., doctor's visit, hospital visits, drugs, etc.), but excludes the costs reimbursed through health insurance. In health economics, OOPME is treated as an important measure of financial risk (Mukherji et al., 2016). In our data, we have three different health insurances: (a) employment insurance, (b) government insurance and (c) other (private) insurance. Note that government insurance includes medicaid insurance for people with low income, and medicare insurance for people with age 65 years or above. Clearly, for an individual if there is no medical expenditure or if all the medical costs are paid by health insurance at a particular wave, then OOPME will be zero for that person at that wave. Thus, our response variable contains excess zeros, and hence we use the proposed two-part model for data analysis. In our data, we have 9–16% of zeros on average in OOPME across the waves. The set of covariates contains body mass index (BMI) as a continuous variable and the presence (yes/no) of eight chronic diseases: hypertension, diabetes, cancer, lung problem, heart problem, stroke, arthritis and psychological problems. We consider BMI, calendar age and three indicator variables for three different health insurances as the covariates with potential time-varying effects on the response. The effects of all the other covariates are considered as fixed.

3.1 Dynamic models and priors

Let $Y_{it}$ be the OOPME for ith subject at Tth wave, for $t = 1, \dots, 8$ . We define the random variables $E_{it}$ and $R_{it}$ as given in Section 2.1. We model the effects of calendar age, BMI and health insurances by Legendre orthogonal polynomial discussed in Section 2.1. We use the following linear mixed model for modelling the latent variable $R_{it}$ :

R_{it} = \sum_{j = 1}^{5} α_{j} (t) x_{ij} (t) + \sum_{j^{'} = 1}^{8} β_{j^{'}} z_{{ij}^{'}} (t) + b_{i} + e_{it},

(3.1)

for $t = 1, \dots, 8$ . Here, $b_{i}$ is the subject-specific random effects, and the residuals $e_{it}$ are assumed to follow normal distribution with mean = 0 and variance = 1.

For non-zero response, however, we consider log transformed data following Mukherji et al. (2016). Typically, medical expenditure data are highly positively skewed and hence a log transformation is useful from the modelling perspective. Thus, our model for the non-zero OOPME is the following:

\log (Y_{i} (t_{il})) = \sum_{j = 1}^{5} γ_{j} (t_{il}) x_{ij} (t_{il}) + \sum_{j^{'} = 1}^{8} δ_{j^{'}} z_{{ij}^{'}} (t_{il}) + c_{i} + ε_{i} (t_{il}),

(3.2)

and we assume that $b_{i}$ and $c_{i}$ jointly follow a bivariate normal distribution as discussed in Section 2.2. The residuals $ε_{i} (t_{il})$ follow $N (0, σ_{ε}^{2})$ .

Table 1:

Sensitivity analysis for the HRS data analysis

Parameter	Prior	Estimate
$a_{1}$	N(0,10)	1.64
–	N(0,100)	1.62
–	N(0,1000)	1.65
$a_{2}$	N(0,10)	0.84
–	N(0,100)	0.86
–	N(0,1000)	0.83
$β_{1}$	N(0,1)	0.065
–	N(0,100)	0.062
$β_{2}$	N(0,1)	0.25
–	N(0,100)	0.23
$δ_{1}$	N(0,1)	0.53
–	N(0,100)	0.54
$δ_{2}$	N(0,1)	0.33
–	N(0,100)	0.31
$σ_{ε}^{2}$	Inverse gamma (1,2.5)	1.88
–	Inverse gamma (2,3.5)	1.85
–	Inverse gamma (0.001,0.001)	1.86

We consider the following prior distributions for model parameters. For each coefficient vector $a_{j}$

( $j = 1, 2, \dots, 5)$ of the LOPs, we consider a multivariate normal prior with mean vector $= 0$ and covariance matrix $Σ = 1 000 I_{k + 1}$ , where $I_{k + 1}$ denotes the identity matrix of order $k + 1$ . For each $β_{j^{'}}$ and $δ_{j^{'}}$ , $j^{'} = 1, 2, \dots, 8$ ; we consider independent normal (0,100) priors. We consider inverse gamma (0.001,0.001) prior for $σ_{ε}^{2}$ . For $2 \times 2$ covariance matrix $Σ$ , we consider an inverse Wishart (IW) prior with the identity matrix of order 2 ( $I_{2}$ ) as scale matrix. Note that most of these are flat priors. Nevertheless, we perform a sensitivity analysis for investigating the effects of the priors on the final estimates. Results for the sensitivity analysis for some of the model parameters are summarized in Table 1, and we note that the priors have minimal effects on the final estimates. Similar results are obtained for the other parameters as well.

3.2 Computational details

First, we check the presence of multicollinearity in our data. Note that we have five covariates with time-varying effects and eight covariates with fixed effects. At each wave, we compute the VIF for each predictor variable, and it turns out the computed VIFs are all below 7. This indicates the absence of multicollinearity in our dataset, and hence we impute the missing covariates separately from a set of univariate models. We impute missing calendar age directly since we have eight waves, and the length of each wave is two years. For the other missing covariates, we use a model-based approach for imputation. We note that the missing covariates should not be imputed for each wave, simply because subjects are of different ages at a particular wave. Hence, we divide the subjects into five ‘age groups’ across all the eight waves: 50–54 (group 1), 55–59 (group 2), 60–64 (group 3), 65–69 (group 4) and ≥70 (group 5). The subjects belonging to the same age group are expected to be similar in terms of physical health and insurance status. Covariates are imputed for each age group. A wider age interval (for example, 50–59) might be inefficient for grouping simply because health conditions change for the older people at a much faster rate. Our choice for five groups in this dataset is subjective, but we consider this reasonable enough for effective imputation of the covariates. Thus, for imputing the missing covariates, we ignore the ‘waves’ but consider the aforementioned ‘age groups’. However, we note that our grouping (based on calendar age) is only for imputing the covariates, and not for any other inference purpose. For each age group, we fit a suitable probability distribution for each covariate based on the available data. For BMI, we fit $N (μ_{1}, σ_{1}^{2})$ distribution and consider weakly informative priors for $μ_{1}$ and $σ_{1}^{2}$ . We consider $N (0, 1000)$ prior for $μ_{1}$ and IG (0.001,0.002) prior for $σ_{1}^{2}$ . For each of the binary covariates for disease, we fit Bernoulli ( $π^{*}$ ) distribution with uniform (0,1) prior on $π^{*}$ . We note that each subject at any particular wave belongs to one and exactly one health insurance category. Hence we fit a multinomial ( $π_{1}, π_{2}, π_{3})$ distribution for three insurance variables. For $(π_{1}, π_{2}, π_{3})$ , we consider a flat Dirichlet prior (each concentration parameter = 1). Missing covariates are imputed from their respective posterior predictive densities as discussed in Section 2.4. Once the missing covariates are imputed, we impute missing responses similar to Section 2.4. The missing non-zero responses are imputed only under ACMV restriction. The other two restrictions (CCMV, NCMV) are not used in our data analysis since in the simulation study (Section 4), we demonstrate the effectiveness of ACMV restriction over the other two identifying restrictions. Once the missing observations are imputed, we use the complete data and fit the joint model given in Equation (2.5). Posterior distribution is sampled using Gibbs sampler. The full conditionals and the computational steps are given in detail in the Appendix. We run 50 000 iterations and discard the first 2 000 iterations as ‘burn-in’. We thin the chains by keeping every 10th iteration. The model parameters are estimated by their respective posterior sample means. The computed scale reduction factors (Brooks and Gelman, 1998) are below 1.1, indicating the convergence of the chains. We consider 1 000 replications of the complete dataset and the aforementioned process is repeated 1,000 times. We compute BIC values for different orders (K) of LOP for modelling $α_{j} (t)$ and $γ_{j} (t)$ in Equations (3.1) and (3.2). We consider $k = 1, 2, 3, 4$ and obtain the smallest BIC for $k = 2$ . The details of the Gibbs sampler are given in the Appendix. All the computations are performed in R (code is available with this article) and the total computational time (including 1 000 replications of the complete dataset) is nearly 29 hours in HPC.

3.3 Results

In Table 2, we show the parameter estimates and the corresponding 95% credible intervals computed from MCMC iterations for the eight predictors with fixed effects on the response. Results are shown for the dynamic models given in the Equations (3.1) and (3.2). We note that for the dynamic model in Equation (3.1), the coefficients corresponding to cancer, heart problem and stroke are significantly different from zero, and the corresponding Bayesian credible intervals do not contain zero. Similar results are obtained for the dynamic model in Equation (3.2). The results are reassuring the fact that serious condition of the patients suffering from cancer, stroke or cardiovascular diseases are associated with frequent hospitalization. The additional cost due to cancer, stroke and heart problem are not, in general, covered by health insurance. We also note that the effects of blood pressure, diabetes and psychological problems are not significant for both the models. The coefficients for lung problem and arthritis are different from zero, but the corresponding 95% credible intervals contain zero. This again accords to our intuition that these diseases are chronic and hence not associated with frequent hospitalization. Also the medical costs for treating these diseases are covered by the health insurance.

Table 2:

Fixed covariate effects for two dynamic models in the HRS data

	First dynamic model		Second dynamic model
Predictor	Coefficient estimate	95% CI	Coefficient estimate	95% CI
Blood pressure	0.065	(−1.89, 0.87)	0.54	(−1.03, 0.88)
Diabetes	0.23	(−1.45, 0.81)	0.31	(−0.98, 1.16)
Cancer	4.64	(2.76, 6.11)	5.08	(3.88, 7.08)
Lung problem	1.33	(−0.65, 2.15)	0.84	(−1.18, 1.25)
Heart problem	3.62	(2.03, 5.14)	2.35	(0.68, 4.23)
Stroke	7.16	(5.91, 9.29)	6.55	(4.87, 8.16)
Arthritis	1.65	(−1.96, 2.84)	2.33	(−0.89, 3.86)
Psychological problem	0.055	(−1.52, 0.81)	0.048	(−0.39, 0.96)

Figure 1:

Time-varying effects of different health insurances in probit model

In Figure 1, we plot the time-varying effects of three insurances for probit model and also plot the corresponding 95% point-wise credible intervals. We note that the curve for employment insurance shows a slow decreasing trend. This indicates that employment insurances cover less for the relatively older people, and hence the effect of this insurance on the probability of zero/non-zero OOPME gradually decreases. The increasing trend of the curve for government insurance after wave 4 reflects that the medicare insurance is more useful than the medicaid insurance. For the other (private) insurance, we observe a sharp decreasing trend after wave 4, indicating that the private insurances are less effective for the relatively older people (age 62 years or above).

Figure 2:

Time-varying effects of calendar age and BMI in probit model

In Figure 2, we plot the time-varying effects of BMI and calendar age. We also plot the corresponding 95% point-wise confidence intervals. Clearly, the effect of the calendar age has a strict increasing trend after wave 3 on the zero-response probability. This is quite intuitive since the older people visit hospital more frequently and hence will result in higher OOPME compared to the relatively younger people. For BMI, the curve shows an increasing trend, indicating that obesity has higher effect on the probability of zero-response for the older people. Since there is a well-known relationship between obesity and many other diseases, this trend accords to our expectation.

Figure 3:

Time-varying effects of different health insurances on non-zero out-of-pocket medical expenditure

Figure 3 shows the effects of the insurances on the non-zero OOPMEs. We observe similar trends for the insurance groups as in Figure 1. The explanations will go along the same line. As mentioned in Min and Agresti (2002), for a two-part model, the inference on the covariates are typically same for both the models. Our approach considers the same set of covariates for both the models and hence get similar inference from those.

Figure 4:

Time-varying effects of calendar age and BMI on non-zero out-of-pocket medical expenditure

In Figure 4, we show the time-varying effects of BMI and the calendar age on the non-zero OOPME, and the corresponding 95% point-wise credible intervals. The trends are similar to the trends observed in Figure 2. Thus, we have the similar statistical inferences for two dynamic model per our expectation. From the estimated covariance matrix $Σ$ , we get $\hat{var} (b_{i}) = 3.89$ , $\hat{var} (c_{i}) = 4.66$ and $\hat{corr} (b_{i}, c_{i}) = 0.74$ . We note that a high value of the correlation between $b_{i}$ and $c_{i}$ indicates that the dynamic models given in Equations (3.1) and (3.2) are highly linked.

4 Simulation study

We perform a simulation study to investigate the operating characteristics of the proposed approach of modelling incomplete zero-inflated longitudinal response. We consider 200 subjects for which a continuous response is measured at 10 different evenly spaced time points. The response for the ith subject ( $i = 1, 2, \dots, 200$ ) at the Tth time point ( $t = 1, 2, \dots, 10$ ) is simulated from the following mixture density: $Y_{it} \sim (1 - p_{it}) δ_{0} + p_{it} N (μ_{it}, σ^{2})$ , where $δ_{0}$ denotes a point mass at zero. The non-zero-probabilities $p_{it}$ are generated from the following probit model:

Φ^{- 1} (p_{it}) = \sum_{j = 1}^{2} α_{j} (t) x_{ij} (t) + \sum_{j^{'} = 1}^{2} β_{j^{'}} z_{{ij}^{'}} (t) + b_{i},

where $β_{1} = 1.4$ , $β_{2} = 3.5$ . We consider $α_{1} (t) = a_{1} + a_{2} t + a_{3} t^{2}$ and $α_{2} (t) = a_{4} + a_{5} t$ , with $a_{1} = 0.87, a_{2} = 2.1, a_{3} = 0.64, a_{4} = 1.34 and a_{5} = 0.94$ . Two covariates $x_{1} (t)$ and $x_{2} (t)$ with time-varying effects are simulated independently from $N (5, t)$ and $N (t, 8)$ , respectively, at each time point T. The covariates $z_{1} (t)$ and $z_{2} (t)$ with fixed effects are simulated from Gamma $(0.5 t, 3)$ and Gamma $(2.5, 0.4 t)$ , respectively, at each T. For non-zero response, the mean function is given in the following equation:

μ_{i} (t_{il}) = \sum_{j = 1}^{2} γ_{j} (t_{il}) x_{ij} (t_{il}) + \sum_{j^{'} = 1}^{2} δ_{j^{'}} z_{{ij}^{'}} (t_{il}) + c_{i},

where we consider $γ_{1} (t) = a_{6} + a_{7} t$ and $γ_{2} (t) = a_{8} + a_{9} t + a_{10} t^{2}$ with $a_{6} = 1.45, a_{7} = 0.8, a_{8} = 1.07, a_{9} = 0.45$ and $a_{10} = 1.32$ . We generate $(b_{i}, c_{i})$ from a bivariate normal density with mean vector = (0,0) and variance $(b_{i}) = 2.5$ , variance $(c_{i}) = 1$ , correlation $(b_{i}, c_{i}) = 0.65$ . Note that two dynamic models for our simulation study are linked by considering correlation between $b_{i}$ and $c_{i}$ .

Once the complete dataset is simulated, we create ‘monotone missingness’ in the data. We choose 30 subjects randomly, and for each of these 30 subjects we randomly select a time point T between 3 and 9 after which the subject is missing. Note that for a missing subject we consider missing response as well as missing covariates. We analyse the generated incomplete data using the method discussed in Section 2. For the imputation of the missing covariates, we use the method proposed in Section 2.4, and we impute the missing responses under three different model restrictions: CCMV, ACMV and NCMV. First, we analyse the complete data with no missing value and estimate model parameters for two dynamic models. We replicate 100 datasets and estimate the model parameters for each dataset separately. Based on replicated data, we compute the average bias and average width of 95% credible intervals, as well as the corresponding coverage probability. Next, we consider the incomplete data generated by the aforementioned method and impute the missing response values under three model restrictions. We consider 100 replications of such incomplete dataset and compute average bias, width of 95% CI and corresponding coverage probability. The results are summarized in Table 3.

Table 3:

Estimated bias, width of 95% CI and coverage probabilities of model parameters for complete data and incomplete data under different imputation methods for the simulation study

	Complete data		NCMV		ACMV		CCMV
Parameter	Bias	CI width(Cov.P)	Bias	CI width(Cov.P)	Bias	CI width(Cov.P)	Bias	CI width(Cov.P)
$a_{1}$	0.11	0.91(0.95)	0.24	1.16(0.96)	0.22	1.08(0.95)	0.56	1.48(0.96)
$a_{2}$	0.10	0.89(0.95)	0.28	1.21(0.95)	0.31	1.25(0.95)	0.63	1.55(0.96)
$a_{3}$	0.09	0.98(0.95)	0.29	1.33(0.95)	0.22	1.28(0.95)	0.59	1.63(0.96)
$a_{4}$	0.11	0.87(0.95)	0.34	1.26(0.96)	0.21	1.32(0.95)	0.62	1.54(0.96)
$a_{5}$	0.13	1.03(0.95)	0.38	1.39(0.96)	0.28	1.22(0.95)	0.65	1.58(0.96)
$a_{6}$	0.12	0.78(0.95)	0.25	0.96(0.95)	0.23	0.84(0.95)	0.51	1.21(0.96)
$a_{7}$	0.11	1.13(0.96)	0.26	1.46(0.96)	0.25	1.39(0.96)	0.39	1.68(0.97)
$a_{8}$	0.10	0.99(0.95)	0.44	1.36(0.95)	0.39	1.07(0.95)	0.55	1.53(0.96)
$a_{9}$	0.11	0.74(0.94)	0.29	1.19(0.96)	0.26	1.03(0.54)	0.53	1.42(0.96)
$a_{10}$	0.08	1.15(0.95)	0.28	1.29(0.95)	0.25	1.24(0.95)	0.56	1.63(0.96)
$β_{1}$	0.13	0.67(0.94)	0.27	1.15(0.95)	0.19	0.93(0.95)	0.46	1.47(0.96)
$β_{2}$	0.11	0.87(0.95)	0.37	1.28(0.95)	0.29	1.23(0.95)	0.66	1.75(0.96)
$δ_{1}$	0.09	0.77(0.94)	0.37	1.13(0.95)	0.39	1.18(0.95)	0.48	1.59(0.96)
$δ_{2}$	0.11	0.69(0.95)	0.27	0.85(0.95)	0.29	0.83(0.95)	0.51	1.28(0.96)

We note that for complete data, our proposed approach estimates the model parameters with very low bias and moderate width of the CIs with desirable coverage probabilities. For the incomplete data almost for all the model parameters, ACMV provides lower bias and shorter CIs. CCMV shows the worst performance. Both CCMV and NCMV result in larger bias and wider CIs compared to ACMV for most of the model parameters. For some parameters, NCMV and ACMV provide similar inference, but CCMV performs worse than the others. We also assess the accuracy of the first dynamic model in predicting a missing zero/non-zero response. For each replicated dataset, we compute the percentage of ‘mismatch’ of the predicted and the actual zero or non-zero response. Based on 100 replicates, we find that our model provides 86% ‘match’ with a standard error of 0.43. Next, we compare the performance of the proposed method to the commonly used Tobit model and the two-part model for our simulated (incomplete) data without imputing the missing data. We consider 100 replicates of the dataset, and for each dataset we fit five competing models: (a) Tobit model for incomplete data, (b) two-part model for incomplete data, (c) proposed model with CCMV imputation, (d) proposed model with ACMV imputation and (e) proposed model with NCMV imputation. In the first two cases missingness is modeled, and in the last three cases missing observations are imputed. We compute three important model selection criteria, for example, (a) mean squared error (MSE), (b) log pseudo-marginal likelihood and (c) deviance information criterion for all the models under consideration. For each replicated dataset, once we fit the models, we compute MSE. Based on 100 replicates, we obtain 100 MSE values for all the five models and compute average MSE (AMSE). A smaller AMSE value indicates a better fit. Second, we use the conditional predictive ordinate (CPO; Gelfand et al., 1992) defined as ${CPO}_{i} = P (y_{i} | y_{- i}) = E_{θ} [P (y_{i} | θ, y_{- i})]$ , where $y_{- i}$ denotes the data excluding $y_{i}$ and $θ$ denotes the set of all model parameters. We compute log pseudo-marginal likelihood (LPML)= $\sum_{i = 1}^{n} \log {\hat{CPO}}_{i}$ . A higher LPML value indicates a better fit.

Finally, we compute conditional deviance information criteria (DIC) proposed in Celeux et al. (2006). This DIC is based on the conditional likelihood $l (y | γ)$ and is given by DIC $= - 4 E + 2 \log l (y | \hat{κ})$ , where $κ$ denotes the vector of random effects and $\hat{κ}$ is the corresponding estimate (posterior mean). A smaller value indicates a better fit.

In Table 4, we show three different measures for goodness of fit. We note that Tobit model and two-part model based on the incomplete data provide poor fit, that is, higher AMSE, lower LPML and higher DIC values. The proposed model with ACMV imputation provides the smallest AMSE and largest LPML. For DIC, the proposed model with NCMV imputation and ACMV imputation are quite comparable. Thus, the effectiveness of the proposed model and ACMV restriction is reassured.

Table 4:

Goodnessof fit measures for different competing models in the simulation study

Model	AMSE	LPML	DIC
Tobit model	48.14	−245.37	64.73
Two-part model	41.08	−223.86	45.19
Proposed model with CCMV	32.94	−188.47	38.25
Proposed model with ACMV	24.18	$-$ 138.26	26.84
Proposed model with NCMV	29.36	−151.43	25.93

note: Values in bold are optimal values.

5 Discussion

Missing observations in the longitudinal study are likely to be informative and hence the analysis of incomplete longitudinal data is challenging. In this article we have addressed the issue of analysing zero-inflated incomplete longitudinal response in a Bayesian framework using a two-part model. Bayesian approach makes the data-augmentation simple and computationally fast because it is based on Gibbs sampler. Our approach is based on the assumption that the covariates are uncorrelated, and hence we impute each covariate separately from their respective posterior predictive distributions. However, in some practical application the covariates can be correlated. In that case, a joint model can be developed for simultaneously imputing all the variables through MCMC iterations. We note that in real applications, longitudinal data can be unbalanced in many different ways. In our current work (both in simulation study and real application), we have considered only the unbalanced case due to monotone missingness. Our estimation method and hence our inference heavily depends on this assumption. Thus our findings should not be generalized for any kind of unbalanced longitudinal dataset and should be interpreted with caution.

Supplementary materials

Supplementary materials for this article (sample data and R code) are available from http://www.statmod.org/smij/archive.html.

Appendix

Some basic notations:

$η = (η_{01}, η_{11}, \dots, η_{r 1}, \dots, η_{0 J}, η_{2 J}, \dots, η_{rJ})^{T}$ ; $η^{*} = (η_{01}^{*}, η_{11}^{*}, \dots, η_{r^{*} 1}^{*}, \dots, η_{0 J}^{*}, η_{1 J}^{*}, \dots, η_{r^{*} J}^{*})^{T}$ ; $β = (β_{1}, β_{2}, β_{3}, \dots, β_{J^{'}})^{T}$ ; $δ = (δ_{1}, δ_{2}, δ_{3}, \dots, δ_{J^{'}})^{T}$ ; $x_{t} (t) = (x_{i 1} (t), x_{i 2} (t), x_{i 3} (t), \dots, x_{iJ} (t))^{T}$ ;

$z_{i} (t) = (z_{i 1} (t), z_{i 2} (t), z_{i 3} (t), \dots, z_{{iJ}^{'}} (t))^{T}$ ; $t^{'} = - 1 + 2 \frac{t - t_{\min}}{t_{\max} - t_{\min}}$ (i.e., $t^{'} \in [- 1, 1]$ );

$P^{(r)} (t^{'}) = (P_{0} (t^{'}), P_{1} (t^{'}), \dots, P_{r} (t^{'}))^{T}$ ; $P^{(r^{*})} (t^{'}) = (P_{0} (t^{'}), P_{1} (t^{'}), \dots, P_{r^{*}} (t^{'}))^{T}$ ;

$x_{i} (t) \otimes P^{(r)} (t^{'}) = (x_{i 1} (t) {P^{(r)} (t^{'})}^{T}, x_{i 2} (t) {P^{(r)} (t^{'})}^{T}, \dots, x_{iJ} (t) {P^{(r)} (t^{'})}^{T})^{T}$ ;

$b = (b_{1}, b_{2}, b_{3}, \dots, b_{n})^{T}$ ; $c = (c_{1}, c_{2}, c_{3}, \dots, c_{n})^{T}$ . Here, $i = 1, 2, 3, \dots, n$ ; $t = 1, 2, 3, \dots, T$ ; $l = 1, 2, 3, \dots, T_{i}$ . Here, $θ ∣ -$ denotes the full conditional distribution of θ, given rest of the parameters, and $μ_{θ}$ is a mean vector of the corresponding parameter θ, $diag (θ_{1}, θ_{2}, \dots, θ_{m})$ means a diagonal matrix of order $m \times m$ in which diagonal entries are $θ_{1}, θ_{2}, \dots, θ_{m}$ , respectively. Note that r and $r^{*}$ may or may not be same. Here, $P_{k} (t^{'})$ denotes Legendre orthogonal polynomial at time point $t^{'}$ . $I_{m}$ is an identity matrix of order $m$ and IG stands for inverse gamma distribution and IW stands for Inverse Wishart distribution.

Dynamic probit model

The model is:

$R_{it} = \sum_{j = 1}^{J} α_{j} (t) x_{ij} (t) + \sum_{j^{'} = 1}^{J^{'}} β_{j^{'}} z_{{ij}^{'}} (t) + b_{i} + e_{it}$ , where $e_{it} iid \sim N (0, 1)$ and $b_{i}$ are independent. Here, $α_{j} (t) = η_{0 j} P_{0} (t^{'}) + η_{1 j} P_{1} (t^{'}) + η_{2 j} P_{2} (t^{'}) + \dots + η_{rj} P_{r} (t^{'}) .$

Rewrite the earlier model: $R_{it} = (x_{i} (t) \otimes P^{(r)} (t^{'}))^{T} η + {z_{i} (t)}^{T} β + b_{i} + e_{it} .$

The priors are the following:

$η \sim N (μ_{η}, σ_{η}^{2} I_{J (r + 1)})$ ; $β \sim N (μ_{β}, σ_{β}^{2} I_{J^{'}}) .$

Dynamic model for non-zero observations

The model is:

$y_{i} (t_{il}) = \sum_{j = 1}^{J} γ_{j} (t_{il}) x_{ij} (t_{il}) + \sum_{j^{'} = 1}^{J^{'}} δ_{j^{'}} z_{{ij}^{'}} (t_{il}) + c_{i} + e_{{it}_{il}}^{*}$ , where $e_{it}^{*} iid \sim N (0, σ_{e^{*}}^{2})$ and $c_{i}$ are independent. Also $(b_{i}, c_{i}) \sim N (0, Λ)$ , where $b_{i}$ comes from the first dynamic model and $c_{i}$ comes from the second dynamic model and $0 = (0, 0)^{T}$ and $Λ$ is a $2 \times 2$ p.d. matrix where $Λ^{- 1} = ((λ_{ij}))_{2 \times 2}$ . Here, $γ_{j} (t_{il}) = η_{0 j}^{*} P_{0} (t_{il}^{'}) + η_{1 j}^{*} P_{1} (t_{il}^{'}) + η_{2 j}^{*} P_{2} (t_{il}^{'}) + \dots + η_{r^{*} j}^{*} P_{r^{*}} (t_{il}^{'}) .$

Rewrite the earlier model: $y_{i} (t_{il}) = (x_{i} (t_{il}) \otimes P^{(r^{*})} (t_{il}^{'}))^{T} η^{*} + {z_{i} (t_{il})}^{T} δ + c_{i} + e_{{it}_{il}}^{*} .$

The priors are the following:

$η^{*} \sim N (μ_{η^{*}}, σ_{η^{*}}^{2} I_{J (r^{*} + 1)})$ ; $δ \sim N (μ_{δ}, σ_{δ}^{2} I_{J^{'}})$ ; $σ_{e^{*}}^{2} \sim IG (α_{e^{*}}, β_{e^{*}})$ ; $Λ \sim IW (ν, Ψ) .$

Full conditional distributions:

$η ∣ - \sim N (V_{η} v_{η}, V_{η})$ , where $v_{η} = \sum_{i = 1}^{n} \sum_{t = 1}^{T} (R_{it} - {z_{i} (t)}^{T} β - b_{i}) (x_{i} (t) \otimes P^{(r)} (t^{'})) + (\frac{1}{σ_{η}^{2}}) μ_{η}$ ; $V_{η}^{- 1} = \sum_{i = 1}^{n} \sum_{t = 1}^{T} (x_{i} (t) \otimes P^{(r)} (t^{'})) (x_{i} (t) \otimes P^{(r)} (t^{'}))^{T} + (\frac{1}{σ_{η}^{2}}) I_{J (r + 1)} .$

$β ∣ - \sim N (V_{β} v_{β}, V_{β})$ , where $v_{β} = \sum_{i = 1}^{n} \sum_{t = 1}^{T} (R_{it} - (x_{i} (t) \otimes P^{(r)} (t^{'}))^{T} η - b_{i}) z_{i} (t) + (\frac{1}{σ_{β}^{2}}) μ_{β}$ ; $V_{β}^{- 1} = \sum_{i = 1}^{n} \sum_{t = 1}^{T} z_{i} (t) {z_{i} (t)}^{T} + (\frac{1}{σ_{β}^{2}}) I_{J^{'}} .$

$b_{i} | - \sim N (V_{b i v b i}, V_{b i}),$ where

$v_{b i} = \sum_{t = 1}^{T} (R_{i t} - {(x_{i} (t) \otimes P^{(r)} (t^{'}))}^{T} η - z_{i} {(t)}^{T} β) - λ_{12 c i}; V_{b i}^{- 1} = (T + λ_{11}) .$

$R_{i t} | - \sim {\begin{matrix} N ({(x_{i} (t) \otimes P^{(r)} (t^{'}))}^{T} η + z_{i} {(t)}^{T} β + b_{i}, 1); left turncated by zero when E_{i t} = 1; \\ N ({(x_{i} (t) \otimes P^{(r)} (t^{'}))}^{T} η + z_{i} {(t)}^{T} β + b_{i}, 1); right turncated by zero when E_{i t} = 0. \end{matrix}$

$η^{*} ∣ - \sim N (V_{η^{*}} v_{η^{*}}, V_{η^{*}})$ , where $v_{η^{*}} = \frac{1}{σ_{e^{*}}^{2}} \sum_{i = 1}^{n} \sum_{l = 1}^{T_{i}} (y_{i} (t_{il}) - {z_{i} (t_{il})}^{T} δ - c_{i}) (x_{i} (t_{il}) \otimes P^{(r^{*})} (t_{il}^{'})) + (\frac{1}{σ_{η^{*}}^{2}}) μ_{η^{*}}$ ; $V_{η^{*}}^{- 1} = \frac{1}{σ_{e^{*}}^{2}} \sum_{i = 1}^{n} \sum_{l = 1}^{T_{i}} (x_{i} (t_{il}) \otimes P^{(r^{*})} (t_{il}^{'})) (x_{i} (t_{il}) \otimes P^{(r^{*})} (t_{il}^{'}))^{T} + (\frac{1}{σ_{η^{*}}^{2}}) I_{J (r^{*} + 1)} .$

$δ ∣ - \sim N (V_{δ} v_{δ}, V_{δ})$ , where $v_{δ} = \frac{1}{σ_{e^{*}}^{2}} \sum_{i = 1}^{n} \sum_{l = 1}^{T_{i}} (y_{i} (t_{il}) - (x_{i} (t_{il}) \otimes P^{(r^{*})} (t_{il}^{'}))^{T} η^{*} - c_{i}) z_{i} (t_{il}) + (\frac{1}{σ_{δ}^{2}}) μ_{δ}$ ; $V_{δ}^{- 1} = \frac{1}{σ_{e^{*}}^{2}} \sum_{i = 1}^{n} \sum_{l = 1}^{T_{i}} z_{i} (t_{il}) {z_{i} (t_{il})}^{T} + (\frac{1}{σ_{δ}^{2}}) I_{J^{'}} .$

$c_{i} ∣ - \sim N (V_{c_{i}} v_{c_{i}}, V_{c_{i}}),$ where

$v_{c_{i}} = \frac{1}{σ_{e^{*}}^{2}} \sum_{l = 1}^{T_{i}} (y_{i} (t_{il}) - (x_{i} (t_{il}) \otimes P^{(r^{*})} (t_{il}^{'}))^{T} η^{*} - {z_{i} (t_{il})}^{T} δ)$ ; $V_{c_{i}}^{- 1} = (\frac{T_{i}}{σ_{e^{*}}^{2}} + λ_{22}) .$

$σ_{e^{*}}^{2} ∣ - \sim IG (\frac{1}{2} \sum_{i = 1}^{n} T_{i} + α_{e^{*}}, \frac{1}{2} \sum_{i = 1}^{n} \sum_{l = 1}^{T_{i}} (y_{i} (t_{il}) - (x_{i} (t_{il}) \otimes P^{(r^{*})} (t_{il}^{'}))^{T} η^{*} - {z_{i} (t_{il})}^{T} δ - c_{i})^{2} + β_{e^{*}}) .$

$Λ ∣ - \sim IW (n + ν, Ψ + \sum_{i = 1}^{n} (b_{i}, c_{i})^{T} (b_{i}, c_{i}))$

Acknowledgments

The authors thank Professor Pulak Ghosh, Indian Institute of Management, Bangalore, India, for sharing the HRS data.

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship and/or publication of this article.

Funding

The authors received no financial support for the research, authorship and/or publication of this article.

References

Aitkin

Rubin

(1985) Estimation and hypothesis testing in finite mixture models. Journal of the Royal Statistical Society, Series B , 47, 67–75.

Albert

Chib

(1993) Bayesian analysis of binary and polychotomous response data. Journal of the American Statistical Association , 88, 669–79.

Alfo

Maruotti

(2010) Two-part regression models for longitudinal zero-inflated count data. Canadian Journal of Statistics , 38, 197–216.

Amemiya

(1984) Tobit models: A survey. Journal of Econometrics , 24, 3–61.

Arias

(2006) United States Life Tables. National Vital Statistics Reports , 58, 1–40.

Bhuyan

, Biswas

, Ghosh

Das

(2018) A Bayesian two-stage regression approach of analysing longitudinal outcomes with endogeneity and incompleteness. Statistical Modelling , doi:~1471082X17747806

Boone

, Stewart- Koster

Kennard

(2012) A hierarchical zero-inflated Poisson regression model for stream fish distribution and abundance. Environmetrics , 23, 207–18.

Brooks

Gelman

(1998) General methods for monitoring convergence of iterative simulations. Journal of Computational and Graphical Statistics , 7, 434–55.

Brown

, Ghosh

Taylor

(2015) Modelling household finances: A Bayesian approach to a multivariate two-part model. Journal of Empirical Finance , 33, 190–207.

10.

Buu

, Li

, Tan

Zucker

(2012) Statistical models for longitudinal zero inflated count data with applications to the substance abuse field. Statistics in Medicine , 31, 4074–86.

11.

Celeux

, Forbes

, Robert

Titterington

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

12.

Chen

(2011) Model selection for zero-inflated regression with missing covariates. Computational Statistics and Data Analysis , 55, 765–73.

13.

Cui

, Zhu

(2006) Functional mapping for genetic control of programmed cell death. Physiological Genomics , 25, 458–69.

14.

Daniels

Hogan

(2008) Missing Data in Longitudinal Studies: Strategies for Bayesian Modeling and Sensitivity Analysis . London: Chapman and Hall (CRC Press).

15.

Das

, Afriyie

Spirko

(2015) A semi-parametric Bayesian model for analysing longitudinal data from multiple related groups. International Journal of Biostatistics , 11, 273–84.

16.

Das

Daniels

(2014) A semiparametric approach to simultaneous covariance estimation for bivariate sparse longitudinal data. Biometrics , 70, 33–43.

17.

Das

, Li

, Sengupta

(2013) A Bayesian semi-parametric model for bivariate sparse longitudinal data. Statistics in Medicine , 32, 3899–3910.

18.

Das

, Pareek

, Brown

Ghosh

(2017) A Semiparametric Bayesian Approach to a New Dynamic Zero-inflated Model. Department of Economics Working Papers, The University of Sheffield, Sheffield, UK.

19.

Diggle

Kenward

(1994) Informative drop-out in longitudinal data analysis (with discussion). Applied Statistics , 43, 49–93.

20.

Dobriansky

, Suzman

Hodes

(2007) Why population aging matters: A global perspective . Bethesda, MD and Washington, DC: National Institute on Aging, National Institutes of Health, US Department of Health and Human Services, US Department of State.

21.

Duan

, Manning

, Morris

Newhouse

(1983) A comparison of alternative models for the demand for medical care (Corr: V2 P413). Journal of Business and Economic Statistics , 1, 115–26.

22.

Farewell

, Long

, Tom

BDM

, Yiu

(2017) Two-part and related regression models for longitudinal data. Annual Review of Statistics and Its Application , 4, 283–315.

23.

Gelfand

, Dey

Chang

(1992) Model determinating using predictive distributions with implementation via sampling-based methods (with discussion). Bayesian Statistics , 4, 147–67.

24.

Ghosh

(2008) Assessing sexual attitudes and behaviors of young women: A joint model with nonlinear time effects, time varying covariates, and dropouts. Journal of the American Statistical Association , 103, 1496–1507.

25.

Hall

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

26.

Hušková

Sen

(1985) On sequentially adaptive asymptotically efficient rank statistics. Sequential Analysis , 4, 125–51.

27.

Laird

Ware

(1982) Random effects model for longitudinal data. Biometrics , 38, 963–74.

28.

Lambert

(1992) Zero-inflated Poisson regression, with an application to defects in manufacturing. Technometrics , 34, 1–14.

29.

Little

(1993) Pattern-mixture models for multivariate incomplete data. Journal of the American Statistical Association , 88, 125–34.

30.

Little

(1995) Modeling the drop-out mechanism in repeated-measures studies. Journal of the American Statistical Association , 90, 1112–21.

31.

Maruotti

(2011) A two-part mixed-effects pattern-mixture model to handle zeroinflation and incompleteness in a longitudinal setting. Biometrical Journal , 53, 716–34.

32.

Meyer

(2000) Random regressions to model phenotypic variation in monthly weights of Australian beef cows. Livestock Production Science , 65, 19–38.

33.

Min

Agresti

(2002) Modeling nonnegative data with clumping at zero: A survey. Journal of Iranian Statistical Society , 1, 7–33.

34.

Molenberghs

, Michiels

, Kenward

Diggle

(1998) Missing data mechanisms and pattern mixture models. Statistica Neerlandica , 52, 153–61.

35.

Mukherji

, Roychoudhury

, Ghosh

Brown

(2016) Estimating health demand for an aging population: A flexible and robust Bayesian joint model. Journal of Applied Econometrics , 31, 1140–1158.

36.

Olsen

Schafer

(2001) A two-part random-effects model for semi-continuous longitudinal data. Journal of the American Statistical Association , 96, 730–45.

37.

Pan

Mackenzie

(2003) On modelling mean-covariance structures in longitudinal studies. Biometrika , 90, 239–44.

38.

Pourahmadi

(1999) Joint mean-covariance models with applications to longitudinal data: Unconstrained parameterisation. Biometrika , 86, 677–90.

39.

Pourahmadi

(2000) Maximum likelihood estimation of generalized linear models for multivariate normal covariance matrix. Biometrika , 87, 425–35.

40.

Rubin

(1976) Inference and missing data. Biometrika , 63, 581–92.

41.

Schafer

Yucel

(2002) Computational strategies for multivariate linear mixedeffects models with missing values. Journal of Computational and Graphical Statistics , 11, 437–57.

42.

Seaman

, Jonathan

White

(2012) Multiple imputation of missing covariates with non-linear effects and interactions: An evaluation of statistical methods. BMC Medical Research Methodology . doi: 10.1186/1471-2288-12-46.

43.

Xie

Paik

(1997) Multiple imputation methods for the missing covariates in generalized estimating equation. Biometrics , 53, 1538–46.

44.

Yang

, Li

, Zucker

Buu

(2016) Two-stage model for time-varying effects of zero-inflated count longitudinal covariates with applications in health behaviour research. Journal of the Royal Statistical Society Series C , 65, 431–44.