DOI: xxx/xxxx

## RESEARCH ARTICLE

## Fast variational Bayesian inference for correlated survival data:

## an application to invasive mechanical ventilation duration analysis

Chengqian Xian 1 | Camila P.E. de Souza 1 | Wenqing He 1 | Felipe F. Rodrigues 1,2 | Renfang Tian 2

1 Department of Statistical and Actuarial Sciences, Western University, London, Ontario, Canada

## Abstract

2 School of Management, Economics, and

Mathematics, King's University College at Western University, London, Ontario, Canada

## Correspondence

Chengqian Xian, Department of Statistical and Actuarial Sciences, Western University, 1151 Richmond Street, London, Ontario, Canada, N6A 5B7. Email: cxian3@uwo.ca

Correlated survival data are prevalent in various clinical settings and have been extensively discussed in literature. One of the most common types of correlated survival data is clustered survival data, where the survival times from individuals in a cluster are associated. Our study is motivated by invasive mechanical ventilation data from different intensive care units (ICUs) in Ontario, Canada, forming multiple clusters. The survival times from patients within the same ICU cluster are correlated. To address this association, we introduce a shared frailty log-logistic accelerated failure time model that accounts for intra-cluster correlation through a cluster-specific random intercept. We present a novel, fast variational Bayes (VB) algorithm for parameter inference and evaluate its performance using simulation studies varying the number of clusters and their sizes. We further compare the performance of our proposed VB algorithm with the h-likelihood method and a Markov Chain Monte Carlo (MCMC) algorithm. The proposed algorithm delivers satisfactory results and demonstrates computational efficiency over the MCMC algorithm. We apply our method to the ICU ventilation data from Ontario to investigate the ICU site random effect on ventilation duration.

## KEYWORDS:

Variational inference, clustered survival data, random effects, ICU ventilation

## 1 INTRODUCTION

Correlated survival data are commonplace in various research contexts and have been extensively studied in the literature 1,2,3 . One such context of correlated survival data is observed in clustered survival data, which are derived from multiple entities such as families or hospitals. Within each cluster, survival data exhibit correlation because of the shared environmental factors, and random effects are often introduced to capture the shared characteristics within the cluster. Conditional on the random effects

(frailty), the survival times within a cluster can be assumed independent, which leads to a shared frailty model to account for cluster-level uncertainty 4,5,6 .

Our research is motivated by data from intensive care units (ICUs) across multiple clinical centers in Ontario, Canada. The ICU data were provided by the Critical Care Information System (CCIS) Ontario database. In Canadian ICUs, invasive mechanical ventilation is prevalent, with approximately one third of patients requiring ventilation during their ICU stay 7 . The duration of ventilation has significant implications for clinical outcomes and is associated with an increased risk of complications 7 . Analysis on the ventilation duration for ICU patients utilizing the associated risk factors, such as patient categories (e.g., medical or surgical), admission diagnoses, and patient severity, can help determine the number of beds with ventilators and therefore, support capacity planning and effective clinical resource management. The severity of patients in ICU is often assessed using the Multiple Organ Dysfunction Score (MODS) 8 that evaluates organ function. Kobara et al. 9 conducted a study on ventilation duration for ICU patients from Ontario using a survival analysis framework. In their study, patient ventilation duration is considered as a time-to-event (survival) outcome and parametric accelerated failure time (AFT) models are applied to predict the ventilation time. If a patient is transferred to another facility without subsequent follow-up information, the ventilation time is considered right censored. Kobara et al. found that the log-logistic AFT model well describes the association between risk factors and patients' ventilation duration in Ontario ICUs. However, Kobara et al. did not consider a possible correlation of ventilation duration times among patients within the same ICU. Previous studies have indicated that patient outcomes within the same ICU site may be correlated, and ignoring this hierarchical structure can result in flawed prediction models 10,11 .

As aforementioned, shared frailty can be used to accommodate cluster-level association among patients within the cluster. Gorfine and Zucker 6 recently provided a comprehensive review of shared frailty methods for complex survival data. To incorporate risk factors in survival data, the Cox proportional hazards (PH) model with a multiplicative shared frailty on the hazard rates has been widely developed 6 . As an alternative to the Cox PH model, AFT models offer an intuitive interpretation of covariate effects on survival time. Lambert et al. 12 developed shared frailty AFT models with different parametric distributions assumed for the survival time and conducted maximum likelihood estimation by integrating out the unobserved frailties. They empirically found that the choice of distribution for the shared frailty is not critical, recommending the normal distribution. Do Ha et al. 13,14 considered a shared frailty log-normal AFT model and used the hierarchical-likelihood (h-likelihood) approach to estimate the model parameters. The h-likelihood method obtains fixed effects estimates by maximizing the h-likelihood, while utilizing the restricted maximum likelihood estimate for the estimation of the variance of the shared frailty 14 . An R package, frailtyHL 15 has been developed to implement the h-likelihood estimation in the shared frailty log-normal AFT model. Building on this, Park and Do Ha 16 introduced a penalized variable selection technique for the shared frailty log-normal AFT model. Additionally, Zhou et al. 17 developed a Markov chain Monte Carlo (MCMC)-based algorithm, called survregBayes , to estimate

the shared frailty AFT model, and an R package, spBayesSurv 18 , is available for its implementation. To our knowledge, there is no work on variational Bayesian inference for shared frailty AFT models.

As an alternative to MCMC methods, which are the gold standard for obtaining posterior distributions under a Bayesian framework, variational inference (VI) has gained popularity due to its favourable results and lower computational cost than MCMC. Recently, several types of VI algorithms have been developed, including mean-field VI 19 , stochastic VI 20 , and blackbox VI 21 . A special case of mean-field VI, called mean-field variational Bayes (VB), arises when the Kullback-Leibler (KL) divergence is utilized to quantify the dissimilarity between exact and approximated posterior distributions. In addition, the approximated posterior distribution, referred to as the variational posterior, is assumed to belong to a mean-field variational family. Under the mean-field VB framework, the solutions to minimizing the KL divergence can be obtained by utilizing the coordinate ascent algorithm 19,22 . Mean-field VB has been widely applied to regression models such as the generalized additive model 23 , nonparametric regression with measurement error 24 , count response semiparametric regression 25 , high-dimensional linear regression 26 , multilevel regression modelling 27 , B-spline regression mixture model for functional data clustering 28 , basis selection for functional data representation 29 , among others. Applications of other types of VI can be found in a comprehensive review of VI from a statistical perspective by Blei et al. 30 .

In this study, we propose a shared frailty log-logistic AFT model to account for the correlation among patient ventilation durations within the ICU sites. We develop a novel and fast mean-field VB algorithm to infer the model parameters. By applying the piece-wise approximation techniques proposed by Xian et al. 31 to avoid intractable calculations, we obtain closed-form posterior distributions. We conduct extensive simulation studies with various numbers of clusters and cluster sizes to evaluate the performance of the proposed method, and compare the performance of our VB algorithm with the h-likelihood method 14 and the MCMC-based algorithm survregBayes 18 . Finally, we apply our methodology to investigate ventilation duration for ICU patients using the same dataset as Kobara et al. 9 . This study was approved by the Research Ethics Review Committee at King's University College at Western University.

The remainder of the paper is organized as follows. Section 2 presents the log-logistic AFT model with a shared frailty under the Bayesian framework. We introduce our proposed VB algorithm in Section 3. In Section 4, simulation studies are conducted to evaluate the performance of our VB algorithm. Section 5 illustrates the application of the proposed method to the ICU ventilation duration data. A discussion is provided in Section 6.

## 2 BAYESIAN LOG-LOGISTIC AFT MODEL WITH A SHARED FRAILTY

Let 𝑇 𝑖𝑗 and 𝐶𝑖𝑗 be the survival and censoring times, respectively, of the 𝑗 𝑡ℎ subject from the 𝑖 𝑡ℎ group (i.e., cluster) in a sample, 𝑖 = 1 , ..., 𝐾 and 𝑗 = 1 , ..., 𝑛 𝑖 . Let 𝑡 𝑖𝑗 = min( 𝑇 𝑖𝑗 , 𝐶 𝑖𝑗 ) and 𝛿 𝑖𝑗 = ✶ ( 𝑇 𝑖𝑗 ≤ 𝐶𝑖𝑗 ) be the subject's observed time and the indicator for

right censoring, respectively. We consider a log-logistic AFT model with shared frailty (a random intercept) specified as follows:

$$\log ( T _ { i j } ) = \gamma _ { i } + X _ { i j } ^ { T } \, \beta + b \epsilon _ { i j }$$

where 𝐗 𝑖𝑗 is a column vector with length 𝑝, 𝑝 ≥ 2 , containing 𝑝 - 1 fixed effects (covariates) and a constant one to incorporate the constant intercept (i.e., 𝐗 𝑖𝑗 = (1 , 𝑥 𝑖𝑗 1 , ..., 𝑥 𝑖𝑗 ( 𝑝 -1) ) 𝑇 ), 𝜷 is the corresponding vector of coefficients for the fixed effects, where 𝛾 𝑖 is a random intercept for the 𝑖 𝑡ℎ cluster and 𝜖 𝑖𝑗 is a random variable following a standard logistic distribution, and 𝑏 is a scale parameter. We further assume that 𝛾 𝑖 iid ∼ 𝑁 (0 , 𝜎 2 𝛾 ) representing the discrepancy between clusters with iid denoting identically and independent distributed. The survival time 𝑇 𝑖𝑗 and censoring time 𝐶𝑖𝑗 are assumed independent given the covariates 𝐗 𝑖𝑗 . Our model follows a structure similar to one presented by Robison (1991) 32 and Nolan et al. 33 for a Gaussian linear mixed effect model.

In our proposed Model (1), we incorporate the unknown and unobserved shared risk through cluster-specific random intercepts and estimate the model parameters, 𝜷 , 𝑏 , and 𝜎 2 𝛾 using a Bayesian framework by further assuming the following independent marginal prior distributions:

- · 𝜷 ∼ 𝑁𝑝 ( 𝝁 0 , 𝜎 2 0 I 𝑝 × 𝑝 ) with precision 𝑣 0 = 1∕ 𝜎 2 0 and I 𝑝 × 𝑝 being a 𝑝 × 𝑝 identity matrix
- · 𝑏 ∼ Inverse-Gamma ( 𝛼 , 𝜔 0 0 )
- · 𝛾 𝑖 | 𝜎 2 𝛾 iid ∼ 𝑁 (0 , 𝜎 2 𝛾 ) , 𝑖 = 1 , ..., 𝐾
- · 𝜎 2 𝛾 ∼ Inverse-Gamma ( 𝜆 , 𝜂 0 0 )

where 𝜇 , 𝑣 0 0 , 𝛼 0 , 𝜔 0 , 𝜆 0 and 𝜂 0 are known hyperparameters (parameters of the prior distributions).

## 3 VARIATIONAL BAYES ALGORITHM

In what follows, we outline our methodology for deriving a mean-field VB algorithm for Model (1). We summarize the resulted VB algorithm in Algorithm 1.

Given the observed data 𝐃 ∶= {( 𝑡 𝑖𝑗 , 𝛿 𝑖𝑗 , 𝐗 𝑖𝑗 ) , 𝑖 = 1 , ..., 𝐾, 𝑗 = 1 , ..., 𝑛 𝑖 } , we denote the complete data log-likelihood by log 𝑝 ( 𝐃 , 𝜷 𝜸 , , 𝑏, 𝜎 2 𝛾 ) , where 𝜸 = ( 𝛾 , ..., 𝛾 1 𝐾 ) . Our objective is to derive a VB algorithm to approximate the exact posterior joint distribution of 𝜷 , 𝑏 , 𝜸 and 𝜎 2 𝛾 given the data 𝐃 by maximizing the evidence lower bound (ELBO) defined as

$$\text{ELBO} ( q ) = \mathbb { E } _ { q } [ \log p ( \mathbf D, \beta, \gamma, b, \sigma _ { \gamma } ^ { 2 } ) ] - \mathbb { E } _ { q } [ \log q ( \beta, b, \gamma, \sigma _ { \gamma } ^ { 2 } ) ]$$

where 𝑞 ( 𝜷 , 𝑏, 𝜸 , 𝜎 2 𝛾 ) is the approximated posterior joint distribution, which is also called the variational density, and the expectation is taken with respect to the variational density 30 . We consider the mean-field variational family which assumes

that 𝑞 ( 𝜷 , 𝑏, 𝜸 , 𝜎 2 𝛾 ) = 𝑞 ( 𝜷 ) 𝑞 𝑏 ( ) 𝑞 𝜎 ( 2 𝛾 ) ∏ 𝐾 𝑖 =1 𝑞 𝛾 ( 𝑖 ) , and apply the coordinate ascent variational inference (CAVI) algorithm 19 to obtain each variational component (e.g., 𝑞 ( 𝜷 ) ) in 𝑞 ( 𝜷 , 𝑏, 𝜸 , 𝜎 2 𝛾 ) ). Under our Model (1), the complete data log-likelihood log 𝑝 ( 𝐃 , 𝜷 𝜸 , , 𝑏, 𝜎 2 𝛾 ) can be obtained by

$$\log p ( \mathbf D, \beta, \gamma, b, \sigma _ { \gamma } ^ { 2 } ) = \log p ( \mathbf D \, | \, \beta, \gamma, b ) + \log p ( \beta ) + \log p ( b ) + \sum _ { i = 1 } ^ { \kappa } \log p ( \gamma _ { i } \, | \sigma _ { \gamma } ^ { 2 } ) + \log p ( \sigma _ { \gamma } ^ { 2 } )$$

