Bayesian varying coefficient mixed-effects joint models with asymmetry and missingness

Abstract

Abstract:

Longitudinal and survival data are often collected from clinical studies. Mixed-effects joint models are commonly used for the analysis of such data. Nevertheless, the following issues may arise in longitudinal survival data analysis: (a) most joint models assume a simple parametric mixed-effects model for longitudinal outcome, which may obscure the important relationship between response and covariates; (b) clinical data often exhibits asymmetry so that symmetric assumption for model errors may lead to biased estimation of parameters; (c) response may be missing and missingness may be informative. There is little work concerning all of these issues simultaneously. We develop a Bayesian varying coefficient mixed-effects joint model with skewness and missingness to study the simultaneous influence of these features. The proposed methods are applied to an AIDS clinical data. Simulation studies are conducted to assess the performance of the method.

Keywords

Bayesian inference Competing risks longitudinal data varying coefficient mixed-effects models missing data proportional hazard models asymmetric distribution survival data

1 Introduction

Joint analysis of longitudinal and survival data has received a great deal of attention in clinical trial studies. Joint modelling and inferential methods have been proposed for analyzing such type of data under various scenarios (Faucett and Thomas, 1996; Song et al., 2002; Brown and Ibrahim, 2003; Elashoff et al., 2008; Hu et al., 2009; Huang et al., 2011; Li et al., 2012). Most of these joint models assume a mixed-effects model for longitudinal process and a proportional hazard (PH) model for competing risks process. There is little work on analyzing the simultaneous impacts induced by asymmetry and non-ignorable missingness which are commonly observed in practice.

In AIDS studies, it is interesting to study the relationship between two important biomarkers: HIV viral load and CD4 counts. Figure 1 illustrates the correlation between viral load and CD4 counts at different time points from a Multicenter AIDS Cohort Study (MACS) (Section 2.1 has the details of this study). It is clear that the relationship between viral load and CD4 counts is varying across the study period. Therefore, a varying coefficient mixed-effects model is more appropriate than a parametric model for the purpose of investigating the changing relationship.

A common assumption for most models is that model errors are symmetrically distributed, which may lack robustness against departures from normality and outliers. Misleading results may occur by making inferences based on symmetric assumptions (Sahu et al., 2003; Ho and Lin, 2010; Lachos et al., 2011). Although a transformation such as logarithm or square root may be adopted for the skewed data which makes feasible the analysis with symmetric assumption, the interpretation of results may be difficult. Sometime, even transformation does not make the validity of symmetric distribution under transformed scale. Figure 2 displays the distribution of residuals based on OLS estimates of the subject-specific regression from 1598 subjects in the MACS study. It is noted that the model errors are highly skewed. Thus, a normality assumption is not quite realistic for this dataset. To model skewed data, it is more appropriate to adopt asymmetric distributions such as skew-t (ST) and skew-normal (SN) distributions. The validity of any statistical inference is based on the assumption that variables that are intended to be measured are actually measured without missingness. However, due to many factors, data are never perfectly collected. Although clinical studies are designed to collect data from each patient at each assessment time, missing observations in the response are often encountered due to many reasons. Missingness may be informative or non-ignorable if it is related to missing values. Without accounting for missingness in analysis may make the inferential results questionable. Many likelihood-based methods have been proposed for handling non-ignorably missing responses. Approaches using selection models include Baker (1995), Fitzmaurice et al. (1996), Molenberghs et al. (1997) and Ibrahim et al. (2001), among others. Approaches based on pattern-mixture models have been taken by Ekholm and Skinner (1998), Fitzmaurice and Laird (2000) and Birmingham and Fitzmaurice (2002), among others. Follmann and Wu (1995) propose a class of shared parameter models for non-ignorable non-response that can be specified as a random effects model for the primary response, combined with a model for the missingness in which the random effects are treated as covariates. They condition on the data that describes missingness and use a conditional model for inference. Rotnitzky et al. (1998) develop a class of estimators for generalized linear models (GLMs) with non-ignorable missing responses that are based on inverse probability weighted estimating equations.

It is unclear how asymmetry and missing observations which are commonly seen in longitudinal-survival data may interact and influence the inferential procedure, especially under the non-parametric framework. We propose a varying coefficient mixed-effects joint model to account for longitudinal survival data with asymmetry and non-ignorable missingness. We develop a Bayesian approach to simultaneously estimate parameters in the joint model based on the ST distribution. Note that normal (N) and skew-normal (SN) distributions are special cases of ST distribution. When degrees of freedom approach infinity, ST distribution is approximate to SN distribution. If the skewness parameter is zero, SN distribution degenerates to normal distribution. Therefore, we built the joint model using ST distributions since it includes important special cases such as SN and normal distributions. We consider multivariate ST and SN distributions (Sahu et al., 2003), which are convenient for Bayesian inference. Appendix 1 briefly discusses these distributions.

Figure 1:

Illustration of varying relationship between viral load and CD4 counts at different times in the MACS. A simple linear regression model was used to quantify their correlation at different study times. Data (’+’) are overlaid with the regression line. The estimates of intercept and slope as well as confidence interval are presented above each plot

The rest of article is organized as follows. In Section 2, we describe the dataset that motivated this research and the joint statistical models that account for skewness and missingness are introduced. The Bayesian approach that estimates the parameters in the joint model is presented in Section 3. In Section 4, the proposed joint models and inferential method are applied to MACS data and analysis results are presented. To validate the proposed model and method, simulation studies are carried out in Section 5. The article is concluded with discussion in Section 6.

2 Motivating data and joint models

2.1 Motivating data

The data that motivates this study was gathered in MACS, an ongoing prospective study of the natural and treated histories of HIV-1 infection in homosexual and bisexual men conducted by multiple sites (Baltimore, Chicago, Pittsburgh and Los Angeles). A total of 6 972 men at risk of HIV infection have been enrolled in three cohorts (first cohort $n = 4 954$ , from 1984 to 1985; second cohort $n = 668$ , from 1987 to 1991; third cohort $n = 1 350$ , from 2001 to 2003). The study participants had baseline and semi-annual follow-up visits. There are 1 598 subjects with more than three measurements on HIV viral load and CD4 counts. Among all patients, 7.3% died from AIDS-related diseases (referred to as AIDS death) and 23% died from other causes (referred to as other death). Roughly, 21% of viral load measurements are missing. Since the original scale has wide ranges, a log₁₀ transformation of viral load and standardized CD4 counts were used in the analysis for the purpose of stabilizing variation of measurement errors.

Figure 1 shows that the relationship between viral load and CD4 counts varies across time during the treatment. Thus, the relationship can not be modelled by a parametrical function. We adopt a varying coefficient mixed-effect joint model in which the non-linear relationship between viral load and CD4 counts is modelled non-parametrially. Further more, Figure 2 shows that the model errors are highly skewed. Therefore, a normality assumption is not quite realistic for this dataset. Alternatively, we adopt an asymmetric ST distribution (Azzalini and Capitanio, 2003; Sahu et al., 2003; Arellano-Valle and Genton, 2005; Azzalini and Genton, 2008; Jara et al., 2008; Ho and Lin, 2010) rather than a symmetric normal distribution.

Figure 2:

Histogram of residual of OLS estimates

Note: The skewed distribution of residuals based on OLS estimates of the subject-specific regression from 6 972 patients in the MACS study 2.2 ST mixed-effects varying-coefficient model

We denote the viral load and CD4 counts for the $i$ th subject at time $t_{ij}$ ( $i = 1, 2, \dots, n$ , $j = 1, 2, \dots, n_{i}$ ) as $y_{i} (t_{ij})$ and $z_{i} (t_{ij})$ respectively, where $n$ is the number of subjects in the study and $n_{i}$ is the number of longitudinal measurements for subject $i$ . We further let $y_{i} = (y_{i} (t_{i 1}), \dots, y_{i} (t_{{in}_{i}}))^{T}$ , $z_{i} = (z_{i} (t_{i 1}), \dots, z_{i} (t_{{in}_{i}}))^{T}$ , $t_{i} = (t_{i 1}, \dots, t_{{in}_{i}})^{T}$ , $ε_{i} = (ε_{i} (t_{i 1}), \dots, ε_{i} (t_{{in}_{i}}))^{T}$ . The mixed-effects varying-coefficient model is as follows:

\begin{matrix} y_{i} = β_{0} (t_{i}) + z_{i} β_{1 i} (t_{i}) + ε_{i} \\ β_{1 i} (t_{i}) = β_{1} (t_{i}) + γ_{i} (t_{i}), \end{matrix}

(2.1)

where

z_{i} * β_{1 i} (t_{i}) = (z_{i 1} β_{1 i} (t_{i 1}), \dots, z_{{in}_{i}} β_{1 i} (t_{{in}_{i}}))^{T}

. The coefficient vectors

β_{0} (t_{i}) = (β_{0} (t_{i 1}),

\dots, β_{0} (t_{{in}_{i}}))^{T}

β_{1 i} (t_{i}) = (β_{1 i} (t_{i 1}), \dots, β_{1 i} (t_{{in}_{i}}))^{T}

represent the time-varying coefficients for intercept and slope, respectively. Both

β_{1} (.)

and

γ_{i} (.)

are unknown smoothing functions.

