Journal of Probability and Statistics

Journal of Probability and Statistics / 2012 / Article
Special Issue

Joint Models and Their Applications

View this Special Issue

Research Article | Open Access

Volume 2012 |Article ID 614102 | 19 pages |

Mixed-Effects Tobit Joint Models for Longitudinal Data with Skewness, Detection Limits, and Measurement Errors

Academic Editor: Lang Wu
Received29 May 2011
Accepted13 Aug 2011
Published18 Oct 2011


Complex longitudinal data are commonly analyzed using nonlinear mixed-effects (NLME) models with a normal distribution. However, a departure from normality may lead to invalid inference and unreasonable parameter estimates. Some covariates may be measured with substantial errors, and the response observations may also be subjected to left-censoring due to a detection limit. Inferential procedures can be complicated dramatically when such data with asymmetric characteristics, left censoring, and measurement errors are analyzed. There is relatively little work concerning all of the three features simultaneously. In this paper, we jointly investigate a skew-t NLME Tobit model for response (with left censoring) process and a skew-t nonparametric mixed-effects model for covariate (with measurement errors) process under a Bayesian framework. A real data example is used to illustrate the proposed methods.

1. Introduction

Modeling of longitudinal data is an active area of biostatistics and statistics research that has received a lot of attention in the recent years. Various statistical modeling and analysis methods have been suggested in the literature for analyzing such data with complex features (Higgins et al. [1], Liu and Wu [2], Wulfsohn and Tsiatis [3], and Wu [4]). However, there is a relatively little work done on simultaneously accounting for skewness, left censoring due to a detection limit (for example, a threshold below which viral loads are not quantifiable) and covariate measurement errors, which are inherent features of longitudinal data. This paper proposes a joint skew-𝑑 NLME Tobit model for a response and measurement errors in covariate by simultaneously accounting for left-censoring and skewness. Thus, the proposed model addresses three important features of longitudinal data such as viral load in an AIDS study.

Firstly, our model relaxes the normality assumption for random errors and random-effects by using flexible skew-normal and skew-𝑑 distributions. It has been documented in the literature that the normality assumption lacks robustness against extreme values, obscures important features of between- and within-subject variations, and leads to biased or misleading results (Huang and Dagne [5], Verbeke and Lesaffre [6], and Sahu et al. [7]). Specially, nonnormal characteristics such as skewness with heavy tails appear very often in virologic responses. For example, Figures 1(a) and 1(b) displays the histograms of repeated viral load (in ln scale) and CD4 cell count measurements for 44 subjects enrolled in an AIDS clinical study (Acosta et al. [8]). For this data set, which is analyzed in this paper, both viral load (even after ln-transformation) and CD4 cell count are highly skewed, and thus a normality assumption may be violated.

Secondly, an outcome of a longitudinal study may be subject to a detection limit because of low sensitivity of current standard assays (Perelson et al. [9]). For example, for a longitudinal AIDS study, designed to collect data on every individual at each assessment, the response (viral load) measurements may be subject to left censoring due to a detection limit of quantification. Figures 1(c) and 1(d) shows the measurements of viral load and CD4 cell count for three randomly selected patients in the study. We can see that for some patients their viral loads are below detection limit (BDL), which is 50 (in copies/mL). When observations fall below the BDL, a common practice is to impute the censored values by either the detection limit or half of the detection limit (Wu [4], Ding and Wu [10], and Davidian and Giltinan [11]). Such ad hoc methods may produce biased results (Hughes [12]). In this paper, instead of arbitrarily imputing the observations below detection limit, we impute them using fully Bayesian predictive distributions based on a Tobit model (Tobin [13]), which is discussed in Section 2.

Thirdly, another feature of a longitudinal data set is the existence of time-varying covariates which suffer from random measurement errors. This is usually the case in a longitudinal AIDS study where CD4 cell counts are often measured with substantial measurement errors. Thus, any statistical inference without considering measurement errors in covariates may result in biased results (Liu and Wu [2], Wu [4], and Huang and Dagne [5]). In this paper, we jointly model measurement errors in covariate process along with the response process. The distributional assumption for the covariate model is a skew-𝑑 distribution which is relatively robust against potential extreme values and heavy tails.

Our research was motivated by the AIDS clinical trial considered by Acosta et al. [8]. In this study, 44 HIV-1-infected patients were treated with a potent artiretroviral regimen. RNA viral load was measured in copies/mL at study days 0, 7, 14, 28, 56, 84, 112, 140, and 168 of followup. Covariates such as CD4 cell counts were also measured throughout the study on similar scheme. In this study, the viral load detectable limit is 50 copies/mL, and there are 107 out of 357 (30 percent) of all viral load measurements that are below the detection limit. Previous studies show that change in viral load may be associated with change in CD4 cell counts. It is important to study the patterns of virological response to treatment in order to make clinical decisions and provide individualized treatments. Since viral load measurements appear to be skewed and censored, and in addition CD4 cell counts are typically measured with substantial errors and skewness, statistical analyses must take all these factors into account.

For longitudinal data, it is not clear how asymmetric nature, left censoring due to BDL, and covariate measurement error may interact and simultaneously influence inferential procedures. It is the objective of this paper to investigate the effects on inference when all of the three typical features exist in the longitudinal data. To achieve our objective, we employ a fairly general framework to accommodate a large class of problems with various features. Accordingly, we explore a flexible class of skew-elliptical (SE) distributions (see the Appendix for details) which include skew-normal (SN) and skew-𝑑 (ST) distributions as special cases for accounting skewness and heavy tails of longitudinal data, extend the Tobit model (Tobin [13]) to treat all left-censored observations as missing values, and investigate nonparametric mixed effects model for covariate measured with error under the framework of joint models. Because the SN distribution is a special case of the ST distribution when the degrees of freedom approach infinity, for the completeness and convenient presentation, we chose ST distributions to develop NLME Tobit joint models (i.e., the ST distribution is assumed for within-subject random errors and between-subject random effects). The skewness in both within-subject random errors and random-effects distributions may jointly contribute to the skewness of response and covariate variables in a longitudinal study, which makes the assumption of normality unrealistic.

The remaining of the paper is structured as follows. In Section 2, we present the joint models with ST distribution and associated Bayesian modeling approach in general forms so that they can be applicable to other scientific fields. In Section 3, we discuss specific joint models for HIV response process with left censoring and CD4 covariate process with measurement error that are used to illustrate the proposed methods using the data set described above and report the analysis results. Finally, the paper concludes with some discussions in Section 4.

2. Joint Models and Bayesian Inferential Methods

2.1. Skew-𝑑 Mixed-Effects Tobit Joint Models

In this section, we present the models and methods in general forms so that our methods may be applicable to other areas of research. An approach we present in this paper treats censored values as realizations of a latent (unobserved) continuous variable that has been left-censored. This idea was popularized by Tobin ([13]) and the resulting model is commonly referred to as the Tobit model. Denote the number of subjects by n and the number of measurements on the ith subject by 𝑛𝑖. Let 𝑦𝑖𝑗=𝑦𝑖(𝑑𝑖𝑗) and 𝑧𝑖𝑗=𝑧𝑖(𝑑𝑖𝑗) be observed response and covariate for individual 𝑖 at time 𝑑𝑖𝑗 (𝑖=1,2,…,𝑛;𝑗=1,2,…,𝑛𝑖) and π‘žπ‘–π‘— denote the latent response variable that would be measured if the assay did not have a lower detectible limit 𝜌. In our case the Tobit model can be formulated as𝑦𝑖𝑗=ξ‚»π‘žπ‘–π‘—ifπ‘žπ‘–π‘—>𝜌,missingifπ‘žπ‘–π‘—β‰€πœŒ,(2.1) where 𝜌 is a nonstochastic BDL, which in our example below is equivalent to ln(50). Note that the value of 𝑦𝑖𝑗 is missing when it is less than or equal to 𝜌.

For the response process with left-censoring, we consider the following NLME model with an ST distribution which incorporates possibly mismeasured time-varying covariates 𝑦𝑖𝑗𝑑=𝑔𝑖𝑗,𝐱𝑖𝑗,πœ·π‘–π‘—ξ€Έ+𝑒𝑖𝑗,πžπ‘–π‘–π‘–π‘‘βˆΌπ‘†π‘‡π‘›π‘–,πœˆπ‘’ξ€·πŸŽ,𝜎2πˆπ‘›π‘–ξ€·πœΉ,πš«π‘’π‘–,πœ·ξ€Έξ€Έπ‘–π‘—ξ€·π‘§=πβˆ—π‘–π‘—,𝜷,𝐛𝑖,π›π‘–π‘–π‘–π‘‘βˆΌπ‘†π‘‡π‘ 3,πœˆπ‘ξ€·πŸŽ,Ξ£b𝛿,πš«π‘,ξ€Έξ€Έ(2.2) where 𝐱𝑖𝑗 is an 𝑠1Γ—1 design vector, 𝑔(β‹…) is a linear or nonlinear known function, 𝐝(β‹…) is an 𝑠1-dimensional vector-valued linear function, πœ·π‘–π‘— is an 𝑠1Γ—1 individual-specific time-dependent parameter vector, 𝜷 is an 𝑠2Γ—1 population parameter vector (𝑠2β‰₯𝑠1); in the model (2.2), we assume that the individual-specific parameters πœ·π‘–π‘— depend on the true (but unobservable) covariate π‘§βˆ—π‘–π‘— rather than the observed covariate 𝑧𝑖𝑗, which may be measured with errors, and we discuss a covariate model (2.3) below.

It is noticed that we assume that an 𝑠3Γ—1 vector of random effects 𝐛𝑖=(𝑏𝑖1,…,𝑏𝑖𝑠3)𝑇 (𝑠3≀𝑠1) follows a multivariate ST distribution with the unrestricted covariance matrix πšΊπ‘, the 𝑠3×𝑠3 skewness diagonal matrix Ξ”(πœΉπ‘)=diag(𝛿𝑏1,…,𝛿𝑏𝑠3), and the degree of freedom πœˆπ‘; the model random error πžπ‘–=(𝑒𝑖1,…,𝑒𝑖𝑛𝑖)𝑇 follows a multivariate ST distribution with the unknown scale parameter 𝜎2, the degree of freedom πœˆπ‘’, and the 𝑛𝑖×𝑛𝑖 skewness diagonal matrix Ξ”(πœΉπ‘’π‘–)=diag(𝛿𝑒𝑖1,…,𝛿𝑒𝑖𝑖𝑛), where the 𝑛𝑖×1 skewness parameter vector πœΉπ‘’π‘–=(𝛿𝑒𝑖1,…,𝛿𝑒𝑖𝑖𝑛)𝑇. In particular, if 𝛿𝑒𝑖1=β‹―=𝛿𝑒𝑖𝑖𝑛=𝛿𝑒, then Ξ”(πœΉπ‘’π‘–)=π›Ώπ‘’πˆπ‘›π‘– and πœΉπ‘’π‘–=π›Ώπ‘’πŸπ‘›π‘– with πŸπ‘›π‘–=(1,…,1)𝑇; this indicates that we are interested in skewness of overall data set and is the case to be used in real data analysis in Section 3.