where

$$\log p ( \mathbf D | \, \beta, \gamma, b ) = - \delta \log b + \sum _ { l = 1 } ^ { K } \sum _ { j = 1 } ^ { n _ { l } } \left [ \delta _ { l j } \frac { y _ { l j } - X _ { l j } ^ { T } \beta - \gamma _ { l } } { b } - ( 1 + \delta _ { l j } ) \log \left \{ 1 + \exp \left ( \frac { y _ { l j } - \mathbf X _ { l j } ^ { T } \beta - \gamma _ { l } } { b } \right ) \right \} \right ] \, \quad \quad ( 4 )$$

with 𝛿 = ∑ 𝐾 𝑖 =1 ∑ 𝑛 𝑖 𝑗 =1 𝛿 𝑖𝑗 being the number of observed uncensored times and 𝑦 𝑖𝑗 = log( 𝑡 𝑖𝑗 ) , 𝑖 = 1 , ..., 𝐾, 𝑗 = 1 , ..., 𝑛 𝑖 , being the log observed survival time.

By maximizing the ELBO, the following solutions are provided by the CAVI algorithm:

$$\log q ^ { * } ( \beta ) \stackrel { + } { \approx } \mathbb { E } _ { - } \beta [ \log p ( \mathbf D \, | \, \beta, \gamma, b ) + \log p ( \beta ) ]$$

$$\log q ^ { * } ( \gamma _ { i } ) \stackrel { + } { \approx } \mathbb { E } _ { - \gamma _ { i } } \left [ \log p ( \mathbf D \, | \, \beta, \gamma, b ) + \sum _ { i = 1 } ^ { K } \log p ( \gamma _ { i } \, | \sigma _ { \gamma } ^ { 2 } ) \right ]$$

$$\log q ^ { * } ( b ) \stackrel { + } { \approx } \mathbb { E } _ { - b } [ \log p ( \mathbf D \, | \, \beta, \gamma, b ) + \log p ( b ) ]$$

$$\log q ^ { * } ( \sigma _ { \gamma } ^ { 2 } ) \stackrel { + } { \approx } \mathbb { E } _ { - \sigma _ { \gamma } ^ { 2 } } \left [ \sum _ { i = 1 } ^ { K } \log p ( \gamma _ { i } \, | \sigma _ { \gamma } ^ { 2 } ) + \log p ( \sigma _ { \gamma } ^ { 2 } ) \right ]$$

where we use + ≈ to denote equality up to a constant additive factor for convenience, and -𝜷 indicates the expectation is taken with respect to the variational density of other latent variables but 𝜷 , same for other solutions. To achieve conjugacy and tractable expectation calculation of log 𝑝 ( 𝐃 | 𝜷 𝜸 , , 𝑏 ) as specified in (4), we apply the proposed method by Xian et al. 31 , piecewise approximations of the function, 𝑓 𝑥 ( ) = log(1+exp( 𝑥 )) , 𝑥 ∈ (-∞ ∞) , , embedded in deriving the update equations for each parameter. As in Nolan et al. 33 and Lee and Wand 34 , we are also interested in the posterior distribution of the random effects, 𝑞 ∗ ( 𝛾 𝑖 ) , 𝑖 = 1 , ..., 𝐾 , in our proposed VB framework.

## 3.1 Update equation for each variational density

The update equations to obtain the optimal variational densities of 𝜷 , 𝛾 𝑖 , 𝑏 and 𝜎 2 𝛾 denoted by 𝑞 ∗ ( 𝜷 ) , 𝑞 ∗ ( 𝛾 𝑖 ) , 𝑞 ∗ ( 𝑏 ) and 𝑞 ∗ ( 𝜎 2 𝛾 ) , respectively, which are the corresponding approximated posterior distributions, are presented as follows within the paper while their derivations are given in the Appendix A. The calculation of expectations in the update equations are given in Section 3.2.

In the update equations, 𝜑𝑖𝑗 , 𝜁 𝑖𝑗 , and 𝜌 𝑖𝑗 represent the piece-wise approximation coefficients proposed by Xian et al. 31 for the 𝑗 𝑡ℎ subject in the 𝑖 𝑡ℎ cluster. Detailed information regarding these coefficients is also provided in the Appendix A.

## (1) Update equation for 𝑞 ∗ ( 𝜷 )

𝑞 ∗ ( 𝜷 ) is an 𝑁𝑝 ( 𝝁 ∗ , Σ ∗ ) where

$$\Sigma ^ { * } = \left [ v _ { 0 } \mathbf I _ { p \times p } + 2 \mathbb { E } _ { q ( b ) } ( \frac { 1 } { b ^ { 2 } } ) \sum _ { i = 1 } ^ { K } \sum _ { j = 1 } ^ { n _ { i } } ( 1 + \delta _ { i j } ) \zeta _ { i j } \mathbf X _ { i j } \mathbf X _ { i j } ^ { T } \right ] ^ { - 1 }$$

and

