Skip to contents

Joint model for Bidirectional survival data using rstanarm

jmstB.Rd

The function fits joint model for survival data with two events. It utilizes the rstanarm package for obtaining the model parameter estimates.

Usage

jmstB(
  dtlong,
  dtsurv,
  longm,
  survm,
  timeVar,
  id,
  nchain = 1,
  refresh = 1000,
  BIGdata = FALSE,
  samplesize = 200
)

Arguments

dtlong

longitudinal data

dtsurv

survival data with two event status along with event time

longm

longitudinal model e.g. list(serBilir~drug * year+(year|id),serBilir ~ drug * year+(year|id))

survm

survival model e.g. list(Surv(years,status2)~drug,Surv(time_2,status_2)~drug+age)

timeVar

time variable

id

ID variable

nchain

number of MCMC chain

refresh

number of refresh sample

BIGdata

logical argument TRUE or FALSE

samplesize

samplesize for bigdata

Value

Estimated model parameters of Joint model with bidirectional survival data

References

Goodrich, B., et al. "rstanarm: Bayesian applied regression modeling via Stan. R package version 2.17. 4." Online< http://mc-stan. org (2018).

Bhattacharjee, A., Rajbongshi, B. K., & Vishwakarma, G. K. (2024). jmBIG: enhancing dynamic risk prediction and personalized medicine through joint modeling of longitudinal and survival data in big routinely collected data. BMC Medical Research Methodology, 24(1), 172.

Author

Atanu Bhattacharjee, Bhrigu Kumar Rajbongshi and Gajendra Kumar Vishwakarma

Examples

 # \donttest{
 ##
library(JMbayes2)
library(rstanarm)
#> Loading required package: Rcpp
#> This is rstanarm version 2.32.1
#> - See https://mc-stan.org/rstanarm/articles/priors for changes to default priors!
#> - Default priors may change, so it's safest to specify priors, even if equivalent to the defaults.
#> - For execution on a local, multicore CPU with excess RAM we recommend calling
#>   options(mc.cores = parallel::detectCores())
#> 
#> Attaching package: 'rstanarm'
#> The following objects are masked from 'package:FastJM':
#> 
#>     fixef, ranef
st_pbcid<-function(){
  new_pbcid<-pbc2.id
  new_pbcid$time_2<-rexp(n=nrow(pbc2.id),1/10)
  cen_time<-runif(nrow(pbc2.id),min(new_pbcid$time_2),max(new_pbcid$time_2))
  status_2<-ifelse(new_pbcid$time_2<cen_time,1,0)
  new_pbcid$status_2<-status_2
  new_pbcid$time_2<-ifelse(new_pbcid$time_2<cen_time,new_pbcid$time_2,cen_time)
  new_pbcid$time_2<-ifelse(new_pbcid$time_2<new_pbcid$years,new_pbcid$years,
                           new_pbcid$time_2)
  new_pbcid
}
new_pbc2id<-st_pbcid()
pbc2$status_2<-rep(new_pbc2id$status_2,times=data.frame(table(pbc2$id))$Freq)
pbc2$time_2<-rep(new_pbc2id$time_2,times=data.frame(table(pbc2$id))$Freq)
pbc2_new<-pbc2[pbc2$id%in%c(1:50),]
new_pbc2id<-new_pbc2id[new_pbc2id$id%in%c(1:50),]
model_jmstBdirect<-jmstB(
  dtlong=pbc2_new,
  dtsurv = new_pbc2id,
  longm=list(serBilir~drug*year+(year|id),albumin~drug+year+(year|id)),
  survm=list(Surv(years,status2)~drug,Surv(time_2,status_2)~drug),
  timeVar="year",
  id='id',
  refresh=400,
  nchain=1)
