Stanford heart transplant survival data
Tommi Härkänen
Plot the figures
Figures of Bite output can be generated running the splus.S file in R e.g. using the source("splus.S") function.
Another option is to drop the pdf(), dev.off() and q() functions from the file and including the source in a code block of an R markdown document.
cmd <- scan("splus.S", what="", sep="\n")
## Drop the function calls for the pdf output:
cmd <- cmd[!grepl("pdf\\(|dev.off()|q()", cmd)]
## Run the plotting commands:
eval(parse(text=paste(cmd, collapse="\n")))## [1] 1 1
## [1] 1 1
## [1] 1 1
Scalar parameters
Regression coefficients
Read the MCMC sample of the regression coefficients from files ending .tmp. Note that the results differ from those in the manual as the seed of the random number generator is not the same.
## Read the output data files and construct a matrix 'betas':
betas <- c()
for (i in dir(pattern="beta_"))
betas <- cbind(betas, scan(i))
dimnames(betas) <- list(1:nrow(betas), sub(".tmp", "", dir(pattern="beta_")))
## Convert the matrix into MCMC object of the coda package:
library(coda)
betas <- mcmc(betas)
## Summary of betas:
summary(betas)##
## Iterations = 1:5000
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 5000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## beta_hla 0.2438 0.4338 0.006136 0.006907
## beta_mism_score 0.5469 0.3304 0.004673 0.007745
## beta_mism -0.4109 0.1742 0.002463 0.005106
## beta_prevs1 -1.7674 1.0335 0.014616 0.015116
## beta_prevs2 -0.9014 0.4647 0.006572 0.007735
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## beta_hla -0.6040 -0.04732 0.2440 0.5430 1.09199
## beta_mism_score -0.1092 0.32276 0.5483 0.7731 1.19652
## beta_mism -0.7501 -0.52673 -0.4110 -0.2942 -0.06875
## beta_prevs1 -4.1563 -2.34183 -1.6475 -1.0496 -0.09648
## beta_prevs2 -1.8748 -1.20138 -0.8865 -0.5860 -0.03227## Summary of hazard ratios (HR):
summary(exp(betas))##
## Iterations = 1:5000
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 5000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## beta_hla 1.4007 0.6270 0.008867 0.009758
## beta_mism_score 1.8240 0.6120 0.008654 0.015234
## beta_mism 0.6732 0.1184 0.001674 0.003494
## beta_prevs1 0.2624 0.2423 0.003427 0.003668
## beta_prevs2 0.4505 0.2095 0.002963 0.003555
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## beta_hla 0.54660 0.95378 1.2764 1.7211 2.9802
## beta_mism_score 0.89659 1.38093 1.7304 2.1664 3.3086
## beta_mism 0.47232 0.59053 0.6630 0.7451 0.9336
## beta_prevs1 0.01567 0.09615 0.1925 0.3501 0.9080
## beta_prevs2 0.15338 0.30078 0.4121 0.5566 0.9682Hazard function values at predefined points
The expectation command can be applied to calculate posterior expectations of different functionals of parameters e.g. hazard functions. The save expectation command saves the values at the MCMC iterations. For example, expectation haz_l : value(lambda0(),t) range t (5, 70, 13); calculates the posterior expectations of \(\lambda_0(t)\) at 13 time points 5, 10, , 65, and the output can be found in file haz_l.0.tmp. The MCMC iterations are saved using the command save expectation haz_l;, and the values are store in file haz_l.0.sav.
## Read the output data files and construct a matrix 'betas':
haz <- NULL
for (i in paste0("haz_", c("l","f"), ".0.sav")) {
tmp <- read.table(i)
colnames(tmp) <- paste0(sub(".0.sav", "", i), "_", 1:ncol(tmp))
if (!is.null(haz))
haz <- cbind(haz, tmp) else haz <- tmp
}
## Convert the matrix into MCMC object of the coda package:
haz <- mcmc(haz)
## Summary of hazard functions:
summary(haz)##
## Iterations = 1:5000
## Thinning interval = 1
## Number of chains = 1
## Sample size per chain = 5000
##
## 1. Empirical mean and standard deviation for each variable,
## plus standard error of the mean:
##
## Mean SD Naive SE Time-series SE
## haz_l_1 1.2108 0.8883 0.012563 0.038038
## haz_l_2 1.3821 0.9971 0.014101 0.033404
## haz_l_3 1.6253 1.3872 0.019617 0.038893
## haz_l_4 1.8656 1.5579 0.022031 0.050622
## haz_l_5 0.9422 0.7988 0.011297 0.013858
## haz_l_6 0.6775 0.3235 0.004575 0.008472
## haz_l_7 0.8837 0.5880 0.008316 0.011886
## haz_l_8 2.1537 0.7384 0.010442 0.022793
## haz_l_9 2.2454 0.6194 0.008760 0.023203
## haz_l_10 2.6925 0.9430 0.013336 0.025520
## haz_l_11 3.4558 1.1776 0.016654 0.037169
## haz_l_12 4.0162 2.0240 0.028624 0.052069
## haz_l_13 4.3118 3.2120 0.045424 0.068258
## haz_f_1 1.5172 0.6970 0.009856 0.033256
## haz_f_2 1.5355 1.0100 0.014283 0.041099
## haz_f_3 0.3287 0.2731 0.003863 0.009047
## haz_f_4 0.3171 0.2376 0.003360 0.008808
## haz_f_5 0.2966 0.1863 0.002635 0.007146
## haz_f_6 0.2859 0.1553 0.002196 0.006373
## haz_f_7 0.2835 0.1473 0.002083 0.006299
## haz_f_8 0.2808 0.1453 0.002055 0.006215
## haz_f_9 0.2771 0.1433 0.002027 0.006152
## haz_f_10 0.2692 0.1397 0.001976 0.006063
##
## 2. Quantiles for each variable:
##
## 2.5% 25% 50% 75% 97.5%
## haz_l_1 0.18681 0.6004 0.9591 1.5961 3.4081
## haz_l_2 0.24282 0.7005 1.1187 1.8294 3.8959
## haz_l_3 0.29118 0.7513 1.2437 2.0785 5.0741
## haz_l_4 0.39141 0.8456 1.3864 2.3400 6.2095
## haz_l_5 0.19972 0.5083 0.7371 1.1049 3.1178
## haz_l_6 0.21786 0.4511 0.6217 0.8427 1.4867
## haz_l_7 0.24713 0.5114 0.7104 1.0306 2.6182
## haz_l_8 0.78439 1.6770 2.1202 2.5864 3.7960
## haz_l_9 1.21538 1.8155 2.1778 2.6234 3.6853
## haz_l_10 1.36445 2.0495 2.5228 3.1274 5.1596
## haz_l_11 1.75572 2.5935 3.2469 4.0892 6.2086
## haz_l_12 1.68075 2.7167 3.5700 4.6546 9.4418
## haz_l_13 0.60246 2.5507 3.5819 5.0914 12.3359
## haz_f_1 0.47327 1.0455 1.3871 1.8940 3.0498
## haz_f_2 0.17476 0.9155 1.3797 1.9874 3.8091
## haz_f_3 0.09097 0.1869 0.2682 0.3761 1.0611
## haz_f_4 0.09348 0.1878 0.2674 0.3710 0.8602
## haz_f_5 0.09243 0.1845 0.2607 0.3617 0.6994
## haz_f_6 0.08970 0.1804 0.2556 0.3543 0.6576
## haz_f_7 0.08970 0.1799 0.2539 0.3534 0.6536
## haz_f_8 0.08943 0.1784 0.2517 0.3512 0.6455
## haz_f_9 0.08724 0.1757 0.2486 0.3476 0.6404
## haz_f_10 0.08025 0.1701 0.2411 0.3377 0.6155