$$\mu ^ { * } = \left [ \left \{ v _ { 0 } \, \mu _ { 0 } ^ { T } + \sum _ { i = 1 } ^ { \kappa } \sum _ { j = 1 } ^ { n _ { i } } \left ( \left ( \mathbb { E } _ { q ( b ) } ( \frac { 1 } { b } ) \left ( - \delta _ { i _ { J } } + ( 1 + \delta _ { i _ { J } } ) \rho _ { i _ { J } } \right ) + 2 \mathbb { E } _ { q ( b ) } ( \frac { 1 } { b ^ { 2 } } ) ( 1 + \delta _ { i _ { J } } ) \zeta _ { i _ { J } } ( y _ { i _ { J } } - \mathbb { E } _ { q ( y _ { i _ { J } } } ( y _ { i _ { J } } ) ) \right ) \mathbb { X } _ { i _ { J } } ^ { T } \right ) \right \} \Sigma ^ { * } \right ] ^ { T }.$$

(2) Update equation for 𝑞 ∗ ( 𝛾 𝑖 )

𝑞 ∗ ( 𝛾 𝑖 ) is an 𝑁𝑙 ( 𝜏 ∗ 𝑖 , 𝜎 2∗ 𝑖 ) where

$$\tau _ { i } ^ { * } = \sigma _ { i } ^ { 2 * } \sum _ { j = 1 } ^ { n _ { i } } \left [ \mathbb { E } _ { q ( b ) } ( \frac { 1 } { b } ) ( - \delta _ { i j } + ( 1 + \delta _ { i j } ) \rho _ { i j } ) + 2 \mathbb { E } _ { q ( b ) } ( \frac { 1 } { b ^ { 2 } } ) ( 1 + \delta _ { i j } ) \zeta _ { i j } ( y _ { i j } - \mathbf X _ { i j } ^ { T } \mathbb { E } _ { q ( \beta ) } \beta ) \right ]$$

and

$$\sigma _ { i } ^ { 2 * } = \left [ \mathbb { E } _ { q ( \sigma _ { \gamma } ^ { 2 } ) } ( \frac { 1 } { \sigma _ { \gamma } ^ { 2 } } ) + 2 \mathbb { E } _ { q ( b ) } ( \frac { 1 } { b ^ { 2 } } ) \sum _ { j = 1 } ^ { n _ { i } } ( 1 + \delta _ { i j } ) \zeta _ { i j } \right ] ^ { - 1 }.$$

## (3) Update equation for 𝑞 ∗ ( 𝑏 )

𝑞 ∗ ( 𝑏 ) is an Inverse-Gamma ( 𝛼 , 𝜔 ∗ ∗ ) where 𝛼 ∗ = 𝛼 0 + 𝛿 and

$$\omega ^ { * } = \omega _ { 0 } - \sum _ { i = 1 } ^ { K } \sum _ { j = 1 } ^ { n _ { i } } \left ( \delta _ { i j } - ( 1 + \delta _ { i j } ) \varphi _ { i j } \right ) \left ( y _ { i j } - { \mathbf X } _ { i j } ^ { T } { \mathbb { E } } _ { q ( \beta ) } ( \beta ) - { \mathbb { E } } _ { q ( \gamma _ { i } ) } \gamma _ { i } \right ).$$

(4) Update equation for 𝑞 ∗ ( 𝜎 2 𝛾 )

𝑞 ∗ ( 𝜎 2 𝛾 ) is an Inverse-Gamma ( 𝜆 , 𝜂 ∗ ∗ ) where 𝜆 ∗ = 𝜆 0 + 𝐾 ∕2 and

$$\eta ^ { * } = \eta _ { 0 } + \frac { 1 } { 2 } \sum _ { i = 1 } ^ { K } \mathbb { E } _ { q ( \gamma _ { i } ) } \gamma _ { i } ^ { 2 }.$$

## 3.2 ELBO calculation

The ELBO under Model (1) is defined in (2) with the complete-data log-likelihood calculated by (3) and note that 𝑞 ( 𝜷 , 𝑏, 𝜸 , 𝜎 2 𝛾 ) = 𝑞 ( 𝜷 ) 𝑞 𝑏 ( ) 𝑞 𝜎 ( 2 𝛾 ) ∏ 𝐾 𝑖 =1 𝑞 𝛾 ( 𝑖 ) . Let

$$d i f f _ { \beta } = \mathbb { E } _ { q } [ \log p ( \beta ) ] - \mathbb { E } _ { q } [ \log q ( \beta ) ]$$

$$d i f f _ { \gamma } = \mathbb { E } _ { q } \left [ \sum _ { i = 1 } ^ { K } \log p ( \gamma _ { i } \left | \sigma _ { \gamma } ^ { 2 } \right ) \right ] - \mathbb { E } _ { q } \left [ \sum _ { i = 1 } ^ { K } \log q ( \gamma _ { i } ) \right ]$$

$$d i f f _ { b } = \mathbb { E } _ { q } [ \log p ( b ) ] - \mathbb { E } _ { q } [ \log q ( b ) ]$$

$$d i f f _ { \sigma _ { \gamma } ^ { 2 } } = \mathbb { E } _ { q } [ \log p ( \sigma _ { \gamma } ^ { 2 } ) ] - \mathbb { E } _ { q } [ \log q ( \sigma _ { \gamma } ^ { 2 } ) ]$$

Then we can calculate the ELBO as follows with proof details provided in the Appendix B.

$$\text{ELBO} ( q ) = \mathbb { E } _ { q } [ \log p ( \mathbf D \, | \, \beta, \gamma, b ) ] + d i f f _ { \beta } + d i f f _ { \gamma } + d i f f _ { b } + d i f f _ { \sigma _ { \gamma } ^ { 2 } }$$

where

$$\mathbb { E } _ { q } [ \log p ( \mathbf D | \, \beta, \gamma, \, b ) ] \stackrel { + } { \approx } - \delta \mathbb { E } _ { q ( b ) } ( \log b ) + \mathbb { E } _ { q ( b ) } ( \frac { 1 } { b } ) \sum _ { i = 1 } ^ { K } \sum _ { j = 1 } ^ { n _ { i } } \left ( \delta _ { i j } - ( 1 + \delta _ { i j } ) \varphi _ { i j } \right ) \left ( y _ { i j } - X _ { i j } ^ { T } \mathbb { E } _ { q ( \beta ) } ( \beta ) - \mathbb { E } _ { q ( y _ { i } ) } y _ { i } \right )$$

$$d i f f _ { \beta } \stackrel { + } { \approx } - \frac { 1 } { 2 } v _ { 0 } [ \text{trace} ( \Sigma ^ { * } ) + ( \mu ^ { * } - \mu _ { 0 } ) ^ { T } ( \mu ^ { * } - \mu _ { 0 } ) ] + \frac { 1 } { 2 } \log ( | \Sigma ^ { * } | ),$$

$$d i f f _ { \gamma } \stackrel { + } { \approx } - \frac { K } { 2 } \mathbb { E } _ { q ( \sigma _ { \gamma } ^ { 2 } ) } ( \log \sigma _ { \gamma } ^ { 2 } ) - \frac { 1 } { 2 } \mathbb { E } _ { q ( \sigma _ { \gamma } ^ { 2 } ) } ( \frac { 1 } { \sigma _ { \gamma } ^ { 2 } } ) \sum _ { i = 1 } ^ { K } \mathbb { E } _ { q ( \gamma _ { i } ) } ( \gamma _ { i } ^ { 2 } ) - \frac { 1 } { 2 } \sum _ { i = 1 } ^ { K } \left ( \log \sigma _ { i } ^ { 2 * } \right ),$$

$$d i f f _ { b } \stackrel { + } { \approx } ( \alpha ^ { * } - \alpha _ { 0 } ) \mathbb { E } _ { q ( b ) } ( \log b ) + ( \omega ^ { * } - \omega _ { 0 } ) \mathbb { E } _ { q ( b ) } \left ( \frac { 1 } { b } \right ) - \alpha ^ { * } \log \omega ^ { * }, \text{ and}$$

$$d i f f _ { \sigma _ { \gamma } ^ { 2 } } \stackrel { + } { \approx } ( \lambda ^ { * } - \lambda _ { 0 } ) \mathbb { E } _ { q ( \sigma _ { \gamma } ^ { 2 } ) } ( \log \sigma _ { \gamma } ^ { 2 } ) + ( \eta ^ { * } - \eta _ { 0 } ) \mathbb { E } _ { q ( \sigma _ { \gamma } ^ { 2 } ) } ( \frac { 1 } { \sigma _ { \gamma } ^ { 2 } } ) - \lambda ^ { * } \log \eta ^ { * }$$

We now present the calculation of the expectations in the update equations and the ELBO calculation. All the expectations are taken with respect to the approximated variational distributions. Since 𝑞 ∗ ( 𝑏 ) is an Inverse-Gamma ( 𝛼 , 𝜔 ∗ ∗ ) , we have

$$\mathbb { E } _ { q ( b ) } \left ( \frac { 1 } { b } \right ) = \frac { \alpha ^ { * } } { \omega ^ { * } },$$

$$\mathbb { E } _ { q ( b ) } \left ( \frac { 1 } { b ^ { 2 } } \right ) = \mathbb { E } _ { q ( b ) } \left [ \left ( \frac { 1 } { b } \right ) ^ { 2 } \right ] = \text{Var} _ { q ( b ) } \left [ \left ( \frac { 1 } { b } \right ) \right ] + \left [ \mathbb { E } _ { q ( b ) } \left ( \frac { 1 } { b } \right ) \right ] ^ { 2 } = \frac { \alpha ^ { * } } { \omega ^ { * 2 } } + \frac { \alpha ^ { * 2 } } { \omega ^ { * 2 } } = \frac { \alpha ^ { * } + \alpha ^ { * 2 } } { \omega ^ { * 2 } }, \text{ and}$$

$$\mathbb { E } _ { q ( b ) } ( \log b ) = \log ( \omega ^ { * } ) - \psi ( \alpha ^ { * } ),$$

where 𝜓 is the digamma function defined as 𝜓 𝑥 ( ) = 𝑑 𝑑𝑥 log Γ( 𝑥 ) .

Similarly,

$$\mathbb { E } _ { q ( \sigma _ { \gamma } ^ { 2 } ) } \left ( \frac { 1 } { \sigma _ { \gamma } ^ { 2 } } \right ) = \frac { \lambda ^ { * } } { \eta ^ { * } } \quad \text{and} \quad \mathbb { E } _ { q ( \sigma _ { \gamma } ^ { 2 } ) } ( \log \sigma _ { \gamma } ^ { 2 } ) = \log ( \eta ^ { * } ) - \psi ( \lambda ^ { * } ).$$

## Algorithm 1 Variational Bayesian inference of correlated survival data using a shared frailty log-logistic AFT model

Result: posterior distributions of 𝜷 , 𝛾 , 𝑖 𝑖 = 1 , ..., 𝐾 , 𝑏 and 𝜎 2 𝛾 , and their parameters: Σ , 𝝁 , 𝜎 2 𝑖 , 𝜏 , 𝛼, 𝜔, 𝜆, 𝜂 𝑖 Initialization : initialize 𝜔 = 𝜔 0 , 𝝁 = 𝝁 0 , 𝜏 𝑖 = 0 and 𝜂 = 𝜂 0 , set 𝑚 = 0 and ELBO = 0 Calculation : obtain 𝛼 by 𝛼 = 𝛼 0 + 𝛿 with 𝛿 = ∑ 𝐾 𝑖 =1 ∑ 𝑛 𝑖 𝑗 =1 𝛿 𝑖𝑗 and 𝜆 by 𝜆 = 𝜆 0 + 𝐾 ∕2

Data: a sample of independent log observed time 𝑦 𝑖𝑗 , their corresponding covariate vectors 𝐗 𝑖𝑗 and the right censoring indicator 𝛿 𝑖𝑗 , 𝑖 = 1 , ..., 𝐾, 𝑗 = 1 , ..., 𝑛 𝑖 for the 𝑗 𝑡ℎ observation from the 𝑖 𝑡ℎ group; values of hyperparameters: 𝝁 0 , 𝑣 0 , 𝛼 0 , 𝜔 0 , 𝜆 0 and 𝜂 0 ; convergence threshold Δ and maximum number of iterations 𝑀

while iteration 𝑚 &lt; 𝑀 and difference of ELBO &gt; Δ do:

end

𝑚 = 𝑚 +1 ; Σ ( 𝑚 ) ← [ 𝑣 0 I 𝑝 × 𝑝 +2 𝔼 𝑞 𝑏 ( ) ( 1 𝑏 2 ) ∑ 𝐾 𝑖 =1 ∑ 𝑛 𝑖 𝑗 =1 (1 + 𝛿 𝑖𝑗 ) 𝜁 𝑖𝑗 𝐗 𝐗 𝑖𝑗 𝑇 𝑖𝑗 ] -1 ; 𝝁 ( 𝑚 ) ← [ { 𝑣 0 𝝁 𝑇 0 + ∑ 𝐾 𝑖 =1 ∑ 𝑛 𝑖 𝑗 =1 ( ( 𝔼 𝑞 𝑏 ( ) ( 1 𝑏 ) ( -𝛿 𝑖𝑗 +(1 + 𝛿 𝑖𝑗 ) 𝜌 𝑖𝑗 ) +2 𝔼 𝑞 𝑏 ( ) ( 1 𝑏 2 ) (1 + 𝛿 𝑖𝑗 ) 𝜁 𝑖𝑗 ( 𝑦 𝑖𝑗 -𝔼 𝑞 𝛾 ( 𝑖 ) ( 𝛾 𝑖 ) ) ) 𝐗 𝑇 𝑖𝑗 )} Σ ( 𝑚 ) ] 𝑇 ; 𝜎 𝑖 2( 𝑚 ) ← [ 𝔼 𝑞 𝜎 ( 2 ) 𝛾 ( 1 𝜎 2 𝛾 ) +2 𝔼 𝑞 𝑏 ( ) ( 1 𝑏 2 ) ∑ 𝑛 𝑖 𝑗 =1 (1 + 𝛿 𝑖𝑗 ) 𝜁 𝑖𝑗 ] -1 , 𝑖 = 1 , ..., 𝐾 ; 𝜏 ( 𝑚 ) 𝑖 ← 𝜎 𝑖 2( 𝑚 ) ∑ 𝑛 𝑖 𝑗 =1 [ 𝔼 𝑞 𝑏 ( ) ( 1 𝑏 ) ( -𝛿 𝑖𝑗 +(1 + 𝛿 𝑖𝑗 ) 𝜌 𝑖𝑗 ) +2 𝔼 𝑞 𝑏 ( ) ( 1 𝑏 2 )(1 + 𝛿 𝑖𝑗 ) 𝜁 𝑖𝑗 ( 𝑦 𝑖𝑗 -𝐗 𝑇 𝑖𝑗 𝔼 𝑞 ( 𝜷 𝜷 ) ) ] , 𝑖 = 1 , ..., 𝐾 ; 𝜔 𝑚 ( ) ← 𝜔 0 -∑ 𝐾 𝑖 =1 ∑ 𝑛 𝑖 𝑗 =1 ( 𝛿 𝑖𝑗 -(1 + 𝛿 𝑖𝑗 ) 𝜑𝑖𝑗 )( 𝑦 𝑖𝑗 -𝐗 𝑇 𝑖𝑗 𝔼 𝑞 ( 𝜷 ) ( 𝜷 ) -𝔼 𝑞 𝛾 ( 𝑖 ) 𝛾 𝑖 ) ; 𝜂 ( 𝑚 ) ← 𝜂 0 + 1 2 ∑ 𝐾 𝑖 =1 𝔼 𝑞 𝛾 ( 𝑖 ) 𝛾 2 𝑖 ; calculate the current ELBO, ELBO ( 𝑚 ) ; calculate the difference of ELBO = ELBO ( 𝑚 ) -ELBO ( 𝑚 -1) ;

## 4 SIMULATION STUDY

## 4.1 Design of simulation

Weconduct a simulation study to evaluate the performance of our proposed VB algorithms across different scenarios by varying the number of clusters and the number of observations within each cluster.

We generate the log of survival time for the 𝑗 𝑡ℎ subject in the 𝑖 𝑡ℎ cluster, log( 𝑇 𝑖𝑗 ) , 𝑖 = 1 , ..., 𝐾 and 𝑗 = 1 , ..., 𝑛 𝑖 as follows:

$$\log T _ { i j } = 0. 5 + \beta _ { 1 } x _ { i j 1 } + \beta _ { 2 } x _ { i j 2 } + \gamma _ { i } + b \epsilon _ { i j },$$

where 𝑥 𝑖𝑗 1 , 𝑥 𝑖𝑗 2 , and 𝜖 𝑖𝑗 are mutually independently generated with 𝑥 𝑖𝑗 1 ∼ 𝑁 (1 , 0 2 2 . ) , 𝑥 𝑖𝑗 2 ∼ Bernoulli (0 5) . and 𝜖 𝑖𝑗 ∼ logistic (0 , 1) . The values of 𝛽 1 , 𝛽 2 and 𝑏 are chosen as 0.2, 0.8 and 0.8, respectively. The random intercept for the 𝑖 𝑡ℎ cluster, 𝛾 𝑖 , is generated from 𝑁 (0 , 𝜎 2 𝛾 ) with 𝜎 2 𝛾 = 1 . The censoring time for the 𝑗 𝑡ℎ subject in the 𝑖 𝑡ℎ cluster, 𝐶𝑖𝑗 , is generated from a uniform distribution, uniform (0 , 𝑑 ) , where 𝑑 is a positive value controlling the percentage of censoring. Then 𝑡 𝑖𝑗 = min( 𝑇 𝑖𝑗 , 𝐶 𝑖𝑗 ) and 𝛿 𝑖𝑗 = ✶ ( 𝑇 𝑖𝑗 ≤ 𝐶𝑖𝑗 ) . Take 𝑑 = 48 to achieve a 15% censoring rate in our simulations. Investigation of the effect of censoring rates can be found in Xian et al. 31 where they show that as the censoring rate increases, an increase in mean squared error (MSE) of estimating parameters in an AFT model using a VB algorithm was observed.

Weexplore various scenarios by varying the number of clusters, 𝐾 , and the number of observations within each cluster, 𝑛 𝑖 = 𝑛 for all 𝑖 = 1 , ..., 𝐾 . Across the experiments, we consider 𝐾 values from the set {15 30 50 80} , , , and 𝑛 values from {5 15 30 50} , , , ,

resulting in a total of 16 unique scenarios. Our objective is to assess the estimation performance of the VB algorithm concerning the variations in 𝐾 and 𝑛 . We consider a prior setting with 𝝁 0 = (0 0 0) , , 𝑇 , 𝑣 0 = 0 1 . , 𝛼 0 = 𝜆 0 = 3 , and 𝜔 0 = 𝜂 0 = 2 , which indicates no strong prior information on the parameters. The ELBO convergence threshold is set as 0 01 . which is the default recommendation 35 , and the maximum number of iterations is 100 .

We conduct 𝑁 = 500 runs (replicates) in each considered scenario and apply our proposed VB algorithm to derive the approximated posterior distribution of each parameter to each run. The mean of each approximated posterior distribution is used as the parameter estimate for the corresponding parameter. The empirical bias and sample standard deviation (SD) as well as the empirical MSE for each estimate are obtained, where

$$\text{MSE} = \frac { \sum _ { i = 1 } ^ { N } ( \theta _ { 0 } - \hat { \theta } _ { i } ) ^ { 2 } } { N }$$

and ̂ 𝜃 𝑖 is the estimate of parameter 𝜃 in the 𝑖 𝑡ℎ replicate, and 𝜃 0 is the true value. In addition, for each parameter of interest, we report the empirical 95% coverage rate (CR) calculated as

$$C R = \frac { \sum _ { i = 1 } ^ { N } I _ { i } } { N }$$

where 𝐼 is the indicator variable which takes 1 if the true parameter value 𝜃 0 falls into the 95% credible interval.

## 4.2 Simulation results

The empirical bias, sample SD, empirical MSE and empirical CR of estimating 𝛽 1 , 𝛽 2 , 𝑏 , and 𝜎 2 𝛾 from 500 replicates are summarized in Table 1. Considering 𝛽 1 , 𝛽 2 , and 𝑏 , we observe that when the number of clusters 𝐾 is fixed, increasing the cluster size 𝑛 from 5 to 50 does not always significantly affect the empirical bias. However, there is a noticeable decrease in the sample SD, leading to a pronounced reduction in the empirical MSE. A similar trend is observed when the cluster size is fixed while increasing 𝐾 . Figure 1 visually illustrates this asymptotic property through boxplots. Our primary focus lies in understanding the impact of 𝐾 and 𝑛 on estimating the variance of the random intercept, denoted as 𝜎 2 𝛾 . In most of the scenarios, when the number of clusters 𝐾 is fixed, we observe that the bias of estimating 𝜎 2 𝛾 decreases as the cluster size increases, while the sample SD remains relatively stable. However, when 𝐾 is increased, both the empirical bias and the sample SD decrease noticeably, except for the case when 𝑛 = 30 as 𝐾 grows from 15 to 30. Particularly noteworthy is the scenario with 𝐾 = 80 and 𝑛 = 50 , where we achieve a less biased and less variable estimation of the variance of the random intercept. In addition, the empirical CRs corresponding to a 95% credible interval for each parameter across all scenarios are close to the nominal level of 95%, with a mean of 94.1% and a standard deviation of 0.014. Furthermore, based on a nominal level of 95% with 500 replicates, the coverage with two standard deviations is 95% ± 2 SD = [93 1% 97 0%] . , . since SD = √ 0 95 × 0 05∕500 ≈ 0 0097 . . . . We observe that the majority of our empirical CRs are within the range of [93 1% 97 0%] . , . .

The left part of Figure 2 presents the computational run time of 500 replicates in minutes against the number of clusters ( 𝐾 ) with different colors for different cluster sizes 𝑛 . Fixing 𝐾 , as 𝑛 increases, the used time increases. The scatter plot in the right of Figure 2 shows the run time against the sample size ( 𝐾 × 𝑛 ). As anticipated, the computational cost rises with increasing sample size.

b

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00177v1-figure-d54593e0-1748-434c-93cd-5ec209f19354.jpg)