#> Fitting a univariate joint model.
#> 
#> Please note the warmup may be much slower than later iterations!
#> 
#> SAMPLING FOR MODEL 'jm' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 0.003416 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 34.16 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 69.702 seconds (Warm-up)
#> Chain 1:                68.198 seconds (Sampling)
#> Chain 1:                137.9 seconds (Total)
#> Chain 1: 
#> Fitting a univariate joint model.
#> 
#> Please note the warmup may be much slower than later iterations!
#> 
#> SAMPLING FOR MODEL 'jm' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 0.000343 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 3.43 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 23.903 seconds (Warm-up)
#> Chain 1:                15.1 seconds (Sampling)
#> Chain 1:                39.003 seconds (Total)
#> Chain 1: 
model_jmstBdirect
#> 
#>   Joint model using rstanarm
#>  ===========================================================
#>  Summary of first joint model with Event 1: 
#>  -------------------------------------------
#>  Longitudinal process:
#>                      Mean StDev   2.5% 97.5% Zvalue Pvalue
#> (Intercept)         3.297 0.885  1.457 5.167  3.724  0.000
#> drugD-penicil      -0.414 1.390 -3.288 2.556 -0.298  0.766
#> year                2.377 0.682  1.012 3.811  3.484  0.000
#> drugD-penicil:year -0.221 1.022 -2.142 1.871 -0.216  0.829
#> sigma               2.192 0.090  2.033 2.365 24.335  0.000
#> mean_PPD            3.509 0.153  3.218 3.799 22.947  0.000
#> 
#>  Survival process for event1: 
#>                         Mean StDev   2.5%  97.5%  Zvalue Pvalue
#> (Intercept)          -3.447 0.188 -3.818 -3.085 -18.374  0.000
#> drugD-penicil         0.072 0.375 -0.649  0.807   0.193  0.847
#> b-splines-coef1      -0.534 1.037 -2.669  1.248  -0.515  0.607
#> b-splines-coef2       0.276 0.871 -1.470  1.949   0.316  0.752
#> b-splines-coef3      -1.985 1.129 -4.309  0.007  -1.759  0.079
#> b-splines-coef4       1.324 1.523 -1.621  4.156   0.869  0.385
#> b-splines-coef5      -0.091 1.637 -3.496  2.923  -0.056  0.956
#> b-splines-coef6      -1.251 1.899 -5.326  1.725  -0.659  0.510
#> Assoc|Long1|etavalue  0.187 0.024  0.140  0.234   7.659  0.000
#> 
#>  Random effect covariance matrix :
#>  Groups Name              Std.Dev. Corr 
#>  id     Long1|(Intercept) 4.9882        
#>         Long1|year        3.8296   0.971
#> 
#>  Summary of second joint model with Event 2: 
#>  -------------------------------------------
#>  Longitudinal process:
#>                 Mean StDev   2.5%  97.5%  Zvalue
#> (Intercept)    3.585 0.083  3.423  3.746  43.326
#> drugD-penicil -0.093 0.145 -0.379  0.182  -0.643
#> year          -0.116 0.014 -0.143 -0.091  -8.079
#> sigma          0.293 0.012  0.270  0.317  23.736
#> mean_PPD       3.353 0.020  3.313  3.393 164.641
#> 
#>  Survival process for event1: 
#>                             Mean StDev   2.5%  97.5%  Zvalue Pvalue
#> Long1|drugD-penicil:year -0.221 1.022 -2.142  1.871  -0.216  0.829
#> (Intercept)              -3.447 0.188 -3.818 -3.085 -18.374  0.000
#> Long1|sigma               2.192 0.090  2.033  2.365  24.335  0.000
#> b-splines-coef1          -0.534 1.037 -2.669  1.248  -0.515  0.607
#> b-splines-coef2           0.276 0.871 -1.470  1.949   0.316  0.752
#> b-splines-coef3          -1.985 1.129 -4.309  0.007  -1.759  0.079
#> b-splines-coef4           1.324 1.523 -1.621  4.156   0.869  0.385
#> b-splines-coef5          -0.091 1.637 -3.496  2.923  -0.056  0.956
#> drugD-penicil             0.072 0.375 -0.649  0.807   0.193  0.847
#> 
#>  Random effect covariance matrix :
#>  Groups Name              Std.Dev. Corr 
#>  id     Long1|(Intercept) 0.457963      
#>         Long1|year        0.060634 0.359
##
 # }