β_{1} (.)

stands for the population smoothing curve while

γ_{i} (.)

represents the random-effects. Altogether,

β_{1 i} (.)

is the smoothing curve for individual

i

. The error term

ε_{i} (.)

and random smoothing function

γ_{i} (.)

are zero mean stochastic processes and are independent from each other. To fit model (2.1), we adopt a regression spline method. The main idea is to approximate

β_{0} (.)

β_{1} (.)

and

γ_{i} (.)

with a linear combination of spline basis functions:

\begin{matrix} β_{0, p} (t) \approx \sum_{k = 1}^{p} ξ_{k} θ_{k} (t) = Θ_{p} (t)^{T} ξ_{p} \\ β_{1, q} (t) \approx \sum_{k = 1}^{q} χ_{k} ϕ_{k} (t) = Φ_{q} (t)^{T} χ_{q}, \\ γ_{i, r} (t) \approx \sum_{k = 1}^{r} b_{ik} ψ_{k} (t) = Ψ_{r} (t)^{T} b_{i} \end{matrix}

(2.2)

where

ξ_{p} = (ξ_{1}, \dots, ξ_{p})^{T}

and

χ_{q} = (χ_{1}, \dots, χ_{q})^{T}

are the fixed-effects and

b_{i} = (b_{i 1}, \dots, b_{ir})^{T}

is regarded as the random-effects. We denote

Θ_{pi} = (Θ_{p} (t_{i 1}), \dots, Θ_{p} (t_{{in}_{i}}))^{T}

Φ_{qi} = (Φ_{q} (t_{i 1}), \dots,

Φ_{q} (t_{{in}_{i}}))^{T}

and

Ψ_{ri} = (Ψ_{r} (t_{i 1}), \dots, Ψ_{r} (t_{{in}_{i}}))^{T}

as basis functions such as cubic B-spline basis. The dimensional parameters

p

q

r

and location of knots for regression splines are based on suggestions in Eubank (1999). Plugging (2.2) into (2.1), we have

\begin{matrix} y_{i} = Θ_{pi} ξ_{p} + z_{i} * Φ_{qi} χ_{q} + z_{i} * Ψ_{ri} b_{i} + ε_{i}, \end{matrix}

(2.3)

where

z_{i} * Φ_{qi} = (z_{i 1} Φ_{q} (t_{i 1}), \dots, z_{{in}_{i}} Φ_{q} (t_{{in}_{i}}))^{T}

and

z_{ij} Φ_{q} (t_{ij})^{T} = (z_{ij} ϕ_{1} (t_{ij}), \dots, z_{ij} ϕ_{q} (t_{ij}))

z_{i} * Ψ_{ri}

follows similar line. We denote

X_{i} = (Θ_{pi}, z_{i} * Φ_{qi})

ζ = (ξ_{p}^{T}, χ_{q}^{T})^{T}

Z_{i} = z_{i} * Ψ_{ri}

. We can write (2.3) as

\begin{matrix} y_{i} = X_{i} ζ + Z_{i} b_{i} + ε_{i}, \end{matrix}

(2.4)

which is a standard linear mixed effects (LME) model if we assume that

X_{i}

and

Z_{i}

are the fixed-effects and random-effects design matrix, respectively;

ζ

and

b_{i}

are the fixed-effects and random-effects parameter vectors, respectively;

b_{i} \sim N (0, Σ)

and

ε_{i} \sim N (0, R_{i})

. However, in practice, the viral load

y_{i}

are most likely not normally distributed. In this case, we may assume

ε_{i} \overset{iid}{\sim} {ST}_{n_{i}, ν} (- J (ν) δ_{e} 1_{n_{i}}, σ^{2} I_{n_{i}}, δ_{e} I_{n_{i}})

, which follows a multivariate ST distribution with degrees of freedom

ν

, unknown scale parameter

σ^{2}

and skewness parameter

δ_{e}

, where

1_{n_{i}} = (1, \dots, 1)^{T}

and

J (ν) = (ν / π)^{1 / 2} [Γ ((ν - 1) / 2) / Γ (ν / 2)]

was set here in order to have a zero mean vector for ST distribution of error term.

2.3 Cause-specific proportional hazard model

In the follow up of a study, a subject may experience multiple types of failure or could be right censored. We assume that there are $K$ types of failure and use indicator $d_{i}$ to represent each type. For example, $d_{i} = 0$ indicates an independent censoring event and $d_{i} = k$ indicates the $k$ th type of failure, where $k = 1, \dots, K$ . We denote $D_{i}$ as the failure time, $C_{i}$ as the censoring time and $T_{i} = min (D_{i}, C_{i})$ as the observed time for subject $i$ . To model the competing risks survival process, we consider a cause-specific hazard model.

λ_{ik} (t) = λ_{0 k} (t) exp (ψ_{k} z_{i} + b_{i}^{T} c_{k}), k = 1, \dots, K,

(2.5)

where the function

λ_{ik} (t)

is the instantaneous hazard rate from cause

k

at time

t

for subject

i

, the coefficient

ψ_{k}

represents the effect of the covariate

z_{i}

to the

k

th type of failure,

c_{k} = (c_{k 1}, \dots, c_{kr})^{T}

is the parameter vector that links the longitudinal and competing risks survival processes. In particular,

c_{k}

measures the strength of association between viral load and each type of failure. Note that

b_{i}

is the random effects coefficients for the spline basis in (2.2) and

b_{i 1}, \dots, b_{ir}

collectively describe the non-linear time effect on the change of viral load in subject

i

relative to the population. This makes the interpretation of the parameters

c_{k}

less straightforward as

c_{k}

is the coefficients associated with

b_{i}

. Intuitively, the interpretation for any

c_{kj}

(

j = 1, \dots, r

) is that for each unit increase of subject-specific

b_{ij}

, the hazard ratio for the

k

th event is exp(

c_{kj}

). The baseline hazard

λ_{0 k} (t)

is approximated by a piecewise step function, that is,

λ_{0 k} (t) = λ_{0 k, j}

for

t

in the interval between

t_{k, j - 1}

and

t_{k, j}

, where

0 < t_{k, 1} < \dots < t_{k, J} < \infty

is the partition of time. For future reference, we denote

λ_{0; k} = (λ_{0 k, 1}, \dots, λ_{0 k, J})

The corresponding survival function for the $k$ th competing terminal event is $S_{k} (t_{i}) = exp [- \int_{0}^{t_{i}} λ_{0 k} (t) exp (ψ_{k} z_{i} + b_{i}^{T} c_{k}) dt]$ . The likelihood of the competing terminal events is

\begin{matrix} L (ϒ) & = & \prod_{k = 1}^{K} {[λ_{0 k} (t_{i}) exp (ψ_{k} z_{i} + b_{i}^{T} c_{k})]^{I (d_{i} = k)} S_{k} (t_{i})}, \end{matrix}

(2.6)

where the parameter vector

ϒ = (ψ_{1}, \dots, ψ_{K}, c_{1}^{T}, \dots, c_{K}^{T})^{T}

2.4 Informative missing data model

In clinical studies, it is very common that some patients may miss scheduled visits due to drug resistance or other reasons. The probability of missing data may be related to the missing values which renders informative or non-ignorable missingness. Statistical inferences ignoring missing data mechanism may lead to biased estimation when missing mechanism is informative or non-ignorable (Little and Rubin, 2002). To achieve reliable estimates of parameters, we consider a missing data model which relates the missingness in response to the missing values:

logit [P (r_{ij} = 1 | η)] = η_{0} + η_{1} y (t_{ij}),

(2.7)

where the missing indicator

r_{ij} = 1

y (t_{ij})

is missing and 0 otherwise. Further more, we denote

r_{i} = (r_{i 1}, \dots, r_{{in}_{i}})^{T}

as the vector of missing indicator for subject

i

η = (η_{0}, η_{1})^{T}

is the unknown parameters. Other missing response models may be considered. However, all models are not testable. Sensitivity analysis may be helpful to check the consistency of these models. For example, we may add a higher-order term of

y (t_{ij})

in model (2.7). Alternatively, covariates such as CD4 counts may be related to the missing probability as well. We tested a few alternative models accounting for the missing data mechanism. The findings suggest that parameter estimation in the joint model is not sensitive to the specification of various missing data models. To achieve parsimony, model (2.7) is employed to account for missing data mechanism.

3 Bayesian inferences

The previously discussed longitudinal process, time-to-event process and missing data process are inherently connected. The computational cost associated with the joint likelihood for the longitudinal data with skewness as well as the competing risks survival data is usually high. Often times, the frequentist's approach fails to converge (Wu, 2002). We propose a fully Bayesian approach for the joint models (2.1), (2.5) and (2.7) which handle longitudinal-competing risk survival data with skew distribution and missingness. The Markov chain Monte Carlo (MCMC) technique is employed to estimate all parameters simultaneously based on the joint likelihood.

Through introduction of a $n_{i} \times 1$ random vectors $w_{e_{i}}$ based on the stochastic representation for ST distributions (see Appendix in detail), we formulate the longitudinal response model (2.4) in association with cause-specific competing risk hazard model (2.5) and missing reponse model (2.7), hierarchically:

