Tutorial 7

require(tidyverse)
library(brms)
library(bayesplot)
library(DT)
library(gtsummary)
library(survival)

3.3 Bayesian parametric survival modelling with Brms

# using default priors;
lung$censored <- ifelse(lung$status == 1, 'right', 'none')

m1a <- brm(time | cens(censored) ~ age+sex,
           data=lung, 
           prior=c(set_prior("normal(0,10)",class="Intercept"),
                   set_prior("normal(0,5)",class="b"),
                   set_prior("cauchy(0,25)", class="shape")),
           family='weibull')

saveRDS(m1a, "data/surv_wb")

m1a <- readRDS("data/surv_wb")
print(m1a)

 Family: weibull 
  Links: mu = log; shape = identity 
Formula: time | cens(censored) ~ age + sex 
   Data: lung (Number of observations: 228) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Regression Coefficients:
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept     6.21      0.48     5.29     7.19 1.00     4242     3136
age          -0.01      0.01    -0.03     0.00 1.00     4287     3167
sex           0.39      0.13     0.15     0.65 1.00     3979     2812

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
shape     1.31      0.08     1.16     1.47 1.00     3708     2799

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

plot(m1a)

prior_summary(m1a)

        prior     class coef group resp dpar nlpar lb ub       source
  normal(0,5)         b                                          user
  normal(0,5)         b  age                             (vectorized)
  normal(0,5)         b  sex                             (vectorized)
 normal(0,10) Intercept                                          user
 cauchy(0,25)     shape                             0            user

s1 <- posterior_samples(m1a)
beta <- data.matrix(s1[,c("b_age","b_sex")])

my.summary <- function(x){
  round(c(Mean=mean(x), Median=median(x), low=quantile(x,0.025), high=quantile(x, 0.975)),2)
}

t(apply(exp(beta),2,my.summary))

      Mean Median low.2.5% high.97.5%
b_age 0.99   0.99     0.97       1.00
b_sex 1.49   1.47     1.16       1.92

# The conditional effect of age on mean survival time 
conditional_effects(m1a, effects='age')

# A posterior predictive check: 
# data are dark points, posterior medians are large points;
pp_check(m1a, type='intervals', subset=TRUE)

m1b <- brm(time | cens(censored) ~ age+sex,
           data=lung, 
           prior=c(set_prior("normal(0,10)",class="Intercept"),
                   set_prior("normal(0,5)",class="b")),
           family='exponential')

saveRDS(m1b, "data/surv_exp")

# using default weakly informative priors;
m1b <- readRDS("data/surv_exp")
print(m1b)

 Family: exponential 
  Links: mu = log 
Formula: time | cens(censored) ~ age + sex 
   Data: lung (Number of observations: 228) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Regression Coefficients:
          Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
Intercept     6.37      0.63     5.14     7.63 1.00     4249     2720
age          -0.02      0.01    -0.03     0.00 1.00     4542     3115
sex           0.48      0.17     0.15     0.81 1.00     3681     2870

Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
and Tail_ESS are effective sample size measures, and Rhat is the potential
scale reduction factor on split chains (at convergence, Rhat = 1).

plot(m1b)

s2 <- posterior_samples(m1b)
beta <- data.matrix(s2[,c("b_age","b_sex")])

my.summary <- function(x){
  round(c(Mean=mean(x), Median=median(x), low=quantile(x,0.025), high=quantile(x, 0.975)),2)
}

t(apply(exp(beta),2,my.summary))

      Mean Median low.2.5% high.97.5%
b_age 0.98   0.98     0.97       1.00
b_sex 1.64   1.62     1.16       2.26

# The conditional effect of age on mean survival time 
conditional_effects(m1a, effects='age')

# A posterior predictive check: 
# data are dark points, posterior medians are large points;
pp_check(m1a, type='intervals', subset=TRUE)

• Comparing between the models using loo

loo1 <- loo(m1a) # better one!
loo2 <- loo(m1b) 
loo_compare(loo1, loo2)

    elpd_diff se_diff
m1a  0.0       0.0   
m1b -7.5       4.1