Covariate models have been investigated extensively in the literature (Higgins et al. [1], Liu and Wu [2], Wu [4], and Carroll et al. [14]). However, those models used the normality assumption for random measurement errors. As we pointed out earlier, this assumption lacks robustness against departures from normality and may also lead to misleading results. In this paper, we extend the covariate models by assuming an ST distribution for the random errors. We adopt a flexible empirical nonparametric mixed-effects model with an ST to quantify the covariate process as follows: 𝑧𝑖𝑗𝑑=𝑀𝑖𝑗+β„Žπ‘–ξ€·π‘‘π‘–π‘—ξ€Έ+πœ–π‘–π‘—ξ€·β‰‘π‘§βˆ—π‘–π‘—+πœ–π‘–π‘—ξ€Έππ‘–π‘–π‘–π‘‘βˆΌπ‘†π‘‡π‘›π‘–,πœˆπœ–ξ€·πŸŽ,𝜏2πˆπ‘›π‘–ξ€·πœΉ,πš«πœ–π‘–,ξ€Έξ€Έ(2.3) where 𝑀(𝑑𝑖𝑗) and β„Žπ‘–(𝑑𝑖𝑗) are unknown nonparametric smooth fixed-effects and random effects functions, respectively, and 𝝐𝑖=(πœ–π‘–1,…,πœ–π‘–π‘›π‘–)𝑇 follows a multivariate ST distribution with degrees of freedom πœˆπœ–, the unknown scale parameter 𝜏2, and the 𝑛𝑖×𝑛𝑖 skewness diagonal matrix Ξ”(πœΉπœ–π‘–)=diag(π›Ώπœ–π‘–1,…,π›Ώπœ–π‘–π‘–π‘›) with 𝑛𝑖×1 skewness parameter vector πœΉπœ–π‘–=(π›Ώπœ–π‘–1,…,π›Ώπœ–π‘–π‘–π‘›)𝑇. In particular, if π›Ώπœ–π‘–1=β‹―=π›Ώπœ–π‘–π‘–π‘›ξ=π›Ώπœ–, then Ξ”(πœΉπœ–π‘–)=π›Ώπœ–πˆπ‘›π‘– and πœΉπœ–π‘–=π›Ώπœ–πŸπ‘›π‘–. π‘§βˆ—π‘–π‘—=𝑀(𝑑𝑖𝑗)+β„Žπ‘–(𝑑𝑖𝑗) are the true but unobservable covariate values at time 𝑑𝑖𝑗. The fixed smooth function 𝑀(𝑑) represents population average of the covariate process, while the random smooth function β„Žπ‘–(𝑑) is introduced to incorporate the large interindividual variation in the covariate process. We assume that β„Žπ‘–(𝑑) is the realization of a zero-mean stochastic process.

Nonparametric mixed-effects model (2.3) is more flexible than parametric mixed-effects models. To fit model (2.3), we apply a regression spline method to 𝑀(𝑑) and β„Žπ‘–(𝑑). The working principle is briefly described as follows and more details can be found in the literature (Davidian and Giltinan [11] and Wu and Zhang [15]). The main idea of regression spline is to approximate 𝑀(𝑑) and β„Žπ‘–(𝑑) by using a linear combination of spline basis functions. For instance, 𝑀(𝑑) and β„Žπ‘–(𝑑) can be approximated by a linear combination of basis functions Ψ𝑝(𝑑)={πœ“0(𝑑),πœ“1(𝑑),...,πœ“π‘βˆ’1(𝑑)}𝑇 and Ξ¦π‘ž(𝑑)={πœ™0(𝑑),πœ™1(𝑑),...,πœ™π‘žβˆ’1(𝑑)}𝑇, respectively. That is, 𝑀(𝑑)β‰ˆπ‘€π‘(𝑑)=π‘βˆ’1𝑙=0π›Όπ‘™πœ“π‘™(𝑑)=πšΏπ‘(𝑑)π‘‡πœΆ,β„Žπ‘–(𝑑)β‰ˆβ„Žπ‘–π‘ž(𝑑)=π‘žβˆ’1𝑙=0π‘Žπ‘–π‘™πœ™π‘™(𝑑)=πš½π‘ž(𝑑)π‘‡πšπ‘–,(2.4) where 𝜢=(𝛼0,…,π›Όπ‘βˆ’1)𝑇 is a 𝑝×1 vector of fixed-effects and πšπ‘–=(π‘Žπ‘–0,…,π‘Žπ‘–,π‘žβˆ’1)𝑇 (π‘žβ‰€π‘) is a π‘žΓ—1 vector of random-effects with πšπ‘–π‘–π‘–π‘‘βˆΌπ‘†π‘‡π‘ž,πœˆπ‘Ž(𝟎,πšΊπ‘Ž,Ξ”(πœΉπ‘Ž)) with the unrestricted covariance matrix πšΊπ‘Ž, the skewness diagonal matrix Ξ”(πœΉπ‘Ž)=diag(π›Ώπ‘Ž1,…,π›Ώπ‘Žπ‘ž), and the degrees of freedom πœˆπ‘Ž. Based on the assumption of β„Žπ‘–(𝑑), we can regard πšπ‘– as 𝑖𝑖𝑑 realizations of a zero-mean random vector. For our model, we consider natural cubic spline bases with the percentile-based knots. To select an optimal degree of regression spline and numbers of knots, that is, optimal sizes of 𝑝 and π‘ž, the Akaike information criterion (AIC) or the Bayesian information criterion (BIC) is often applied (Davidian and Giltinan [11] and Wu and Zhang [15]). Replacing 𝑀(𝑑) and β„Žπ‘–(𝑑) by their approximations 𝑀𝑝(𝑑) and β„Žπ‘–π‘ž(𝑑), we can approximate model (2.3) by the following linear mixed-effects (LME) model:π‘§π‘–π‘—β‰ˆπšΏπ‘ξ€·π‘‘π‘–π‘—ξ€Έπ‘‡πœΆ+πš½π‘žξ€·π‘‘π‘–π‘—ξ€Έπ‘‡πšπ‘–+πœ–π‘–π‘—β‰ˆπ‘§βˆ—π‘–π‘—+πœ–π‘–π‘—,ππ‘–π‘–π‘–π‘‘βˆΌπ‘†π‘‡π‘›π‘–,πœˆπœ–ξ€·πŸŽ,𝜏2πˆπ‘›π‘–ξ€·πœΉ,πš«πœ–π‘–.ξ€Έξ€Έ(2.5)

2.2. Simultaneous Bayesian Inference

In a longitudinal study, such as the AIDS study described previously, the longitudinal response and covariate processes are usually connected physically or biologically. Statistical inference based on the commonly used two-step method may be undesirable since it fails to take the covariate estimation into account (Higgins et al. [1]). Although a simultaneous inference method based on a joint likelihood for the covariate and response data may be favorable, the computation associated with the joint likelihood inference in joint models of longitudinal data can be extremely intensive and may lead to convergence problems and in some cases it can even be computationally infeasible (Liu and Wu [2] and Wu [4]). Here we propose a simultaneous Bayesian inference method based on MCMC procedure for longitudinal data of response with left censoring and covariate with measurement error. The Bayesian joint modeling approach may pave a way to alleviate the computational burdens and to overcome convergence problems.