\begin{matrix} y_{i} | z_{i}, b_{i}, w_{e_{i}}; ζ, σ^{2}, δ_{e}, ν \sim t_{n_{i}, ν + n_{i}} (X_{i} ζ + Z_{i} b_{i} + δ_{e} [w_{e_{i}} - J (ν + n_{i}) 1_{n_{i}}], ω_{i} σ^{2} I_{n_{i}}), \\ w_{e_{i}} \sim t_{n_{i}, ν} (0, I_{n_{i}}) I (w_{e_{i}} > 0), b_{i} \sim N (0, Σ), t_{i} \sim F_{c} (t_{i} | ϒ, λ_{0; k}) . r_{ij} \sim Bernoulli (P_{ij}), \end{matrix}

(3.1)

where

ω_{i} = (ν + w_{e_{i}}^{T} w_{e_{i}}) / (ν + n_{i})

t_{n_{i}, ν} (μ, A)

denotes the

n_{i}

-variate

t

distribution with parameters

μ

A

and degrees of freedom

ν

I (w > 0)

is an indicator function and

w = | ς |

with

ς \sim t_{n_{i}, ν} (0, I_{n_{i}})

and

P_{ij} = P (r_{ij} = 1 | η)

F_{c} (\cdot)

is the distribution function related to the competing risks survival model (2.5). The full conditional distributions and corresponding sampling schemes are listed in Appendix B.

In the Bayesian framework, we specify the models for the observed data and the prior distributions for the unknown model parameters. Subsequently, based on the posterior distributions, we make statistical inference for the unknown parameters. We denote $θ = {ζ, ϒ, η, λ_{0; k}, σ^{2}, Σ, ν, δ_{e}}$ as the collection of unknown population parameters in models (2.4), (2.5) and (2.7). The prior distributions for $θ$ are as follows:

\begin{matrix} ζ \sim N (τ_{1}, Λ_{1}), σ^{2} \sim IG (γ_{σ 1}, γ_{σ 2}), Σ \sim IW (Ω, ρ), δ_{e} \sim N (0, γ), \\ ν \sim Exp (ν_{0}) I (ν > 3), ϒ \sim N (τ_{2}, Λ_{2}), λ_{0; k} \sim Γ (γ_{k 1}, γ_{k 2}), η \sim N (τ_{3}, Λ_{3}), \end{matrix}

(3.2)

where the mutually independent Normal (

N

), Inverse Gamma (

IG

), Inverse Wishart (

IW

), Gamma (

Γ

) and Exponential (

Exp

) prior distributions are chosen to facilitate computations. The hyper-parameter matrices

Λ_{1}

Λ_{2}

Λ_{3}

and

Ω

are assumed to be diagonal for implementation convenience. The exponential prior for

ν

is truncated to lie above 3 to make variance of ST distribution well defined.

We denote $f (\cdot)$ , $f (\cdot | \cdot)$ and $π (\cdot)$ as a generic density function, a conditional density function and a prior density function, respectively. If we assume that $ζ, ϒ, η, λ_{0; k}, σ^{2}, Σ, ν, δ_{e}$ are independent of each other, then $π (θ) = π (ζ) π (ϒ) π (η) π (λ_{0; k}) π (σ^{2}) π (Σ) π (ν) π (δ_{e})$ . As a result, the joint posterior density of $θ$ is

\begin{matrix} f (θ | D) & \propto & {\prod_{i}^{n} \int \int f (y_{i} | z_{i}, b_{i}, w_{e_{i}}; ζ, σ^{2}, δ_{e}, ν) f (w_{e_{i}} | w_{e_{i}} > 0) f (b_{i}) f (r_{i} | η) F_{c} (t_{i} | ϒ, λ_{0; k}) \\ d b_{i} d y_{mis, i}} π (θ), \end{matrix}

(3.3)

Generally, the integrals in equation (3.3) are of high dimension and do not have closed form. It is often not accurate to approximate the integrals numerically. Therefore, it is prohibitive to directly calculate the posterior distribution of

θ

based on the observed data. Alternatively, the MCMC procedure is powerful for drawing samples from posterior distributions, based on equation (3.3), by employing Metropolis–Hastings (M–H) algorithm along with the the Gibbs sampler.

4 Application

4.1 Model specification

The dataset that motivates this study was described in Section (2.1). As shown in Figure 1, the relationship between viral load and CD4 counts is varying during the treatment. In addition, the viral load in log $_{10}$ scale is highly skewed (Figure 2). Therefore, it is critical to consider the skew distribution, such as ST, for the varying coefficient mixed-effects model. From the biological point of view, the longitudinal process and competing risk time-to-event process are inherently connected. Ignoring either one may lead to severe bias on parameter estimation. The proposed joint models take into consideration all of these factors. Thus, we believe it will outperform other potential models. Towards this end, we first compare three statistical models with different specification of random errors for the longitudinal response model.

Model N: A joint model with the independent multivariate normal distribution of random errors for response model (2.4).

Model SN: A joint model with the independent multivariate SN distribution of random errors for response model (2.4).

Model ST: A joint model with the independent multivariate ST distribution of random errors for response model (2.4).

By comparing Models N through ST, we test if the skew distribution assumed for error terms works better than that of symmetric distribution.

To approximate the time-varying coefficients in model (2.4), we use the same basis functions for $β_{0} (t)$ , $β_{1} (t)$ and $γ_{i} (t)$ . The smoothing parameters $p$ , $q$ and $r$ are determined by Deviance Information Criterion (DIC) and the location of knots is selected at the quantiles of the data (Eubank, 1999). For MACS dataset, the optimal smoothing parameters are set at 3 for $p$ , $q$ and $r$ . To stabilize the variance and reach quick convergence, we use the log $_{10}$ transformed viral load and standardized CD4 counts.

In terms of MACS data, we are interested in two cause-specific time-to-events: HIV death and other death. Thus, the competing risk hazard models for the two risks are:

\begin{matrix} λ_{i 1} (t) = λ_{01} (t) exp (ψ_{1} z_{i}^{*} + b_{i}^{T} c_{1}) and \\ λ_{i 2} (t) = λ_{02} (t) exp (ψ_{2} z_{i}^{*} + b_{i}^{T} c_{2}) . \end{matrix}

(4.1)

We take following prior distribution for the parameters in the joint models: (a) fixed-effects are taken to be independent normal distributions $N (0, 100)$ for each component of the population parameter vectors $ζ$ and $ϒ$ and $η$ ; (b) for the scale parameters $σ^{2}$ , we assume an inverse gamma prior distribution, $IG (0.01, 0.01)$ so that the distribution has mean 1 and variance 100; (c) the degree of freedom parameters $ν$ follows truncated exponential distribution with $ν_{0} = 0.5$ ; (d) for the skewness parameters $δ_{e}$ , we choose independent normal distribution $N (0, 100)$ ; (e) inverse Wishart distributions $IW (Ω, ρ)$ is taken as priors for the variance–covariance matrice of the random-effects $Σ$ , with covariance matrices $Ω = diag (0.01, 0.01, 0.01)$ and $ρ = 3$ , respectively; and (f) for each of piecewise baseline hazard $λ_{0 k, j}$ , gamma distributions $Γ (0.1, 0.1)$ are employed.

The proposed method is implemented with WinBUGS software (Lunn et al., 2000). The MCMC scheme iterates between the Gibbs sampler and the M–H algorithm for drawing samples from posterior distributions of parameters in the joint models. Based on MCMC samples, we are able to draw statistical inference for interested parameters. Specifically, we are interested in the posterior means (PM) and quantiles. Convergence of MCMC algorithm is assessed by standard tools in WinBUGS software such as trace plots and Gelman–Rubin diagnostics (Gelman and Rubin, 1992). We will run multiple chains. After convergency, we keep sampling for certain period and extract samples after every few iterations. The collected samples will be used to make inferences for the posterior distribution of interested parameters.

When convergence was achieved, for each of the three chains, after an initial number of 50 000 burn-in iterations, every 50th MCMC sample is retained from the next 50 000. Thus, we obtain 3 000 samples of targeted posterior distributions of the unknown parameters for statistical inference.

Figure 3:

Gelman–Rubin diagnostic plot based on Model ST for representative parameters from three Markov chains. The middle and bottom curves below the dashed line which indicates value one, stand for the pooled posterior variance and average within-sample variance separately, and the top curve represents their ratio

Figure 3 shows the dynamic Gelman–Rubin diagnostics based on model ST as obtained from the WinBUGS output for the representative parameters where the three curves are given: the middle and bottom curves below the dashed horizontal line (indicating value one) represent the pooled posterior variance and average within-sample variance, respectively, and the top curve above the dashed horizontal line represents their ratio. It is seen that the ratio tends to 1, and both pooled posterior variance and average within-sample variance settle down as the number of iterations increases, indicating that the algorithm has reached convergence.

4.2 Model comparison

