Bivariate Model Example

In BGPhazard: Markov Beta and Gamma Processes for Modeling Hazard Rates

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.width=8,
  fig.height=6,
  fig.align = "center"
)

library(BGPhazard)
library(dplyr)
library(ggplot2)

We will use the built-in dataset KIDNEY to show how the bivariate model functions work. All the functions for the bivariate model start with the letters BSB, which stand for Bayesian Semiparametric Bivariate.

KIDNEY

Initial setup

First, we use the BSBInit function to create the necessary data structure that we have to feed the Gibbs Sampler. We can skim the data structure with the summary and print methods.

bsb_init <- BSBInit(
  KIDNEY,
  alpha = 0.001,
  beta = 0.001,
  c = 1000,
  part_len = 10,
  seed = 42
  )

summary(bsb_init)

Our data consists of 38 individuals with two failure times each. For the first failure time t1 we have six censored observations, while for the second failure time we have twelve. The model will use sex as a predictor variable.

Gibbs Sampler

To obtain the posterior samples, we use the function BSBHaz. We run 100 iterations with a burn-in period of 10. The number of simulations is low in order to reduce the complexity of building this vignette. In practice, you should see how many iterations the model needs to reach convergence.

samples <- BSBHaz(
  bsb_init,
  iter = 100,
  burn_in = 10,
  gamma_d = 0.6, 
  theta_d = 0.3, 
  seed = 42
)

print(samples)

The print method shows that we only kept the last 90 iterations as posterior simulations.

Summaries

Tables

We can get posterior sample summaries with the function get_summaries. This function returns the posterior mean and a 0.95 probability interval for all the model parameters. Additionally, it returns the acceptance rate for variables sampled using the Metropolis-Hastings algorithm.

BSBSumm(samples, "omega1")
BSBSumm(samples, "lambda1")

It is important to notice that lambda1 and lambda2 are the estimated hazard rates for the baseline hazards $h_0$. They do not include the effect of predictor variables. The same applies for the survival function estimates s1 and s2.

Plots

We can get two summary plots: estimated hazard rates and estimated survival functions.

Baseline hazards

BSBPlotSumm(samples, "lambda1")
BSBPlotSumm(samples, "lambda2")

Survival functions

BSBPlotSumm(samples, "s1")
BSBPlotSumm(samples, "s2")

You can also get diagnostic plots for the simulated variables. Choose the type of plot with the argument type.

BSBPlotDiag(samples, "omega1", type = "traceplot")
BSBPlotDiag(samples, "omega1", type = "ergodic_means")

BGPhazard documentation built on Sept. 3, 2023, 5:09 p.m.