We assume that πšπ‘–, 𝐛𝑖, 𝝐𝑖, and πžπ‘– are mutually independent of each other. Following Sahu et al. [7] and properties of ST distribution, in order to specify the models (2.5) and (2.2) for MCMC computation, it can be shown that by introducing four random variable vectors 𝐰𝑒𝑖=(𝑀𝑒𝑖1,…,𝑀𝑒𝑖𝑖𝑛)𝑇,π°πœ–π‘–=(π‘€πœ–π‘–1,…,π‘€πœ–π‘–π‘–π‘›)𝑇,𝐰𝑏𝑖=(𝑀𝑏𝑖1,…,𝑀𝑏3𝑖𝑠)𝑇 and π°π‘Žπ‘–=(π‘€π‘Žπ‘–1,…,π‘€π‘Žπ‘–π‘ž)𝑇 and four random variables πœ‰π‘’π‘–, πœ‰πœ–π‘–, πœ‰π‘π‘–, and πœ‰π‘Žπ‘– (𝑖=1,…,𝑛) based on the stochastic representation for the ST distribution (see the Appendix for details), 𝑧𝑖𝑗 and 𝑦𝑖𝑗 can be hierarchically formulated as π‘¦π‘–π‘—βˆ£π›π‘–,𝑀𝑒𝑖𝑗,πœ‰π‘’π‘–;𝜷,𝜎2,π›Ώπ‘’π‘–π‘—ξ‚€π‘”ξ€·π‘‘βˆΌπ‘π‘–π‘—,𝐱𝑖𝑗𝑧,πβˆ—π‘–π‘—,𝜷,𝐛𝑖+𝛿𝑒𝑖𝑗𝑀𝑒𝑖𝑗,πœ‰π‘’βˆ’1π‘–πœŽ2,π‘€π‘’π‘–π‘—ξ‚€π‘€βˆΌπ‘(0,1)𝐼𝑒𝑖𝑗>0,πœ‰π‘’π‘–βˆ£πœˆπ‘’ξ‚€πœˆβˆΌπΊπ‘’2,πœˆπ‘’2,π›π‘–βˆ£π°π‘π‘–,πœ‰π‘π‘–;πšΊπ‘,πœΉπ‘βˆΌπ‘π‘ 3ξ‚€πš«ξ€·πœΉπ‘ξ€Έπ°π‘π‘–,πœ‰π‘βˆ’1π‘–πšΊπ‘ξ‚,π°π‘π‘–βˆΌπ‘π‘ 3ξ€·πŸŽ,πˆπ‘ 3𝐼𝐰𝑏𝑖>𝟎,πœ‰π‘π‘–βˆ£πœˆπ‘ξ‚€πœˆβˆΌπΊπ‘2,πœˆπ‘2,π‘§π‘–π‘—βˆ£πšπ‘–,π‘€πœ–π‘–π‘—,πœ‰πœ–π‘–;𝜢,𝜏2,π›Ώπœ–π‘–π‘—ξ‚€π‘§βˆΌπ‘βˆ—π‘–π‘—+π›Ώπœ–π‘–π‘—π‘€πœ–π‘–π‘—,πœ‰πœ–βˆ’1π‘–πœ2,π‘€πœ–π‘–π‘—ξ‚€π‘€βˆΌπ‘(0,1)πΌπœ–π‘–π‘—ξ‚>0,πœ‰πœ–π‘–βˆ£πœˆπœ–ξ‚€πœˆβˆΌπΊπœ–2,πœˆπœ–2,πšπ‘–βˆ£π°π‘Žπ‘–,πœ‰π‘Žπ‘–;πšΊπ‘Ž,πœΉπ‘ŽβˆΌπ‘π‘žξ€·πš«ξ€·πœΉπ‘Žξ€Έπ°π‘Žπ‘–,πœ‰π‘Žβˆ’1π‘–πšΊπ‘Žξ€Έ,π°π‘Žπ‘–βˆΌπ‘π‘žξ€·πŸŽ,πˆπ‘žξ€ΈπΌξ€·π°π‘Žπ‘–ξ€Έ>𝟎,πœ‰π‘Žπ‘–βˆ£πœˆπ‘Žξ‚€πœˆβˆΌπΊπ‘Ž2,πœˆπ‘Ž2,(2.6) where 𝐺(β‹…) is a gamma distribution, 𝐼(𝑀𝑒𝑖𝑗>0) is an indicator function, and π‘€π‘’π‘–π‘—βˆΌπ‘(0,1) truncated in the space 𝑀𝑒𝑖𝑗>0 (standard half-normal distribution); π‘€πœ–π‘–π‘—,π°π‘Žπ‘–, and 𝐰𝑏𝑖 can be defined similarly. π‘§βˆ—π‘–π‘— is viewed as the true but unobservable covariate values at time 𝑑𝑖𝑗. It is noted that, as discussed in the Appendix, the hierarchical model with the ST distribution (2.6) can be reduced to the following three special cases: (i) a model with skew-normal (SN) distribution as πœˆπ‘’,πœˆπœ–,πœˆπ‘,πœˆπ‘Žβ†’βˆž and πœ‰π‘’π‘–, πœ‰πœ–π‘–, πœ‰π‘π‘– and πœ‰π‘Žπ‘–β†’1 with probability 1 (i.e., the four corresponding distributional specifications are omitted in (2.6)); (ii) a model with standard 𝑑-distribution as π›Ώπœ–π‘–π‘—=𝛿𝑒𝑖𝑗=0, πœΉπ‘=πœΉπ‘Ž=𝟎, and thus the four distributional specifications of π‘€πœ–π‘–π‘—, 𝑀𝑒𝑖𝑗, π°π‘Žπ‘–, and 𝐰𝑏𝑖 are omitted in (2.6); (iii) a model with standard normal distribution as πœˆπœ–,πœˆπ‘’,πœˆπ‘Ž,πœˆπ‘β†’βˆž and π›Ώπœ–π‘–π‘—=𝛿𝑒𝑖𝑗=0 and πœΉπ‘=πœΉπ‘Ž=𝟎; in this case, the eight corresponding distributional specifications are omitted in (2.6).

Let 𝜽={𝜢,𝜷,𝜏2,𝜎2,πšΊπ‘Ž,πšΊπ‘,πœˆπœ–,πœˆπ‘’,πœˆπ‘Ž,πœˆπ‘,πœΉπ‘Ž,πœΉπ‘,πœΉπœ–π‘–,πœΉπ‘’π‘–;𝑖=1,…,𝑛} be the collection of unknown parameters in models (2.2) and (2.5). To complete the Bayesian formulation, we need to specify prior distributions for unknown parameters in the models (2.2) and (2.5) as follows:πœΆβˆΌπ‘π‘ξ€·πœΆ0,𝚲1ξ€Έ,𝜏2ξ€·πœ”βˆΌIG1,πœ”2ξ€Έ,πšΊπ‘Žξ€·π›€βˆΌIW1,𝜌1ξ€Έ,πœΉπœ–π‘–βˆΌπ‘π‘›π‘–ξ€·πŸŽ,πšͺ1ξ€Έ,πœ·βˆΌπ‘π‘ 2ξ€·πœ·0,𝚲2ξ€Έ,𝜎2ξ€·πœ”βˆΌIG3,πœ”4ξ€Έ,πšΊπ‘ξ€·π›€βˆΌIW2,𝜌2ξ€Έ,πœΉπ‘’π‘–βˆΌπ‘π‘›π‘–ξ€·πŸŽ,πšͺ2ξ€Έ,πœˆπœ–ξ€·πœˆβˆΌπΊπœ–0,πœˆπœ–1ξ€ΈπΌξ€·πœˆπœ–ξ€Έ>3,πœˆπ‘’ξ€·πœˆβˆΌπΊπ‘’0,πœˆπ‘’1ξ€ΈπΌξ€·πœˆπ‘’ξ€Έ>3,πœˆπ‘Žξ€·πœˆβˆΌπΊπ‘Ž0,πœˆπ‘Ž1ξ€ΈπΌξ€·πœˆπ‘Žξ€Έ,𝜈>3π‘ξ€·πœˆβˆΌπΊπ‘0,πœˆπ‘1ξ€ΈπΌξ€·πœˆπ‘ξ€Έ>3,πœΉπ‘ŽβˆΌπ‘π‘žξ€·πŸŽ,πšͺ3ξ€Έ,πœΉπ‘βˆΌπ‘π‘ 3ξ€·πŸŽ,πšͺ4ξ€Έ,(2.7) where the mutually independent Inverse Gamma (IG), Normal (𝑁), Gamma (𝐺), and Inverse Wishart (IW) prior distributions are chosen to facilitate computations (Pinheiro and Bates [16]). The hyperparameter matrices Ξ›1,Ξ›2,Ξ©1,Ξ©2,Ξ“1,Ξ“2,Ξ“3, and Ξ“4 can be assumed to be diagonal for convenient implementation.

Let 𝑓(β‹…βˆ£β‹…),𝐹(β‹…βˆ£β‹…) and πœ‹(β‹…) denote a probability density function (pdf), cumulative density function (cdf), and prior density function, respectively. Conditional on the random variables and some unknown parameters, a detectable measurement 𝑦𝑖𝑗 contributes 𝑓(π‘¦π‘–π‘—βˆ£π›π‘–,𝑀𝑒𝑖𝑗,𝑒𝑒𝑖), whereas a nondetectable measurement contributes 𝐹(πœŒβˆ£π›π‘–,𝑀𝑒𝑖𝑗,𝑒𝑒𝑖)≑𝑃(𝑦𝑖𝑗<πœŒβˆ£π›π‘–,𝑀𝑒𝑖𝑗,𝑒𝑒𝑖) in the likelihood. We assume that 𝜢,𝜷,𝜏2,𝜎2,πšΊπ‘Ž,πšΊπ‘,πœˆπœ–,πœˆπ‘’,πœΉπœ–π‘–,πœΉπ‘’π‘–(𝑖=1,…,𝑛) are independent of each other, that is, πœ‹(𝜽)=πœ‹(𝜢)πœ‹(𝜷)πœ‹(𝜏2)πœ‹(𝜎2)πœ‹(πšΊπ‘Ž)πœ‹(πšΊπ‘)πœ‹(πœˆπœ–)πœ‹(πœˆπ‘’)πœ‹(πœˆπ‘Ž)πœ‹(πœˆπ‘)πœ‹(πœΉπ‘Ž)πœ‹(πœΉπ‘)βˆπ‘–πœ‹(πœΉπœ–π‘–)πœ‹(πœΉπ‘’π‘–). After we specify the models for the observed data and the prior distributions for the unknown model parameters, we can make statistical inference for the parameters based on their posterior distributions under the Bayesian framework. Letting 𝐲𝑖=(𝑦𝑖1,…,𝑦𝑖𝑛𝑖)𝑇 and 𝐳𝑖=(𝑧𝑖1,…,𝑧𝑖𝑛𝑖)𝑇, the joint posterior density of 𝜽 based on the observed data can be given by 𝑓(𝜽∣data)βˆπ‘›ξ‘π‘–=1πΏξ€œξ€œπ²π‘–πΏπ³π‘–πΏπšπ‘–πΏπ›π‘–π‘‘πšπ‘–π‘‘π›π‘–ξƒ°πœ‹(πœƒ),(2.8) where 𝐿𝐲𝑖=βˆπ‘›π‘–π‘—=1𝑓(π‘¦π‘–π‘—βˆ£π›π‘–,𝑀𝑒𝑖𝑗,πœ‰π‘’π‘–)1βˆ’π‘π‘–π‘—πΉ(πœŒβˆ£π›π‘–,𝑀𝑒𝑖𝑗,πœ‰π‘’π‘–)𝑐𝑖𝑗𝑓(π‘€π‘’π‘–π‘—βˆ£π‘€π‘’π‘–π‘—>0)𝑓(πœ‰π‘’π‘–) is the likelihood for the observed response data, 𝑐𝑖𝑗 is the censoring indicator such that 𝑦𝑖𝑗 is observed if 𝑐𝑖𝑗=0, and 𝑦𝑖𝑗 is left-censored if 𝑐𝑖𝑗=1, that is, 𝑦𝑖𝑗=π‘žπ‘–π‘— if 𝑐𝑖𝑗=0, and 𝑦𝑖𝑗 is treated as missing if 𝑐𝑖𝑗=1, and 𝐿𝐳𝑖=βˆπ‘›π‘–π‘—=1𝑓(π‘§π‘–π‘—βˆ£πšπ‘–,π‘€πœ–π‘–π‘—,πœ‰πœ–π‘–)𝑓(π‘€πœ–π‘–π‘—βˆ£π‘€πœ–π‘–π‘—>0)𝑓(πœ‰πœ–π‘–) is the likelihood for the observed covariate data {𝐳𝑖,𝑖=1,…,𝑛},𝐿𝐛𝑖=𝑓(π›π‘–βˆ£π°π‘π‘–,πœ‰π‘π‘–)𝑓(π°π‘π‘–βˆ£π°π‘π‘–>𝟎)𝑓(πœ‰π‘π‘–), and πΏπšπ‘–=𝑓(πšπ‘–βˆ£π°π‘Žπ‘–,πœ‰π‘Žπ‘–)𝑓(π°π‘Žπ‘–βˆ£π°π‘Žπ‘–>𝟎)𝑓(πœ‰π‘Žπ‘–).