Tables 1 and 2 present the population PM estimates, standard deviation (SI) as well as associated 95% credible intervals CI for fixed-effects, scale and skewness parameters in the three models. We observe quite distinct estimates among three models. In particular, the estimated variance parameter $σ^{2}$ is fairly large in Model N (1.63) comparing to the values in Models ST (0.18) and SN (0.39). This finding is expected since Model N does not account for skewness and heavy tails in the data. The estimates of skewness parameter $δ$ are $-$ 0.16 with 95% CI ( $-$ 0.26, $-$ 0.04) and $-$ 1.07 with 95% CI ( $-$ 1.14, $-$ 1.03) in Models SN and ST, separately. This finding suggests that there is a significantly left skewness in the data which confirms the fact that the viral load data is left skewed (Figure 1). Therefore, it is critical to incorporate skewness in the modelling process. The estimates of $c_{11}, c_{12}, c_{13}, c_{21}, c_{22} and c_{23}$ are mostly significant in all models. Since they are parameters that link the longitudinal and survival processes, the non-zero estimates demonstrate the strong association between the two processes.

Figure 4 examines the goodness-of-fit for the three models. Obviously, comparing to Model N which assumes normal distribution for random errors, Models ST and SN that assume asymmetric distribution for random errors provide a better fit to the observed data. In addition, Q–Q plots–of residuals show fewer outliers from Models ST/SN than Model N. Moreover, Model ST fits data better than Model SN since it accounts for heavy tails. Hence, accounting for skewness and heavy tails in data is essential for better model performance.

Figure 4:

Diagnostics of model fitting for the three models. Top panel: fitted versus observed values of log $_{10}$ viral load; bottom panel: Q–Q plots for residuals in three models

To select the best model that fits the data, we use a Bayesian model selection criterion: DIC introduced by Spiegelhalter et al. (2002). In addition, we compute and compare the expected predictive deviance (EPD) for each model. EPD is defined as $EPD = E \{\sum_{i, j} (y_{rep, ij} - y_{obs, ij})^{2}\}$ where the predictive value $y_{rep, ij}$ is a replicate of the observed $y_{obs, ij}$ and the expectation is taken over the posterior distribution of the model parameters $θ$ (Gelman et al., 2003). The model has the least discrepancy between the observed values and the predictive values is the best one. The DICs for the three models are 27154, 11383 and 8117, respectively. The best selected model is Model ST, followed by Models SN and N. According to EPDs, Model ST is still the best selected model (0.16), followed by Models SN (0.51) and N (1.85). The results are consistent with the goodness of fit diagnostics in Figure 4 in which Model ST fits the data better than others. Based on these results, we can conclude that accounting for skewness in the joint model is essential in order to achieve more reliable results when data exhibit skewness.

Table 1:

A summary of estimated PM for the population parameters in longitudinal model (2.4), cause-specific hazard model (4.1), missing data model (2.7) and the corresponding SD, lower limit ( $L_{CI}$ ) and upper limit ( $U_{CI}$ ) of $95 %$ equal-tail CI

Model		ξ₁	ξ₂	ξ₃	χ₁	χ₂	χ₃	ψ₁	ψ₂	η₀	η₁
N	PM	8.36	3.03	4.03	-2.95	-3.9	-1.28	-2.72	-1.33	12.63	-6.0
	$L_{CI}$	4.3	2.71	2.85	-3.79	-4.13	-2.34	-4.84	-2.47	10.22	-7.78
	$U_{CI}$	12.27	3.35	5.14	-2.1	-3.65	-0.34	-0.72	-0.24	15.01	-4.12
	SD	1.94	0.17	0.61	0.43	0.12	0.53	1.04	0.57	1.21	0.86

SN	PM	8.12	4.4	-1.24	-3.5	-1.39	-0.87	-2.3	-1.48	6.44	-5.3
	$L_{CI}$	6.38	4.34	-1.45	-4.01	-1.49	-1.39	-4.4	-2.59	2.96	-8.3
	$U_{CI}$	9.87	4.47	-1.03	-2.97	-1.29	-0.44	-0.45	-0.44	9.42	-3.78
	SD	0.84	0.03	0.11	0.26	0.05	0.25	1.0	0.56	2.14	1.33

ST	PM	7.13	4.29	-0.89	-1.79	-0.37	-0.64	-2.22	-1.56	15.04	-3.52
	$L_{CI}$	5.74	4.23	-1.1	-2.15	-0.44	-0.99	-4.38	-2.73	8.28	-6.6
	$U_{CI}$	8.46	4.35	-0.68	-1.45	-0.29	-0.31	-0.29	-0.42	22.9	-0.43
	SD	0.66	0.03	0.11	0.16	0.04	0.17	1.05	0.59	3.52	1.57

T	PM	8.46	3.54	5.22	-3.08	-5.41	-1.48	-2.68	-1.56	11.53	-2.03
	$L_{CI}$	4.26	3.26	4.02	-3.95	-5.66	-2.55	-4.76	-2.73	10.46	-1.22
	$U_{CI}$	12.65	3.82	6.42	-2.21	-5.17	-0.4	-0.67	-0.42	12.6	-2.97
	SD	2.1	0.14	0.6	0.46	0.19	0.58	1.01	0.57	0.6	0.51

S	PM	7.28	5.42	-0.98	-2.46	-0.27	-0.58	-2.01	-1.23	18.24	-3.70
	$L_{CI}$	6.23	5.37	-1.17	-2.75	-0.33	-0.87	-3.25	-2.09	13.83	-5.53
	$U_{CI}$	8.33	5.46	-0.78	-2.18	-0.21	-0.29	-0.76	-0.36	22.66	-1.88
	SD	0.53	0.02	0.10	0.14	0.03	0.14	0.62	0.43	2.21	0.91

NM	PM	4.19	-1.94	4.82	-0.54	-0.59	-2.05	-2.63	-1.61	–	–
	$L_{CI}$	4.11	-2.12	3.45	-0.66	-1.02	-3.53	-4.74	-2.74	–	–
	$U_{CI}$	4.25	-1.75	6.16	-0.42	-0.18	-0.47	-0.68	-0.53	–	–
	SD	0.03	0.1	0.68	0.06	0.21	0.78	1.04	0.56	–	–

Table 2:

The estimated PM for parameters of scale, skewness and the corresponding SD and lower limit ( $L_{CI}$ ) and upper limit ( $U_{CI}$ ) of $95 %$ equal-tail CI as well as DIC and EPD values for MACS study

,Model		$σ^{2}$	$δ_{e}$	$c_{11}$	$c_{12}$	$c_{13}$	$c_{21}$	$c_{22}$	$c_{23}$	DIC	EPD
N	PM	1.63	–	1.89	-0.31	-0.52	-0.34	-0.32	0.23	27 154	1.85-
ampersand 1.74	$L_{CI}$	1.14	–	1.36	-0.65	-0.89	-1.64	-0.71	-0.13
	$U_{CI}$	2.01	–	2.37	0.02	-0.19	-1.08	-0.07	0.63
	SD	0.26	–	0.24	0.16	0.17	0.14	0.18	0.19

SN	PM	0.39	-0.16	1.14	-0.38	-0.69	-0.75	-0.24	0.52	11 383	0.51
ampersand 2.84	$L_{CI}$	0.33	-0.26	0.71	-0.75	-0.94	-1.03	-0.44	0.31
ampersand 2.17	$U_{CI}$	0.47	-0.04	1.86	-0.02	-0.41	-0.45	-0.01	0.79
ampersand 3.52	SD	0.04	0.06	0.22	0.17	0.13	0.13	0.11	0.12
ampersand 0.32
ST	PM	0.18	-1.07	0.96	-0.24	-0.35	-0.45	-0.11	0.81	8 117	0.16
ampersand 2.51	$L_{CI}$	0.15	-1.14	0.51	-0.45	-0.53	-0.82	-0.19	0.54
ampersand 2.01	$U_{CI}$	0.22	-1.01	1.34	-0.05	-0.15	-0.08	-0.04	1.13
ampersand 3.03	SD	0.02	0.03	0.23	0.1	0.09	0.18	0.03	0.14
ampersand 0.25
T	PM	1.22	–	1.88	-0.17	-1.42	-0.4	-1.2	0.01	24 378	1.59
	$L_{CI}$	0.57	–	1.32	-0.5	-1.88	-0.58	-1.48	-0.12
	$U_{CI}$	1.87	–	2.44	0.16	-0.96	-0.22	-0.92	0.13
	SD	0.32	–	0.28	0.17	0.23	0.09	0.14	0.06

S	PM	0.17	-1.18	–	–	–	–	–	–	8 426	0.24
ampersand 2.51	$L_{CI}$	0.15	-1.24	–	–	–	–	–	–
ampersand 2.01	$U_{CI}$	0.2	-1.13	–	–	–	–	–	–
ampersand 3.03	SD	0.01	0.03	–	–	–	–	–	–
ampersand 0.25
NM	PM	0.29	-1.44	1.69	-0.45	0.04	-0.3	-0.28	0.4	12 166	0.54
ampersand 2.51	$L_{CI}$	0.25	-1.52	0.84	-0.67	-0.03	-0.37	-0.56	-0.19
ampersand 2.01	$U_{CI}$	0.36	-1.36	2.55	-0.04	0.09	-0.21	-0.01	0.11
ampersand 3.03	SD	0.03	0.04	0.42	0.19	0.03	0.04	0.15	0.28
ampersand 0.25

4.3 Results based on Model ST

