Bayesian inference on longitudinal-survival data with multiple features

Abstract

The modeling of longitudinal and survival data is an active research area. Most of researches focus on improving the estimating efficiency but ignore many data features frequently encountered in practice. In this article, we develop a joint model that concurrently accounting for longitudinal-survival data with multiple features. Specifically, our joint model handles skewness, limit of detection, missingness and measurement errors in covariates which are typical observed in the collection of longitudinal-survival data from many studies. We employ a Bayesian approach for making inference on the joint model. The proposed model and method are applied to an AIDS study. A few alternative models under different conditions are compared. Some interesting results are reported. Simulation studies are conducted to assess the performance of the proposed methods.

Appendix: Multivariate skew distributions

Different versions of multivariate skew distributions have been introduced in the literature (Sahu et al. 2003; Azzalini and Capitanio 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 skew-normal (SN) and skew-t (ST) distributions as special cases. A k-dimensional random vector ${\varvec{Y}}$ follows a k-variate skew-elliptical (SE) distribution if its probability density function (pdf) is given by

$$\begin{aligned} f({\varvec{y}}|{\varvec{\mu }},{\varvec{\Sigma }},{\varvec{\Delta }};m^{(k)}_{\nu })= 2^k f({\varvec{y}}|{\varvec{\mu }}, {\varvec{A}};m^{(k)}_{\nu })P({\varvec{V}}>\mathbf 0 ), \end{aligned}$$

(14)

where ${\varvec{A}}={\varvec{\Sigma }}+{\varvec{\Delta }}^2$, ${\varvec{\mu }}$ is a location parameter vector, ${\varvec{\Sigma }}$ is a $k \times k$ positive (diagonal) covariance matrix, ${\varvec{\Delta }}=\text {diag}(\delta _1, \delta _2,\ldots , \delta _k)$ is a $k \times k$ skewness matrix with the skewness parameter vector ${\varvec{\delta }}=(\delta _1,\delta _2,\ldots ,\delta _k)^T$; ${\varvec{V}}$ follows the elliptical distribution $El({\varvec{\Delta }}A^{-1}({\varvec{y}}-{\varvec{\mu }}), {\varvec{I}}_{k}-{\varvec{\Delta }}A^{-1}{\varvec{\Delta }}; m^{(k)}_{\nu })$ and the density generator function $m^{(k)}_{\nu }(\zeta )=\frac{\Gamma (k/2)}{\pi ^{k/2}}\frac{m_{\nu }(\zeta )}{\int _0^{\infty }r^{k/2-1}m_{\nu }(\zeta )dr}$, with $m_{\nu }(\zeta )$ being a function such that $\int _0^{\infty }r^{k/2-1}m_{\nu }(\zeta )dr$ exists. The function $m_{\nu }(\zeta )$ provides the kernel of the original elliptical density and may depend on the parameter $\nu$. This SE distribution is denoted by $SE({\varvec{\mu }},{\varvec{\Sigma }},{\varvec{\Delta }};m^{(k)})$. Two examples of $m_{\nu }(\zeta )$, leading to important special cases used throughout the paper, are $m_{\nu }(\zeta )=\exp (-\zeta /2)$ and $m_{\nu }(\zeta )=(1+\zeta /\nu )^{-(\nu +k)/2}$, where $\nu >0$. These two expressions lead to the multivariate SN and ST distributions, respectively. In the latter case, $\nu$ 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 Reference Sahu et al. (2003). Assume a k-dimensional random vector ${\varvec{Y}}$ follows a k variate ST distribution with location vector ${\varvec{\mu }}$, $k \times k$ positive (diagonal) covariance matrix ${\varvec{\Sigma }}$ and $k \times k$ skewness matrix ${\varvec{\Delta }}=\text {diag}(\delta _1, \delta _2,\ldots , \delta _k)$ or the degrees of freedom $\nu$.

A k-dimensional random vector ${\varvec{Y}}$ follows an m-variate ST distribution if its probability density function (pdf) is given by

$$\begin{aligned} f({\varvec{y}}|{\varvec{\mu }},{\varvec{\Sigma }},{\varvec{\Delta }},\nu )= 2^k t_{k,\nu }({\varvec{y}}|{\varvec{\mu }}, {\varvec{A}})P({\varvec{V}}>\mathbf 0 ), \end{aligned}$$

(15)

we denote the k-variate t distribution with parameters ${\varvec{\mu }}$, ${\varvec{A}}$ and degrees of freedom $\nu$ by $t_{k,\nu }({\varvec{\mu }}, {\varvec{A}})$ and the corresponding pdf by $t_{k,\nu }({\varvec{y}}|{\varvec{\mu }}, {\varvec{A}})$ henceforth, ${\varvec{V}}$ follows the t distribution $t_{k,\nu +k}$. We denote this distribution by $ST_{k,\nu }({\varvec{\mu }},{\varvec{\Sigma }},{\varvec{\Delta }})$. In particular, when ${\varvec{\Sigma }}=\sigma ^2 {\varvec{I}}_k$ and ${\varvec{\Delta }}=\delta {\varvec{I}}_k$, the Eq. (15) simplifies to

$$\begin{aligned} f({\varvec{y}}|{\varvec{\mu }},\sigma ^2,\delta ,\nu )= 2^k (\sigma ^2+\delta ^2)^{-k/2}\frac{\Gamma ((\nu +k)/2)}{\Gamma (\nu /2)(\nu \pi )^{k/2}} \left\{ 1+\frac{({\varvec{y}}-{\varvec{\mu }})^T({\varvec{y}}-{\varvec{\mu }})}{\nu (\sigma ^2+\delta ^2)} \right\} ^{-(\nu +k)/2} \\ \times \, T_{k,\nu +k} \left[ \left\{ \frac{\nu +(\sigma ^2+\delta ^2)^{-1}({\varvec{y}}-{\varvec{\mu }})^T({\varvec{y}}-{\varvec{\mu }})}{\nu +k}\right\} ^{-1/2}\frac{\delta ({\varvec{y}}-{\varvec{\mu }})}{\sigma \sqrt{\sigma ^2+\delta ^2}}\right] , \end{aligned}$$

where $T_{k,\nu +k}(\cdot )$ denotes the cumulative distribution function (cdf) of $t_{k,\nu +k}(\mathbf 0 ,{\varvec{I}}_k)$. However, unlike in the SN distribution, the ST density can not be written as the product of univariate ST densities. Here ${\varvec{Y}}$ is dependent but uncorrelated. It is noted that when ${\varvec{\delta }}=\mathbf 0$, the ST distribution reduces to usual the t-distribution. It can be shown that the mean and covariance matrix of the ST distribution $ST_{k,\nu }({\varvec{\mu }},\sigma ^2 {\varvec{I}}_k,{\varvec{\Delta }})$ are given by

$$\begin{aligned} E({\varvec{Y}})= & {} {\varvec{\mu }}+(\nu /\pi )^{1/2}\frac{\Gamma ((\nu -1)/2)}{\Gamma (\nu /2)}{\varvec{\delta }}, \nonumber \\ \text {cov}({\varvec{Y}})= & {} \left[ \sigma ^2{\varvec{I}}_k+{\varvec{\Delta }}^2({\varvec{\delta }})\right] \frac{\nu }{\nu -2}-\frac{\nu }{\pi }\left[ \frac{\Gamma \{(\nu -1)/2\}}{\Gamma (\nu /2)}\right] ^2{\varvec{\Delta }}^2({\varvec{\delta }}). \end{aligned}$$

(16)

In order to have a zero mean, we should assume the location parameter ${\varvec{\mu }}=-(\nu /\pi )^{1/2}$ $\frac{\Gamma ((\nu -1)/2)}{\Gamma (\nu /2)}{\varvec{\delta }}$. According to Lemma 1 of Azzalini and Capitanio (2003), if ${\varvec{Y}}$ follows $ST_{k,\nu }({\varvec{\mu }},{\varvec{\Sigma }},{\varvec{\Delta }})$, it can be represented by

$$\begin{aligned} {\varvec{Y}}={\varvec{\mu }}+\zeta ^{-1/2}{\varvec{X}}\end{aligned}$$

(17)

where $\zeta$ follows a Gamma distribution $\Gamma (\nu /2,\nu /2)$, which is independent of ${\varvec{X}}$, and ${\varvec{X}}$ follows a k-dimensional skew-normal (SN) distribution, denoted by $SN_k(\mathbf 0 ,{\varvec{\Sigma }},{\varvec{\Delta }})$. It follows from (17) that ${\varvec{Y}}|\zeta \sim SN_k({\varvec{\mu }}, \zeta ^{-1}{\varvec{\Sigma }},\zeta ^{-1/2}{\varvec{\Delta }})$. Following studies by Azzalini and Genton (2008), the SN distribution of ${\varvec{Y}}$, conditional on $\zeta$, has a convenient stochastic representation as

$$\begin{aligned} {\varvec{Y}}={\varvec{\mu }}+\zeta ^{-1/2}{\varvec{\Delta }}|{\varvec{X}}_0|+\zeta ^{-1/2}{\varvec{\Sigma }}^{1/2}{\varvec{X}}, \end{aligned}$$

(18)

where ${\varvec{X}}_0$ and ${\varvec{X}}$ are two independent $N_k(\mathbf 0 ,{\varvec{I}}_k)$ random vectors. Note that the expression (18) provides a convenience device for random number generation and for implementation purpose. Let ${\varvec{w}}=\zeta ^{-1/2}|{\varvec{X}}_0|$; then ${\varvec{w}}$, conditional on $\zeta$, follows a k-dimensional normal distribution $N_k(\mathbf 0 , \zeta ^{-1}{\varvec{I}}_k)$ truncated in the space ${\varvec{w}}>\mathbf 0$ (i.e., the half-normal distribution). Thus, following (Jara et al. 2008), a hierarchical representation of (18) is given by

$$\begin{aligned} {\varvec{Y}}|{\varvec{w}},\zeta \sim N_k({\varvec{\mu }}+{\varvec{\Delta }}{\varvec{w}}, \zeta ^{-1}{\varvec{\Sigma }}),\; {\varvec{w}}|\zeta \sim N_k(\mathbf 0 ,\zeta ^{-1}{\varvec{I}}_k){\varvec{I}}({\varvec{w}}>\mathbf 0 ),\; \zeta \sim \Gamma (\nu /2, \nu /2), \end{aligned}$$

(19)

Note that the ST distribution presented in (19) can be reduced to the following three special cases: (i) as $\nu \rightarrow \infty$ and $\zeta \rightarrow 1$ with probability 1 (i.e., the last distributional specification is omitted), then the hierarchical expression (19) becomes an SN distribution $SN_k({\varvec{\mu }},{\varvec{\Sigma }},{\varvec{\Delta }})$; (ii) as ${\varvec{\Delta }}=\mathbf 0$, then the hierarchical expression (19) is a standard multivariate t-distribution; (iii) as $\nu \rightarrow \infty$, $\zeta \rightarrow 1$ with probability 1, and ${\varvec{\Delta }}=\mathbf 0$, then the hierarchical expression (19) reverts to a standard multivariate normal distribution.