In general, the integrals in (2.8) are of high dimension and do not have closed form solutions. Therefore, it is prohibitive to directly calculate the posterior distribution of 𝜽 based on the observed data. As an alternative, MCMC procedures can be used to sample based on (2.8) using the Gibbs sampler along with the Metropolis-Hasting (M-H) algorithm. An important advantage of the above representations based on the hierarchical models (2.6) and (2.7) is that they can be very easily implemented using the freely available WinBUGS software (Lunn et al. [17]) and that the computational effort is equivalent to the one necessary to fit the normal version of the model. Note that when using WinBUGS to implement our modeling approach, it is not necessary to explicitly specify the full conditional distributions. Thus we omit those here to save space.

3. Data Analysis

3.1. Specification of Models

We now analyze the data set described in Section 1 based on the proposed method. Among the 44 eligible patients, the number of viral load measurements for each patient varies from 4 to 9 measurements. As is evident from Figures 1(c) and 1(d), the interpatient variations in viral load appear to be large and these variations appear to change over time. Previous studies suggest that the interpatient variation in viral load may be partially explained by time-varying CD4 cell count (Wu [4] and Huang et al. [18]).

Models for covariate processes are needed in order to incorporate measurement errors in covariates. CD4 cell counts often have nonnegligible measurement errors, and ignoring these errors can lead to severely misleading results in a statistical inference (Carroll et al. [14]). In A5055 study, roughly 10 per cent of the CD4 measurement times are inconsistent with the viral load measurement times. Consequently, CD4 measurements may be missed at viral load measurement times mainly due to a different CD4 measurement scheme as designed in the study (e.g., CD4 measurements were missed at day 7 as displayed in Figures 1(c) and 1(d)). There seem to be no particular patterns for the missingness. Thus we assume that the missing data in CD4 are missing at random (MAR) in the sense of Rubin [19], so that the missing data mechanism can be ignored in the analysis. With CD4 measures collected over time from the AIDS study, we may model the CD4 process to partially address the measurement errors (Wu [4]). However, the CD4 trajectories are often complicated, and there is no well-established model for the CD4 process. We, thus, model the CD4 process empirically using a nonparametric mixed-effects model, which is flexible and works well for complex longitudinal data. We use linear combinations of natural cubic splines with percentile-based knots to approximate 𝑀(𝑑) and β„Žπ‘–(𝑑). Following the study in (Liu and Wu [2]), we set πœ“0(𝑑)=πœ™0(𝑑)=1 and take the same natural cubic splines in the approximations (2.4) with π‘žβ‰€π‘ (in order to limit the dimension of random-effects). The values of 𝑝 and π‘ž are determined based on the AIC/BIC criteria. The AIC/BIC values are evaluated for various models with (𝑝,π‘ž)={(1,1),(2,1),(2,2),(3,1),(3,2),(3,3)} which was found that the model with (𝑝,π‘ž)=(3,3) has the smallest AIC/BIC values being 703.6/744.4. We thus adopted the following ST nonparametric mixed-effects CD4 covariate model:𝑧𝑖𝑗=𝛼0+π‘Žπ‘–0ξ€Έ+𝛼1+π‘Žπ‘–1ξ€Έπœ“1𝑑𝑖𝑗+𝛼2+π‘Žπ‘–2ξ€Έπœ“2𝑑𝑖𝑗+πœ–π‘–π‘—ξ€·β‰‘π‘§βˆ—π‘–π‘—+πœ–π‘–π‘—ξ€Έ,(3.1) where 𝑧𝑖𝑗 is the observed CD4 value at time 𝑑𝑖𝑗, πœ“1(β‹…) and πœ“2(β‹…) are two basis functions given in Section 2.1 and taking the same natural cubic splines for πœ™(β‹…), 𝜢=(𝛼0,𝛼1,𝛼2)𝑇 is a vector of population parameters (fixed-effects), πšπ‘–=(π‘Žπ‘–0,π‘Žπ‘–1,π‘Žπ‘–2)𝑇 is a vector of random-effects, and 𝝐𝑖=(πœ–π‘–1,…,πœ–π‘–π‘›π‘–)π‘‡βˆΌπ‘†π‘‡π‘›π‘–,πœˆπœ–(𝟎,𝜏2πˆπ‘›π‘–,π›Ώπœ–πˆπ‘›π‘–). In addition, in order to avoid too small or large estimates which may be unstable, we standardize the time-varying covariate CD4 cell counts (each CD4 value is subtracted by mean 375.46 and divided by standard deviation 228.57) and rescale the original time (in days) so that the time scale is between 0 and 1.

For the initial stage of viral decay after treatment, a biologically reasonable viral load model can be formulated by the uniexponential form (Ho et al. [20]), 𝑉(𝑑)=𝑉(0)exp(βˆ’πœ†π‘‘), where 𝑉(𝑑) is the total virus at time 𝑑 and πœ† is the rate of change in viral load. To model the complete viral load trajectory, one possible extension of the model given above is to allow πœ† to vary over time. A simple determinant for time-varying πœ† is the linear function πœ†(𝑑)=π‘Ž+𝑏𝑑. For HIV viral dynamic models, it is typical to take ln-transformation of the viral load in order to stabilize the variance and to speed up estimation algorithm (Ding and Wu [10]). After ln-transformation of 𝑉(𝑑), substituting πœ† by the linear function πœ†(𝑑)=π‘Ž+𝑏𝑑, we obtain the following quadratic linear mixed-effects model:𝑦𝑖𝑗=𝛽𝑖0+𝛽𝑖𝑗1𝑑𝑖𝑗+𝛽𝑖𝑗2𝑑2𝑖𝑗+𝑒𝑖𝑗,(3.2) where 𝑦𝑖𝑗=ln(𝑉𝑖(𝑑𝑖𝑗)), parameter 𝛽𝑖0 represents the initial viral load in ln scale, and parameters 𝛽𝑖𝑗1 and 𝛽𝑖𝑗2 incorporate change in viral decay rate over time, with πœ†π‘–π‘—β‰‘βˆ’(𝛽𝑖𝑗1+𝛽𝑖𝑗2𝑑𝑖𝑗) being the time-varying exponential decay rate. πžπ‘–=(𝑒𝑖1,…,𝑒𝑖𝑛𝑖)π‘‡βˆΌST𝑛𝑖,πœˆπ‘’(𝟎,𝜎2πˆπ‘›π‘–,π›Ώπ‘’πˆπ‘›π‘–);πœ·π‘–π‘—=(𝛽𝑖𝑗0,𝛽𝑖𝑗1,𝛽𝑖𝑗2)𝑇 is a vector of individual parameters for the 𝑖th subject at time 𝑑𝑖𝑗.

Since CD4 cell counts are measured with errors, we assume that the individual-specific and time-varying parameters πœ·π‘–π‘— are related to the summary of true CD4 values π‘§βˆ—π‘–π‘— which may be interpreted as the β€œregularized” CD4 covariate value. As discussed by Wu [21], to determine whether CD4 values influence the dynamic parameters πœ·π‘–π‘—, AIC/BIC criteria are used again as guidance (Pinheiro and Bates [16]) to find the following model 𝛽𝑖0=𝛽1+𝑏𝑖1,𝛽𝑖𝑗1=𝛽2+𝛽3π‘§βˆ—π‘–π‘—+𝑏𝑖2,𝛽𝑖𝑗2=𝛽4+𝛽5π‘§βˆ—π‘–π‘—+𝑏𝑖3,(3.3) where 𝐛𝑖=(𝑏𝑖1,…,𝑏𝑖3)𝑇 is individual random-effect, and 𝜷=(𝛽1,𝛽2,…,𝛽5)𝑇 is a vector of population parameters. The model (3.3) indicates that the current (regularized) CD4 values π‘§βˆ—π‘–π‘— rather than the past (observed) CD4 values 𝑧𝑖𝑗 are most predictive of the change in viral load at time 𝑑𝑖𝑗. One possible explanation is that, since CD4 measurements for each individual are often sparse, the current CD4 value may be the best summary of immediate past CD4 values, while the early CD4 values may not be very predictive of the current change in viral load.

3.2. Model Implementation

In this section, we analyze the AIDS data set described in Section 1 to illustrate the proposed joint modeling method (denoted by JM) based on the joint models (3.2) in conjunction with the covariate model (3.1) and the corresponding specifications of prior distributions. As shown in Figures 1(a) and 1(b), the histograms of viral load in ln scale and CD4 cell count clearly indicate their asymmetric nature and it seems logical to fit the joint model with a skew distribution to the data set. Along with this consideration, the following statistical models with different distributions of both model errors and random-effects for both the response model (3.2) and the covariate model (3.1) are employed to compare their performance.(i)SN Model: πžπ‘–, 𝝐𝑖, 𝐛𝑖, and πšπ‘– follow an SN-distribution.(ii)ST Model: πžπ‘–, 𝝐𝑖, 𝐛𝑖, and πšπ‘– follow an ST-distribution.(iii)𝐍 Model: πžπ‘–, 𝝐𝑖, 𝐛𝑖, and πšπ‘– follow a normal (𝑁) distribution.

We investigate the following three scenarios. First, since a normal distribution is a special case of an SN distribution when skewness parameter is zero, while the ST distribution reduces to the SN distribution when the degree of freedom approaches infinity, we investigate how an asymmetric (SN or ST) distribution contributes to modeling results and parameter estimation in comparison with a symmetric (normal) distribution. Second, we estimate the model parameters by using the β€œnaive” method (denoted by NV), which ignores measurement errors in CD4, and missing responses are imputed by the half (i.e., ln(25)) of the BDL. That is, the β€œnaive” method only uses the observed CD4 values 𝑧𝑖𝑗 rather than true (unobservable) CD4 values π‘§βˆ—π‘–π‘— in the response model (3.2) and the missing data in the Tobit model (2.1) is imputed by ln(25). We use it as a comparison to the JM proposed in Section 2. This comparison attempts to investigate how the measurement errors in CD4 and missing data in viral load together contribute to modeling results. Third, when covariates are measured with errors, a common approach is the so-called two-step (TS) method (Higgins et al. [1]): the first step estimates the β€œtrue” covariate values based on the covariate model (3.1); at the second step the covariate in the response model (2.6) is substituted by the estimate from the first step. Thus we use the two-step (TS) method to assess the performance of the JM method.

