Load packages
library(tidyverse)
library(rstanarm)
Load datasets
## multiple-row per patient longitudinal biomarker information
data(pbcLong, package = "survival")
## Warning in data(pbcLong, package = "survival"): data set 'pbcLong' not found
pbcLong <- as_data_frame(pbcLong)
pbcLong
## # A tibble: 304 x 8
## id age sex trt year logBili albumin platelet
## * <int> <dbl> <fct> <int> <dbl> <dbl> <dbl> <int>
## 1 1 58.8 f 1 0 2.67 2.6 190
## 2 1 58.8 f 1 0.526 3.06 2.94 183
## 3 2 56.4 f 1 0 0.0953 4.14 221
## 4 2 56.4 f 1 0.498 -0.223 3.6 188
## 5 2 56.4 f 1 0.999 0 3.55 161
## 6 2 56.4 f 1 2.10 0.642 3.92 122
## 7 2 56.4 f 1 4.90 0.956 3.32 135
## 8 2 56.4 f 1 5.89 1.28 2.92 100
## 9 2 56.4 f 1 6.89 1.44 2.73 103
## 10 2 56.4 f 1 7.89 1.28 2.8 113
## # ... with 294 more rows
## single-row per patient survival information
data(pbcSurv, package = "survival")
## Warning in data(pbcSurv, package = "survival"): data set 'pbcSurv' not found
pbcSurv <- as_data_frame(pbcSurv)
pbcSurv
## # A tibble: 40 x 7
## id age sex trt futimeYears status death
## * <int> <dbl> <fct> <int> <dbl> <int> <dbl>
## 1 1 58.8 f 1 1.10 2 1
## 2 2 56.4 f 1 14.2 0 0
## 3 3 70.1 m 1 2.77 2 1
## 4 4 54.7 f 1 5.27 2 1
## 5 5 38.1 f 0 4.12 1 0
## 6 6 66.3 f 0 6.85 2 1
## 7 7 55.5 f 0 6.85 0 0
## 8 8 53.1 f 0 6.75 2 1
## 9 9 42.5 f 1 6.57 2 1
## 10 10 70.6 f 0 0.140 2 1
## # ... with 30 more rows
Univariate joint model (current value association structure)
mod1 <- stan_jm(
## Longitudinal submodel
formulaLong = logBili ~ sex + trt + year + (year | id),
dataLong = pbcLong,
## Event submodel
formulaEvent = survival::Surv(futimeYears, death) ~ sex + trt,
dataEvent = pbcSurv,
## Name of the variable in ‘dataLong’ which represents time.
time_var = "year",
## MCMC specifications
chains = 1, refresh = 2000, seed = 12345)
## Fitting a univariate joint model.
##
## Please note the warmup may be much slower than later iterations!
##
## SAMPLING FOR MODEL 'jm' NOW (CHAIN 1).
##
## Gradient evaluation took 0.000492 seconds
## 1000 transitions using 10 leapfrog steps per transition would take 4.92 seconds.
## Adjust your expectations accordingly!
##
##
## Iteration: 1 / 2000 [ 0%] (Warmup)
## Iteration: 1001 / 2000 [ 50%] (Sampling)
## Iteration: 2000 / 2000 [100%] (Sampling)
##
## Elapsed Time: 21.6344 seconds (Warm-up)
## 22.0055 seconds (Sampling)
## 43.6399 seconds (Total)
Priors
prior_summary(mod1)
## Priors for model 'mod1'
## ------
## Long1|Intercept (after predictors centered)
## ~ normal(location = 0, scale = 10)
## **adjusted scale = 10.86
##
## Long1|Coefficients
## ~ normal(location = [0,0,0], scale = [2.5,2.5,2.5])
## **adjusted scale = [2.72,2.72,0.76]
##
## Long1|Auxiliary (sigma)
## ~ half-cauchy(location = 0, scale = 5)
## **adjusted scale = 5.43
##
## Event|Coefficients
## ~ normal(location = [0,0], scale = [2.5,2.5])
##
## Event|Auxiliary (spline-coefficients)
## ~ cauchy(location = [0,0,0,...], scale = [20,20,20,...])
##
## Association parameters
## ~ normal(location = 0, scale = 2.5)
## **adjusted scale = 2.20
##
## Covariance
## ~ lkj(reg. = 1, conc. = NULL, shape = NULL, scale = 10)
## ------
## See help('prior_summary.stanreg') for more details
print(mod1)
## stan_jm
## formula (Long1): logBili ~ sex + trt + year + (year | id)
## family (Long1): gaussian [identity]
## formula (Event): survival::Surv(futimeYears, death) ~ sex + trt
## baseline hazard: bs
## assoc: etavalue (Long1)
## ------
##
## Longitudinal submodel: logBili
## Median MAD_SD
## (Intercept) 0.268 0.592
## sexf 0.537 0.542
## trt -0.146 0.374
## year 0.208 0.044
## sigma 0.354 0.017
##
## Event submodel:
## Median MAD_SD exp(Median)
## (Intercept) -3.163 0.664 0.042
## sexf -0.334 0.661 0.716
## trt -0.740 0.459 0.477
## Long1|etavalue 1.351 0.225 3.860
## b-splines-coef1 -0.920 1.034 NA
## b-splines-coef2 0.502 0.853 NA
## b-splines-coef3 -1.814 1.163 NA
## b-splines-coef4 0.377 1.557 NA
## b-splines-coef5 -0.094 1.664 NA
## b-splines-coef6 -0.823 1.770 NA
##
## Group-level error terms:
## Groups Name Std.Dev. Corr
## id Long1|(Intercept) 1.2884
## Long1|year 0.1847 0.49
## Num. levels: id 40
##
## Sample avg. posterior predictive distribution
## of longitudinal outcomes:
## Median MAD_SD
## Long1|mean_PPD 0.586 0.030
##
## ------
## For info on the priors used see help('prior_summary.stanreg').
summary(mod1)
##
## Model Info:
##
## function: stan_jm
## formula (Long1): logBili ~ sex + trt + year + (year | id)
## family (Long1): gaussian [identity]
## formula (Event): survival::Surv(futimeYears, death) ~ sex + trt
## baseline hazard: bs
## assoc: etavalue (Long1)
## algorithm: sampling
## priors: see help('prior_summary')
## sample: 1000 (posterior sample size)
## num obs: 304 (Long1)
## num subjects: 40
## num events: 29 (72.5%)
## groups: id (40)
## runtime: 0.8 mins
##
## Estimates:
## mean sd 2.5% 25% 50% 75% 97.5%
## Long1|(Intercept) 0.251 0.609 -0.988 -0.135 0.268 0.662 1.465
## Long1|sexf 0.538 0.581 -0.528 0.139 0.537 0.881 1.745
## Long1|trt -0.153 0.405 -0.956 -0.410 -0.146 0.092 0.663
## Long1|year 0.211 0.042 0.140 0.181 0.208 0.240 0.298
## Long1|sigma 0.355 0.016 0.326 0.343 0.354 0.366 0.387
## Long1|mean_PPD 0.587 0.029 0.530 0.567 0.586 0.607 0.644
## Event|(Intercept) -3.194 0.646 -4.600 -3.616 -3.163 -2.725 -2.113
## Event|sexf -0.312 0.642 -1.486 -0.780 -0.334 0.111 1.030
## Event|trt -0.749 0.458 -1.680 -1.042 -0.740 -0.426 0.103
## Event|b-splines-coef1 -0.996 1.014 -2.985 -1.693 -0.920 -0.287 0.852
## Event|b-splines-coef2 0.477 0.886 -1.327 -0.080 0.502 1.071 2.101
## Event|b-splines-coef3 -1.832 1.159 -4.123 -2.611 -1.814 -1.061 0.486
## Event|b-splines-coef4 0.356 1.623 -3.019 -0.667 0.377 1.430 3.386
## Event|b-splines-coef5 -0.093 1.637 -3.369 -1.176 -0.094 1.074 3.013
## Event|b-splines-coef6 -1.055 1.875 -5.390 -2.132 -0.823 0.244 1.890
## Assoc|Long1|etavalue 1.363 0.232 0.934 1.206 1.351 1.504 1.855
## Sigma[id:Long1|(Intercept),Long1|(Intercept)] 1.660 0.435 1.004 1.364 1.588 1.891 2.680
## Sigma[id:Long1|year,Long1|(Intercept)] 0.117 0.066 0.004 0.071 0.109 0.155 0.266
## Sigma[id:Long1|year,Long1|year] 0.034 0.016 0.014 0.023 0.030 0.041 0.077
## log-posterior -328.214 10.133 -349.416 -334.756 -328.072 -321.189 -309.559
##
## Diagnostics:
## mcse Rhat n_eff
## Long1|(Intercept) 0.037 1.006 277
## Long1|sexf 0.029 1.000 405
## Long1|trt 0.030 1.013 182
## Long1|year 0.003 0.999 233
## Long1|sigma 0.001 1.000 712
## Long1|mean_PPD 0.001 0.999 1000
## Event|(Intercept) 0.020 0.999 1000
## Event|sexf 0.020 0.999 1000
## Event|trt 0.014 0.999 1000
## Event|b-splines-coef1 0.032 0.999 1000
## Event|b-splines-coef2 0.028 0.999 1000
## Event|b-splines-coef3 0.045 0.999 652
## Event|b-splines-coef4 0.062 1.000 690
## Event|b-splines-coef5 0.060 0.999 742
## Event|b-splines-coef6 0.065 0.999 834
## Assoc|Long1|etavalue 0.007 0.999 1000
## Sigma[id:Long1|(Intercept),Long1|(Intercept)] 0.028 1.007 243
## Sigma[id:Long1|year,Long1|(Intercept)] 0.004 1.002 268
## Sigma[id:Long1|year,Long1|year] 0.001 1.008 357
## log-posterior 0.860 1.031 139
##
## For each parameter, mcse is Monte Carlo standard error, n_eff is a crude measure of effective sample size, and Rhat is the potential scale reduction factor on split chains (at convergence Rhat=1).
plot(mod1)