=

b

=

FIGURE 1 Boxplots of parameter estimates using posterior means from 500 replicates under various scenarios with different number of clusters 𝐾 and cluster sizes 𝑛 based on our proposed VB algorithm. The horizontal dashed line on each plot represents the true value of the corresponding parameter used when generating the data.

Furthermore, among those considered 16 scenarios, we select three scenarios ( 𝐾 = 30 with 𝑛 = 5 , 𝐾 = 50 with 𝑛 = 15 , and 𝐾 = 80 with 𝑛 = 30 ) and conduct a comparative analysis of estimation results obtained from the proposed VB algorithm with those from two alternative methods which we introduced in Section 1: the h-likelihood method proposed by Do Ha et al. 14 and the MCMC-based survregbayes developed by Zhou et al. 17 . Summaries from the MCMC-based survregbayes algorithm are derived from one Markov chain, subsampled every 5 iterates to achieve a final chain size of 2,000 after a burn-in period of 5,000 iterates 17 . We further refer to our VB method and the h-likelihood method as survregVBfrailty and survregHL , respectively. Numerical estimation results are presented in Table 2, where we focus on the comparison of estimation for 𝛽 1 , 𝛽 2 and 𝜎 2 𝛾 since the MCMC-based survregbayes does not directly return the estimate for 𝑏 . In the scenario with 𝐾 = 30 and 𝑛 = 5 , the proposed survregVBfrailty exhibits a larger empirical bias compared to both the survregHL and survregbayes methods. However, as 𝐾 and 𝑛 increase, the discrepancy in empirical bias among the three methods becomes negligible. In terms of the sample SD, on average, across the three parameters 𝛽 1 , 𝛽 2 , and 𝜎 2 𝛾 , survregVBfrailty yields a 16.5% smaller sample SD when 𝐾 = 30 and 𝑛 = 5 ,

a 7.6% smaller sample SD when 𝐾 = 50 and 𝑛 = 15 , and a 3.7% smaller sample SD when 𝐾 = 80 and 𝑛 = 30 , compared to survregHL . Compared with the survregBayes algorithm, we observe a decrease in sample SD of 7.8% in the 𝐾 = 30 , 𝑛 = 5 scenario, mainly due to 𝛽 1 , and only 2.0% in the 𝐾 = 50 , 𝑛 = 15 scenario. However, in the scenario with 𝐾 = 80 , 𝑛 = 30 , no average difference in sample SD is observed between VB and MCMC. Therefore, the survregVBfrailty algorithm generally yields a smaller SD, resulting in a lower MSE compared to the other two methods in most cases. On average, across the three parameters 𝛽 1 , 𝛽 2 , and 𝜎 2 𝛾 , survregVBfrailty achieves a 27.2% lower MSE when 𝐾 = 30 , 𝑛 = 5 , a 12.0% lower MSE when 𝐾 = 50 , 𝑛 = 15 , and a 13.4% lower MSE when 𝐾 = 80 , 𝑛 = 30 , compared to survregHL . Compared with survregBayes , we observe a reduction in MSE of 12.4% in the 𝐾 = 30 , 𝑛 = 5 scenario, mainly because of 𝛽 1 , and only 3.0% in the 𝐾 = 50 , 𝑛 = 15 scenario. However, in the scenario with 𝐾 = 80 , 𝑛 = 30 , no difference in MSE is observed between VB and MCMC. In addition, we compare the posterior densities for 𝛽 1 , 𝛽 2 , and 𝜎 2 𝛾 obtained from both survregVBfrailty and the MCMC-based survregBayes using a simulated dataset with 𝐾 = 50 and 𝑛 = 15 . These densities present a strong consistency between the VB and MCMC algorithms as illustrated in Figure 3, which is also reflected in the results presented in Table 2.

7.5

time used in minutes

n

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00177v1-figure-5063e1e3-26b9-418c-b8d5-d7114ce28ae0.jpg)

15

30

50

80

number of clusters K

FIGURE 2 Left : the run time in minutes used for 500 replicates under various scenarios with different number of clusters 𝐾 and cluster sizes 𝑛 based on our proposed VB algorithm. Right : the run time in minutes used for 500 replicates in different sample sizes based on VB.

To assess the computational efficiency of the proposed VB algorithm, we compare the run times in minutes for 500 replicates in these three scenario across survregVBfrailty , survregHL , and survregbayes methods and present the results in Table 3. In Table 3, we observe that in the scenario with a small sample size ( 𝐾 = 30 and 𝑛 = 5 ), there is minimal difference in the run times between the VB and h-likelihood methods. However, as both 𝐾 and 𝑛 increase, the h-likelihood method requires progressively longer computation times to obtain estimates. Particularly in the scenario with 𝐾 = 80 and 𝑛 = 30 , the computational time is over 10 times longer than that required by the VB algorithm. Comparatively, the VB algorithm demonstrates significantly higher computational efficiency when contrasted with the MCMC-based survregbayes algorithm. Notably, in the scenario with

5.0

2.5

5

15

30

50

FIGURE 3 A comparison of posterior densities of 𝛽 1 , 𝛽 2 and 𝜎 2 𝛾 obtained from our VB method survregVBfrailty and MCMCbased survregbayes from a simulated data set in the scenario with 𝐾 = 50 and 𝑛 = 15 .

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00177v1-figure-1c3be911-d56c-4a3d-a814-4ebf8f068ecf.jpg)

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00177v1-figure-22e67f4c-43a7-4071-b582-7d7b746e4f08.jpg)

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00177v1-figure-efcf5851-11f6-4688-a5a8-d341ffde1d38.jpg)

𝐾 = 80 and 𝑛 = 30 , the survregVBfrailty algorithm runs approximately 150 times faster than the survregbayes algorithm. All algorithms were implemented within R version 4.3.2 and simulations were conducted on a computer operating the Mac OS X platform, equipped with a 4.05 GHz CPU and 8 GB RAM.

## TABLE 1 Numerical estimation results (point estimate using the mean of the corresponding posterior distribution) including the empirical Bias, sample SD, empirical MSE and coverage rate (CR) for parameters in each scenario with different number of clusters 𝐾 and cluster sizes 𝑛 from the proposed VB algorithm.

| CR   | 91.4 93.1 95.5             | 94.2              | 92.6 91.1 92.4 95.2                | 93.2                                | 92.4 93.3 93.4 94.4 92.2 92.8 94.4     |
|------|----------------------------|-------------------|------------------------------------|-------------------------------------|----------------------------------------|
| MSE  | 0.107 0.089 0.087          | 0.100 0.084       | 0.062 0.063 0.052 0.059            | 0.042 0.037 0.040 0.041 0.026 0.026 | 0.025                                  |
| SD   | 0.312 0.287 0.290 0.299    | 0.274             | 0.242 0.248 0.228 0.235 0.203      | 0.191 0.201 0.197                   | 0.160 0.161 0.157                      |
| Bias | -0.101 -0.081 0.003 -0.040 | -0.092            | -0.057 -0.036 -0.026 -0.062 -0.033 | -0.027 0.016 -0.048                 | -0.023 0.003 0.014                     |
| CR   | 94.5 94.8 96.0 94.7        | 94.4              | 96.1 95.6 95.8 95.1                | 94.2 95.6 95.7 95.2                 | 95.5 95.6 95.5                         |
| MSE  | 0.009 0.003 0.001 0.001    | 0.005             | 0.001 0.001 < 0 . 001 0.003 0.001  | < 0 . 001 < 0 . 001 0.002 0.001     | < 0 . 001 < 0 . 001                    |
| SD   | 0.095 0.052 0.036 0.029    | 0.071             | 0.038 0.026 0.020 0.051            | 0.029 0.020 0.016 0.041             | 0.023 0.015 0.013                      |
| Bias | 0.005 0.006 -0.006 0.003   | 0.015             | 0.008 0.005 0.003                  | 0.014 0.018 0.003                   | 0.003                                  |
| CR   |                            |                   |                                    | 0.002                               | 91.3 0.014 92.4 0.008 95.1 -0.003 95.5 |
|      | 90.0 94.0 92.4             | 94.8 92.8         | 93.0 95.7 95.0 91.4                | 94.8 94.4 95.6                      |                                        |
| MSE  | 0.124 0.037 0.019 0.011    | 0.066             | 0.020 0.008 0.005 0.041            | 0.011 0.006 0.003                   | 0.026 0.008 0.003 0.002                |
| SD   | 0.352 0.192 0.139          | 0.104 0.257       | 0.140 0.088 0.074                  | 0.199 0.104 0.074 0.057             | 0.161 0.090 0.058 0.043                |
| Bias | -0.026 0.007 -0.017 -0.001 | -0.015            | -0.018 -0.004 -0.005               | -0.033 -0.001 -0.006 -0.004         | -0.005 -0.005 -0.005 -0.002            |
| CR   | 94.2 93.2 93.2             | 95.6 92.8         | 93.0 94.6 93.0 94.3                | 95.2 92.3 94.4 93.2                 | 93.4 95.8 95.4                         |
| MSE  | 0.587 0.239 0.125          |                   | 0.129 0.057 0.037                  | 0.147                               |                                        |
|      | 0.764 0.489 0.354 0.254    | 0.065 0.589 0.360 | 0.346 0.239 0.192                  | 0.220 0.064 0.037 0.020             | 0.223 0.050 0.141 0.020 0.108 0.012    |
| SD   |                            |                   | -0.007                             | 0.470 0.253 0.192 0.143 0.384       | -0.008 0.008                           |
| Bias | 0.066 0.018 0.006 0.009    | 0.004             | 0.002 0.012                        | 0.003 -0.006 -0.002                 | -0.005                                 |
| 𝑛    | 5 15 30 50                 | 5 15              | 30 50 5                            | -0.011 15 30 50                     | 5 0.013 15 30 50                       |
| 𝐾    | 15                         |                   | 30                                 | 50                                  | 80                                     |

TABLE 2 A comparison of numerical estimation results including the empirical Bias, sample SD and MSE, from our VB method survregVBfrailty , the h-likelihood method survregHL and MCMC-based survregbayes method in each scenario ( 𝐾 : number of clusters, 𝑛 : number of observations in each cluster).

|               |           | survregVBfrailty   | survregVBfrailty   | survregVBfrailty   | survregHL   | survregHL   | survregHL   | survregbayes   | survregbayes   | survregbayes   |
|---------------|-----------|--------------------|--------------------|--------------------|-------------|-------------|-------------|----------------|----------------|----------------|
| Scenario      | Parameter | Bias               | SD                 | MSE                | Bias        | SD          | MSE         | Bias           | SD             | MSE            |
| = 30 , 𝑛 = 5  | 𝛽 1       | 0.004              | 0.589              | 0.346              | 0.001       | 0.664       | 0.440       | 0.001          | 0.664          | 0.440          |
| = 30 , 𝑛 = 5  | 𝛽 2       | -0.015             | 0.257              | 0.066              | -0.002      | 0.270       | 0.073       | -0.002         | 0.266          | 0.071          |
| = 30 , 𝑛 = 5  | 𝜎 2 𝛾     | -0.092             | 0.274              | 0.084              | 0.001       | 0.412       | 0.170       | -0.042         | 0.300          | 0.092          |
| = 50 , 𝑛 = 15 | 𝛽 1       | 0.003              | 0.253              | 0.064              | 0.001       | 0.268       | 0.072       | -0.005         | 0.259          | 0.067          |
| = 50 , 𝑛 = 15 | 𝛽 2       | -0.001             | 0.104              | 0.011              | 0.007       | 0.107       | 0.011       | -0.004         | 0.104          | 0.011          |
| = 50 , 𝑛 = 15 | 𝜎 2 𝛾     | -0.033             | 0.203              | 0.042              | 0.011       | 0.237       | 0.056       | -0.001         | 0.211          | 0.044          |
| = 80 , 𝑛 = 30 | 𝛽 1       | 0.008              | 0.141              | 0.020              | 0.009       | 0.143       | 0.021       | -0.009         | 0.142          | 0.020          |
| = 80 , 𝑛 = 30 | 𝛽 2       | -0.005             | 0.058              | 0.003              | 0.001       | 0.060       | 0.004       | 0.001          | 0.058          | 0.003          |
| = 80 , 𝑛 = 30 | 𝜎 2 𝛾     | 0.003              | 0.161              | 0.026              | 0.007       | 0.172       | 0.029       | 0.001          | 0.160          | 0.026          |

TABLE 3 Times in minutes for 500 replicates from our VB algorithm survregVBfrailty , the h-likelihood method survregHL and the MCMC-based survregbayes algorithm, respectively, under each scenario.

| Scenario        |   survregVBfrailty |   survregHL |   survregbayes |
|-----------------|--------------------|-------------|----------------|
| 𝐾 = 30 , 𝑛 = 5  |               0.61 |        0.67 |          47.65 |
| 𝐾 = 50 , 𝑛 = 15 |               1.7  |        9.59 |         232.31 |
| 𝐾 = 80 , 𝑛 = 30 |               8.47 |      124.53 |        1302.8  |

## 5 APPLICATION TO VENTILATION DURATION ANALYSIS

In this section, we apply the shared frailty log-logistic AFT model with the proposed VB algorithm to the ICU data as we descried in the Introduction section and conduct a retrospective study on the ventilation duration time. We aim to investigate the ICU site-specific random effect on patient's ventilation duration. We extend the work by Kobara et al. 9 by incorporating the uncertainty from the ICU sites via a group-specific random intercept under a Bayesian analysis framework.