The progress in the Bayesian posterior computation due to MCMC procedures has made it possible to fit increasingly complex statistical models (Lunn et al. [17] and Huang et al. [18]). To choose the best model among candidate models, it has become more important to develop efficient model selection criteria. A recent publication by Spiegelhalter et al. [22] suggested a generalization of AIC called deviance information criterion (DIC). Since the structure of DIC allows for an automatic computation in WinBUGS, we use DIC to compare the models in this paper. As with other model selection criteria, we caution that DIC is not intended for identification of the β€œcorrect” model, but rather merely as a method of comparing a collection of alternative formulations. In our models with different distribution specifications for model errors, DIC can be used to find out how assumption of a skew-normal distribution contributes to virologic response in comparison with that of a normal distribution and how the proposed joint modeling approach influences parameter estimation compared with the β€œnaive” method and imputation method.

To carry out the Bayesian inference, we need to specify the values of the hyperparameters in the prior distributions. In the Bayesian approach, we only need to specify the priors at the population level. The values of the hyperparameters were mostly chosen from previous studies in the literature (Liu and Wu [2], Huang and Dagne [5], Sahu et al. [7], Wu [21], and among others). We take weakly informative prior distribution for the parameters in the models. In particular, (i) fixed-effects were taken to be independent normal distribution 𝑁(0,100) for each component of the population parameter vectors 𝜢 and 𝜷. (ii) For the scale parameters 𝜏2 and 𝜎2, we assume a limiting noninformative inverse gamma prior distribution, IG(0.01,0.01) so that the distribution has mean 1 and variance 100. (iii) The priors for the variance-covariance matrices of the random-effects πšΊπ‘Ž and πšΊπ‘ are taken to be inverse Wishart distributions IW(Ξ©1,𝜌1) and IW(Ξ©2,𝜌2) with covariance matrices Ξ©1=Ξ©2=diag(0.01,0.01,0.01) and degrees of freedom 𝜌1=𝜌2=5, respectively. (iv) The degrees of freedom parameters πœˆπœ–, πœˆπ‘’, πœˆπ‘Ž, and πœˆπ‘ follow a truncated gamma distribution with two hyperparameter values being 1 and 0.1, respectively. (v) For each of the skewness parameters 𝛿𝑒, π›Ώπœ–, π›Ώπ‘Žπ‘˜, and π›Ώπ‘π‘˜ (π‘˜=1,2,3), we choose independent normal distribution 𝑁(0,100), where we assume that πœΉπ‘’π‘–=π›Ώπ‘’πŸπ‘›π‘– and πœΉπœ–π‘–=π›Ώπœ–πŸπ‘›π‘– to indicate that we are interested in skewness of overall viral load data and overall CD4 cell count data. The MCMC sampler was implemented using WinBUGS software, and the program codes are available from authors on request. The convergence of MCMC implementation was assessed using standard tools (such as trace plots which are not shown here to save space) within WinBUGS, and convergence was achieved after initial 50,000 burn-in iterations. After convergence diagnostics was done, one long chain of 200,000 iterations, retaining every 20th, was run to obtain 10,000 samples for Bayesian inference. Next, we report analysis results of the three scenarios proposed above.

3.3. Comparison of Joint Modeling Results

The population posterior mean (PM), the corresponding standard deviation (SD), and 95% credible interval for fixed-effects parameters based on the three models (SN, ST, and 𝑁) for JM method are presented in the upper part of Table 1. The significant findings are presented as follows. (i) For the response model (3.2), where the most substantively interesting parameters are (𝛽2,𝛽3,𝛽4,𝛽5), the estimates of 𝛽2 and 𝛽4, the linear coefficient and quadratic coefficient of time, respectively, under the three models, are significant since the 95% credible intervals do not contain zero. Among the coefficients of the true CD4 covariate (𝛽3,𝛽5) in model (3.3), the posterior means of 𝛽5 are significantly different from zero for all the three models under JM method. Moreover, the posterior mean values for 𝛽5 are quite different between models SN (βˆ’4.76), ST (βˆ’6.31), and 𝑁 (βˆ’6.26), implying that the posterior means may be substantially biased if model distribution ignores skewness. We will see later that SN gives better fit than either ST or 𝑁. In addition, for the scale parameter 𝜎2, the posterior mean value (2.63) in 𝑁 model is much larger than that of any other corresponding posterior means in SN and ST models. (ii) For parameter estimates of the CD4 covariate model (3.1), the posterior means of intercept 𝛼0 and coefficient 𝛼1 based on SN and ST models are significant, while the posterior mean of 𝛼2 turns out to be nonsignificant under all the three models. For the scale parameter 𝜏2 of the covariate model, the posterior mean value (0.13) is the largest under 𝑁 model. This is expected since the model based on ordinary normal distribution does not account for skewness and heaviness in tails for the type of data analyzed here.

MethodModel 𝛼 0 𝛼 1 𝛼 2 𝛽 1 𝛽 2 𝛽 3 𝛽 4 𝛽 5 𝜏 2 𝜎 2 DIC

JM SNPM βˆ’0.95 0.15 βˆ’ 0 . 2 3 5 . 6 2 βˆ’ 1 4 . 6 βˆ’ 2 . 3 4 1 1 . 7 βˆ’ 4 . 7 6 0 . 0 7 0 . 1 4 863.0
𝐿 C I βˆ’ 1 . 5 8 0.06 βˆ’ 1 5 . 2 4 . 1 7 βˆ’ 2 2 . 1 βˆ’ 5 . 1 4 4 . 5 2 βˆ’ 9 . 9 2 0 . 0 4 0 . 0 1
π‘ˆ C I βˆ’ 0 . 0 1 0 . 9 0 1 4 . 8 7 . 5 9 βˆ’ 8 . 1 4 1 . 4 4 2 1 . 7 βˆ’ 0 . 6 2 0 . 1 1 0 . 6 4
SD 0 . 4 7 0 . 3 7 7 . 6 3 0 . 9 6 3 . 9 8 1 . 6 5 5 . 2 5 2 . 3 4 0 . 0 2 0 . 1 8
STPM βˆ’ 0 . 9 4 0 . 3 4 βˆ’ 0 . 3 1 5 . 8 4 βˆ’ 1 2 . 0 βˆ’ 1 . 2 0 8 . 1 2 βˆ’ 6 . 3 1 0 . 0 4 0 . 2 1 9 8 5 . 6
𝐿 C I βˆ’ 1 . 4 1 0 . 1 8 βˆ’ 1 4 . 1 4 . 1 5 βˆ’ 1 6 . 5 βˆ’ 5 . 7 2 2 . 2 0 βˆ’ 1 2 . 6 0 . 0 2 0 . 0 1
π‘ˆ C I βˆ’ 0 . 0 6 0 . 8 8 1 3 . 4 8 . 0 2 βˆ’ 7 . 7 2 2 . 7 2 1 9 . 2 βˆ’ 1 . 4 1 0 . 0 5 0 . 8 6
SD 0 . 3 5 0.26 7 . 0 9 1 . 1 0 2 . 2 2 2 . 2 2 4 . 1 4 2 . 7 7 0 . 0 1 0 . 2 6
𝑁 PM βˆ’ 0 . 2 1 0 . 4 5 βˆ’ 2 . 8 7 7 . 7 4 βˆ’ 1 5 . 4 βˆ’ 0 . 8 0 1 3 . 6 βˆ’ 6 . 2 6 0 . 1 3 2 . 6 3 1242.3
𝐿 C I βˆ’ 0 . 4 6 0 . 2 2 βˆ’ 1 5 . 9 7 . 2 0 βˆ’ 1 8 . 3 βˆ’ 4 . 1 6 9 . 9 7 βˆ’ 1 1 . 7 0 . 1 1 2 . 0 6
π‘ˆ C I 0 . 0 4 0 . 6 8 9 . 9 0 8 . 2 9 βˆ’ 1 2 . 6 2 . 5 3 1 7 . 2 βˆ’ 1 . 4 3 0 . 1 6 3 . 3 5
SD 0 . 1 3 0 . 1 2 6 . 5 4 0 . 2 8 1 . 4 8 1 . 7 3 1 . 8 5 2 . 6 1 0 . 0 1 0 . 3 3

NVSNPMβ€”β€”β€” 5 . 0 3 βˆ’ 1 1 . 1 0 . 5 8 6 . 8 3 βˆ’2.10 β€” 0 . 1 0 1083.5
𝐿 C I β€”β€”β€” 3 . 8 2 βˆ’ 1 3 . 6 βˆ’ 0 . 9 4 4 . 5 2 βˆ’4.18 β€” 0 . 0 1
π‘ˆ C I β€”β€”β€” 6 . 5 9 βˆ’ 8 . 7 3 2 . 0 8 9 . 1 8 0 . 0 7 β€” 0 . 3 5
SDβ€”β€”β€” 0 . 7 5 1 . 2 3 0 . 7 7 1 . 1 9 1 . 0 4 β€” 0 . 0 9
TS SN PM βˆ’ 0 . 9 9 0 . 1 9 2 . 7 1 5 . 9 1 βˆ’ 1 4 . 4 βˆ’ 1 . 2 4 8 . 4 7 βˆ’ 5 . 9 0 0 . 0 9 0 . 1 4 1 0 2 3 . 8
𝐿 C I βˆ’ 1 . 5 8 βˆ’ 0 . 4 3 βˆ’ 1 2 . 1 4 . 1 2 βˆ’ 2 2 . 1 βˆ’ 5 . 0 1 1 . 8 3 βˆ’ 1 0 . 6 0 . 0 5 0 . 0 1
π‘ˆ C I 0 . 0 7 0 . 9 0 1 7 . 1 7 . 7 2 βˆ’ 8 . 5 0 2 . 1 6 2 1 . 2 βˆ’ 0 . 8 0 0 . 1 4 0 . 6 5
SD 0 . 4 2 0 . 3 6 7 . 5 4 1 . 0 5 3 . 8 8 1 . 7 9 5 . 1 4 2 . 5 2 0 . 0 2 0 . 1 8

