#' Performs the DIC-tau_g procedure and returns the posterior quantities of the optimal model.
#'
#' Performs the DIC-tau_g procedure by first running the function SCRSELECTRUN with 60 percent burnin, which performs SVSS on two disperse starting values for beta1,beta2,beta3. Afterwards, the function DICTAUG is used to extract the DIC values for unique models visited by the grid search and the optimal model is determined as the one with the lowest DIC which is the most parsimonius. After the optimal model is determined, one final MCMC is performed to obtain posterior beta1,beta2 and beta3 quantities for this model, returning summary values for each hazard.
#'
#' @importFrom graphics par plot
#' @importFrom stats dgamma dnorm dpois rgamma rnorm runif median quantile
#' @importFrom utils write.table
#' @importFrom mvtnorm rmvnorm dmvnorm
#' @param Y1 Vector Containing non-terminal event times (or censoring time due to death/censoring)
#' @param I1 Vector Containing non-terminal event indicators (1 if non-terminal event for a patient, 0 otherwise)
#' @param Y2 Vector Containing Terminal Event times (or censoring)
#' @param I2 Vector Containing Terminal event indicators (1 if a patients experiences a non-ternminal event, 0 if censored)
#' @param X Matrix of Patient Covariates. The last inc will be left out of variable selection.
#' @param hyperparameters List containing 29 hyperparameters and four starting values. In order they are: psi-the swap rate of the SVSS algorithm.
#' c-parameter involved in Sigma matrix for selection. z1a, z1b, z2a, z2b, z3a, z3b - beta hyper parameters on probability of inclusion for each of the three hazard functions.
#' a1,b1,a2,b2,a3,b3- hyperparameters on sigma_lambda_1, sigma_lambda_2, and sigma_lambda_3.
#' clam1, clam2, clam3 - spatial dependency of baseline hazard (between 0 and 1) for the three hazard functions.
#' Alpha1, Alpha2, Alpha3 - The parameter for the number of split points in hazards 1,2 and 3 (must be whole number).
#' J1max, J2max, J3max - Maximum number of split points allowed (must be whole number).
#' J1, J2, J3- Starting number of split points. w, psi1- hyperparameters on theta^{-1}. cep=Tuning Parameter for theta^{-1} sampler.
#' epstart-Starting value for theta^{-1}. cl1,cl2,cl3-Tuning parameters for log baseline hazard height sampler.
#' @param c sparsity parameter involved in Sigma matrix for selection. This should be the same c as that used in the hyperparameters vector.
#' @param BSVSS Number of iterations to perform during the SVSS procedure. 100,000 is a reccomended value to achieve convergence.
#' @param BDIC Number of iterations to perform during the DIC-tau_g grid search. 10,000 is a reccomended value to achieve convergence in a reasonable amount of time.
#' @param inc Number of variables left out of selection.
#' @param Path Where to save posterior coefficient samples for the optimal model.
#' @return Returns the optimal model determined by the DIC-Tau_g procedure and it's DIC along with summaries of these posterior quantities. Additionally, this function saves these posterior samples to a desired path.
#'
#' @references
#' [1] Lee, K. H., Haneuse, S., Schrag, D. and Dominici, F. (2015), Bayesian semi-parametric analysis of semi-competing risks data: investigating hospital readmission after a pancreatic cancer diagnosis. Journal of the Royal Statistical Society: Series C (Applied Statistics), 64: 253-273. doi: 10.1111/rssc.12078
#' [2] Chapple, A.C., Vannucci, M., Thall, P.F., Lin, S.(2017), Bayesian Variable selection for a semi-competing risks model with three hazard functions. Journal of Computational Statistics & Data Analysis, Volume 112, August 2017, Pages 170-185
#' [3] https://adventuresinstatistics.wordpress.com/2017/04/10/package-scrselect-using-returnmodel/
#'
#' @examples
#' ####Randomly Generate Semicompeting Risks Data
#' set.seed(1)
#' ####Generates random patient time, indicator and covariates.
#' n=100
#' Y1=runif(n,0,100)
#' I1=rbinom(n,1,.5)
#' Y2=Y1
#' I2=I1
#' for(i in 1:n){if(I1[i]==0){Y2[i]=Y1[i]}else{Y2[i]=Y1[i]+runif(1,0,100)}}
#' I2=rbinom(n,1,.5)
#' library(mvtnorm)
#' X=rmvnorm(n,rep(0,7),diag(7))
#' ####Read in Hyperparameters
#' ##Swap Rate
#' psi=.5
#' c=5
#' ###Eta Beta function probabilities
#' z1a=.4
#' z1b=1.6
#' z2a=.4
#' z2b=1.6
#' z3a=.4
#' z3b=1.6
#' ####Hierarchical lam params
#' ###Sigma^2 lambda_g hyperparameters
#' a1=.7
#' b1=.7
#' a2=a1
#' b2=b1
#' a3=a1
#' b3=b1
#' ##Spacing dependence c in [0,1]
#' clam1=1
#' clam2=1
#' clam3=1
#' #####NumSplit
#' alpha1=3
#' alpha2=3
#' alpha3=3
#' J1max=10
#' J2max=10
#' J3max=10
#' ####Split Point Starting Value ###
#' J1=3
#' J2=3
#' J3=3
#' ###epsilon starting values/hyperparameters###
#' w=.7
#' psi1=.7
#' cep=2.4
#' #############
#' epstart=1.5
#' cl1=.25
#' cl2=.25
#' cl3=.25
#' ###Beta Starting Values
#' hyper1=c(psi,c,z1a,z1b,z2a,z2b,z3a,z3b,a1,b1,a2,b2,a3,b3,clam1,clam2,clam3)
#' hyper2=c(alpha1,alpha2,alpha3,J1max,J2max,J3max,J1,J2,J3,w,psi1,cep,epstart,cl1,cl2,cl3)
#' hyper=c(hyper1,hyper2)
#' ###Number of iterations and output location
#' BSVSS=10
#' BDIC=4
#'Path=tempdir()
#' ###Number of variables to exclude from selection and burnin percent
#'inc=2
#' ReturnModel(Y1,I1,Y2,I2,X,hyper,inc,c,BSVSS,BDIC,Path)
#' @export
ReturnModel=function(Y1,I1,Y2,I2,X,hyperparameters,inc,c,BSVSS,BDIC,Path){
cat("Running SVSS MCMC
")
B=BSVSS
burn=.6
beta1start=rep(1,ncol(X))
beta2start=beta1start
beta3start=beta1start
z1=SCRSELECTRETURN(Y1,I1,Y2,I2,X,hyperparameters,beta1start,beta2start,beta3start,B,inc,burn)
cat("
Chain 1 Finished
")
beta1start=c(rep(0,ncol(X)-inc),rep(-1,inc))
beta2start=beta1start
beta3start=beta1start
z2=SCRSELECTRETURN(Y1,I1,Y2,I2,X,hyperparameters,beta1start,beta2start,beta3start,B,inc,burn)
cat("
Chain 2 Finished
")
PCT1=colMeans(rbind(z1[[1]],z2[[1]]))
PCT2=colMeans(rbind(z1[[2]],z2[[2]]))
PCT3=colMeans(rbind(z1[[3]],z2[[3]]))
gam=colMeans(rbind(z1[[4]],z2[[4]]))
J1=c(z1[[5]],z2[[5]])
J2=c(z1[[6]],z2[[6]])
J3=c(z1[[7]],z2[[7]])
Mode <- function(x) {
ux <- unique(x)
ux[which.max(tabulate(match(x, ux)))]
}
J1=Mode(J1)
J2=Mode(J2)
J3=Mode(J3)
s1=rbind(z1[[8]],z2[[8]])
s2=rbind(z1[[9]],z2[[9]])
s3=rbind(z1[[10]],z2[[10]])
lam1=rbind(z1[[11]],z2[[11]])
lam2=rbind(z1[[12]],z2[[12]])
lam3=rbind(z1[[13]],z2[[13]])
s=matrix(rep(NA,length(s1)),nrow=nrow(s1))
l=matrix(rep(NA,length(lam1)),nrow=nrow(lam1))
for(b in 1:nrow(s)){
if(is.na(s1[b,(J1+2)])==FALSE && is.na(s1[b,(J1+3)])==TRUE){
s[b,1:(J1+2)]=s1[b,1:(J1+2)]
l[b,1:(J1+1)]=lam1[b,1:(J1+1)]
}
}
s1=colMeans(s,na.rm=TRUE)[1:(J1+2)]
lam1=colMeans(l,na.rm=TRUE)[1:(J1+1)]
s=matrix(rep(NA,length(s2)),nrow=nrow(s2))
l=matrix(rep(NA,length(lam2)),nrow=nrow(lam2))
for(b in 1:nrow(s)){
if(is.na(s2[b,(J2+2)])==FALSE && is.na(s2[b,(J2+3)])==TRUE){
s[b,1:(J2+2)]=s2[b,1:(J2+2)]
l[b,1:(J2+1)]=lam2[b,1:(J2+1)]
}
}
s2=colMeans(s,na.rm=TRUE)[1:(J2+2)]
lam2=colMeans(l,na.rm=TRUE)[1:(J2+1)]
s=matrix(rep(NA,length(s3)),nrow=nrow(s3))
l=matrix(rep(NA,length(lam3)),nrow=nrow(lam3))
for(b in 1:nrow(s)){
if(is.na(s3[b,(J3+2)])==FALSE && is.na(s3[b,(J3+3)])==TRUE){
s[b,1:(J3+2)]=s3[b,1:(J3+2)]
l[b,1:(J3+1)]=lam3[b,1:(J3+1)]
}
}
s3=colMeans(s,na.rm=TRUE)[1:(J3+2)]
lam3=colMeans(l,na.rm=TRUE)[1:(J3+1)]
Z=list(PCT1,PCT2,PCT3,gam,s1,s2,s3,lam1,lam2,lam3)
cat("
Marginal Posterior Probabilities of Inclusion
")
PCT1=Z[[1]]
PCT2=Z[[2]]
PCT3=Z[[3]]
cat("
Hazard 1: Non-terminal Event
")
print(PCT1)
cat("
Hazard 2: Death Before Non-terminal Event
")
print(PCT2)
cat("
Hazard 3: Death After Non-terminal Event
")
print(PCT3)
gam=Z[[4]]
s1=Z[[5]]
s2=Z[[6]]
s3=Z[[7]]
lam1=Z[[8]]
lam2=Z[[9]]
lam3=Z[[10]]
B=BDIC
cat("Baseline Hazard Quantities used for DIC-TAUG
")
cat("Haz1: S, lam")
print(s1)
print(lam1)
cat("Haz2: S, lam")
print(s2)
print(lam2)
cat("Haz3: S, lam")
print(s3)
print(lam3)
cat("
Starting DIC-tau_g Procedure
")
X1=DICTAUG(PCT1,PCT2,PCT3,X,Y1,Y2,I1,I2,s1,lam1,s2,lam2,s3,lam3,gam,c,B,inc)
Inc=inc
ext=function(X){
m1=rep(NA,length(X))
for(b in 1:length(X)){
if(is.null(dim(X[[b]]))==FALSE){
m1[b]=suppressWarnings(min(as.numeric(X[[b]]),na.rm=TRUE))
}
}
return(m1)
}
Vec=ext(X1)
VEC1=sort(unique(Vec))
VEC1=VEC1[!is.na(VEC1)]
## How many different TAU1 candidates are there?
##Min DIC
G1=min(Vec,na.rm=TRUE)
t1=1
for(b in 2:length(VEC1)){
if(VEC1[b]<(G1+1)){
t1=t1+1
}
}
tau1=rep(NA,t1)
tau12=rep(NA,t1)
tau13=rep(NA,t1)
for(k in 1:t1){
for(b in 1:length(Vec)){
if(is.na(Vec[b])==FALSE){
if(VEC1[k]==Vec[b]){
tau1[k]=b
}
}
}
}
minDIC=rep(NA,t1)
for( k in 1:t1){
G2=as.numeric(min(X1[[tau1[k]]],na.rm=TRUE))
NumCand=sum(X1[[tau1[k]]]<(G2+1))
H=suppressWarnings(as.numeric(sort(unique(X1[[tau1[k]]]))))
H=H[!is.na(H)]
H=H[1:NumCand]
j1=1
tau2=rep(NA,NumCand)
tau3=rep(NA,NumCand)
for(j in 1:NumCand){
for(l in 1:18){
for(m in 1:18){
if(X1[[tau1[k]]][l,m]==H[j]){
tau2[j]=l/20
tau3[j]=m/20
}
}
}
}
if(NumCand==1){
tau12[k]=tau2[1]
tau13[k]=tau3[1]
minDIC[k]=min(H[1])
}else{
sumINC=rep(NA,NumCand)
for(j in 1:NumCand){
sumINC[j]=sum(PCT2>tau2[j])+sum(PCT3>tau3[j])
}
minINC = min(sumINC)
H1=sumINC==minINC
spot1=rep(NA,length(sumINC))
for(j in 1:NumCand){
if(H1[j]==TRUE){
spot1[j]=j
}
}
spot1=spot1[!is.na(spot1)]
minDIC[k]=min(H[spot1])
for(j in 1:length(spot1)){
if(H[spot1[j]]==minDIC[k]){
spot2=spot1[j]
}
}
tau12[k]=tau2[spot2]
tau13[k]=tau3[spot2]
}
}
taumark=tau1
tau1=tau1/20
sumINC=rep(NA,length(tau1))
for(j in 1:t1){
sumINC[j]=sum(PCT1>tau1[j])+sum(PCT2>tau12[j])+sum(PCT3>tau13[j])
}
minINC = min(sumINC)
H1=sumINC==minINC
spot1=rep(NA,length(sumINC))
for(j in 1:t1){
if(H1[j]==TRUE){
spot1[j]=j
}
}
spot1=spot1[!is.na(spot1)]
tauG=c(tau1[spot1],tau12[spot1],tau13[spot1])
cat("Grid Search Finished
")
cat("DIC of best model: ",minDIC)
cat("
Tau_g value for optimum model
")
cat(tauG)
cat("
Variables Included
")
cat("
Hazard 1: Non-terminal Event
")
E1=PCT1>tauG[1]
print(E1)
cat("
Hazard 2: Death Before Non-terminal Event
")
E2=PCT2>tauG[2]
print(E2)
cat("
Hazard 3: Death After Non-terminal Event
")
E3=PCT3>tauG[3]
print(E3)
cat("
Performing MCMC for Final Model
")
COV=X
eta1=rep(0,length(E1))
eta2=rep(0,length(E2))
eta3=rep(0,length(E3))
for(i in 1:(ncol(COV)-inc)){
if(E1[i]==1){
eta1[i]=i
}
}
for(i in 1:(ncol(COV)-inc)){
if(E2[i]==1){
eta2[i]=i
}
}
for(i in 1:(ncol(COV)-inc)){
if(E3[i]==1){
eta3[i]=i
}
}
eta1=eta1[!(eta1==0)]
eta2=eta2[!(eta2==0)]
eta3=eta3[!(eta3==0)]
B=BSVSS
if(inc>0){
Include=(ncol(COV)-Inc+1):ncol(COV)
###Covariate Matrices
COV1=as.matrix(COV[,c(eta1,Include)])
COV2=as.matrix(COV[,c(eta2,Include)])
COV3=as.matrix(COV[,c(eta3,Include)])
}else{
COV1=as.matrix(COV[,eta1])
COV2=as.matrix(COV[,eta2])
COV3=as.matrix(COV[,eta3])
}
p1=ncol(COV1)
p2=ncol(COV2)
p3=ncol(COV3)
### Sets up Storage Matrices ##
##Hundred Thousand
###Beta/Eta###
beta1=matrix(rep(0,B*(p1)),nrow=B)
beta2=matrix(rep(0,B*(p2)),nrow=B)
beta3=matrix(rep(0,B*(p3)),nrow=B)
n=length(Y1)
Like=rep(0,B)
Indcond1=matrix(rep(0,p1*B),nrow=B)
Indcond2=matrix(rep(0,p2*B),nrow=B)
Indcond3=matrix(rep(0,p3*B),nrow=B)
Sigma1=c*solve(t(COV1)%*%COV1)
Sigma2=c*solve(t(COV2)%*%COV2)
Sigma3=c*solve(t(COV3)%*%COV3)
Indmix1=rep(0,B)
Indmix2=rep(0,B)
Indmix3=rep(0,B)
G1=length(s1)-1
G2=length(s2)-1
G3=length(s3)-1
LK1=function(Y1,Y2,I1,I2,Beta1){
LOGBH=0
et1=COV1%*%Beta1
LOGBH=LOGBH+sum(I1*et1)
for(k in 1:G1){
Del=pmax(0,pmin(Y1,s1[k+1])-s1[k])
LOGBH=LOGBH-sum(gam*Del*exp(lam1[k])*exp(et1))
}
return(LOGBH)
}
##
LK2=function(Y1,Y2,I1,I2,Beta2){
LOGBH=0
et1=COV2%*%Beta2
LOGBH=LOGBH+sum(I2*(1-I1)*et1)
Y=Y1
Y[I1==0]=Y2[I1==0]
for(k in 1:G2){
Del=pmax(0,pmin(Y,s2[k+1])-s2[k])
LOGBH=LOGBH-sum(gam*Del*exp(lam2[k])*exp(et1))
}
return(LOGBH)
}
###
##
LK3=function(Y1,Y2,I1,I2,Beta3){
LOGBH=0
et1=COV3%*%Beta3
LOGBH=LOGBH+sum(I2*(I1)*et1)
for(k in 1:G3){
Del=pmax(0,pmin(Y2[I1==1]-Y1[I1==1],s3[k+1])-s3[k])
LOGBH=LOGBH-sum(gam[I1==1]*Del*exp(lam3[k])*exp(et1[I1==1]))
}
return(LOGBH)
}
iter=0
for(b in 2:B){
iter="haz1"
##Print iteration
beta1[b,]=beta1[b-1,]
if(p1>1){
for(m in 1:p1){
V1 = Sigma1[m,m]
V2 = as.matrix(Sigma1[-m,-m])
V12 = as.matrix(Sigma1[m,-m])
thetab=beta1[b,]
thetano = as.matrix(thetab[-m])
meannew = t(V12)%*%solve(V2)%*%thetano
varnew = sqrt(V1 - t(V12)%*%solve(V2)%*%V12)
##################
beta1[b,m]=rnorm(1,meannew,varnew)
#beta1[b,m]=beta1[b-1,m] + runif(1,-clb,clb)
dn=log(dnorm(beta1[b,m],meannew,varnew))
###density old
do=log(dnorm(thetab[m],meannew,varnew))
Likeo=LK1(Y1,Y2,I1,I2,thetab)
Liken=LK1(Y1,Y2,I1,I2,beta1[b,])
alpha=Liken-Likeo+dn-do
U=log(runif(1,0,1))
if(U>alpha){
Indcond1[b,m]=0
beta1[b,]=thetab
}else{
Indcond1[b,m]=1
}
}
}else{
for(m in 1:p1){
thetab=beta1[b,]
meannew = 0
varnew = sqrt(Sigma1[m,m])
##################
beta1[b,m]=rnorm(1,meannew,varnew)
#beta1[b,m]=beta1[b-1,m] + runif(1,-clb,clb)
dn=log(dnorm(beta1[b,m],meannew,varnew))
###density old
do=log(dnorm(thetab[m],meannew,varnew))
Likeo=LK1(Y1,Y2,I1,I2,thetab)
Liken=LK1(Y1,Y2,I1,I2,beta1[b,])
alpha=Liken-Likeo+dn-do
U=log(runif(1,0,1))
if(U>alpha){
Indcond1[b,m]=0
beta1[b,]=thetab
}else{
Indcond1[b,m]=1
}
}
}
iter="haz2"
beta2[b,]=beta2[b-1,]
if(p2>1){
for(m in 1:p2){
V1 = Sigma2[m,m]
V2 = as.matrix(Sigma2[-m,-m])
V12 = as.matrix(Sigma2[m,-m])
thetab=beta2[b,]
thetano = as.matrix(thetab[-m])
meannew = t(V12)%*%solve(V2)%*%thetano
varnew = sqrt(V1 - t(V12)%*%solve(V2)%*%V12)
##################
beta2[b,m]=rnorm(1,meannew,varnew)
dn=log(dnorm(beta2[b,m],meannew,varnew))
###density old
do=log(dnorm(thetab[m],meannew,varnew))
Likeo=LK2(Y1,Y2,I1,I2,thetab)
Liken=LK2(Y1,Y2,I1,I2,beta2[b,])
alpha=Liken-Likeo+dn-do
U=log(runif(1,0,1))
if(U>alpha){
Indcond2[b,m]=0
beta2[b,]=thetab
}else{
Indcond2[b,m]=1
}
}
}else{
for(m in 1:p2){
thetab=beta2[b,]
meannew = 0
varnew = sqrt(Sigma2[m,m])
##################
beta2[b,m]=rnorm(1,meannew,varnew)
dn=log(dnorm(beta2[b,m],meannew,varnew))
###density old
do=log(dnorm(thetab[m],meannew,varnew))
Likeo=LK2(Y1,Y2,I1,I2,thetab)
Liken=LK2(Y1,Y2,I1,I2,beta2[b,])
alpha=Liken-Likeo+dn-do
U=log(runif(1,0,1))
if(U>alpha){
Indcond2[b,m]=0
beta2[b,]=thetab
}else{
Indcond2[b,m]=1
}
}
}
iter="haz3"
##Print iteration
beta3[b,]=beta3[b-1,]
if(p3>1){
for(m in 1:p3){
V1 = Sigma3[m,m]
V2 = as.matrix(Sigma3[-m,-m])
V12 = as.matrix(Sigma3[m,-m])
thetab=beta3[b,]
thetano = as.matrix(thetab[-m])
meannew = t(V12)%*%solve(V2)%*%thetano
varnew = sqrt(V1 - t(V12)%*%solve(V2)%*%V12)
##################
beta3[b,m]=rnorm(1,meannew,varnew)
dn=log(dnorm(beta3[b,m],meannew,varnew))
###density old
do=log(dnorm(thetab[m],meannew,varnew))
Likeo=LK3(Y1,Y2,I1,I2,thetab)
Liken=LK3(Y1,Y2,I1,I2,beta3[b,])
alpha=Liken-Likeo+dn-do
U=log(runif(1,0,1))
if(U>alpha){
Indcond3[b,m]=0
beta3[b,]=thetab
}else{
Indcond3[b,m]=1
}
}
}else{
for(m in 1:p3){
thetab=beta3[b,]
meannew = 0
varnew = sqrt(Sigma3[m,m])
##################
beta3[b,m]=rnorm(1,meannew,varnew)
dn=log(dnorm(beta3[b,m],meannew,varnew))
###density old
do=log(dnorm(thetab[m],meannew,varnew))
Likeo=LK3(Y1,Y2,I1,I2,thetab)
Liken=LK3(Y1,Y2,I1,I2,beta3[b,])
alpha=Liken-Likeo+dn-do
U=log(runif(1,0,1))
if(U>alpha){
Indcond3[b,m]=0
beta3[b,]=thetab
}else{
Indcond3[b,m]=1
}
}
}
Like[b]=LK3(Y1,Y2,I1,I2,beta3[b,])+LK2(Y1,Y2,I1,I2,beta2[b,])+LK1(Y1,Y2,I1,I2,beta1[b,])
###End
}
cat("MCMC is finished, returning posterior quantities
")
##burnin
beta1=beta1[(B*burn):B,]
beta2=beta2[(B*burn):B,]
beta3=beta3[(B*burn):B,]
##Re-parse out vector and return 0 for non-included variables.
BETA1=matrix(rep(0,nrow(beta1)*ncol(COV)),nrow=nrow(beta1))
BETA2=matrix(rep(0,nrow(beta1)*ncol(COV)),nrow=nrow(beta1))
BETA3=matrix(rep(0,nrow(beta1)*ncol(COV)),nrow=nrow(beta1))
if(inc>0){
BETA1[,c(eta1,Include)]=beta1
BETA2[,c(eta2,Include)]=beta2
BETA3[,c(eta3,Include)]=beta3
}else{
BETA1[,eta1]=beta1
BETA2[,eta2]=beta2
BETA3[,eta3]=beta3
}
Path1= paste0(Path,"/beta1.txt")
write.table(BETA1, Path1, sep="\t")
Path1= paste0(Path,"/beta2.txt")
write.table(BETA2, Path1, sep="\t")
Path1= paste0(Path,"/beta3.txt")
write.table(BETA3, Path1, sep="\t")
haz1=colMeans(beta1>0)
haz2=colMeans(beta2>0)
haz3=colMeans(beta3>0)
HAZ1=rep(NA,ncol(COV))
HAZ2=HAZ1
HAZ3=HAZ1
if(inc>0){
HAZ1[c(eta1,Include)]=haz1
HAZ2[c(eta2,Include)]=haz2
HAZ3[c(eta3,Include)]=haz3
}else{
HAZ1[eta1]=haz1
HAZ2[eta2]=haz2
HAZ3[eta3]=haz3
}
cat("
Posterior Probability of Increasing Hazard
")
cat("Hazard 1: Non-terminal Event
")
print(HAZ1)
cat("
Hazard 2: Death Before Non-terminal Event
")
print(HAZ2)
cat("
Hazard 3: Death After Non-terminal Event
")
print(HAZ3)
cat("
Posterior Median Hazard Ratio and Credible Interval
")
med1=exp(apply(beta1,2,median))
med2=exp(apply(beta2,2,median))
med3=exp(apply(beta3,2,median))
ql1=exp(apply(beta1,2,quantile,probs=.025))
ql2=exp(apply(beta2,2,quantile,probs=.025))
ql3=exp(apply(beta3,2,quantile,probs=.025))
qu1=exp(apply(beta1,2,quantile,probs=.975))
qu2=exp(apply(beta2,2,quantile,probs=.975))
qu3=exp(apply(beta3,2,quantile,probs=.975))
MED1=rep(NA,ncol(COV))
MED2=MED1
MED3=MED1
QL1=MED1
QL2=MED1
QL3=MED1
QU1=MED1
QU2=MED1
QU3=MED1
if(inc>0){
MED1[c(eta1,Include)]=med1
MED2[c(eta2,Include)]=med2
MED3[c(eta3,Include)]=med3
QL1[c(eta1,Include)]=ql1
QL2[c(eta2,Include)]=ql2
QL3[c(eta3,Include)]=ql3
QU1[c(eta1,Include)]=qu1
QU2[c(eta2,Include)]=qu2
QU3[c(eta3,Include)]=qu3
}else{
MED1[eta1]=med1
MED2[eta2]=med2
MED3[eta3]=med3
QL1[eta1]=ql1
QL2[eta2]=ql2
QL3[eta3]=ql3
QU1[eta1]=qu1
QU2[eta2]=qu2
QU3[eta3]=qu3
}
cat("
Hazard 1: Non-terminal Event
")
print(MED1)
print(QL1)
print(QU1)
cat("
Hazard 2: Death Before Non-terminal Event
")
print(MED2)
print(QL2)
print(QU2)
cat("
Hazard 3: Death After Non-terminal Event
")
print(MED3)
print(QL3)
print(QU3)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.