The CCIS ICU data were collected between July 2015 and December 2016 and contained 49,467 patients receiving invasive mechanical ventilation upon arrival to ICU. About 3% of these patients were discharged/transferred to a Complex Continuing Care Facility, other hospitals, the Level 3 Unit, and Outside the ICU while still on a ventilator, and therefore, their ventilation time were considered as right-censored data 9 . The data were from 66 ICU sites (centers) and each site has a unique site code, for example, 3970. In the CCIS dataset, the average number of ventilated patients per ICU site over the study period of about 1.5 years is 749.5. We consider the significant covariates investigated by Kobara et al. 9 as fixed effects which are admission

source (e.g., from operation rooms), admission diagnosis, patient type (medical or surgical), scheduled admission (yes or no), scheduled surgery (yes or no), referring physician specialty, other interventions (yes or no), central venous line (CVL, yes or no), arterial line (AL, yes or no), intra-cranial pressure monitor (IPM, yes or no), extracorporeal membrane oxygen (EMO, yes or no), intra-aortic balloon pump (IABP, yes or no), age group (18-39, 40-80 or above 80 years of age), pre-LOS (no more than 1 day, between 2 and 7 days, or no less than 7 days), and the MODS score (none with MODS ≤ 1 , minimal with 1 - 4 scores, mild with 4 - 8 scores, moderate with 8 - 12 scores or, severe with scores &gt; 12 ).

To apply our proposed VB algorithm in real data analysis, we consider the same weak prior setting in our simulation study described in Section 4: 𝝁 0 = 0 , 𝑣 0 = 0 1 . , 𝛼 0 = 𝜆 0 = 3 , and 𝜔 0 = 𝜂 0 = 2 . Table 4 displays the estimated regression coefficients from the fitted shared frailty log-logistic AFT model. For comparative analysis, we also employed the h-likelihood method survregHL and the MCMC-based survregbayes to model the data incorporating shared frailty. In the absence of shared frailty, we used the likelihood-based survreg from the survival package in R, and the VB method survregVB proposed by Xian et al. 31 . Furthermore, for the Bayesian methods ( survregVBfrailty , survregVB , and survregbayes ), we report the 95% credible intervals, whereas for the likelihood-based methods, we provide the 95% confidence intervals. In Table 4, we first observe that there is no significant difference in the estimated regression coefficients between survregVBfrailty and survregbayes . However, some of the estimated coefficients from survregHL , such as scheduled surgery, CVL, AL, IPM and EMO, are different from those based on survregVBfrailty and survregbayes , which is further discussed in Section 6. In comparison to the models with frailty, the survregVB and survreg methods, which do not account for frailty, exhibit some differences in estimating certain coefficients. However, they maintain a strong overall consistency in estimating the regression coefficients with the results from the frailty models. This consistency is expected, as the inclusion of frailty does not affect the fixed effects, implying that the regression lines from models with or without frailty should be parallel. To illustrate the difference in interval estimation between models with and without frailty, we calculated the mean overlap percentage of the two 95% credible intervals from VB with or without frailty. Using the 95% credible intervals from VB without frailty ( survregVB ) as references, we found a 77.3% overlap on average. This means that if we fit a model using VB without considering the ICU site as a shared frailty, about 77% of the credible interval will fall within the corresponding credible interval obtained from the model that includes the shared frailty.

In what follows, we summarize the fixed effects on ventilation duration based on estimates from our proposed VB algorithm, survregVBfrailty . Regarding the admission source, patients arriving from the emergency department (ED), operating room (OR), or ward have ventilation times that are (1 - exp(-0 134))100% = 12 54% . . shorter (credible interval (CI): [9.24%, 12.72%]), 27.89% shorter (CI: [25.62%, 30.02%]), and 13.06% shorter (CI: [9.70%, 16.31%]), respectively, compared to patients admitted from a downstream unit (baseline). Conversely, patients admitted from home or another hospital have ventilation durations that are 10.63% (CI: [-10.40%, 23.49%]) and 10.30% (CI: [7.04%, 13.54%]) longer, respectively. Since the CI for patients admitted

from home includes zero, we conclude that admission from home is not a statistically significant factor. Additionally, patients from other sources, such as from outside the province, have ventilation durations that are 12.98% longer (CI: [5.13%, 21.29%]).

Compared to cardiovascular patients, those with gastrointestinal, neurological, and trauma diagnoses experience significantly longer ventilation times. Specifically, ventilation duration increases by 29.69% (CI: [25.36%, 34.18%]) for gastrointestinal patients, 33.11% (CI: [29.05%, 37.44%]) for neurological patients, and 77.71% (CI: [69.72%, 86.08%]) for trauma patients. Patient categories (surgical or medical) do not show significant differences in ventilation time. However, scheduled ICU admissions or surgeries are important factors. Patients with a scheduled ICU admission have a 24.95% longer ventilation time (CI: [21.89%, 27.82%]), while those with a scheduled surgery have a 13.32% shorter ventilation time (CI: [9.70%, 16.81%]) compared to patients without scheduled admissions or surgeries. The referral physician service is also a significant risk factor. Compared with medical referrals, surgical referrals result in a 5.82% shorter ventilation time (CI: [8.97%, 2.57%]), while respirology referrals result in a 16.77% longer ventilation time (CI: [11.29%, 22.51%]).

Specific treatment interventions upon ICU arrival significantly impact ventilation duration. Patients receiving a CVL, AL, IPM, EMO, or IABP have increased ventilation times by 20.20% (CI: [17.82%, 22.51%]), 21.90% (CI: [19.36%, 24.48%]), 70.40% (CI: [60.16%, 81.48%]), 148.93% (CI: [127.50%, 172.10%]), and 40.92% (CI: [32.45%, 49.78%]), respectively. Age is also a significant factor. Compared to patients aged 18 to 39, those aged 40 to 80 and those over 80 have longer ventilation periods, increasing by 14.22% (CI: [11.29%, 17.23%]) and 7.79% (CI: [4.50%, 11.18%]), respectively.

Both pre-ICU LOS and patients' severity scores, as measured by MODS, are significant risk factors for longer ventilation duration. Longer pre-ICU LOS correlates with longer ventilation times. Patients with a pre-ICU LOS of 2-7 days or ≥ 7 days experience increases in ventilation duration by 3.87% (CI: [1.71%, 6.08%]) and 12.98% (CI: [10.52%, 15.60%]), respectively, compared to those with ≤ 1 day pre-ICU LOS. Compared to patients with a MODS score ≤ 1 , those with minimal (1-4), mild (4-8), moderate (8-12), or severe ( &gt; 12 ) MODS scores experience increases in ventilation duration by 10.63% (CI: [7.90%, 13.54%]), 23.49% (CI: [20.56%, 26.49%]), 30.34% (CI: [26.49%, 34.31%]), and 14.22% (CI: [7.47%, 21.53%]), respectively.

Wenowturn to the analysis of variance in terms of the variations from individual patients and from ICU sites. Using the mean of the posterior distribution for 𝜎 2 𝛾 , the estimated variance ( 𝜎 2 𝛾 ) for the shared ICU site effect from survregVBfrailty is 0.1, with a square root of 0.316, which indicates that the average spread of the ventilation time among ICU sites is exp(0 316) . , or 1.372. The estimated scale parameter 𝑏 is 0.444, resulting in an estimated variance of the logarithms of individual patient ventilation times of 0 444 2 . × 𝜋 2 ∕3 , or 0.649, where 𝜋 2 ∕3 is the variance of a standard log-logistic distribution. To measure the strength of the correlation between patients within the same ICU site, the intra-class correlation coefficient between the logarithms of ventilation time can be estimated by 0 1∕(0 1 + 0 649) ≈ 0 134 . . . . .

In the left of Figure 4, we observe a strong consistency in estimating the variance of the ICU-site specific random effect between survregVBfrailty and the MCMC-based survregbayes . We further visualize the estimated ICU-site random effects with

the 95% credible intervals based on their posterior distributions obtained from survregVBfrailty , as shown in the right of Figure 4, where the random effects are ranked from smallest to largest. Wider credible intervals correspond to ICUs with larger numbers of patients. As discussed in Lambert et al. 12 , if the intervals overlap, there are no significant center random effects. We observe that some of the intervals do not overlap, indicating that these ICU sites perform differently in terms of patient ventilation duration.

For reference on computational efficiency, the run times in minutes for each method are as follows: survregVBfrailty took 1.45 minutes, h-likelihood took 106.18 minutes, MCMC-based survregbayes took 267.04 minutes, survreg took 0.02 minutes, and survregVB took 0.13 minutes.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00177v1-figure-ae761fae-7347-4f46-b6ab-8ba12f9da7c0.jpg)

Estimate

FIGURE 4 Left : Posterior distribution of the variance of the random intercept from VB and MCMC. Right : Estimated ICU site specific random effects with their 95% credible interval from the proposed VB algorithm. The random effects have been ranked in an increasing order.