To assess the goodness-of-fit of the proposed JM method, the diagnosis plots for the SN, ST, and 𝑁 models comparing the residuals and the fitted values (Figures 2(a)–2(c)) and the observed values versus the fitted values (Figures 2(d)–2(f)). The distribution of the residuals for SN model looks tighter than those for either ST model or 𝑁 model, showing a better fit. Similar results are observed by looking at the plots in Figures 2(d)–2(f). The plot for SN model has most of the points close the line showing a strong agreement between the observed and the fitted values. Clearly, it can be seen from the plots that 𝑁 model, which ignores skewness, does not fit the data very well as compared to either SN model or ST model. Note that the horizontal line designates the below detection limit (BDL), which is at ln(50). The recorded observations less than BDL are not accurate and, therefore, have not been used in the analysis, but instead they were treated as missing and predicted values are obtained. These predicted values are plotted against the recorded observations below detection limit as shown in the lower-row plots. In general, from the model fitting results, both SN and ST models provide a reasonably good fit to the observed data even though SN model is slightly better than ST model.

In order to further investigate whether SN model under JM method can provide better fit to the data than ST model, the DIC values are obtained and found to be 863.0 for SN model and 985.6 for ST model. The DIC value for SN model is smaller than that of ST model, confirming that SN model is better than ST model in fitting the proposed joint model. As mentioned before, it is hard sometimes to tell which model is β€œcorrect” but which one fits data better. The model which fits the data better may be more appealing in order to describe the mechanism of HIV infection and CD4 changing process. Thus, based on the DIC criterion, the results indicate that SN model is relatively better than either ST model or 𝑁 model. These findings are consistent with those displayed in the goodness-of-fit in Figure 2 indicating that SN model outperforms both ST model and 𝑁 model. In summary, our results suggest that it is very important to assume an SN distribution for the response Tobit model and the CD4 covariate model in order to achieve reliable results, in particular if the data exhibit skewness, but not heaviness in the tails. Along with these observations, next we provide detailed fitting results and interpretations based on the SN Model.

3.4. Estimation Results Based on SN Model

For comparison, we used the β€œnaive” (NV) method to estimate the model parameters presented in the lower part of Table 1 where the raw (observed) CD4 values 𝑧𝑖𝑗 rather than the true (unobserved) CD4 values π‘§βˆ—π‘–π‘— are substituted in the response model (3.3). It can be seen that there are important differences in the posterior means for the parameters 𝛽3 and 𝛽5, which are coefficients of CD4 covariate. These posterior means are ̂𝛽3=0.58 and ̂𝛽5=βˆ’2.10 for the NV method, and ̂𝛽3=βˆ’2.34 and ̂𝛽5=βˆ’4.76 for the JM method. The NV method may produce biased posterior means and may substantially overestimate the covariate CD4 effect. The estimated standard deviations (SD) for the CD4 effect (𝛽3 and 𝛽5) using the JM method are 1.65 and 2.34, which are approximately twice as large as those (0.77 and 1.04) using the NV method, respectively, probably because the JM method incorporates the variation from fitting the CD4 process. The differences of the NV estimates and the JM estimates suggest that the estimated parameters may be substantially biased if measurement errors in CD4 covariate are ignored. We also obtained DIC value of 1083.5 for the NV method, while the DIC value for the JM method is 863.0. We can see from the estimated DIC values that the JM approach provides a better fit to the data in comparison with the NV method. Thus it is important to take the measurement errors into account when covariates are measured with errors.

Comparing the JM method against the two-step (TS) method from the lower part of Table 1, we can see that the TS estimates and the JM estimates are somewhat different. In particular, there are important differences in the posterior means for the parameters 𝛽4 and 𝛽5 which is directly associated CD4 covariate. For the parameter 𝛽5, the posterior means are βˆ’4.76 (95%CI=(βˆ’9.92,βˆ’0.62)) and βˆ’5.90 (95%CI=(βˆ’10.60,βˆ’0.80)) for the JM and TS methods, respectively. The TS method slightly underestimates the effect of CD4 covariate.

The estimated results based on the JM method for SN model in Table 2 presents the estimated skewness parameters, and the only significant skewness parameters are those for the response model errors and CD4 covariate model errors, but not random-effects. These are ̂𝛿𝑒=2.34 (95%CI=(1.93,2.73)) and Μ‚π›Ώπœ–=0.41 (95%CI=(0.25,0.54)) for viral load and CD4 cell count, respectively. They are significantly positive confirming the right-skewed viral load and CD4 cell count as was depicted in Figure 1. Thus, the results suggest that accounting for significant skewness, when the data exhibit skewness, provides a better model fit to the data and gives more accurate estimates to the parameters.

MethodModel 𝛿 πœ– 𝛿 𝑒 𝛿 π‘Ž 1 𝛿 π‘Ž 2 𝛿 π‘Ž 3 𝛿 𝑏 1 𝛿 𝑏 2 𝛿 𝑏 3 𝜈 πœ– 𝜈 𝑒 𝜈 π‘Ž 𝜈 𝑏

JM SNPM 0 . 4 1 2 . 3 4 0 . 5 8 0 . 3 4 0 . 2 6 0 . 5 2 βˆ’ 1 . 8 1 2 . 5 8 β€” β€”β€”β€”
𝐿 C I 0 . 2 5 1 . 9 3 βˆ’ 0 . 6 2 βˆ’ 0 . 5 8 βˆ’ 1 6 . 1 βˆ’ 1 . 9 0 βˆ’ 1 0 . 1 βˆ’ 1 0 . 7 β€” β€”β€”β€”
π‘ˆ C I 0 . 5 4 2 . 7 3 1 . 4 1 1 . 0 9 1 7 . 0 2 . 3 1 7 . 1 2 1 1 . 0 β€” β€”β€”β€”
SD 0 . 0 7 0 . 2 0 0 . 6 0 0 . 4 7 8 . 5 7 1 . 1 8 5 . 3 0 6 . 9 3 β€” β€”β€”β€”
STPM 0 . 0 5 2 . 2 6 0 . 8 7 0 . 0 4 βˆ’ 0 . 5 6 0 . 3 2 βˆ’ 4 . 7 0 6 . 4 5 3 . 3 2 1 0 . 2 1 4 . 0 1 4 . 4
𝐿 C I βˆ’0.14 1 . 5 9 βˆ’ 0 . 3 5 βˆ’ 0 . 6 0 βˆ’ 1 6 . 5 βˆ’ 2 . 2 3 βˆ’ 1 0 . 1 βˆ’ 8 . 2 3 3 . 0 1 3 . 0 7 3 . 5 2 3 . 5 2
π‘ˆ C I 0 . 2 5 2 . 7 0 1 . 5 4 0 . 6 5 1 4 . 7 2 . 3 7 1 . 9 1 1 2 . 3 4 . 1 8 3 5 . 2 4 1 . 1 4 1 . 9
SD 0 . 1 1 0 . 3 3 0 . 5 0 0 . 3 2 8 . 4 1 1 . 3 1 2 . 9 3 4 . 8 8 0 . 3 2 8 . 9 8 1 0 . 3 1 0 . 3

NV SNPM β€” 2 . 2 4 β€” β€”β€” 0 . 8 0 0 . 1 5 5 . 5 3 β€” β€”β€” β€”
𝐿 C I β€” 1 . 9 5 β€” β€”β€” βˆ’ 1 . 0 5 βˆ’1.71 3 . 6 2 β€” β€”β€” β€”
π‘ˆ C I β€” 2 . 5 5 β€” β€”β€” 2 . 2 5 2 . 3 0 7 . 7 4 β€” β€”β€” β€”
SD β€” 0 . 1 5 β€” β€”β€” 0 . 9 2 1 . 0 0 1 . 0 6 β€” β€”β€”β€”
TS SNPM 0 . 1 6 2 . 4 4 0 . 8 9 0 . 2 8 3 . 0 6 0 . 0 4 βˆ’ 0 . 9 4 5 . 1 8 β€” β€”β€” β€”
𝐿 C I βˆ’ 0 . 3 9 2 . 0 7 βˆ’0.48 βˆ’ 0 . 5 8 βˆ’ 1 1 . 7 βˆ’ 2 . 2 0 βˆ’8.53 βˆ’ 1 2 . 4 β€” β€”β€” β€”
π‘ˆ C I 0 . 5 1 2 . 7 9 1 . 5 5 1 . 0 4 2 1 . 0 2 . 2 3 7 . 4 9 1 2 . 2 β€” β€”β€” β€”
SD 0 . 2 9 0 . 1 8 0 . 5 0 0 . 4 5 8 . 3 0 1 . 3 1 4 . 7 5 6 . 6 7 β€” β€”β€” β€”

In summary, the results indicate that the SN model under the JM method is a better suited model for viral loads and CD4 covariate with measurement errors. Looking now at the estimated population initial stage of viral decay after treatment bases on the JM method, we get Μ‚πœ†(𝑑)=βˆ’(βˆ’14.6βˆ’2.34π‘§βˆ—(𝑑)+11.7π‘‘βˆ’4.76π‘§βˆ—(𝑑)𝑑), where π‘§βˆ—(𝑑) is the standardized true CD4 value at time 𝑑 which may be interpreted as the β€œregularized” covariate value. Thus, the population viral load process may be approximated by ̂𝑉(𝑑)=exp[5.62βˆ’πœ†(𝑑)𝑑]. Since the viral decay rate (πœ†(𝑑)) is significantly associated with the true CD4 values (due to statistically significant estimate of 𝛽5), this suggests that the viral load change 𝑉(𝑑) may be significantly associated with the true CD4 process. Note that, although the true association described above may be complicated, the simple approximation considered here may provide a reasonable guidance and suggest a further research.

4. Discussion

Attempts to jointly fit the viral load data and CD4 cell counts with measurement errors are compromised by left censoring in viral load response due to detection limits. We addressed this problem using Bayesian nonlinear mixed-effects Tobit models with skew distributions. The models were fitted based on the assumption that the viral dynamic model (2.2) continues to hold for those unobserved left-censored viral loads. This assumption may be reasonable since the dynamic model considered here is a natural extension of a biologically justified model (Ding and Wu [10]). Even though left censoring effects are the focus of this paper, right-censoring (ceiling) effects can also be dealt with in very similar ways. It is therefore important for researchers to pay attention to censoring effects in a longitudinal data analysis, and Bayesian Tobit models with skew distributions make best use of both censored and uncensored data information.