Since Model ST is superior to Models SN/N, we carry out further analysis based on Model ST. To assess how the missing data in response influences modelling outcome, we fit another model (Model NM) in which the missing data mechanism is not taken into account in the joint model. It is observed based on Tables 1 and 2 that there are important differences in the estimates of parameters between Models NM and ST. The estimated DIC/EPD values are larger for Model NM (12 166/0.54) than Model ST (8 117/0.16). In addition, the estimates of parameters $η$ in the informative missing data model (2.7) based on Model ST are significantly different from zero. This finding demonstrates that missingness of response is strongly associated with the missing viral load values. Thus, it is important to account for the informative missing viral load data when collected data are incompletely measured.

Figure 5:

The population estimating curve of varying-coefficients $β_{0} (t)$ , $β_{1} (t)$ for longitudinal model (2.1) based on the best selected Model ST. The estimates (solid curve) along with the 95% CI (dashed curves) are presented

The primary objective is to investigate the time-varying relationship between viral load and CD4 counts. The smoothing estimate of population varying-coefficients ( $β_{0} (t)$ and $β_{1} (t)$ in model (2.1)) against time is plotted in Figure 5 along with 95% CI. We see that there is a strong negative relationship between viral load and CD4 counts. This is because CD4s are one of the major sources that eradicate HIV virus. Thus, the negative relationship reflects the biological mechanism. We further note that the relationship becomes stronger (more negative) during the treatment period. The increment of negative relationship ( $β_{1} (t)$ ) is almost linear before the first 12 years but it attenuates afterwards. The 95% CI is narrower at the beginning of the treatment but gradually becomes wider at the later stage. This is because there is more dropout at the later stage than the earlier period so that more uncertainty evolves alongside the parameter estimates.

Figure 6:

Illustration of four represented individual estimating curve (dashed lines) of $β_{1 i} (t)$ based on Model ST. Population curves (solid lines) are shown for comparison

The mixed-effects modelling approach makes it possible for the parameter estimation at both the population and individual level, simultaneously, since it accounts for both within and between subjects variation. The estimates at the individual level may not follow the pattern of the population level exactly if the between subject variation is large. Figure 6 presents a few patterns for the time-varying relationship between viral load and CD4 counts at the individual level. For comparison purpose, the population estimating curve ( $β_{1} (t)$ , solid line) is overlayed with the individual estimating curve ( $β_{1 i} (t)$ , dotted line). We observe that some of the individuals, such as subject 69, follow a similar pattern as the population curve, that is, the negative relationship between viral load and CD4 counts varies across time in a similar fashion as the population. Nevertheless, some individuals, for instance subject 186, demonstrate stronger negative relationship as reflected by the steeper slope than the population when the treatment progresses but the strength is weaken later. On the other hand, there are also subjects who display opposite trend as compared to the population, for example, subject 4 who shows a gradually increasing curve suggesting a reduced negative relationship between viral load and CD4 counts. Another pattern typically observed in the study is the flat curves, such as subject 122. For these subjects, the estimates of $β_{1} (t)$ do not vary across time in general. Based on these observations, we see that there are much variation in the varying-coefficients between individuals. The mixed-effects varying-coefficient joint model characterizes the relationship between the viral load and CD4 counts at both the individual and population levels.

5 Simulation

We conduct simulation studies to evaluate the performance of the proposed joint models. The design of simulated data is similar to the joint models used for the real data. Specifically, the longitudinal outcomes are simulated from the mixed-effects model (2.4). The measurement time $t_{ij}$ are same to those in the real data. Similar to the real data, we simulate two competing risks according to the cause-specific PH model (4.1). To generate event time data, we use constant baseline hazard of 0.1 and 0.2 for two risks, respectively. An exponential distribution with mean equal to 0.1 is used to generate censoring time. The informative missing data model (2.7) is used to generate missing data in response. We set the true values of model parameters as those obtained from real data analysis for Model ST which are listed in Tables 1 and 2. We simulate data for 300 subjects in the simulation study. One advantage of the proposed joint model is the capability of handling skewness exhibited in data. To simulate data with skewness, we generate the random errors for model (2.4) as follows: first, data are sampled from a gamma distribution Γ(3,1); second, the samples are subtracted by 3 which yields skewed data with mean 0 and variance 3. Under the settings described above, we simulate 200 datasets.

Table 3:

Simulation results based on 200 simulated dataset. The true values, bias (%) and coverage rate (%) of 95% CI are presented for each parameter

,		Model N		Model NM		Model SN		Model ST
Parameter	True value	Bias	Coverage	Bias	Coverage	Bias	Coverage	Bias	Coverage
$α_{1}$	0.54	-5.55	93.0	-4.03	91.0	-2.96	91.5	-1.84	93.5
$α_{2}$	-4.52	-8.62	88.5	-8.00	92.0	-5.47	95.0	-3.38	96.5
$α_{3}$	3.82	-10.91	89.0	-8.13	93.0	-4.90	92.0	-3.44	97.5
$ξ_{1}$	4.41	-5.28	94.0	-2.99	95.0	-2.77	96.0	-1.62	97.5
$ξ_{2}$	-1.23	6.15	94.0	3.41	92.5	2.57	96.0	2.02	94.5
$ξ_{3}$	2.82	14.70	92.0	9.44	93.5	9.14	92.0	4.87	96.0
$χ_{1}$	-0.45	-3.84	89.5	-3.05	95.5	-1.95	93.5	-1.49	97.5
$χ_{2}$	-0.44	14.12	94.5	10.81	93.0	6.12	95.0	4.42	99.0
$χ_{3}$	-1.4	4.10	91.0	2.64	91.0	2.30	94.5	1.56	97.5
$ψ_{1}$	-2.46	-5.21	91.5	-3.80	96.0	-2.77	91.5	-1.82	96.0
$ψ_{2}$	-1.56	11.72	89.0	8.22	93.5	6.59	94.5	4.36	95.5
$c_{11}$	0.51	6.57	90.0	3.80	92.0	3.68	96.0	2.14	96.5
$c_{12}$	1.28	11.21	93.0	6.89	94.5	6.16	92.5	3.80	97.0
$c_{13}$	-0.1	14.50	95.0	9.08	96.0	7.43	93.0	4.74	94.5
$c_{21}$	0.21	11.79	89.0	9.16	94.5	7.00	94.5	4.52	96.0
$c_{22}$	-0.41	10.11	90.5	6.59	92.0	5.99	95.0	3.69	94.5
$c_{23}$	0.42	-2.88	94.5	-2.29	90.0	-2.10	92.0	-1.11	95.0

For each simulated dataset, we adopted similar models and MCMC sampling schemes that were used for the real data analysis. Priors with large variances are used in the Bayesian inference. We fit each of the Models ST, SN, N and NM to each simulated dataset and compare their performances. Table 3 presents the simulation results. As can be seen, Models ST (bias ranges from $-$ 3.44% to 4.87%) performs best with estimates close to true values, followed by Model SN (bias ranges from $-$ 5.47% to 9.14%) and Model NM (bias ranges from $-$ 8.13% to 10.81%). However, Model N produces larger bias (ranges from $-$ 10.91% to 14.7%). In summary, the simulation results show that it is important to account for both skewness and missingness present in the data.

6 Discussion

We developed a Bayesian varying coefficient mixed-effects joint models to deal with longitudinal survival data with asymmetry and missingness. This type of joint modelling approach is useful in many application areas, producing accurate estimation of parameters while accounting for skewness and incompleteness in longitudinal survival data. Although we illustrated the method with an AIDS study, the basic idea of the proposed model and method has broader applications whenever technical conditions are met. The proposed joint modelling approach using skew distributions proves that there is potential gain of efficiency and accuracy in parameter estimation when symmetric distribution does not hold. In particular, we assessed the contribution of asymmetric distributions for model errors to model building and parameter estimation, comparing to a symmetric distribution. The findings suggest that to achieve reliable results for skewed data, it is critical to assume a skew distribution in the joint model. We further investigated how missing data affect modelling outcomes. We found that the informative missing data play an important role in parameter estimation. When data are not completely measured, it is crucial to adopt an appropriate missing data model to avoid biased estimation. The proposed model and method is useful for evaluating the relationship between viral load and CD4 count in AIDS studies which are two important indices for treatment assessment.

Note that the varying coefficient mixed-effects models are not typical non-parametric models. Formally, a non-parametric model contains infinite-dimensional parameters and, therefore, not all varying coefficient mixed-effects model are formally non-parametric. As a matter of fact, the model considered in this article, based on a finite number of spline bases, is in fact parametric.

Table 4:

Sensitivity analysis for different sets of hyper-parameters based on Model ST

,Set		$ξ_{1}$	$ξ_{2}$	$ξ_{3}$	$χ_{1}$	$χ_{2}$	$χ_{3}$	$ψ_{1}$	$ψ_{2}$	$η_{0}$	$η_{1}$	$σ^{2}$	$δ_{e}$
1	PM	7.13	4.29	-0.89	-1.79	-0.37	-0.64	-2.22	-1.56	15.04	-3.52	0.18	-1.07
	$L_{CI}$	5.74	4.23	-1.1	-2.15	-0.44	-0.99	-4.38	-2.73	8.28	-6.6	0.15	-1.14
	$U_{CI}$	8.46	4.35	-0.68	-1.45	-0.29	-0.31	-0.29	-0.42	22.9	-0.43	0.22	-1.01
	SD	0.66	0.03	0.11	0.16	0.04	0.17	1.05	0.59	3.52	1.57	0.02	0.03