![Image](https://ac.marcodsn.me/arxiv-markdown-images/2408.00177v1-figure-70f270e8-df10-43d6-bb70-3df8aed3afd5.jpg)

## 6 DISCUSSION

In this paper, we have proposed a fast variational Bayesian algorithm, called survregVBfrailty , for statistical inference using a shared frailty log-logistic AFT model, which can be applied to analyze clustered survival data. We demonstrated that our proposed survregVBfrailty algorithm achieves satisfactory estimation performance through simulation studies under various scenarios with different numbers of clusters and cluster sizes. A strong consistency in the estimated posterior distributions was observed between VB and MCMC methods in both the simulation study and the application to ICU ventilation data. Moreover, the survregVBfrailty algorithm is significantly more computationally efficient, running over 150 times faster than the MCMCbased survregbayes algorithm.

The h-likelihood method proposed by Do Ha et al. 13 for analyzing clustered survival data in a log-normal AFT model has been shown to be robust against violations of the normality assumption (e.g., extreme value distribution) for the logarithm of

survival time. In our simulation study, we investigated the performance of the h-likelihood method when the survival data, given the covariates in a specific cluster, follow a log-logistic distribution. We compared its estimation results to those obtained using VB and MCMC algorithms. We found that the h-likelihood method, survregHL , results in higher mean squared errors (MSEs) compared to VB and MCMC algorithms. Specifically, survregHL produced a 49.1% higher MSE than VB and a 41.2% higher MSEthan MCMC for estimating the variance of the random intercept in our simulation study. Additionally, in our application to ICU ventilation duration analysis, we observed significant differences in some regression coefficients between h-likelihood and VB or MCMC methods. Therefore, our proposed shared frailty log-logistic AFT model using variational Bayes can be viewed as a better approach to the shared frailty log-normal AFT model using h-likelihood.

Our application of the proposed method to the CCIS ICU data for ventilation duration analysis reveals a moderate correlation of 0.134 among patients within the same ICU site. By incorporating the ICU site-specific random effect as an unknown shared frailty, we further validate the significant risk factors including the patient severity score MODS identified in the study by Kobara et al. 9 based on the 95% credible intervals obtained from the posterior distributions of each regression coefficient. We observe a similar trend regarding the effect of MODS on ventilation duration as reported by Kobara et al. 9 . Specifically, when MODS reaches a severe level, the increase in ventilation duration becomes less pronounced. This may be due to patient mortality, which warrants further investigation in future studies. We demonstrate that different ICU sites have varying effects on patient ventilation duration, as observed through differences in the estimated ICU site-specific random effects. Our research provides valuable insights and practical implications for clinical practice and resource management. Specifically, ICUs with smaller estimated random intercepts tend to have shorter overall ventilation durations compared to other ICUs. Understanding the patient characteristics and ventilation practices within these ICUs may help improve clinical performance. For ICUs with larger estimated random intercepts, equipping more ventilators could enhance the overall efficiency of the ICU ventilation procedures.

We present several open problems and directions for future work related to our proposed methodology. As discussed, we assume that the survival time from a specific cluster follows a log-logistic distribution. The robustness of the proposed shared frailty log-logistic AFT model to other distributions remains unknown and warrants further study and comparison with the shared frailty log-normal AFT model. Additionally, extending the current VB algorithm to accommodate other survival distributions, such as the Weibull distribution, could be an interesting avenue of research. In our current framework, we account for clusterlevel uncertainty using a random intercept. This approach could be extended to a more general model that includes cluster-level covariates (e.g., the ICU type, general or specialized), resulting in a mixed-effects model. Such a general model can be used to better assess differences in performance across ICU sites regarding the duration of invasive mechanical ventilation. Another potential area for development is the integration of variable selection techniques within the VB algorithm, see Park and Do Ha 16 as a reference.

## ACKNOWLEDGEMENTS

This research is supported by the Natural Sciences and Engineering Research Council of Canada (NSERC). We extend our gratitude to Dr. Yawo M. Kobara for facilitating access to the CCIS data set used in our real data analysis.

## DATA AVAILABILITY STATEMENT

The ICU data analyzed in Section 5 were from Critical Care Information System (CCIS) Ontario database. Under the authors' data use agreement, they are unable to share the files directly.

## ORCID

Chengqian Xian: https://orcid.org/0000-0002-6855-3221

Camila P.E. de Souza: https://orcid.org/0000-0002-4465-5792

Wenqing He: https://orcid.org/0000-0002-8913-9273

Felipe F. Rodrigues: https://orcid.org/0000-0001-9370-5998

Renfang Tian: https://orcid.org/0000-0001-8176-2469

## References

- 1. Luo S, Yi M, Huang X, Hunt KK. A Bayesian model for misclassified binary outcomes and correlated survival data with applications to breast cancer. Statistics in medicine 2013; 32(13): 2320-2334.
- 2. Honerkamp-Smith G, Xu R. Three measures of explained variation for correlated survival data under the proportional hazards mixed-effects model. Statistics in medicine 2016; 35(23): 4153-4165.
- 3. Liu XR, Pawitan Y, Clements MS. Generalized survival models for correlated time-to-event data. Statistics in medicine 2017; 36(29): 4743-4762.
- 4. Hougaard P. Frailty models for survival data. Lifetime data analysis 1995; 1: 255-273.
- 5. Hanagal DD. Modeling survival data using frailty models . Springer . 2011.
- 6. Gorfine M, Zucker DM. Shared frailty methods for complex survival data: a review of recent advances. Annual Review of Statistics and Its Application 2023; 10: 51-73.

- 7. Canadian Institute for Health Information . Care in Canadian ICUs. tech. rep., Canadian Institute for Health Information; Canada: 2020.
- 8. Marshall JC, Cook DJ, Christou NV, Bernard GR, Sprung CL, Sibbald WJ. Multiple organ dysfunction score: a reliable descriptor of a complex clinical outcome. Critical care medicine 1995; 23(10): 1638-1652.
- 9. Kobara YM, Wismer M, Rodrigues FF, de Souza CPE. Invasive mechanical ventilation duration prediction using survival analysis. International Journal of Healthcare Management 2023: 1-11.
- 10. Burgess Jr JF, Christiansen CL, Michalak SE, Morris CN. Medical profiling: improving standards and risk adjustments using hierarchical models. Journal of health economics 2000; 19(3): 291-309.
- 11. Glance LG, Dick AW, Osler TM, Mukamel D. Using hierarchical modeling to measure ICU quality. Intensive Care Medicine 2003; 29(12): 2223-2229. doi: 10.1007/s00134-003-1959-9
- 12. Lambert P, Collett D, Kimber A, Johnson R. Parametric accelerated failure time models with random effects and an application to kidney transplant survival. Statistics in medicine 2004; 23(20): 3177-3192.
- 13. Do Ha I, Lee Y, Song JK. Hierarchical-likelihood approach for mixed linear models with censored data. Lifetime data analysis 2002; 8: 163-176.
- 14. Do Ha I, Jeong JH, Lee Y. Statistical modelling of survival data with random effects. Statistics for Biology and Health 2017.
- 15. Do Ha I, Noh M, Lee Y. frailtyHL: A Package for Fitting Frailty Models with H-likelihood.. R J. 2012; 4(2): 28.
- 16. Park E, Do Ha I. Penalized variable selection for accelerated failure time models with random effects. Statistics in medicine 2019; 38(5): 878-892.
- 17. Zhou H, Hanson T, Zhang J. Generalized accelerated failure time spatial frailty model for arbitrarily censored data. Lifetime data analysis 2017; 23: 495-515.
- 18. Zhou H, Hanson T, Zhang J. spBayesSurv: Fitting Bayesian Spatial Survival Models Using R. Journal of Statistical Software 2020; 92(9): 1-33. doi: 10.18637/jss.v092.i09
- 19. Bishop C. Pattern Recognition and Machine Learning . Springer . 2006.
- 20. Hoffman MD, Blei DM, Wang C, Paisley J. Stochastic variational inference. Journal of Machine Learning Research 2013; 14: 1303-1347.
- 21. Ranganath R, Gerrish S, Blei D. Black box variational inference. In: PMLR. ; 2014: 814-822.

- 22. Jordan MI, Ghahramani Z, Jaakkola T, Saul L. Introduction to variational methods for graphical models. Machine Learning 1999; 37: 183-233.
- 23. Neville SE, Palmer M, Wand M. Generalized extreme value additive model analysis via mean field variational bayes. Australian &amp; New Zealand Journal of Statistics 2011; 53(3): 305-330.
- 24. Pham TH, Ormerod JT, Wand M. Mean field variational Bayesian inference for nonparametric regression with measurement error. Computational Statistics &amp; Data Analysis 2013; 68: 375-387. doi: https://doi.org/10.1016/j.csda.2013.07.014
- 25. Luts J, Wand MP. Variational inference for count response semiparametric regression. Bayesian Analysis 2015; 10(4): 991 - 1023. doi: 10.1214/14-BA932
- 26. Ray K, Szabó B. Variational Bayes for high-dimensional linear regression with sparse priors. Journal of the American Statistical Association 2022; 117(539): 1270-1281. doi: 10.1080/01621459.2020.1847121
- 27. Lee CYY, Wand MP. Streamlined mean field variational Bayes for longitudinal and multilevel data analysis. Biometrical Journal 2016; 58(4): 868-895.
- 28. Xian C, de Souza CPE, Jewell J, Dias R. Clustering functional data via variational inference. Advances in Data Analysis and Classification 2024: 1-50.
- 29. da Cruz AC, de Souza CPE, Sousa PHTO. Fast Bayesian Basis Selection for Functional Data Representation with Correlated Errors. arXiv 2024.
- 30. Blei DM, Kucukelbir A, McAuliffe JD. Variational inference: A review for statisticians. Journal of the American Statistical Association 2017; 112(518): 859-877.
- 31. Xian C, de Souza CPE, He W, Rodrigues FF, Tian R. Variational Bayesian analysis of survival data using a log-logistic accelerated failure time model. Statistics and Computing 2024; 34(2): 67.
- 32. Robinson GK. That BLUP is a good thing: the estimation of random effects. Statistical science 1991: 15-32.
- 33. Nolan TH, Menictas M, Wand MP. Streamlined variational inference with higher level random effects. Journal of Machine Learning Research 2020; 21(157): 1-62.
- 34. Lee CYY, Wand MP. Variational methods for fitting complex Bayesian mixed effects models to health data. Statistics in medicine 2016; 35(2): 165-188.
- 35. Yes, but did it Work?: Evaluating variational inference . 80 of Proceedings of Machine Learning Research ; PMLR: 2018.

## APPENDIX A. UPDATE EQUATIONS

## (1) Update equation for 𝑞 ∗ ( 𝜷 )

$$\log q ^ { * } ( \beta ) \stackrel { + } { \approx } \mathbb { E } _ { - \beta } [ \log p ( \mathbf D \, | \, \beta, \gamma, b ) ] + \mathbb { E } _ { - \beta } [ \log p ( \beta ) ]$$

where

$$& \mathbb { E } _ { - \beta } [ \log p ( \mathbf D | \, \beta, \gamma, \, b ) ] \\ = & \mathbb { E } _ { - \beta } \left [ - \delta \log b + \sum _ { i = 1 } ^ { K } \sum _ { j = 1 } ^ { n _ { i } } \left \{ \delta _ { i j } \frac { y _ { i j } - \mathbf X _ { i j } ^ { T } \beta - \gamma _ { i } } { b } - ( 1 + \delta _ { i j } ) \log \left \{ 1 + \exp \left ( \frac { y _ { i j } - \mathbf X _ { i j } ^ { T } \beta - \gamma _ { i } } { b } \right ) \right \} \right \} \right ]$$

with 𝛿 = ∑ 𝐾 𝑖 =1 ∑ 𝑛 𝑖 𝑗 =1 𝛿 𝑖𝑗 being the number of observed survival times, and

$$\mathbb { E } _ { - \beta } [ \log p ( \beta ) ] & \stackrel { + } { \approx } \frac { p } { 2 } \log v _ { 0 } - \frac { 1 } { 2 } v _ { 0 } ( \beta - \mu _ { 0 } ) ^ { T } ( \beta - \mu _ { 0 } ) \\ & \stackrel { + } { \approx } - \frac { 1 } { 2 } v _ { 0 } \left [ \beta ^ { T } \beta - 2 \mu _ { 0 } ^ { T } \beta \right ] = v _ { 0 } \, \mu _ { 0 } ^ { T } \beta - \frac { 1 } { 2 } v _ { 0 } \beta ^ { T } \beta$$

In Equation (6), we use the quadratic piece-wise approximation proposed by Xian et al to approximate 1 + exp ( 𝑦 𝑖𝑗 -𝐗 𝑇 𝑖𝑗 𝜷 -𝛾 𝑖 𝑏 ) :

$$1 + \exp { ( \frac { y _ { i j } - \mathbf X _ { i j } ^ { T } \beta - \gamma _ { i } } { b } ) } \approx \rho _ { i j } \frac { y _ { i j } - \mathbf X _ { i j } ^ { T } \beta - \gamma _ { i } } { b } + \zeta _ { i j } ( \frac { y _ { i j } - \mathbf X _ { i j } ^ { T } \beta - \gamma _ { i } } { b } ) ^ { 2 }$$

where 𝜌 𝑖𝑗 ∶= 0 𝜈 𝑖𝑗 1 ×0 1696 . 𝜈 𝑖𝑗 2 ×0 5 . 𝜈 𝑖𝑗 3 ×0 8303 . 𝜈 𝑖𝑗 4 ×1 1-∑ 4 𝑘 =1 𝜈 𝑖𝑗𝑘 and 𝜁 𝑖𝑗 ∶= 0 𝜈 𝑖𝑗 1 ×0 0189 . 𝜈 𝑖𝑗 2 ×0 1138 . 𝜈 𝑖𝑗 3 ×0 0190 . 𝜈 𝑖𝑗 4 ×0 1-∑ 4 𝑘 =1 𝜈 𝑖𝑗𝑘 with 𝜈 𝑖𝑗 1 = ✶ ( 𝑦 𝑖𝑗 -𝐗 𝑇 𝑖𝑗 𝜷 -𝛾 𝑖 𝑏 ≤ -5 ) , 𝜈 𝑖𝑗 2 = ✶ ( -5 &lt; 𝑦 𝑖𝑗 -𝐗 𝑇 𝑖𝑗 𝜷 -𝛾 𝑖 𝑏 ≤ -1 7 . ) , 𝜈 𝑖𝑗 3 = ✶ ( -1 7 . &lt; 𝑦 𝑖𝑗 -𝐗 𝑇 𝑖𝑗 𝜷 -𝛾 𝑖 𝑏 ≤ 1 7 . ) , and 𝜈 𝑖𝑗 4 = ✶ ( 1 7 . &lt; 𝑦 𝑖𝑗 -𝐗 𝑇 𝑖𝑗 𝜷 -𝛾 𝑖 𝑏 ≤ 5 ) . Then, we obtain,

$$& \nu ^ { + } & \quad \log q ^ { * } ( \beta ) & \stackrel { + } { \approx } & \left [ v _ { 0 } \mu _ { 0 } ^ { T } + \left ( \sum _ { i = 1 } ^ { K } \sum _ { j = 1 } ^ { n _ { i } } \left ( \mathbb { E } _ { q ( b ) } ( \frac { 1 } { b } ) \left ( - \delta _ { i j } + ( 1 + \delta _ { i j } ) \rho _ { i j } \right ) + 2 \mathbb { E } _ { q ( b ) } ( \frac { 1 } { b ^ { 2 } } ) ( 1 + \delta _ { i j } ) \zeta _ { i j } ( y _ { i j } - \mathbb { E } _ { q ( i j } ) ( y _ { i j } ) \right ) \mathbf X _ { i j } ^ { T } \right ) \right ] \beta \\ & \quad - \frac { 1 } { 2 } \beta ^ { T } \left ( v _ { 0 } \mathbf 1 + 2 \mathbb { E } _ { q ( b ) } ( \frac { 1 } { b ^ { 2 } } ) \sum _ { i = 1 } ^ { K } \sum _ { j = 1 } ^ { n _ { i } } ( 1 + \delta _ { i j } ) \zeta _ { i j } X _ { i j } X _ { i j } ^ { T } \right ) \beta \\ & \quad -$$

Let and

$$\mu ^ { * } \, \colon = \left [ \left \{ v _ { 0 } \, \mu _ { 0 } ^ { T } + \sum _ { i = 1 } ^ { K } \sum _ { j = 1 } ^ { n _ { i } } \left ( \mathbb { E } _ { q ( b ) } ( \frac { 1 } { b } ) \left ( - \delta _ { i j } + ( 1 + \delta _ { i j } ) \rho _ { i j } \right ) + 2 \mathbb { E } _ { q ( b ) } ( \frac { 1 } { b ^ { 2 } } ) ( 1 + \delta _ { i j } ) \zeta _ { i j } ( y _ { i j } - \mathbb { E } _ { q ( y _ { i j } } ( \gamma _ { i j } ) ) \right ) \mathbb { X } _ { i j } ^ { T } \right \} \Sigma ^ { T } \right ] ^ { T }$$

Then, 𝑞 ∗ ( 𝜷 ) is 𝑁𝑝 ( 𝝁 ∗ , Σ ∗ ) .

$$\Sigma ^ { * } \coloneqq \left [ v _ { 0 } \mathbf I + 2 \mathbb { E } _ { q ( b ) } ( \frac { 1 } { b ^ { 2 } } ) \sum _ { i = 1 } ^ { K } \sum _ { j = 1 } ^ { n _ { i } } ( 1 + \delta _ { i j } ) \zeta _ { i j } \mathbf X _ { i j } \mathbf X _ { i j } ^ { T } \right ] ^ { - 1 }$$

## (2) Update equation for 𝑞 ∗ ( 𝛾 𝑖 )

$$\log q ^ { * } ( \gamma _ { i } ) \stackrel { + } { \approx } \mathbb { E } _ { - \gamma _ { i } } \left [ \log p ( \mathbf D \, | \, \beta, \gamma, b ) \right ] + \mathbb { E } _ { - \gamma _ { i } } \left [ \sum _ { i = 1 } ^ { K } \log p ( \gamma _ { i } \, | \sigma _ { \gamma } ^ { 2 } ) \right ]$$

where

$$\mathbb { E } _ { - \gamma _ { i } } \left [ \log p ( \mathbf D \, | \, \beta, \gamma, b ) \right ] \\ \overset { \ast } { \sim } \mathbb { E } _ { - \gamma _ { i } } \left [ \sum _ { j = 1 } ^ { n _ { i } } \left \{ \delta _ { i j } \frac { y _ { i j } - \mathbf X _ { i j } ^ { T } \beta - \gamma _ { i } } { b } - ( 1 + \delta _ { i j } ) \log \left ( 1 + \exp \left ( \frac { y _ { i j } - \mathbf X _ { i j } ^ { T } \beta - \gamma _ { i } } { b } \right ) \right ) \right \} \right ]$$

and

$$\mathbb { E } _ { - \gamma _ { i } } \left [ \sum$$

We again apply the quadratic approximation in (7) to (8), and obtain

$$\begin{array} {$$

Let and

$$\tau _ { i } ^ { * } = \sigma _ { i } ^ { 2$$

Then, 𝑞 ∗ ( 𝛾 𝑖 ) is 𝑁𝑙 ( 𝜏 ∗ 𝑖 , 𝜎 2∗ 𝑖 ) .

(3) Update equation for 𝑞 ∗ ( 𝑏 )

$$\log q ^ { * } ( b ) \stackrel { + } { \approx } \$$

where

$$& \mathbb { E } _ { - b } [ \log p ( \mathbf D |$$

𝔼 -𝑏 [log 𝑝 𝑏 ( )] ≈ -( + 𝛼 0 +1)log 𝑏 -𝜔 0 𝑏 since we assume the prior 𝑏 ∼ Inverse-Gamma ( 𝛼 , 𝜔 0 0 )

and

$$\sigma _ { i } ^ { 2 * } = \left [ \mathbb { E } _ { q ( \sigma _ { \gamma } ^ { 2 } ) } ( \frac { 1 } { \sigma _ { \gamma } ^ { 2 } } ) + 2 \mathbb { E } _ { q ( b ) } ( \frac { 1 } { b ^ { 2 } } ) \sum _ { j = 1 } ^ { n _ { i } } ( 1 + \delta _ { i j } ) \zeta _ { i j } \right ] ^ { - 1 }$$

In Equation (9), we apply the linear piece-wise approximation proposed by Xian et al. to approximate 1 + exp ( 𝑦 𝑖𝑗 -𝐗 𝑇 𝑖𝑗 𝜷 -𝛾 𝑖 𝑏 ) to achieve conjugacy:

$$1 + \exp { ( \frac { y _ { i j } - { \mathbf X } _ { i j } ^ { T } \beta - \gamma _ { i } } { b } ) } \approx \varphi _ { i j } \frac { y _ { i j } - { \mathbf X } _ { i j } ^ { T } \beta - \gamma _ { i } } { b }.$$

where 𝜑𝑖𝑗 = 0 𝜂 𝑖𝑗 1 ×0 0426 . 𝜂 𝑖𝑗 2 ×0 3052 . 𝜂 𝑖𝑗 3 ×0 6950 . 𝜂 𝑖𝑗 4 × 0 9574 . 𝜂 𝑖𝑗 5 ×1 1-∑ 5 𝑘 =1 𝜂 𝑖𝑗𝑘 with 𝜂 𝑖𝑗 1 = ✶ ( 𝑦 𝑖𝑗 -𝐗 𝑇 𝑖𝑗 𝜷 -𝛾 𝑖 𝑏 ≤ -5 ) , 𝜂 𝑖𝑗 2 = ✶ ( -5 &lt; 𝑦 𝑖𝑗 -𝐗 𝑇 𝑖𝑗 𝜷 -𝛾 𝑖 𝑏 ≤ -1 701 . ) , 𝜂 𝑖𝑗 3 = ✶ ( - 1 701 . &lt; 𝑦 𝑖𝑗 -𝐗 𝑇 𝑖𝑗 𝜷 -𝛾 𝑖 𝑏 ≤ 0 ) , 𝜂 𝑖𝑗 4 = ✶ ( 0 &lt; 𝑦 𝑖𝑗 -𝐗 𝑇 𝑖𝑗 𝜷 -𝛾 𝑖 𝑏 ≤ 1 702 . ) , and 𝜂 𝑖𝑗 5 = ✶ ( 1 702 . &lt; 𝑦 𝑖𝑗 -𝐗 𝑇 𝑖𝑗 𝜷 -𝛾 𝑖 𝑏 ≤ 5 ) .

Then, we obtain

$$\log q ^ { * } ( b ) \stackrel { * } { \approx } - ( \alpha _ { 0 } + \delta + 1 ) \log b - \frac { 1 } { b } \left [ \omega _ { 0 } - \sum _ { i = 1 } ^ { K } \sum _ { j = 1 } ^ { n _ { i } } \left ( \delta _ { i j } - ( 1 + \delta _ { i j } ) \varphi _ { i j } \right ) \left ( y _ { i j } - \mathbf X _ { i j } ^ { T } \mathbb { E } _ { q ( \beta ) } ( \beta ) - \mathbb { E } _ { q ( \gamma _ { i j } ) } \gamma _ { i } \right ) \right ]$$

Let 𝛼 ∗ = 𝛼 0 + 𝛿 and

$$\omega ^ { * } = \omega _ { 0 } - \sum _ { i = 1 } ^ { K } \sum _ { j = 1 } ^ { n _ { i } } \left ( \delta _ { i j } - ( 1 + \delta _ { i j } ) \varphi _ { i j } \right ) \left ( y _ { i j } - \mathbf X _ { i j } ^ { T } \mathbb { E } _ { q } ( \beta ) ( \beta ) - \mathbb { E } _ { q ( \gamma _ { i } ) } \gamma _ { i } \right ),$$

we have 𝑞 ∗ ( 𝑏 ) is Inverse-Gamma ( 𝛼 , 𝜔 ∗ ∗ ) .

## (4) Update equation for 𝑞 ∗ ( 𝜎 2 𝛾 )

$$\log q ^ { * } ( \sigma _ { \gamma } ^ { 2 } ) \stackrel { + } { \approx } \mathbb { E } _ { - \sigma _ { \gamma } ^ { 2 } } \left [ \sum _ { i = 1 } ^ { K } \log p ( \gamma _ { i } \left | \sigma _ { \gamma } ^ { 2 } \right ) \right ] + \mathbb { E } _ { - \sigma _ { \gamma } ^ { 2 } } \left [ \log p ( \sigma _ { \gamma } ^ { 2 } ) \right ]$$

where

$$\mathbb { E } _ { - \sigma _ { \gamma } ^ { 2 } } \Big { [ } \sum _ { i = 1 } ^ { K } \log p ( \gamma _ { i } \, | \sigma _ { \gamma } ^ { 2 } ) \Big { ] } & = \mathbb { E } _ { - \sigma _ { \gamma } ^ { 2 } } \Big { [ } \sum _ { i = 1 } ^ { K } \log \Big { ( } \frac { 1 } { \sqrt { 2 \pi \sigma _ { \gamma } ^ { 2 } } } \exp ( - \frac { 1 } { 2 \sigma _ { \gamma } ^ { 2 } } \gamma _ { i } ^ { 2 } ) \Big { ) } \Big { ] } \\ & \quad \stackrel { + } { \approx } - \frac { K } { 2 } \log \sigma _ { \gamma } ^ { 2 } - \frac { 1 } { 2 \sigma _ { \gamma } ^ { 2 } } \sum _ { i = 1 } ^ { K } \mathbb { E } _ { q ( \gamma _ { i } ) } \gamma _ { i } ^ { 2 } \\. \, 1 \quad. \, \quad. \quad. \quad. \quad.$$

and 𝔼 - 2 𝜎 𝛾 [ log 𝑝 𝜎 ( 2 𝛾 ) ] = -( 𝜆 0 +1)log 𝜎 2 𝛾 -𝜂 0 ∕ 𝜎 2 𝛾 .

Therefore,

$$\log q ^ { * } ( \sigma _ { \gamma } ^ { 2 } ) \stackrel { + } { \approx } - ( \lambda _ { 0 } + \frac { K } { 2 } + 1 ) \log \sigma _ { \gamma } ^ { 2 } - [ \eta _ { 0 } + \frac { 1 } { 2 } \sum _ { i = 1 } ^ { K } \mathbb { E } _ { q ( \gamma _ { i } ) } \gamma _ { i } ^ { 2 } ] / \sigma _ { \gamma } ^ { 2 }$$

Let 𝜆 ∗ = 𝜆 0 + 𝐾 ∕2 and

$$\eta ^ { * } = \eta _ { 0 } + \frac { 1 } { 2 } \sum _ { i = 1 } ^ { K } \mathbb { E } _ { q ( \gamma _ { i } ) } \gamma _ { i } ^ { 2 },$$

𝑞 ∗ ( 𝜎 2 𝛾 ) is an Inverse-Gamma ( 𝜆 , 𝜂 ∗ ∗ ) .

## APPENDIX B. ELBO CALCULATION

We now present details of calculating diff 𝜷 + diff 𝜸 + diff 𝑏 + diff 𝜎 2 𝛾 in the ELBO defined in (5).

$$\text{diff} _ { \beta } & = \mathbb { E } _ { q } [ \log p ( \beta ) ] - \mathbb { E } _ { q } [ \log q ( \beta ) ] \\ & \stackrel { + } { \approx } \mathbb { E } _ { q } [ - \frac { 1 } { 2 } v _ { 0 } ( \beta - \mu _ { 0 } ) ^ { T } ( \beta - \mu _ { 0 } ) ] - \mathbb { E } _ { q } [ - \frac { 1 } { 2 } \log ( | \Sigma ^ { * } | ) - \frac { 1 } { 2 } ( \beta - \mu ^ { * } ) ^ { T } ( \Sigma ^ { * } ) ^ { - 1 } ( \beta - \mu ^ { * } ) ] \\ & \stackrel { + } { \approx } - \frac { 1 } { 2 } v _ { 0 } [ \text{trace} ( \Sigma ^ { * } ) + ( \mu ^ { * } - \mu _ { 0 } ) ^ { T } ( \mu ^ { * } - \mu _ { 0 } ) ] + \frac { 1 } { 2 } \log ( | \Sigma ^ { * } | ).$$

Note that

$$\mathbb { E } _ { q } [ \frac { 1 } { 2 } ( \beta - \mu ^ { * } ) ^ { T } ( \Sigma ^ { * } ) ^ { - 1 } ( \beta - \mu ^ { * } ) ] = \frac { p } { 2 },$$

which is always a constant at each iteration and therefore we ignore it.

$$\text{when is always a constant an aracian aracian aracian aracian.r} \\ \text{diff} \gamma & = \mathbb { E } _ { q } [ \sum _ { i = 1 } ^ { K } \log p ( \gamma _ { i } | \sigma _ { \gamma } ^ { 2 } ) ] - \mathbb { E } _ { q } [ \sum _ { i = 1 } ^ { K } \log q ( \gamma _ { i } ) ] \\ & \stackrel { + } { \sim } - \frac { K } { 2 } \mathbb { E } _ { q ( \sigma _ { \gamma } ^ { 2 } ) } ( \log \sigma _ { \gamma } ^ { 2 } ) - \frac { 1 } { 2 } \mathbb { E } _ { q ( \sigma _ { \gamma } ^ { 2 } ) } ( \frac { 1 } { \sigma _ { \gamma } ^ { 2 } ) } ) \sum _ { i = 1 } ^ { K } \mathbb { E } _ { q ( \gamma _ { i } ) } ( \gamma _ { i } ^ { 2 } ) - \frac { 1 } { 2 } \sum _ { i = 1 } ^ { K } \left [ \log \sigma _ { i } ^ { 2 * } - \frac { 1 } { \sigma _ { i } ^ { 2 * } } ( \mathbb { E } _ { q ( \gamma _ { i } ) } ( \gamma _ { i } ^ { 2 } ) - \mu _ { i } ^ { * 2 } \right ) \right ] \\ & \stackrel { + } { \sim } - \frac { K } { 2 } \mathbb { E } _ { q ( \sigma _ { \gamma } ^ { 2 } ) } ( \log \sigma _ { \gamma } ^ { 2 } ) - \frac { 1 } { 2 } \mathbb { E } _ { q ( \sigma _ { \gamma } ^ { 2 } ) } ( \frac { 1 } { \sigma _ { \gamma } ^ { 2 } ) } ) \sum _ { i = 1 } ^ { K } \mathbb { E } _ { q ( \gamma _ { i } ) } ( \gamma _ { i } ^ { 2 } ) - \frac { 1 } { 2 } \sum _ { i = 1 } ^ { K } \left ( \log \sigma _ { i } ^ { 2 * } \right ) \\ \text{Since} \left ( \mathbb { F } _ { 1 } \. \left ( \gamma _ { i } ^ { 2 } \right ) - u ^ { * 2 } \right ) / \sigma ^ { 2 * } = 1 \, \text{for} \, i = 1 \, \text{end} \, K \text{ and therefore}.$$

Since ( 𝔼 𝑞 𝛾 ( 𝑖 ) ( 𝛾 2 𝑖 ) -𝜇 𝑖 ∗2 ) ∕ 𝜎 2∗ 𝑖 = 1 for 𝑖 = 1 , ...𝐾 , and therefore,

$$\left ( \sum _ { i = 1 } ^ { K } \left ( \mathbb { E } _ { q ( \gamma _ { i } ) } ( \gamma _ { i } ^ { 2 } ) - \mu _ { i } ^ { * 2 } \right ) \right ) / ( 2 \sigma _ { i } ^ { 2 * } ) = K / 2$$

which is a constant, and we can ignore it.

$$\text{diff} _ { b } & = \mathbb { E } _ { q } [ \log p ( b ) ] - \mathbb { E } _ { q } [ \log q ( b ) ] \\ & \stackrel { * } { \sim } \mathbb { E } _ { q } \left [ - ( \alpha _ { 0 } + 1 ) \log b - \frac { \omega _ { 0 } } { b } \right ] \\ & \quad - \mathbb { E } _ { q } \left [ \alpha ^ { * } \log \omega ^ { * } - \log ( \Gamma ( \alpha ^ { * } ) ) - ( \alpha ^ { * } + 1 ) \log b - \frac { \omega ^ { * } } { b } \right ] \\ & = ( \alpha ^ { * } - \alpha _ { 0 } ) \mathbb { E } _ { q ( b ) } ( \log b ) + ( \omega ^ { * } - \omega _ { 0 } ) \mathbb { E } _ { q ( b ) } ( \frac { 1 } { b } ) - \alpha ^ { * } \log \omega ^ { * }.$$

Since 𝛼 ∗ does not change at each iteration, we remove log(Γ( 𝛼 ∗ )) in the calculation of the ELBO.

$$d i f f _ { \sigma _ { \gamma } ^ { 2 } } & = \mathbb { E } _ { q } [ \log p ( \sigma _ { \gamma } ^ { 2 } ) ] - \mathbb { E } _ { q } [ \log q ( \sigma _ { \gamma } ^ { 2 } ) ] \\ & \stackrel { * } { \approx } \mathbb { E } _ { q } \left [ - ( \lambda _ { 0 } + 1 ) \log \sigma _ { \gamma } ^ { 2 } - \frac { \eta _ { 0 } } { \sigma _ { \gamma } ^ { 2 } } \right ] \\ & \quad - \mathbb { E } _ { q } \left [ \lambda ^ { * } \log \eta ^ { * } - \log ( \Gamma ( \lambda ^ { * } ) ) - ( \lambda ^ { * } + 1 ) \log \sigma _ { \gamma } ^ { 2 } - \frac { \eta ^ { * } } { \sigma _ { \gamma } ^ { 2 } } \right ] \\ & = ( \lambda ^ { * } - \lambda _ { 0 } ) \mathbb { E } _ { q ( \sigma _ { \gamma } ^ { 2 } ) } ( \log \sigma _ { \gamma } ^ { 2 } ) + ( \eta ^ { * } - \eta _ { 0 } ) \mathbb { E } _ { q ( \sigma _ { \gamma } ^ { 2 } ) } ( \frac { 1 } { \sigma _ { \gamma } ^ { 2 } } ) - \lambda ^ { * } \log \eta ^ { * }. \\ \text{change at each iteration. we remove lose} ( \Gamma ( \lambda ^ { * } ) ) \text{ in the calculation of the FI.Bo.}$$

Since 𝜆 ∗ does not change at each iteration, we remove log(Γ( 𝜆 ∗ )) in the calculation of the ELBO.

## APPENDIX C. REAL DATA ANALYSIS

TABLE 4 Results for ICU ventilation duration analysis: Posterior means (Mean) with 95% credible intervals (95% Cred. Int.) from the VB algorithms and MCMCbased survregbayes , respectively. Point estimates (Est.) with 95% confidence interval (95% Conf. Int.) from the h-likelihood survregHL and survreg methods.

Without frailty

## With frailty

## survregVBfrailty survregHL survregbayes survreg survregVB Mean 95% Cred. Int. Est. 95% Conf. Int. Mean 95% Cred. Int. Est. 95% Conf. Int. Mean 95% Cred. Int.

-

[-0.171, -0.097]

[-0.013, 0.197]

-

[0.039, 0.119]

-0.134

[-0.375, -0.293]

0.092

-

[0.049, 0.252]

0.079

[-0.173, -0.098]

[-0.175, -0.099]

-0.334

[-0.011, 0.201]

-

0.151

[0.037, 0.117]

-0.135

-0.137

[-0.374, -0.292]

0.095

-

[0.043, 0.248]

0.077

[-0.172, -0.103]

[-0.174, -0.097]

-0.333

[-0.014, 0.211]

-

0.146

[0.064, 0.128]

-0.134

-0.135

[-0.342, -0.281]

0.096

-

[0.039, 0.191]

0.097

[-0.181, -0.086]

[-0.178, -0.101]

-0.312

[-0.016, 0.191]

-

0.115

[0.033, 0.135]

-0.134

-0.140

[-0.372, -0.266]

0.088

-

[-0.062, 0.186]

0.084

[-0.171, -0.097]

[-0.191, -0.094]

-0.319

[-0.009, 0.211]

-

0.062

[0.068, 0.127]

-0.134

-0.142

[-0.357, -0.296]

0.101

Downstream (ref)

[0.050, 0.193]

0.098

ED

[-0.178, -0.102]

-0.327

Home

0.122

Hospital

-0.140

OR

Other

Ward

Diagnosis

Admission Source

| [0.230, 0.299] [0.270, 0.335]           | [0.363, 0.410] [0.543, 0.635]   |
|-----------------------------------------|---------------------------------|
| 0.265 0.303                             | 0.386 0.589                     |
| [0.231, 0.300] [0.272, 0.336]           | [0.363, 0.410] [0.539, 0.631]   |
| 0.265 0.304 0.386                       | 0.585                           |
| - [0.226, 0.294] [0.250, 0.317] [0.349, | 0.394] [0.530, 0.622]           |
| - 0.260 0.280                           | 0.370 0.574                     |
| 0.310] 0.329]                           | 0.397] 0.586]                   |
| - 0.265 0.286                           | 0.366 0.527                     |
| 0.294] 0.318]                           |                                 |
| [0.255,                                 | 0.402] 0.621]                   |
|                                         | Other Trauma                    |
| - [0.226,                               | [0.356,                         |
|                                         | [0.529,                         |
| - 0.260 0.286                           |                                 |
|                                         | 0.379 0.575                     |
| Cardiovascular (ref) Gastrointestinal   |                                 |
| Neurological                            | Neurological                    |

[-0.064, 0.015]

-0.024

[-0.061, 0.017]

-0.022

-

[-0.058, 0.013]

-

-0.023

-

[-0.031, 0.069]

-

-0.019

-

[-0.060, 0.017]

-

-0.021

Medical (ref)

Patient type

[-0.317, -0.237]

-0.277

[-0.329, -0.249]

-0.289

-

[-0.329, -0.246]

-

-0.285

-

[-0.336, -0.226]

-

-0.281

-

[-0.326, -0.247]

-

Surgical

-0.287

No (ref)

Scheduled Admission

[-0.186, -0.102]

-0.144

[-0.183, -0.099]

-0.141

-

[-0.188, -0.104]

-

-0.146

-

[-0.218, -0.106]

-

-0.162

-

[-0.184, -0.102]

-

Yes

-0.143

No (ref)

Scheduled Surgery

Yes

Referral Medical (ref) ------Other -0.028 [-0.058, 0.001] -0.012 [-0.052, 0.027] -0.034 [-0.065, -0.004] -0.051 [-0.080, -0.021] -0.052 [-0.081, -0.022] Respirology 0.155 [0.107, 0.203] 0.112 [0.047, 0.177] 0.156 [0.108, 0.205] 0.165 [0.117, 0.214] 0.168 [0.120, 0.216] Surgical -0.060 [-0.094, -0.026] -0.052 [-0.098, -0.005] -0.061 [-0.099, -0.024] -0.048 [-0.083, -0.013] -0.049 [-0.085, -0.014]

## TABLE 4 Continued Results for ICU ventilation duration analysis: Posterior means (Mean) with 95% credible intervals (95% Cred. Int.) from the VB algorithms and MCMC-based survregbayes , respectively. Point estimates (Est.) with 95% confidence interval (95% Conf. Int.) from the h-likelihood survregHL and survreg methods.

Without frailty

With frailty survregVB

survreg survregbayes

survregHL

survregVBfrailty

95% Cred. Int.

Mean

95% Conf. Int.

Est.

95% Cred. Int.

Mean

95% Conf. Int.

Est.

95% Cred. Int.

[0.046, 0.078]

0.062

[0.047, 0.079]

0.063

-

[0.075, 0.105]

-

0.089

-

[0.069, 0.112]

-

0.091

-

[0.073, 0.104]

Mean

-

0.088

No (ref)

Other Interventions

[0.181, 0.220]

0.201

[0.178, 0.218]

0.198

-

[0.157, 0.198]

-

0.178

-

[0.137, 0.188]

-

0.163

-

[0.164, 0.203]

-

Yes

0.184

No (ref)

CVL

[0.181, 0.223]

0.202

[0.178, 0.220]

0.199

-

[0.172, 0.213]

-

0.193

-

[0.149, 0.205]

-

0.177

-

[0.177, 0.219]

-

Yes

0.198

No (ref)

AL

[0.479, 0.604]

0.542

[0.481, 0.606]

0.544

-

[0.473, 0.599]

-

0.535

-

[0.401, 0.552]

-

0.476

-

[0.471, 0.596]

-

Yes

0.533

No (ref)

IPM

[0.904, 1.085]

0.994

[0.902, 1.083]

0.993

-

[0.817, 0.996]

-

0.908

-

[0.691, 0.906]

-

0.798

-

[0.822, 1.001]

-

Yes

0.912

No (ref)

EMO

[0.254, 0.379]

0.317

[0.261, 0.387]

0.324

-

[0.274, 0.398]

-

0.340

-

[0.216, 0.377]

-

0.297

-

[0.281, 0.404]

-

Yes

0.343

No (ref)

IABP

[0.089, 0.143]

[0.027, 0.090]

0.116

0.058

[0.026, 0.089]

[0.089, 0.142]

0.115

0.057

-

[0.102, 0.158]

[0.040, 0.110]

-

0.130

0.077

-

[0.093, 0.162]

[0.022, 0.105]

-

0.128

0.063

-

[0.107, 0.159]

[0.044, 0.106]

-

Yes

0.133

0.075

18-39 (ref)

40-80

80

&gt;

Pre-LOS

Age

≤ 1 day (ref) ------2-7 days 0.038 [0.017, 0.059] 0.040 [0.010, 0.069] 0.038 [0.016, 0.060] 0.040 [0.019, 0.062] 0.040 [0.019, 0.062] ≥ 7 days 0.122 [0.100, 0.145] 0.140 [0.109, 0.171] 0.121 [0.100, 0.143] 0.130 [0.107, 0.153] 0.129 [0.106, 0.152]

TABLE 4 Continued Results for ICU ventilation duration analysis: Posterior means (Mean) with 95% credible intervals (95% Cred. Int.) from the VB algorithms and MCMC-based survregbayes , respectively. Point estimates (Est.) with 95% confidence interval (95% Conf. Int.) from the h-likelihood survregHL and survreg methods.

Without frailty

With frailty survregVB

survreg survregbayes

survregHL

## survregVBfrailty

95% Cred. Int.

Mean

95% Conf. Int.

Est.

95% Cred. Int.

Mean

95% Conf. Int.

Est.

95% Cred. Int.

[0.056, 0.108]

[0.156, 0.205]

[0.215, 0.276]

0.082

[0.071, 0.193]

0.180

0.245

[0.057, 0.109]

0.132

[0.158, 0.207]

[0.217, 0.278]

0.083

[0.075, 0.198]

0.182

0.248

-

[0.078, 0.130]

0.137

[0.186, 0.233]

-

[0.231, 0.293]

0.103

[0.070, 0.189]

0.210

-

0.262

[0.076, 0.148]

0.130

[0.191, 0.264]

-

[0.247, 0.332]

0.112

[0.122, 0.274]

0.228

-

0.290

[0.076, 0.127]

Mean

0.198

[0.187, 0.235]

-

[0.235, 0.295]

0.101

[0.072, 0.195]

0.211

None (ref)

0.265

Minimal

0.133

Mild

Moderate

Severe

MODS

30
Chengqian Xian Camila P. E. de Souza Wenqing He Felipe F. Rodrigues Renfang Tian
Fast variational Bayesian inference for correlated survival data: an application to invasive mechanical ventilation duration analysis
Correlated survival data are prevalent in various clinical settings and have
been extensively discussed in literature. One of the most common types of
correlated survival data is clustered survival data, where the survival times
from individuals in a cluster are associated. Our study is motivated by
invasive mechanical ventilation data from different intensive care units (ICUs)
in Ontario, Canada, forming multiple clusters. The survival times from patients
within the same ICU cluster are correlated. To address this association, we
introduce a shared frailty log-logistic accelerated failure time model that
accounts for intra-cluster correlation through a cluster-specific random
intercept. We present a novel, fast variational Bayes (VB) algorithm for
parameter inference and evaluate its performance using simulation studies
varying the number of clusters and their sizes. We further compare the
performance of our proposed VB algorithm with the h-likelihood method and a
Markov Chain Monte Carlo (MCMC) algorithm. The proposed algorithm delivers
satisfactory results and demonstrates computational efficiency over the MCMC
algorithm. We apply our method to the ICU ventilation data from Ontario to
investigate the ICU site random effect on ventilation duration.