Our results suggest that both ST (skew-𝑑) and SN (skew-normal) models show superiority to the 𝑁 (normal) model. Our results also indicate that the JM method outperformed the NV and TS methods in the sense that it produces more accurate parameter estimates. The JM method is quite general and so can be applied to other application areas, allowing accurate inferences of parameters while adjusting for skewness, left-censoring, and measurement errors. In short, skew distributions show potentials to gain efficiency and accuracy in estimating certain parameters when the normality assumption does not hold in the data.

The proposed NLME Tobit joint model with skew distributions can be easily fitted using MCMC procedure by using the WinBUGS package that is available publicly and has a computational cost similar to the normal version of the model due to the features of its hierarchically stochastic representations. Implementation via MCMC makes it straightforward to compare the proposed models and methods with various scenarios for real data analysis in comparison with symmetric distributions and asymmetric distributions for model errors. This makes our approach quite powerful and also accessible to practicing statisticians in the fields. In order to examine the sensitivity of parameter estimates to the prior distributions and initial values, we also conducted a limited sensitivity analysis using different values of hyperparameters of prior distributions and different initial values (data not shown). The results of the sensitivity analysis showed that the estimated dynamic parameters were not sensitive to changes of both priors and initial values. Thus, the final results are reasonable and robust, and the conclusions of our analysis remain unchanged (see Huang et al. [18] for more details).

The methods of this paper may be extended to accommodate various subpopulations of patients whose viral decay trajectories after treatment may differ. In addition, the purpose of this paper is to demonstrate the proposed models and methods with various scenarios for real data analysis for comparing asymmetric distributions for model errors to a symmetric distribution, although a limited simulation study might have been conducted to evaluate our results from different model specifications and the corresponding methods. However, since this paper investigated many different scenarios-based models and methods with real data analysis, the complex natures considered, especially skew distributions involved, will pose some challenges for such a simulation study which requires additional efforts, and it is beyond the purpose of this paper. We are currently investigating these related problems and will report the findings in the near future.


A. Multivariate Skew-𝑑 and Skew Normal Distributions

Different versions of the multivariate skew-elliptical (SE) distributions have been proposed and used in the literature (Sahu et al. [7], Azzalini and Capitanio [23], Jara et al. [24], Arellano-Valle et al. [25], and among others). We adopt a class of multivariate SE distributions proposed by Sahu et al. [7], which is obtained by using transformation and conditioning and contains multivariate skew-𝑑 (ST) and skew-normal (SN) distributions as special cases. A π‘˜-dimensional random vector 𝐘 follows a π‘˜-variate SE distribution if its probability density function (pdf) is given by π‘“ξ‚€π²βˆ£π,𝚺,𝚫(𝜹);π‘šπœˆ(π‘˜)=2π‘˜π‘“ξ‚€π²βˆ£π,𝐀;π‘šπœˆ(π‘˜)𝑃(𝐕>𝟎),(A.1) where 𝐀=𝚺+Ξ”2(𝜹), 𝝁 is a location parameter vector, 𝚺 is a covariance matrix, Ξ”(𝜹) is a skewness diagonal matrix with the skewness parameter vector 𝜹=(𝛿1,𝛿2,…,π›Ώπ‘˜)𝑇; 𝐕 follows the elliptical distribution 𝐸𝑙(Ξ”(𝜹)π΄βˆ’1(π²βˆ’π),πˆπ‘˜βˆ’Ξ”(𝜹)π΄βˆ’1Ξ”(𝜹);π‘šπœˆ(π‘˜)) and the density generator function π‘šπœˆ(π‘˜)(𝑒)=(Ξ“(π‘˜/2)/πœ‹π‘˜/2)(π‘šπœˆβˆ«(𝑒)/∞0π‘Ÿπ‘˜/2βˆ’1π‘šπœˆ(𝑒)π‘‘π‘Ÿ), with π‘šπœˆ(𝑒) being a function such that ∫∞0π‘Ÿπ‘˜/2βˆ’1π‘šπœˆ(𝑒)π‘‘π‘Ÿ exists. The function π‘šπœˆ(𝑒) provides the kernel of the original elliptical density and may depend on the parameter 𝜈. We denote this SE distribution by SE(𝝁,𝚺,Ξ”(𝜹);π‘š(π‘˜)). Two examples of π‘šπœˆ(𝑒), leading to important special cases used throughout the paper, are π‘šπœˆ(𝑒)=exp(βˆ’π‘’/2) and π‘šπœˆ(𝑒)=(𝑒/𝜈)βˆ’(𝜈+π‘˜)/2, where 𝜈>0. These two expressions lead to the multivariate ST and SN distributions, respectively. In the latter case, 𝜈 corresponds to the degree of freedom parameter.

Since the SN distribution is a special case of the ST distribution when the degree of freedom approaches infinity, for completeness, this section is started by discussing the multivariate ST distribution that will be used in defining the ST joint models considered in this paper. For detailed discussions on properties and differences among various versions of ST and SN distributions, see the references above. We consider a multivariate ST distribution introduced by Sahu et al. [7], which is suitable for a Bayesian inference since it is built using conditional method and is defined below.

An π‘˜-dimensional random vector 𝐘 follows an π‘˜-variate ST distribution if its pdf is given by 𝑓(𝐲∣𝝁,𝚺,𝚫(𝜹),𝜈)=2π‘˜π‘‘π‘˜,𝜈(𝐲∣𝝁,𝐀)𝑃(𝐕>𝟎).(A.2) We denote the π‘˜-variate 𝑑 distribution with parameters 𝝁,𝐀 and degrees of freedom 𝜈 by π‘‘π‘˜,𝜈(𝝁,𝐀) and the corresponding pdf by π‘‘π‘˜,𝜈(𝐲∣𝝁,𝐀) henceforth, 𝐕 follows the 𝑑 distribution π‘‘π‘˜,𝜈+π‘˜. We denote this distribution by STπ‘˜,𝜈(𝝁,𝚺,Ξ”(𝜹)). In particular, when 𝚺=𝜎2πˆπ‘˜ and Ξ”(𝜹)=π›Ώπˆπ‘˜, (A.2) simplifies to π‘“ξ€·π²βˆ£π,𝜎2ξ€Έ,𝛿,𝜈=2π‘˜ξ€·πœŽ2+𝛿2ξ€Έβˆ’π‘˜/2Ξ“((𝜈+π‘š)/2)Ξ“(𝜈/2)(πœˆπœ‹)π‘˜/2ξƒ―1+(π²βˆ’π)𝑇(π²βˆ’π)πœˆξ€·πœŽ2+𝛿2ξ€Έξƒ°βˆ’(𝜈+π‘˜)/2Γ—π‘‡π‘˜,𝜈+π‘˜βŽ‘βŽ’βŽ’βŽ£ξƒ―ξ€·πœŽπœˆ+2+𝛿2ξ€Έβˆ’1(π²βˆ’π)𝑇(π²βˆ’π)ξƒ°πœˆ+π‘˜βˆ’1/2𝛿(π²βˆ’π)𝜎√𝜎2+𝛿2⎀βŽ₯βŽ₯⎦,(A.3) where π‘‡π‘˜,𝜈+π‘˜(β‹…) denotes the cumulative distribution function (cdf) of π‘‘π‘˜,𝜈+π‘˜(𝟎,πˆπ‘˜). However, unlike in the SN distribution to be discussed below, the ST density cannot be written as the product of univariate ST densities. Here 𝐘 are dependent but uncorrelated.

The mean and covariance matrix of the ST distribution π‘†π‘‡π‘˜,𝜈(𝝁,𝜎2πˆπ‘˜,Ξ”(𝜹)) are given by ξ‚€πœˆπΈ(𝐘)=𝝁+πœ‹ξ‚1/2Ξ“((πœˆβˆ’1)/2)ξ€ΊπœŽΞ“(𝜈/2)𝜹,cov(𝐘)=2πˆπ‘˜+𝚫2ξ€»πœˆ(𝜹)βˆ’πœˆπœˆβˆ’2πœ‹ξ‚ΈΞ“{(πœˆβˆ’1)/2}ξ‚ΉΞ“(𝜈/2)2𝚫2(𝛿).(A.4)

According to Lemma 1 of Azzalini and Capitanio [23], if 𝐘 follows STπ‘˜,𝜈(𝝁,𝚺,Ξ”(𝜹)), it can be represented by 𝐘=𝝁+πœ‰βˆ’1/2𝐗,(A.5) where πœ‰ follows a Gamma distribution Ξ“(𝜈/2,𝜈/2), which is independent of 𝐗, and 𝐗 follows a π‘˜-dimensional skew-normal (SN) distribution, denoted by SNπ‘˜(𝟎,𝚺,Ξ”(𝜹)). It follows from (A.5) that π˜βˆ£πœ‰βˆΌSNπ‘˜(𝝁,𝚺/πœ‰,Ξ”(𝜹)). By Proposition 1 of Arellano-Valle et al. [25], the SN distribution of 𝐘 conditional on πœ‰ has a convenient stochastic representation as follows: ||π‘‹π˜=𝝁+𝚫(𝜹)0||+πœ‰βˆ’1/2𝚺1/2𝑋1,(A.6) where 𝑋0 and 𝑋1 are two independent π‘π‘˜(𝟎,πˆπ‘˜) random vectors. Note that the expression (A.6) provides a convenience device for random number generation and for implementation purpose. Let 𝐰=|𝑋0|; then 𝐰 follows a π‘˜-dimensional standard normal distribution π‘π‘˜(𝟎,πˆπ‘š) truncated in the space𝐰>𝟎 (i.e., the standard half-normal distribution). Thus, following Sahu et al. [7], a hierarchical representation of (A.6) is given by 𝐘∣𝐰,πœ‰βˆΌπ‘π‘˜ξ€·π+𝚫(𝜹)𝐰,πœ‰βˆ’1πšΊξ€Έ,π°βˆΌπ‘π‘˜ξ€·πŸŽ,πˆπ‘˜ξ€Έξ‚€πœˆπˆ(𝐰>𝟎),πœ‰βˆΌπΊ2,𝜈2,(A.7) where 𝐺(β‹…) is a gamma distribution. Note that the ST distribution presented in (A.7) can be reduced to the following three special cases: (i) as πœˆβ†’βˆž and πœ‰β†’1 with probability 1 (i.e., the last distributional specification is omitted), then the hierarchical expression (A.7) becomes an SN distribution SNπ‘˜(𝝁,𝚺,Ξ”(𝜹)); (ii) as Ξ”(𝜹)=𝟎, then the hierarchical expression (A.7) is a standard multivariate 𝑑-distribution; (iii) as πœˆβ†’βˆž, πœ‰β†’1 with probability 1, and Ξ”(𝜹)=𝟎, then the hierarchical expression (A.7) is a standard multivariate normal distribution.

