Bivariate Model Example

library(BGPhazard)
library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
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
#> # A tibble: 38 x 6
#>       id   sex    t1    t2 delta1 delta2
#>    <dbl> <dbl> <dbl> <dbl>  <dbl>  <dbl>
#>  1     1     1     8    16      1      1
#>  2     2     0    23    13      1      0
#>  3     3     1    22    28      1      1
#>  4     4     0   447   318      1      1
#>  5     5     1    30    12      1      1
#>  6     6     0    24   245      1      1
#>  7     7     1     7     9      1      1
#>  8     8     0   511    30      1      1
#>  9     9     0    53   196      1      1
#> 10    10     1    15   154      1      1
#> # … with 28 more rows

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)
#>  
#>  Individuals: 38
#>  Time partition intervals: 57
#>  Censored t1: 6
#>  Censored t2: 12
#>  Predictors: TRUE     sex

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)
#>  
#>  Samples: 90
#>  Individuals: 38
#>  Time partition intervals: 57
#>  Predictors: TRUE

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")
#>    Individual      Mean Prob. Low 95% Prob. High 95% Acceptance Rate
#> 1           1 3.5161998    2.93559241      3.9604323      0.04494382
#> 2           2 1.1065573    0.18571302      1.6881673      0.06741573
#> 3           3 0.7416358    0.36037618      1.5592198      0.15730337
#> 4           4 1.2766251    0.47604867      2.3175307      0.03370787
#> 5           5 0.6756267    0.38536523      1.3612876      0.11235955
#> 6           6 0.5794627    0.13703276      0.7149068      0.08988764
#> 7           7 0.3970680    0.14758059      0.7535655      0.10112360
#> 8           8 1.4426608    1.34200170      1.5006374      0.02247191
#> 9           9 0.8038096    0.46468569      1.4004952      0.10112360
#> 10         10 2.7630778    2.05615796      3.6564398      0.06741573
#> 11         11 0.8593926    0.40157655      3.5427683      0.12359551
#> 12         12 2.8803886    2.89006674      2.8900667      0.01123596
#> 13         13 0.7382331    0.44092970      0.9362285      0.06741573
#> 14         14 1.2081150    0.86413885      1.4607704      0.07865169
#> 15         15 1.7902202    1.29562801      2.2319309      0.06741573
#> 16         16 1.3627014    0.82697995      1.9902757      0.04494382
#> 17         17 0.5523012    0.33161142      1.3102429      0.03370787
#> 18         18 0.5630820    0.49259116      0.7187072      0.04494382
#> 19         19 1.7544980    0.17037680      2.8154057      0.08988764
#> 20         20 2.0120744    1.97106409      2.0230490      0.01123596
#> 21         21 8.2350556    1.85218433     12.9039642      0.11235955
#> 22         22 0.6949245    0.48424296      0.8675818      0.06741573
#> 23         23 0.2214919    0.08392286      0.6664909      0.15730337
#> 24         24 0.2812495    0.08070932      1.2778070      0.13483146
#> 25         25 1.2771252    0.67825416      2.5472605      0.11235955
#> 26         26 2.8058518    1.83041665      3.3339332      0.06741573
#> 27         27 0.8531273    0.47196104      1.7233703      0.06741573
#> 28         28 0.8204746    0.52322421      1.4459134      0.08988764
#> 29         29 0.9169544    0.30734051      1.4592335      0.10112360
#> 30         30 0.5382379    0.22762565      0.7643319      0.07865169
#> 31         31 0.8629524    0.09194932      2.5330043      0.08988764
#> 32         32 0.7209464    0.19416043      1.3089953      0.17977528
#> 33         33 2.1644123    1.65823576      2.8291761      0.04494382
#> 34         34 1.3729098    0.29452080      2.4280549      0.08988764
#> 35         35 1.7693823    1.09815920      2.5212324      0.10112360
#> 36         36 1.6086777    0.78312506      1.8237267      0.05617978
#> 37         37 0.9144461    0.43560334      1.5250703      0.11235955
#> 38         38 0.7551478    0.54353895      1.1552375      0.05617978
BSBSumm(samples, "lambda1")
#>    Interval         Mean Prob. Low 95% Prob. High 95%
#> 1         1 0.0025857030  1.095838e-03    0.004529922
#> 2         2 0.0028201306  1.072372e-03    0.006209417
#> 3         3 0.0028701368  1.028741e-03    0.004577725
#> 4         4 0.0014738311  4.456613e-04    0.002929883
#> 5         5 0.0010384908  1.513226e-04    0.002427628
#> 6         6 0.0015676796  2.653196e-04    0.003885577
#> 7         7 0.0016694291  2.945030e-04    0.003660394
#> 8         8 0.0009113779  1.018708e-04    0.002452031
#> 9         9 0.0009840101  1.035397e-04    0.002962389
#> 10       10 0.0014811821  9.997731e-05    0.003422372
#> 11       11 0.0010404547  4.245591e-05    0.002902399
#> 12       12 0.0015906226  7.036417e-05    0.004392124
#> 13       13 0.0021560652  8.016327e-05    0.006407150
#> 14       14 0.0019912459  1.563375e-05    0.004534501
#> 15       15 0.0018077036  3.339335e-05    0.004162587
#> 16       16 0.0022137172  5.657934e-05    0.006073974
#> 17       17 0.0013511934  1.000000e-05    0.003463942
#> 18       18 0.0013227631  1.000000e-05    0.003636501
#> 19       19 0.0017753832  2.399019e-05    0.004081210
#> 20       20 0.0013563878  1.000000e-05    0.003769219
#> 21       21 0.0010545479  1.000000e-05    0.003959654
#> 22       22 0.0013524100  1.000000e-05    0.003525046
#> 23       23 0.0013050529  1.000000e-05    0.004016113
#> 24       24 0.0013235438  1.000000e-05    0.004210845
#> 25       25 0.0011104413  1.000000e-05    0.003055708
#> 26       26 0.0007253316  1.000000e-05    0.002132804
#> 27       27 0.0009401048  1.000000e-05    0.002904711
#> 28       28 0.0010058183  1.000000e-05    0.002763355
#> 29       29 0.0012358063  1.000000e-05    0.004230714
#> 30       30 0.0018967682  1.828473e-05    0.004760684
#> 31       31 0.0018014437  5.402194e-05    0.004791687
#> 32       32 0.0011479481  8.625434e-05    0.003700070
#> 33       33 0.0011435303  2.343359e-04    0.002767294
#> 34       34 0.0012261302  7.771174e-05    0.003636052
#> 35       35 0.0013919237  1.270542e-04    0.004347896
#> 36       36 0.0012904054  1.361895e-04    0.003190328
#> 37       37 0.0015285924  1.837256e-04    0.004401775
#> 38       38 0.0015326806  1.507929e-04    0.004137058
#> 39       39 0.0019105030  5.312515e-04    0.004687678
#> 40       40 0.0026009009  3.822327e-04    0.006309135
#> 41       41 0.0033056801  9.513998e-04    0.007796441
#> 42       42 0.0039823749  1.055664e-03    0.008110874
#> 43       43 0.0042108617  1.282347e-03    0.007878363
#> 44       44 0.0045044175  1.930756e-03    0.007491707
#> 45       45 0.0052750300  2.019230e-03    0.009657346
#> 46       46 0.0049449366  1.368544e-03    0.009301134
#> 47       47 0.0045744088  2.421738e-04    0.011400391
#> 48       48 0.0047951245  7.088428e-04    0.011686684
#> 49       49 0.0053951413  1.780167e-03    0.011833267
#> 50       50 0.0070323069  3.014840e-03    0.011168967
#> 51       51 0.0088950441  2.929981e-03    0.015553042
#> 52       52 0.0109052593  4.808901e-03    0.018417408
#> 53       53 0.0121661018  4.616973e-03    0.026259576
#> 54       54 0.0145331873  3.989288e-03    0.029947030
#> 55       55 0.0154093520  4.827066e-03    0.034733410
#> 56       56 0.0139422058  2.577616e-03    0.030471030
#> 57       57 0.0124603843  3.500266e-03    0.024895142

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")