2	PM	7.16	4.28	-0.9	-1.8	-0.37	-0.64	-2.23	-1.55	15.15	-3.54	0.18	-1.07
	$L_{CI}$	5.97	4.22	-1.11	-2.13	-0.46	-0.94	-4.61	-2.76	7.08	-6.84	0.14	-1.15
	$U_{CI}$	8.35	4.34	-0.67	-1.46	-0.28	-0.34	-0.22	-0.33	23.23	-0.23	0.22	-1
	SD	0.6	0.03	0.11	0.17	0.05	0.15	1.09	0.61	4.04	1.65	0.02	0.04

3	PM	7.13	4.26	-0.89	-1.78	-0.37	-0.65	-2.24	-1.55	15.04	-3.54	0.18	-1.06
	$L_{CI}$	5.85	4.2	-1.09	-2.07	-0.45	-0.97	-4.72	-2.77	8.08	-6.44	0.15	-1.12
	$U_{CI}$	8.41	4.33	-0.7	-1.5	-0.29	-0.32	-0.24	-0.38	22.0	-0.65	0.22	-1
	SD	0.64	0.03	0.1	0.14	0.04	0.16	1.14	0.63	3.48	1.45	0.02	0.03

We assume an inverse gamma prior distribution, $IG (γ, γ)$ , for the scale parameters $σ^{2}$ . A difficulty of this prior distribution is that in the limit of $γ$ $\to$ 0, it yields an improper posterior density, and thus $γ$ must be set to a reasonable value. For datasets in which low values of $σ^{2}$ are possible, inferences become very sensitive to $γ$ in this model, and the prior distribution hardly looks non-informative (Gelman, 2006). Fortunately, in our data, the estimate of $σ^{2}$ is not very small. Thus, the use of inverse gamma prior is not an issue.

In Bayesian analysis, it is critical to perform sensitivity analysis to see if the posterior estimates change significantly when priors are different. Towards this end, we carried out sensitivity analysis by employing a few sets of different values for the hyper-parameters in (3.2) and re-run the MCMC sampling scheme. Basically, we increase the hyper-parameter values adopted in Section 4.1 by 10 and 100 times for Model ST and re-run the analysis. The parameter estimates based on different set of hyper-parameter values are presented in Table 4. For comparison purpose, the results based on original set of hyper-parameter values in Section 4.1 are denoted as ‘Set 1’. The findings show that, except for some numerical fluctuation, the parameter estimating results are essentially same based on the different set of hyper-parameter values. Thus, we are confident that the obtained results are robust against hyper-parameter values. In model (2.1), we adopted the regression spline basis to represent the unknown smoothing function. There are a lot of alternative ways for approximating the unknown function such as local polynomial kernel and smoothing splines. It is interesting to compare the modelling results based on various methods.

The proposed model and method rely on the conditional PH assumption on the survival sub model, which is rather strong and could lead to wrong inferences when the PH assumption does not hold. Our model can be extended to handle clustered data. Clustered data arise frequently from multi-site clinical trials in which each site can be viewed as a cluster, from studies with families data, or from studies with recurrent events for each subject. A simple method is to add an additional random effect in our joint model to adjust for heterogeneity across the clusters. Our computational MCMC algorithm can be adopted without much modification.

According to Associate Editor's suggestion, we test a few alternative models: Model T assume a student $t$ distribution for the model errors. Based on ST distribution for random errors, we fit a separate model (Model S) by setting $c_{k}$ to 0 in the cause-specific hazard model (2.5). In this way, the longitudinal and survival processes are separated. The parameter estimations are presented in Tables 1 and 2. We find that Model ST outperforms Models T and S based on the comparison of DIC/EPD (8 117/0.16 vs. 24 378/1.59 and 8 426/0.24). Therefore, it is important to assume asymmetric distribution for random errors and the joint modelling of longitudinal and survival processes is necessary to achieve better model fitting results.

7 Appendix A. Multivariate skew distributions

Different versions of multivariate skew distributions have been introduced in the literature (Azzalini and Capitanio, 2003; Sahu et al., 2003; Azzalini and Genton, 2008; Jara et al., 2008). A new class of distributions by introducing skewness in multivariate elliptically distributions were developed in publication (Sahu et al., 2003). The class, which is obtained by using transformation and conditioning, contains many standard families including the multivariate SN and skew- $t$ ST distributions as special cases. A $k$ -dimensional random vector $Y$ follows a $k$ -variate skew-elliptical (SE) distribution if its probability density function (pdf) is given by

\begin{matrix} f (y | μ, Σ, Δ; m_{ν}^{(k)}) = 2^{k} f (y | μ, A; m_{ν}^{(k)}) P (V > 0), \end{matrix}

(A.1)

where

A = Σ + Δ^{2}

μ

is a location parameter vector,

Σ

is a

k \times k

positive (diagonal) covariance matrix,

Δ = diag (δ_{1}, δ_{2}, \dots, δ_{k})

is a

k \times k

skewness matrix with the skewness parameter vector

δ = (δ_{1}, δ_{2}, \dots, δ_{k})^{T}

;

V

follows the elliptical distribution

El (Δ A^{- 1} (y - μ), I_{k} - Δ A^{- 1} Δ; m_{ν}^{(k)})

and the density generator function

m_{ν}^{(k)} (ζ) = \frac{Γ (k / 2)}{π^{k / 2}} \frac{m_{ν} (ζ)}{\int_{0}^{\infty} r^{k / 2 - 1} m_{ν} (ζ) dr}

, with

m_{ν} (ζ)

being a function such that

\int_{0}^{\infty} r^{k / 2 - 1} m_{ν} (ζ) dr

exists. The function

m_{ν} (ζ)

provides the kernel of the original elliptical density and may depend on the parameter

ν

. This SE distribution is denoted by

SE (μ, Σ, Δ; m^{(k)})

. Two examples of

m_{ν} (ζ)

, leading to important special cases used throughout the paper, are

m_{ν} (ζ) = exp (- ζ / 2)

and

m_{ν} (ζ) = (1 + ζ / ν)^{- (ν + k) / 2}

, where

ν > 0

. These two expressions lead to the multivariate SN and ST distributions, respectively. In the latter case,

ν

corresponds to the degrees of freedom parameter.

As we know, a normal distribution is a special case of an SN distribution when the skewness parameter is zero, and the ST distribution reduces to the SN distribution when degrees of freedom are large. For completeness, this Appendix briefly summarizes the multivariate ST distribution introduced by (Sahu et al., 2003) to be suitable for a Bayesian inference since it is built using the conditional method. For detailed discussions on properties of ST distribution, see Sahu et al. (2003). Assume that a $k$ -dimensional random vector $Y$ follows a $k$ variate ST distribution with location vector $μ$ , $k \times k$ positive (diagonal) covariance matrix $Σ$ and $k \times k$ skewness matrix $Δ = diag (δ_{1}, δ_{2}, \dots, δ_{k})$ or the degrees of freedom $ν$ .

A $k$ -dimensional random vector $Y$ follows an $m$ -variate ST distribution if its pdf is given by

\begin{matrix} f (y | μ, Σ, Δ, ν) = 2^{k} t_{k, ν} (y | μ, A) P (V > 0), \end{matrix}

(A.2)

we denote the

k

-variate

t

distribution with parameters

μ

A

and degrees of freedom

ν

t_{k, ν} (μ, A)

and the corresponding pdf by

t_{k, ν} (y | μ, A)

henceforth,

V

follows the

t

distribution

t_{k, ν + k}

. We denote this distribution by

{ST}_{k, ν} (μ, Σ, Δ)

. In particular, when

Σ = σ^{2} I_{k}

and

Δ = δ I_{k}

, the equation (A.2) simplifies to

\begin{matrix} f (y | μ, σ^{2}, δ, ν) = 2^{k} (σ^{2} + δ^{2})^{- k / 2} \frac{Γ ((ν + k) / 2)}{Γ (ν / 2) (ν π)^{k / 2}} {\{1 + \frac{(y - μ)^{T} (y - μ)}{ν (σ^{2} + δ^{2})}\}}^{- (ν + k) / 2} \\ \times T_{k, ν + k} [{\{\frac{ν + (σ^{2} + δ^{2})^{- 1} (y - μ)^{T} (y - μ)}{ν + k}\}}^{- 1 / 2} \frac{δ (y - μ)}{σ \sqrt{σ^{2} + δ^{2}}}], \end{matrix}

(A.3)

where

T_{k, ν + k} (\cdot)

denotes the cumulative distribution function (cdf) of

t_{k, ν + k} (0, I_{k})

. However, unlike in the SN distribution, the ST density cannot be written as the product of univariate ST densities. Here

Y

is dependent but uncorrelated. It is noted that when

δ = 0

, the ST distribution reduces to the usual

t

-distribution. It can be shown that the mean and covariance matrix of the ST distribution

{ST}_{k, ν} (μ, σ^{2} I_{k}, Δ)

are given by