Specifically, if a π‘˜-dimensional random vector 𝐘 follows a π‘˜-variate SN distribution, then (A.2)–(A.4) revert to the following forms, respectively: 𝑓(𝐲∣𝝁,𝚺,𝚫(𝜹))=2π‘˜||𝐀||βˆ’1/2πœ™π‘˜ξ€½π€βˆ’1/2ξ€Ύ(π²βˆ’π)𝑃(𝐕>𝟎),(A.8) where π•βˆΌπ‘π‘˜{Ξ”(𝜹)π€βˆ’1(π²βˆ’π),πˆπ‘˜βˆ’Ξ”(𝜹)π€βˆ’1Ξ”(𝜹)}, and πœ™π‘˜(β‹…) is the pdf of π‘π‘˜(𝟎,πˆπ‘˜). We denote the above distribution by SNπ‘˜(𝝁,𝚺,Ξ”(𝜹)). An appealing feature of (A.8) is that it gives independent marginal when 𝚺=diag(𝜎21,𝜎22,…,𝜎2π‘˜). The pdf (A.8) thus simplifies to 𝑓(𝐲∣𝝁,𝚺,𝚫(𝜹))=π‘˜ξ‘π‘–=1⎑⎒⎒⎒⎣2ξ”πœŽ2𝑖+𝛿2π‘–πœ™βŽ§βŽͺ⎨βŽͺβŽ©π‘¦π‘–βˆ’πœ‡π‘–ξ”πœŽ2𝑖+𝛿2π‘–βŽ«βŽͺ⎬βŽͺ⎭Φ⎧βŽͺ⎨βŽͺβŽ©π›Ώπ‘–πœŽπ‘–π‘¦π‘–βˆ’πœ‡π‘–ξ”πœŽ2𝑖+𝛿2π‘–βŽ«βŽͺ⎬βŽͺ⎭⎀βŽ₯βŽ₯βŽ₯⎦,(A.9) where πœ™(β‹…) and Ξ¦(β‹…) are the pdf and cdf of the standard normal distribution, respectively. The mean and covariance matrix are given by √𝐸(𝐘)=𝝁+2/πœ‹πœΉ,cov(𝐘)=𝚺+(1βˆ’2/πœ‹)Ξ”(𝜹)2. It is noted that when 𝜹=𝟎, the SN distribution reduces to usual normal distribution.


The authors are grateful to the Guest Editor and three reviewers for their insightful comments and suggestions that led to a marked improvement of the paper. They gratefully acknowledge A5055 study investigators for allowing them to use the clinical data from their study. This research was partially supported by NIAID/NIH Grant R03 AI080338 and MSP/NSA Grant H98230-09-1-0053 to Y. Huang.


  1. M. Higgins, M. Davidian, and D. M. Gilitinan, β€œA two-step approach to measurement error in time-dependent covariates in nonlinear mixed-effects models, with application to IGF-I pharmacokinetics,” Journal of the American Statistical Association, vol. 92, pp. 436–448, 1997. View at: Google Scholar
  2. W. Liu and L. Wu, β€œSimultaneous inference for semiparametric nonlinear mixed-effects models with covariate measurement errors and missing responses,” Biometrics, vol. 63, no. 2, pp. 342–350, 2007. View at: Publisher Site | Google Scholar
  3. M. S. Wulfsohn and A. A. Tsiatis, β€œA joint model for survival and longitudinal data measured with error,” Biometrics, vol. 53, no. 1, pp. 330–339, 1997. View at: Publisher Site | Google Scholar | Zentralblatt MATH
  4. L. Wu, β€œA joint model for nonlinear mixed-effects models with censoring and covariants measured with error, with application to AIDS studies,” Journal of the American Statistical Association, vol. 97, no. 460, pp. 955–964, 2002. View at: Publisher Site | Google Scholar | Zentralblatt MATH
  5. Y. Huang and G. Dagne, β€œA Bayesian approach to joint mixed-effects models with a skew-normal distribution and measurement errors in covariates,” Biometrics, vol. 67, no. 1, pp. 260–269, 2011. View at: Publisher Site | Google Scholar
  6. G. Verbeke and E. Lesaffre, β€œA linear mixed-effects model with heterogeneity in the random-effects population,” Journal of the American Statistical Association, vol. 91, no. 433, pp. 217–221, 1996. View at: Google Scholar
  7. S. K. Sahu, D. K. Dey, and M. D. Branco, β€œA new class of multivariate skew distributions with applications to Bayesian regression models,” The Canadian Journal of Statistics, vol. 31, no. 2, pp. 129–150, 2003. View at: Publisher Site | Google Scholar | Zentralblatt MATH
  8. E. P. Acosta, H. Wu, S. M. Hammer et al., β€œComparison of two indinavir/ritonavir regimens in the treatment of HIV-infected individuals,” Journal of Acquired Immune Deficiency Syndromes, vol. 37, no. 3, pp. 1358–1366, 2004. View at: Publisher Site | Google Scholar
  9. A. S. Perelson, P. Essunger, Y. Cao et al., β€œDecay characteristics of HIV-1-infected compartments during combination therapy,” Nature, vol. 387, no. 6629, pp. 188–191, 1997. View at: Publisher Site | Google Scholar
  10. A. A. Ding and H. Wu, β€œRelationships between antiviral treatment effects and biphasic viral decay rates in modeling HIV dynamics,” Mathematical Biosciences, vol. 160, no. 1, pp. 63–82, 1999. View at: Publisher Site | Google Scholar
  11. M. Davidian and D. M. Giltinan, Nonlinear Models for Repeated Measurement Data, Chapman & Hall, London, UK, 1995.
  12. J. P. Hughes, β€œMixed effects models with censored data with application to HIV RNA levels,” Biometrics, vol. 55, no. 2, pp. 625–629, 1999. View at: Google Scholar
  13. J. Tobin, β€œEstimation of relationships for limited dependent variables,” Econometrica, vol. 26, pp. 24–36, 1958. View at: Google Scholar | Zentralblatt MATH
  14. R. J. Carroll, D. Ruppert, L. A. Stefanski, and C. M. Crainiceanu, Measurement Error in Nonlinear Models: A Modern Perspective, Chapman & Hall, Boca Raton, Fla, USA, 2nd edition, 2006. View at: Publisher Site
  15. H. Wu and J. T. Zhang, β€œThe study of long-term HIV dynamics using semi-parametric non-linear mixed-effects models,” Statistics in Medicine, vol. 21, no. 23, pp. 3655–3675, 2002. View at: Publisher Site | Google Scholar
  16. J. C. Pinheiro and M. D. Bates, Mixed-Effects Models in S and S-PLUS, Springer, New York, NY, USA, 2000.
  17. D. J. Lunn, A. Thomas, N. Best, and D. Spiegelhalter, β€œWinBUGSβ€”a Bayesian modelling framework: concepts, structure, and extensibility,” Statistics and Computing, vol. 10, no. 4, pp. 325–337, 2000. View at: Google Scholar
  18. Y. Huang, D. Liu, and H. Wu, β€œHierarchical Bayesian methods for estimation of parameters in a longitudinal HIV dynamic system,” Biometrics, vol. 62, no. 2, pp. 413–423, 2006. View at: Publisher Site | Google Scholar | Zentralblatt MATH
  19. D. B. Rubin, β€œInference and missing data,” Biometrika, vol. 63, no. 3, pp. 581–592, 1976. View at: Publisher Site | Google Scholar | Zentralblatt MATH
  20. D. D. Ho, A. U. Neumann, A. S. Perelson, W. Chen, J. M. Leonard, and M. Markowitz, β€œRapid turnover of plasma virions and CD4 lymphocytes in HIV-1 infection,” Nature, vol. 373, no. 6510, pp. 123–126, 1995. View at: Publisher Site | Google Scholar
  21. L. Wu, β€œSimultaneous inference for longitudinal data with detection limits and covariates measured with errors, with application to AIDS studies,” Statistics in Medicine, vol. 23, no. 11, pp. 1715–1731, 2004. View at: Publisher Site | Google Scholar
  22. D. J. Spiegelhalter, N. G. Best, B. P. Carlin, and A. Van der Linde, β€œBayesian measures of model complexity and fit,” Journal of the Royal Statistical Society. Series B, vol. 64, no. 4, pp. 583–639, 2002. View at: Publisher Site | Google Scholar | Zentralblatt MATH
  23. A. Azzalini and A. Capitanio, β€œDistributions generated by perturbation of symmetry with emphasis on a multivariate skew t-distribution,” Journal of the Royal Statistical Society. Series B, vol. 65, no. 2, pp. 367–389, 2003. View at: Publisher Site | Google Scholar
  24. A. Jara, F. Quintana, and E. San MartΓ­n, β€œLinear mixed models with skew-elliptical distributions: a Bayesian approach,” Computational Statistics & Data Analysis, vol. 52, no. 11, pp. 5033–5045, 2008. View at: Publisher Site | Google Scholar
  25. R. B. Arellano-Valle, H. Bolfarine, and V. H. Lachos, β€œBayesian inference for skew-normal linear mixed models,” Journal of Applied Statistics, vol. 34, no. 5-6, pp. 663–682, 2007. View at: Publisher Site | Google Scholar

Copyright Β© 2012 Getachew A. Dagne and Yangxin Huang. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

More related articles

2078Β Views | 753Β Downloads | 5Β Citations
 PDF  Download Citation  Citation
 Download other formatsMore
 Order printed copiesOrder

Related articles

We are committed to sharing findings related to COVID-19 as quickly and safely as possible. Any author submitting a COVID-19 paper should notify us at to ensure their research is fast-tracked and made available on a preprint server as soon as possible. We will be providing unlimited waivers of publication charges for accepted articles related to COVID-19. Sign up here as a reviewer to help fast-track new submissions.