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 × 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
#> # ℹ 28 more rows
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.
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.
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
.
We can get two summary plots: estimated hazard rates and estimated survival functions.
Baseline hazards
Survival functions
You can also get diagnostic plots for the simulated variables. Choose
the type of plot with the argument type
.