\begin{matrix} E (Y) & = & μ + (ν / π)^{1 / 2} \frac{Γ ((ν - 1) / 2)}{Γ (ν / 2)} δ, \\ cov (Y) & = & [σ^{2} I_{k} + Δ^{2} (δ)] \frac{ν}{ν - 2} - \frac{ν}{π} {[\frac{Γ {(ν - 1) / 2}}{Γ (ν / 2)}]}^{2} Δ^{2} (δ) . \end{matrix}

(A.4)

In order to have a zero mean, we should assume the location parameter

μ = - (ν / π)^{1 / 2}

\frac{Γ ((ν - 1) / 2)}{Γ (ν / 2)} δ

. According to Lemma 1 of (Azzalini and Capitanio, 2003), if

Y

follows

{ST}_{k, ν} (μ, Σ, Δ)

, it can be represented by

Y = μ + ζ^{- 1 / 2} X

(A.5)

where

ζ

follows a Gamma distribution

Γ (ν / 2, ν / 2)

, which is independent of

X

, and

X

follows a

k

-dimensional SN distribution, denoted by

{SN}_{k} (0, Σ, Δ)

. It follows from equation (A.5) that

Y | ζ \sim {SN}_{k} (μ, ζ^{- 1} Σ, ζ^{- 1 / 2} Δ)

. Following studies by (Azzalini and Genton, 2008), the SN distribution of

Y

, conditional on

ζ

, has a convenient stochastic representation as

Y = μ + ζ^{- 1 / 2} Δ | X_{0} | + ζ^{- 1 / 2} Σ^{1 / 2} X,

(A.6)

where

X_{0}

and

X

are two independent

N_{k} (0, I_{k})

random vectors. Note that the expression (A.6) provides a convenience device for random number generation and for implementation purpose. Let

w = ζ^{- 1 / 2} | X_{0} |

; then

w

, conditional on

ζ

, follows a

k

-dimensional normal distribution

N_{k} (0, ζ^{- 1} I_{k})

truncated in the space

w > 0

(i.e., the half-normal distribution). Thus, following (Jara et al., 2008), a hierarchical representation of equation (A.6) is given by

Y | w, ζ \sim N_{k} (μ + Δ w, ζ^{- 1} Σ), w | ζ \sim N_{k} (0, ζ^{- 1} I_{k}) I (w > 0), ζ \sim Γ (ν / 2, ν / 2),

(A.7)

Note that the ST distribution presented in (A.7) can be reduced to the following three special cases: (a) as $ν \to \infty$ and $ζ \to 1$ with probability 1 (i.e., the last distributional specification is omitted), then the hierarchical expression (A.7) becomes an SN distribution ${SN}_{k} (μ, Σ, Δ)$ ; (b) as $Δ = 0$ , then the hierarchical expression (A.7) is a standard multivariate $t$ -distribution; (c) as $ν \to \infty$ , $ζ \to 1$ with probability 1, and $Δ = 0$ , then the hierarchical expression (A.7) reverts to a standard multivariate normal distribution.

8 Appendix B. Full conditional distributions and MCMC implementation

B.1 Full conditional distributions

The full conditional distribution of the coefficients of the mixed-effects submodel (2.4) $ζ$ is

\begin{matrix} P (ζ | .) \propto \prod_{i = 1}^{n} exp {- 0.5 (ζ - τ_{1})^{T} Λ_{1}^{- 1} (ζ - τ_{1})} \times \\ {\{1 + \frac{(y_{i} - X_{i} ζ - Z_{i} b_{i} - δ_{e} [w_{e_{i}} - J (ν + n_{i}) 1_{n_{i}}])^{T} (y_{i} - X_{i} ζ - Z_{i} b_{i} - δ_{e} [w_{e_{i}} - J (ν + n_{i}) 1_{n_{i}}])}{(ν + n_{i}) (ω_{i} σ^{2})}\}}^{- (ν + 2 n_{i}) / 2} . \end{matrix}

(B.1)

The full conditional distribution of the random effects of the mixed-effects submodel (2.4) $b_{i}$ is

\begin{matrix} P (b_{i} | .) \propto exp {- 0.5 b_{i}^{T} Σ^{- 1} b_{i}} \times \\ {\{1 + \frac{(y_{i} - X_{i} ζ - Z_{i} b_{i} - δ_{e} [w_{e_{i}} - J (ν + n_{i}) 1_{n_{i}}])^{T} (y_{i} - X_{i} ζ - Z_{i} b_{i} - δ_{e} [w_{e_{i}} - J (ν + n_{i}) 1_{n_{i}}])}{(ν + n_{i}) (ω_{i} σ^{2})}\}}^{- (ν + 2 n_{i}) / 2} \times \\ \prod_{k = 1}^{K} [\sum_{j = 1}^{J} λ_{0 k, j} I (t_{k, j - 1} < T_{i} \leq t_{k, j}) exp {ψ_{k} z_{i} + b_{i}^{T} c_{k}}]^{I (d_{i} = k)} exp {- \prod_{k = 1}^{K} [exp {ψ_{k} z_{i} + b_{i}^{T} c_{k}}] \\ \sum_{j = 1}^{J} λ_{0 k, j} T_{ij}}, \end{matrix}

(B.2)

where

T_{ij} = \min (T_{i}, t_{k, j}) - \min (T_{i}, t_{k, j - 1})

denote the intervals for the piecewise constant baseline hazard.

The full conditional distribution of the variance of the mixed-effects submodel (2.4) $σ^{2}$ is

\begin{matrix} P (σ^{2} | .) \propto \prod_{i = 1}^{n} σ^{2 (- γ_{σ 1} - 1)} exp (- γ_{σ 2} / σ^{2}) (ω_{i} σ^{2})^{- n_{i} / 2} \times \\ {\{1 + \frac{(y_{i} - X_{i} ζ - Z_{i} b_{i} - δ_{e} [w_{e_{i}} - J (ν + n_{i}) 1_{n_{i}}])^{T} (y_{i} - X_{i} ζ - Z_{i} b_{i} - δ_{e} [w_{e_{i}} - J (ν + n_{i}) 1_{n_{i}}])}{(ν + n_{i}) (ω_{i} σ^{2})}\}}^{- (ν + 2 n_{i}) / 2} . \end{matrix}

(B.3)

The full conditional distribution of the covariance–variance of the random effects of the mixed-effects submodel (2.4) $Σ$ is

\begin{matrix} P (Σ | .) & \propto & \prod_{i = 1}^{n} Σ^{- \frac{1}{2}} exp {- \frac{1}{2} b_{i}^{T} Σ^{- 1} b_{i}} Σ^{- \frac{n + ρ + 1}{2}} exp \{- \frac{1}{2} tr ({Ω Σ}^{- 1})\} \\ \propto & \prod_{i = 1}^{n} Σ^{- \frac{n + ρ + 2}{2}} exp {- \frac{1}{2} tr (b_{i}^{T} b_{i} Σ^{- 1} + {Ω Σ}^{- 1})} . \end{matrix}

(B.4)

Thus,

[Σ | .] \sim IW (b_{i}^{T} b_{i} + Ω, ρ + 1)

The full conditional distribution of the skewness parameter of the mixed-effects submodel (2.4) $δ_{e}$ is

\begin{matrix} P (δ_{e} | .) \propto \prod_{i = 1}^{n} exp (- \frac{1}{2} \frac{δ_{e}^{2}}{γ}) \times \\ {\{1 + \frac{(y_{i} - X_{i} ζ - Z_{i} b_{i} - δ_{e} [w_{e_{i}} - J (ν + n_{i}) 1_{n_{i}}])^{T} (y_{i} - X_{i} ζ - Z_{i} b_{i} - δ_{e} [w_{e_{i}} - J (ν + n_{i}) 1_{n_{i}}])}{(ν + n_{i}) (ω_{i} σ^{2})}\}}^{- (ν + 2 n_{i}) / 2} . \end{matrix}

(B.5)

The full conditional distribution of the degree of freedom of the mixed-effects submodel (2.4) $ν$ is

\begin{matrix} P (ν | .) \propto \prod_{i = 1}^{n} \frac{Γ [(ν + 2 n_{i}) / 2]}{Γ [(ν + n_{i}) / 2] [(ν + n_{i}) π]^{n_{i} / 2}} (ω_{i} σ^{2})^{- n_{i} / 2} \times \\ {\{1 + \frac{(y_{i} - X_{i} ζ - Z_{i} b_{i} - δ_{e} [w_{e_{i}} - J (ν + n_{i}) 1_{n_{i}}])^{T} (y_{i} - X_{i} ζ - Z_{i} b_{i} - δ_{e} [w_{e_{i}} - J (ν + n_{i}) 1_{n_{i}}])}{(ν + n_{i}) (ω_{i} σ^{2})}\}}^{- (ν + 2 n_{i}) / 2} \\ ν_{0} exp (- ν_{0} ν) I (ν > 3) . \end{matrix}

(B.6)

The full conditional distribution of coefficients of the survival submodel (2.5) $ϒ$ is

\begin{matrix} P (ϒ | .) \propto \prod_{i = 1}^{n} exp {- \frac{1}{2} (ϒ - τ_{2})^{T} Λ_{2}^{- 1} (ϒ - τ_{2})} \times \\ \prod_{k = 1}^{K} [\sum_{j = 1}^{J} λ_{0 k, j} I (t_{k, j - 1} < T_{i} \leq t_{k, j}) exp {ψ_{k} z_{i} + b_{i}^{T} c_{k}}]^{I (d_{i} = k)} exp {- \prod_{k = 1}^{K} [exp {ψ_{k} z_{i} + b_{i}^{T} c_{k}}] \\ \sum_{j = 1}^{J} λ_{0 k, j} T_{ij}} . \end{matrix}

(B.7)

The full conditional distribution of the baseline hazards of the survival submodel (2.5) $λ_{0 k, j}$ is

\begin{matrix} P (λ_{0 k, j} | .) \propto \prod_{i = 1}^{n} \prod_{k = 1}^{K} λ_{0 k, j}^{γ_{k 1} - 1} exp (- γ_{k 2} λ_{0 k, j}) [\sum_{j = 1}^{J} λ_{0 k, j} I (t_{k, j - 1} < T_{i} \leq t_{k, j}) \times \\ exp {ψ_{k} z_{i} + b_{i}^{T} c_{k}}]^{I (d_{i} = k)} exp {- \prod_{k = 1}^{K} [exp {ψ_{k} z_{i} + b_{i}^{T} c_{k}}] \sum_{j = 1}^{J} λ_{0 k, j} T_{ij}} . \end{matrix}

(B.8)

The full conditional distribution of the coefficients of the missing data submodel (2.7) $η$ is

\begin{matrix} P (η | .) \propto \prod_{i = 1}^{n} \prod_{j = 1}^{n_{i}} [\frac{e^{η_{0} + η_{1} y (t_{ij})}}{1 + e^{n_{0} + n_{1} y (t_{ij})}}]^{r_{ij}} [\frac{1}{1 + e^{η_{0} + η_{1} y (t_{ij})}}]^{1 - r_{ij}} exp {- \frac{1}{2} (η - τ_{3})^{T} Λ_{3}^{- 1} (η - τ_{3})} . \end{matrix}

(B.9)

B.2 MCMC implementation

Gibbs sampling was combined with M-H sampling for posterior. Particularly, for the parameters Σ, Gibbs sampling was applied since the full conditional distribution is standard. For the other parameters, the M–H algorithm was applied.

References

Arellano-Valle

Genton

(2005) On fundamental skew distributions. Journal of Multivariate Analysis , 96, 93–116.

Azzalini

Capitanio

(2003) On fundamental skew distributions. Journal of Royal Statistical Society, Series B , 65, 367–89.

Azzalini

Genton

(2008) Robust likelihood methods based on the skew-t and related distributions. International Statistical Review , 76, 106–29.

Baker

(1995) Marginal regression for repeated binary data with outcome subject to non-ignorable non-response. Biometrics , 51, 1042–52.

Birmingham

Fitzmaurice

(2002) A pattern-mixture model for longitudinal binary responses with nonignorable non- response. Biometrics , 58, 989–96.

Brown

Ibrahim

(2003) A Bayesian semiparametric joint hierarchical model for longitudinal and survival data. Biometrics , 59, 221–28.

Ekholm

Skinner

(1998) The muscatine children's obesity data reanalysed using pattern mixture models. Applied Statistics , 47, 251–63.

Elashoff

(2008) A joint model for longitudinal measurements and survival data in the presence of multiple failure types. Biometrics , 64, 762–71. doi 10.1111/j.1541-0420.2007.00952.x

Eubank

(1999) Nonparametric Regression and Spline Smoothing . New York: Dekker.

10.

Faucett

Thomas

(1996) Simultaneously modelling censored survival data and repeatedly measured covariates: A Gibbs sampling approach. Statistics in Medicine , 15, 1663–85.

11.

Fitzmaurice

Laird

(2000) Generalized linear mixture models for handling non-ignorable dropouts in longitudinal studies. Biostatistics , 1, 141–56.

12.

Fitzmaurice

Laird

Zahner

GEP

(1996) Multivariate logistic models for incomplete binary responses. Journal of the American Statistical Association , 91, 99–108.

13.

Follmann

(1995) An approximate generalized linear model with random effects for informative missing data. Biometrics , 51, 151–68.

14.

Gelman

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

15.

Gelman

Rubin

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

16.

Gelman

Carlin

Stern

Rubin

(2003) Bayesian Data Analysis . London: Chapman and Hall.

17.

Lin

(2010) Robust linear mixed models using the skew-t distribution with application to schizophrenia data. Biometrical Journal , 52, 449–69.

18.

(2009) A Bayesian approach to joint analysis of longitudinal measurements and competing risks failure time data. Statistics in Medicine , 28, 1601–19.

19.

Huang

Elashoff

Pan

(2011) A general joint model for longitudinal measurements and competing risks survival data with heterogeneous random effects. Lifetime Data Analysis , 17, 80–100.

20.

Ibrahim

Chen

M-H

Lipsitz

(2001) Missing responses in generalised linear mixed models when the missing data mechanism is nonignorable. Biometrika , 88, 551–64.

21.

Jara

Quintana

Martin

(2008) Linear mixed models with skew-elliptical distributions: A bayesian approach. Computational Statistics and Data Analysis , 52, 5033–45.

22.

Lachos

Bandyopadhyay

Dey

(2011) Linear and nonlinear mixed-effects models for censored hiv viral loads using normal/independent distributions. Biometrics , 67, 1594–1604.

23.

Elashoff

Tseng

(2012) Joint analysis of bivariate longitudinal ordinal outcomes and competing risks survival times with nonparametric distributions for random effects. Statistics in Medicine , 31, 1707–21.

24.

Little

Rubin

(2002) Statistical Analysis with Missing Data , 2nd edition. Hoboken, New Jersey: John Wiley & Sons.

25.

Lunn

Thomas

Best

Spiegelhalter

(2000) WinBUGS—A Bayesian modelling framework: Concepts, structure, and extensibility. Statistics and Computing , 10, 325–37.

26.

Molenberghs

Kenward

Lesaffre

(1997) The analysis of longitudinal ordinal data with nonrandom drop-out. Biometrika , 84, 33–44.

27.

Rotnitzky

Robins

Scharfstein

(1998) Semiparametric regression for repeated outcomes with nonignorable nonresponse. Journal of the American Statistical Association , 93, 1321–39.

28.

Sahu

Dey

Branco

(2003) A new class of multivariate skew distributions with applications to bayesian regression models. The Canadian Journal of Statistics , 31, 129–50.

29.

Song

Davidian

Tsiatis

(2002) A semiparametric likelihood approach to joint modeling of longitudinal and time-to-event data. Biometrics , 58, 742–53.

30.

Spiegelhalter

Best

Carlin

Van der Linde

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

31.

(2002) A joint model for nonlinear mixed-effects models with censoring and covariates measured with error. Journal of the American Statistical Association , 97, 955–64.

Bayesian varying coefficient mixed-effects joint models with asymmetry and missingness

Abstract

Abstract:

Keywords

1 Introduction

Figure 1:

2.1 Motivating data

Figure 2:

Histogram of residual of OLS estimates

4.1 Model specification

Figure 4:

Diagnostics of model fitting for the three models. Top panel: fitted versus observed values of log 10 viral load; bottom panel: Q–Q plots for residuals in three models

A summary of estimated PM for the population parameters in longitudinal model (2.4), cause-specific hazard model (4.1), missing data model (2.7) and the corresponding SD, lower limit ( L CI ) and upper limit ( U CI ) of 95 % equal-tail CI

The estimated PM for parameters of scale, skewness and the corresponding SD and lower limit ( L CI ) and upper limit ( U CI ) of 95 % equal-tail CI as well as DIC and EPD values for MACS study

Figure 5:

The population estimating curve of varying-coefficients β 0 ( t ) , β 1 ( t ) for longitudinal model (2.1) based on the best selected Model ST. The estimates (solid curve) along with the 95% CI (dashed curves) are presented

Illustration of four represented individual estimating curve (dashed lines) of β 1 i ( t ) based on Model ST. Population curves (solid lines) are shown for comparison

Table 3:

Simulation results based on 200 simulated dataset. The true values, bias (%) and coverage rate (%) of 95% CI are presented for each parameter

Table 4:

Sensitivity analysis for different sets of hyper-parameters based on Model ST

B.1 Full conditional distributions

B.2 MCMC implementation

References

Diagnostics of model fitting for the three models. Top panel: fitted versus observed values of log $_{10}$ viral load; bottom panel: Q–Q plots for residuals in three models

A summary of estimated PM for the population parameters in longitudinal model (2.4), cause-specific hazard model (4.1), missing data model (2.7) and the corresponding SD, lower limit ( $L_{CI}$ ) and upper limit ( $U_{CI}$ ) of $95 %$ equal-tail CI

The estimated PM for parameters of scale, skewness and the corresponding SD and lower limit ( $L_{CI}$ ) and upper limit ( $U_{CI}$ ) of $95 %$ equal-tail CI as well as DIC and EPD values for MACS study

The population estimating curve of varying-coefficients $β_{0} (t)$ , $β_{1} (t)$ for longitudinal model (2.1) based on the best selected Model ST. The estimates (solid curve) along with the 95% CI (dashed curves) are presented

Illustration of four represented individual estimating curve (dashed lines) of $β_{1 i} (t)$ based on Model ST. Population curves (solid lines) are shown for comparison