Nothing
####
## PRINT METHOD
####
##
print.Freq_HReg <- function(x, digits=3, alpha=0.05, ...)
{
conf.level = alpha
obj <- x
##
logEst <- obj$estimate
logSE <- sqrt(diag(obj$Finv))
value <- cbind(logEst, logSE, logEst - abs(qnorm(conf.level/2, 0, 1))*logSE, logEst + abs(qnorm(conf.level/2, 0, 1))*logSE)
##
dimnames(value) <- list(obj$myLabels, c( "Estimate", "SE", "LL", "UL"))
##
if(obj$class[2] == "Surv")
{
##
cat("\nAnalysis of independent univariate time-to-event data \n")
##
#cat("\nBaseline hazard function components:\n")
#cat(round(value[c(1:2),], digits=digits))
##
if(length(obj$myLabels) >= 3)
{
cat("Confidence level: ", conf.level, "\n", sep = "")
cat("\nRegression coefficients:\n")
cat(round(value[-c(1:2),], digits=digits))
}
}
##
if(obj$class[2] == "ID")
{
##
cat("\nAnalysis of independent semi-competing risks data \n")
cat(obj$class[5], "assumption for h3\n")
##
#cat("\nBaseline hazard function components:\n")
#cat(round(value[c(1:6),], digits=digits))
##
value_theta <- matrix(exp(value[7,]), ncol = 4)
dimnames(value_theta) <- list("", c( "Estimate", "SE", "LL", "UL"))
value_theta[1,2] <- value[7,2] * exp(value[7,1])
cat("Confidence level: ", conf.level, "\n", sep = "")
cat("\nVariance of frailties, theta:\n")
if(obj$frailty == TRUE) print(round(value_theta, digits=digits))
if(obj$frailty == FALSE) cat("NA")
##
if(sum(obj$nP) != 0)
{
cat("\nRegression coefficients:\n")
if(obj$frailty == TRUE) print(round(value[-c(1:7),], digits=digits))
if(obj$frailty == FALSE) print(round(value[-c(1:6),], digits=digits))
}
cat("\nNote: Covariates are arranged in order of transition number, 1->3.\n")
}
##
invisible()
}
print.Bayes_HReg <- function(x, digits=3, alpha=0.05, ...)
{
conf.level = alpha
nChain = x$setup$nChain
if(x$class[2] == "ID")
{
if(x$class[3] == "Cor")
{
##
cat("\nAnalysis of cluster-correlated semi-competing risks data \n")
}
if(x$class[3] == "Ind")
{
##
cat("\nAnalysis of independent semi-competing risks data \n")
}
##
cat(x$setup$model, "assumption for h3\n")
}
if(x$class[2] == "Surv")
{
if(x$class[3] == "Cor")
{
##
cat("\nAnalysis of cluster-correlated univariate time-to-event data \n")
}
if(x$class[3] == "Ind")
{
##
cat("\nAnalysis of independent univariate time-to-event data \n")
}
}
##
cat("\nNumber of chains: ", nChain,"\n")
##
cat("Number of scans: ", x$setup$numReps,"\n")
##
cat("Thinning: ", x$setup$thin,"\n")
##
cat("Percentage of burnin: ", x$setup$burninPerc*100, "%\n", sep = "")
# convergence diagnostics
if(nChain > 1){
cat("\n######\n")
cat("Potential Scale Reduction Factor\n")
if(x$class[2] == "ID")
{
theta <- x$chain1$theta.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
theta <- cbind(theta, x[[nam]]$theta.p)
}
psrftheta <- matrix(calcPSR(theta), 1, 1)
dimnames(psrftheta) <- list("", "")
cat("\nVariance of frailties, theta:")
print(round(psrftheta, digits=digits))
beta.names <- unique(c(x$chain1$covNames1, x$chain1$covNames2, x$chain1$covNames3))
nP <- length(beta.names)
output <- matrix(NA, nrow=nP, ncol=3)
dimnames(output) <- list(beta.names, c("beta1", "beta2", "beta3"))
if(length(x$chain1$beta1.p) != 0){
#beta1
p1 = dim(x$chain1$beta1.p)[2]
psrfBeta1 <- rep(NA, p1)
for(j in 1:p1){
#namPara = paste("beta_", j, sep = "")
beta1 <- x$chain1$beta1[,j]
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
beta1 <- cbind(beta1, x[[nam]]$beta1[,j])
}
psrfBeta1[j] <- calcPSR(beta1)
}
for(i in 1:nP)
{
for(k in 1:p1) if(x$chain1$covNames1[k] == beta.names[i]) output[i,1] <- psrfBeta1[k]
}
}
if(length(x$chain1$beta2.p) != 0){
#beta2
p2 = dim(x$chain1$beta2.p)[2]
psrfBeta2 <- rep(NA, p2)
for(j in 1:p2){
#namPara = paste("beta_", j, sep = "")
beta2 <- x$chain1$beta2[,j]
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
beta2 <- cbind(beta2, x[[nam]]$beta2[,j])
}
psrfBeta2[j] <- calcPSR(beta2)
}
for(i in 1:nP)
{
for(k in 1:p2) if(x$chain1$covNames2[k] == beta.names[i]) output[i,2] <- psrfBeta2[k]
}
}
if(length(x$chain1$beta3.p) != 0){
#beta3
p3 = dim(x$chain1$beta3.p)[2]
psrfBeta3 <- rep(NA, p3)
for(j in 1:p3){
#namPara = paste("beta_", j, sep = "")
beta3 <- x$chain1$beta3[,j]
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
beta3 <- cbind(beta3, x[[nam]]$beta3[,j])
}
psrfBeta3[j] <- calcPSR(beta3)
}
for(i in 1:nP)
{
for(k in 1:p3) if(x$chain1$covNames3[k] == beta.names[i]) output[i,3] <- psrfBeta3[k]
}
}
if(nP > 0)
{
cat("\nRegression coefficients:\n")
output.coef <- output
print(round(output.coef, digits=digits))
}
##
cat("\nBaseline hazard function components:\n")
if(x$class[4] == "WB")
{
##
# alpha
alpha <- x$chain1$alpha1.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
alpha <- cbind(alpha, x[[nam]]$alpha1.p)
}
psrfAlpha1 <- calcPSR(alpha)
alpha <- x$chain1$alpha2.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
alpha <- cbind(alpha, x[[nam]]$alpha2.p)
}
psrfAlpha2 <- calcPSR(alpha)
alpha <- x$chain1$alpha3.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
alpha <- cbind(alpha, x[[nam]]$alpha3.p)
}
psrfAlpha3 <- calcPSR(alpha)
# kappa
kappa <- x$chain1$kappa1.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
kappa <- cbind(kappa, x[[nam]]$kappa1.p)
}
psrfKappa1 <- calcPSR(kappa)
kappa <- x$chain1$kappa2.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
kappa <- cbind(kappa, x[[nam]]$kappa2.p)
}
psrfKappa2 <- calcPSR(kappa)
kappa <- x$chain1$kappa3.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
kappa <- cbind(kappa, x[[nam]]$kappa3.p)
}
psrfKappa3 <- calcPSR(kappa)
bh_WB <- matrix(c(psrfKappa1, psrfKappa2, psrfKappa3, psrfAlpha1, psrfAlpha2, psrfAlpha3), 2, 3, byrow = T)
dimnames(bh_WB) <- list(c("kappa", "alpha"), c("h1", "h2", "h3"))
print(round(bh_WB, digits=digits))
}
if(x$class[4] == "PEM")
{
##
ntime1 = length(x$chain1$time_lambda1)
ntime2 = length(x$chain1$time_lambda2)
ntime3 = length(x$chain1$time_lambda3)
# lambda's
psrfLam <- rep(NA, ntime1)
for(j in 1:ntime1){
lambda1 <- x$chain1$lambda1.fin[,j]
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
lambda1 <- cbind(lambda1, x[[nam]]$lambda1.fin[,j])
}
psrfLam[j] <- calcPSR(lambda1)
}
cat("\nlambda1: summary statistics", "\n")
print(round(summary(psrfLam), digits=digits))
psrfLam <- rep(NA, ntime2)
for(j in 1:ntime2){
lambda2 <- x$chain1$lambda2.fin[,j]
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
lambda2 <- cbind(lambda2, x[[nam]]$lambda2.fin[,j])
}
psrfLam[j] <- calcPSR(lambda2)
}
cat("\nlambda2: summary statistics", "\n")
print(round(summary(psrfLam), digits=digits))
psrfLam <- rep(NA, ntime3)
for(j in 1:ntime3){
lambda3 <- x$chain1$lambda3.fin[,j]
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
lambda3 <- cbind(lambda3, x[[nam]]$lambda3.fin[,j])
}
psrfLam[j] <- calcPSR(lambda3)
}
cat("\nlambda3: summary statistics", "\n")
print(round(summary(psrfLam), digits=digits))
# mu_lam
mu <- x$chain1$mu_lam1.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
mu <- cbind(mu, x[[nam]]$mu_lam1.p)
}
psrfMu1 <- calcPSR(mu)
mu <- x$chain1$mu_lam2.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
mu <- cbind(mu, x[[nam]]$mu_lam2.p)
}
psrfMu2 <- calcPSR(mu)
mu <- x$chain1$mu_lam3.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
mu <- cbind(mu, x[[nam]]$mu_lam3.p)
}
psrfMu3 <- calcPSR(mu)
# sigSq_lam
sig <- x$chain1$sigSq_lam1.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
sig <- cbind(sig, x[[nam]]$sigSq_lam1.p)
}
psrfSig1 <- calcPSR(sig)
sig <- x$chain1$sigSq_lam2.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
sig <- cbind(sig, x[[nam]]$sigSq_lam2.p)
}
psrfSig2 <- calcPSR(sig)
sig <- x$chain1$sigSq_lam3.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
sig <- cbind(sig, x[[nam]]$sigSq_lam3.p)
}
psrfSig3 <- calcPSR(sig)
# J
J <- x$chain1$K1.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
J <- cbind(J, x[[nam]]$K1.p)
}
psrfJ1 <- calcPSR(J)
J <- x$chain1$K2.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
J <- cbind(J, x[[nam]]$K2.p)
}
psrfJ2 <- calcPSR(J)
J <- x$chain1$K3.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
J <- cbind(J, x[[nam]]$K3.p)
}
psrfJ3 <- calcPSR(J)
bh_PEM <- matrix(c(psrfMu1, psrfMu2, psrfMu3, psrfSig1, psrfSig2, psrfSig3, psrfJ1, psrfJ2, psrfJ3), 3, 3, byrow = T)
dimnames(bh_PEM) <- list(c("mu", "sigmaSq", "K"), c("h1", "h2", "h3"))
cat("\n")
print(round(bh_PEM, digits=digits))
}
}
if(x$class[2] == "Surv")
{
beta.names <- c(x$chain1$covNames)
nP <- length(beta.names)
output <- matrix(NA, nrow=nP, ncol=1)
dimnames(output) <- list(beta.names, c("beta"))
if(length(x$chain1$beta.p) != 0){
#beta
p = dim(x$chain1$beta.p)[2]
psrfBeta <- rep(NA, p)
for(j in 1:p){
beta <- x$chain1$beta[,j]
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
beta <- cbind(beta, x[[nam]]$beta[,j])
}
psrfBeta[j] <- calcPSR(beta)
}
for(i in 1:nP)
{
for(k in 1:p) if(x$chain1$covNames[k] == beta.names[i]) output[i,1] <- psrfBeta[k]
}
}
if(nP > 0)
{
cat("\nRegression coefficients:\n")
output.coef <- output
print(round(output.coef, digits=digits))
}
if(x$class[4] == "WB")
{
##
# alpha
alpha <- x$chain1$alpha.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
alpha <- cbind(alpha, x[[nam]]$alpha.p)
}
psrfAlpha <- calcPSR(alpha)
# kappa
kappa <- x$chain1$kappa.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
kappa <- cbind(kappa, x[[nam]]$kappa.p)
}
psrfKappa <- calcPSR(kappa)
bh_WB <- matrix(c(psrfKappa, psrfAlpha), 2, 1, byrow = T)
dimnames(bh_WB) <- list(c("kappa", "alpha"), c("h"))
print(round(bh_WB, digits=digits))
}
if(x$class[4] == "PEM")
{
##
ntime = length(x$chain1$time_lambda)
# lambda
psrfLam <- rep(NA, ntime)
for(j in 1:ntime){
namPara = paste("beta_", j, sep = "")
lambda <- x$chain1$lambda.fin[,j]
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
lambda <- cbind(lambda, x[[nam]]$lambda.fin[,j])
}
psrfLam[j] <- calcPSR(lambda)
}
cat("\n lambda: summary statistics", "\n")
print(round(summary(psrfLam), 2))
# mu_lam
mu <- x$chain1$mu_lam.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
mu <- cbind(mu, x[[nam]]$mu_lam.p)
}
psrfMu <- calcPSR(mu)
# sigSq_lam
sig <- x$chain1$sigSq_lam.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
sig <- cbind(sig, x[[nam]]$sigSq_lam.p)
}
psrfSig <- calcPSR(sig)
# J
J <- x$chain1$K.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
J <- cbind(J, x[[nam]]$K.p)
}
psrfJ <- calcPSR(J)
bh_PEM <- matrix(c(psrfMu, psrfSig, psrfJ), 3, 1, byrow = T)
dimnames(bh_PEM) <- list(c("mu", "sigmaSq", "K"), c("h"))
cat("\n")
print(round(bh_PEM, digits=digits))
}
}
}
else if(nChain == 1)
{
cat("Potential scale reduction factor cannot be calculated. \n")
cat("The number of chains must be larger than 1. \n")
}
cat("\n######\n")
cat("Estimates\n")
cat("Credibility level: ", conf.level, "\n", sep = "")
if(x$class[2] == "ID")
{
##
cat("\nVariance of frailties, theta:\n")
theta.p <- x$chain1$theta.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
theta.p <- rbind(theta.p, x[[nam]]$theta.p)
}
}
theta.pMed <- apply(theta.p, 2, median)
theta.pSd <- apply(theta.p, 2, sd)
theta.pUb <- apply(theta.p, 2, quantile, prob = (1-conf.level/2))
theta.pLb <- apply(theta.p, 2, quantile, prob = conf.level/2)
tbl <- matrix(NA, 1, 4)
dimnames(tbl) <- list("", c( "Estimate", "SD", "LL", "UL"))
tbl[,1] <- theta.pMed
tbl[,2] <- theta.pSd
tbl[,3] <- theta.pLb
tbl[,4] <- theta.pUb
print(round(tbl, digits=digits))
##
tbl_beta <- NULL
if(length(x$chain1$beta1.p) != 0){
p1 = dim(x$chain1$beta1.p)[2]
beta.p <- x$chain1$beta1.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
beta.p <- rbind(beta.p, x[[nam]]$beta1.p)
}
}
beta.pMed <- apply(beta.p, 2, median)
beta.pSd <- apply(beta.p, 2, sd)
beta.pUb <- apply(beta.p, 2, quantile, prob = (1-conf.level/2))
beta.pLb <- apply(beta.p, 2, quantile, prob = conf.level/2)
tbl1 <- matrix(NA, p1, 4)
rownames(tbl1) <- x$chain1$covNames1
tbl1[,1] <- beta.pMed
tbl1[,2] <- beta.pSd
tbl1[,3] <- exp(beta.pLb)
tbl1[,4] <- exp(beta.pUb)
tbl_beta <- tbl1
}
if(length(x$chain1$beta2.p) != 0){
p2 = dim(x$chain1$beta2.p)[2]
beta.p <- x$chain1$beta2.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
beta.p <- rbind(beta.p, x[[nam]]$beta2.p)
}
}
beta.pMed <- apply(beta.p, 2, median)
beta.pSd <- apply(beta.p, 2, sd)
beta.pUb <- apply(beta.p, 2, quantile, prob = (1-conf.level/2))
beta.pLb <- apply(beta.p, 2, quantile, prob = conf.level/2)
tbl2 <- matrix(NA, p2, 4)
rownames(tbl2) <- x$chain1$covNames2
tbl2[,1] <- beta.pMed
tbl2[,2] <- beta.pSd
tbl2[,3] <- exp(beta.pLb)
tbl2[,4] <- exp(beta.pUb)
tbl_beta <- rbind(tbl_beta, tbl2)
}
if(length(x$chain1$beta3.p) != 0){
p3 = dim(x$chain1$beta3.p)[2]
beta.p <- x$chain1$beta3.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
beta.p <- rbind(beta.p, x[[nam]]$beta3.p)
}
}
beta.pMed <- apply(beta.p, 2, median)
beta.pSd <- apply(beta.p, 2, sd)
beta.pUb <- apply(beta.p, 2, quantile, prob = (1-conf.level/2))
beta.pLb <- apply(beta.p, 2, quantile, prob = conf.level/2)
tbl3 <- matrix(NA, p3, 4)
rownames(tbl3) <- x$chain1$covNames3
tbl3[,1] <- beta.pMed
tbl3[,2] <- beta.pSd
tbl3[,3] <- exp(beta.pLb)
tbl3[,4] <- exp(beta.pUb)
tbl_beta <- rbind(tbl_beta, tbl3)
}
}
if(x$class[2] == "Surv")
{
tbl_beta <- NULL
if(length(x$chain1$beta.p) != 0){
p = dim(x$chain1$beta.p)[2]
beta.p <- x$chain1$beta.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
beta.p <- rbind(beta.p, x[[nam]]$beta.p)
}
}
beta.pMed <- apply(beta.p, 2, median)
beta.pSd <- apply(beta.p, 2, sd)
beta.pUb <- apply(beta.p, 2, quantile, prob = (1-conf.level/2))
beta.pLb <- apply(beta.p, 2, quantile, prob = conf.level/2)
tbl_beta <- matrix(NA, p, 4)
rownames(tbl_beta) <- x$chain1$covNames
tbl_beta[,1] <- beta.pMed
tbl_beta[,2] <- beta.pSd
tbl_beta[,3] <- exp(beta.pLb)
tbl_beta[,4] <- exp(beta.pUb)
}
}
if(!is.null(tbl_beta))
{
cat("\nRegression coefficients:\n")
colnames(tbl_beta) <- c( "Estimate", "SD", "LL", "UL")
print(round(tbl_beta, digits=digits))
if(x$class[2] == "ID")
{
cat("\nNote: Covariates are arranged in order of transition number, 1->3.\n")
}
}
invisible()
}
print.Bayes_AFT <- function(x, digits=3, alpha=0.05, ...)
{
conf.level = alpha
nChain = x$setup$nChain
if(x$class[2] == "ID")
{
if(x$class[3] == "Cor")
{
##
cat("\nAnalysis of cluster-correlated semi-competing risks data \n")
}
if(x$class[3] == "Ind")
{
##
cat("\nAnalysis of independent semi-competing risks data \n")
}
}
if(x$class[2] == "Surv")
{
if(x$class[3] == "Cor")
{
##
cat("\nAnalysis of cluster-correlated univariate time-to-event data \n")
}
if(x$class[3] == "Ind")
{
##
cat("\nAnalysis of independent univariate time-to-event data \n")
}
}
##
cat("\nNumber of chains: ", nChain,"\n")
##
cat("Number of scans: ", x$setup$numReps,"\n")
##
cat("Thinning: ", x$setup$thin,"\n")
##
cat("Percentage of burnin: ", x$setup$burninPerc*100, "%\n", sep = "")
# convergence diagnostics
if(nChain > 1){
cat("\n######\n")
cat("Potential Scale Reduction Factor\n")
if(x$class[2] == "ID")
{
theta <- x$chain1$theta.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
theta <- cbind(theta, x[[nam]]$theta.p)
}
psrftheta <- matrix(calcPSR(theta), 1, 1)
dimnames(psrftheta) <- list("", "")
cat("\nVariance of frailties, theta:")
print(round(psrftheta, digits=digits))
beta.names <- unique(c(x$chain1$covNames1, x$chain1$covNames2, x$chain1$covNames3))
nP <- length(beta.names)
output <- matrix(NA, nrow=nP, ncol=3)
dimnames(output) <- list(beta.names, c("beta1", "beta2", "beta3"))
if(length(x$chain1$beta1.p) != 0){
#beta1
p1 = dim(x$chain1$beta1.p)[2]
psrfBeta1 <- rep(NA, p1)
for(j in 1:p1){
#namPara = paste("beta_", j, sep = "")
beta1 <- x$chain1$beta1[,j]
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
beta1 <- cbind(beta1, x[[nam]]$beta1[,j])
}
psrfBeta1[j] <- calcPSR(beta1)
}
for(i in 1:nP)
{
for(k in 1:p1) if(x$chain1$covNames1[k] == beta.names[i]) output[i,1] <- psrfBeta1[k]
}
}
if(length(x$chain1$beta2.p) != 0){
#beta2
p2 = dim(x$chain1$beta2.p)[2]
psrfBeta2 <- rep(NA, p2)
for(j in 1:p2){
#namPara = paste("beta_", j, sep = "")
beta2 <- x$chain1$beta2[,j]
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
beta2 <- cbind(beta2, x[[nam]]$beta2[,j])
}
psrfBeta2[j] <- calcPSR(beta2)
}
for(i in 1:nP)
{
for(k in 1:p2) if(x$chain1$covNames2[k] == beta.names[i]) output[i,2] <- psrfBeta2[k]
}
}
if(length(x$chain1$beta3.p) != 0){
#beta3
p3 = dim(x$chain1$beta3.p)[2]
psrfBeta3 <- rep(NA, p3)
for(j in 1:p3){
#namPara = paste("beta_", j, sep = "")
beta3 <- x$chain1$beta3[,j]
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
beta3 <- cbind(beta3, x[[nam]]$beta3[,j])
}
psrfBeta3[j] <- calcPSR(beta3)
}
for(i in 1:nP)
{
for(k in 1:p3) if(x$chain1$covNames3[k] == beta.names[i]) output[i,3] <- psrfBeta3[k]
}
}
if(nP > 0)
{
cat("\nRegression coefficients:\n")
output.coef <- output
print(round(output.coef, digits=digits))
}
##
cat("\nBaseline survival function components:\n")
if(x$class[4] == "LN")
{
##
# mu
mu <- x$chain1$mu1.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
mu <- cbind(mu, x[[nam]]$mu1.p)
}
psrfmu1 <- calcPSR(mu)
mu <- x$chain1$mu2.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
mu <- cbind(mu, x[[nam]]$mu2.p)
}
psrfmu2 <- calcPSR(mu)
mu <- x$chain1$mu3.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
mu <- cbind(mu, x[[nam]]$mu3.p)
}
psrfmu3 <- calcPSR(mu)
# sigSq
sigSq <- x$chain1$sigSq1.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
sigSq <- cbind(sigSq, x[[nam]]$sigSq1.p)
}
psrfSigSq1 <- calcPSR(sigSq)
sigSq <- x$chain1$sigSq2.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
sigSq <- cbind(sigSq, x[[nam]]$sigSq2.p)
}
psrfSigSq2 <- calcPSR(sigSq)
sigSq <- x$chain1$sigSq3.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
sigSq <- cbind(sigSq, x[[nam]]$sigSq3.p)
}
psrfSigSq3 <- calcPSR(sigSq)
bh_LN <- matrix(c(psrfmu1, psrfmu2, psrfmu3, psrfSigSq1, psrfSigSq2, psrfSigSq3), 2, 3, byrow = T)
dimnames(bh_LN) <- list(c("mu", "sigmaSq"), c("g=1", "g=2", "g=3"))
print(round(bh_LN, digits=digits))
}
if(x$class[4] == "DPM")
{
##
# tau
tau <- x$chain1$tau1.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
tau <- cbind(tau, x[[nam]]$tau1.p)
}
psrfTau1 <- calcPSR(tau)
tau <- x$chain1$tau2.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
tau <- cbind(tau, x[[nam]]$tau2.p)
}
psrfTau2 <- calcPSR(tau)
tau <- x$chain1$tau3.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
tau <- cbind(tau, x[[nam]]$tau3.p)
}
psrfTau3 <- calcPSR(tau)
bh_DPM <- matrix(c(psrfTau1, psrfTau2, psrfTau3), 1, 3, byrow = T)
dimnames(bh_DPM) <- list(c("tau"), c("g=1", "g=2", "g=3"))
cat("\n")
print(round(bh_DPM, digits=digits))
}
}
if(x$class[2] == "Surv")
{
beta.names <- c(x$chain1$covNames)
nP <- length(beta.names)
output <- matrix(NA, nrow=nP, ncol=1)
dimnames(output) <- list(beta.names, c("beta"))
if(length(x$chain1$beta.p) != 0){
#beta
p = dim(x$chain1$beta.p)[2]
psrfBeta <- rep(NA, p)
for(j in 1:p){
beta <- x$chain1$beta.p[,j]
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
beta <- cbind(beta, x[[nam]]$beta.p[,j])
}
psrfBeta[j] <- calcPSR(beta)
}
for(i in 1:nP)
{
for(k in 1:p) if(x$chain1$covNames[k] == beta.names[i]) output[i,1] <- psrfBeta[k]
}
}
if(nP > 0)
{
cat("\nRegression coefficients:\n")
output.coef <- output
print(round(output.coef, digits=digits))
}
if(x$class[4] == "LN")
{
##
# mu
mu <- x$chain1$mu.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
mu <- cbind(mu, x[[nam]]$mu.p)
}
psrfmu <- calcPSR(mu)
# sigSq
sigSq <- x$chain1$sigSq.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
sigSq <- cbind(sigSq, x[[nam]]$sigSq.p)
}
psrfSigSq <- calcPSR(sigSq)
bh_LN <- matrix(c(psrfmu, psrfSigSq), 2, 1, byrow = T)
dimnames(bh_LN) <- list(c("mu", "sigmaSq"), c("h"))
print(round(bh_LN, digits=digits))
}
if(x$class[4] == "DPM")
{
##
# tau
tau <- x$chain1$tau.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
tau <- cbind(tau, x[[nam]]$tau.p)
}
psrfTau <- calcPSR(tau)
bh_DPM <- matrix(c(psrfTau), 1, 1, byrow = T)
dimnames(bh_DPM) <- list(c("tau"), c("h"))
cat("\n")
print(round(bh_DPM, digits=digits))
}
}
}
else if(nChain == 1)
{
cat("Potential scale reduction factor cannot be calculated. \n")
cat("The number of chains must be larger than 1. \n")
}
cat("\n######\n")
cat("Estimates\n")
cat("Credibility level: ", conf.level, "\n", sep = "")
if(x$class[2] == "ID")
{
##
cat("\nVariance of frailties, theta:\n")
theta.p <- x$chain1$theta.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
theta.p <- rbind(theta.p, x[[nam]]$theta.p)
}
}
theta.pMed <- apply(theta.p, 2, median)
theta.pSd <- apply(theta.p, 2, sd)
theta.pUb <- apply(theta.p, 2, quantile, prob = (1-conf.level/2))
theta.pLb <- apply(theta.p, 2, quantile, prob = conf.level/2)
tbl <- matrix(NA, 1, 4)
dimnames(tbl) <- list("", c( "Estimate", "SD", "LL", "UL"))
tbl[,1] <- theta.pMed
tbl[,2] <- theta.pSd
tbl[,3] <- theta.pLb
tbl[,4] <- theta.pUb
print(round(tbl, digits=digits))
##
tbl_beta <- NULL
if(length(x$chain1$beta1.p) != 0){
p1 = dim(x$chain1$beta1.p)[2]
beta.p <- x$chain1$beta1.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
beta.p <- rbind(beta.p, x[[nam]]$beta1.p)
}
}
beta.pMed <- apply(beta.p, 2, median)
beta.pSd <- apply(beta.p, 2, sd)
beta.pUb <- apply(beta.p, 2, quantile, prob = (1-conf.level/2))
beta.pLb <- apply(beta.p, 2, quantile, prob = conf.level/2)
tbl1 <- matrix(NA, p1, 4)
rownames(tbl1) <- x$chain1$covNames1
tbl1[,1] <- beta.pMed
tbl1[,2] <- beta.pSd
tbl1[,3] <- exp(beta.pLb)
tbl1[,4] <- exp(beta.pUb)
tbl_beta <- tbl1
}
if(length(x$chain1$beta2.p) != 0){
p2 = dim(x$chain1$beta2.p)[2]
beta.p <- x$chain1$beta2.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
beta.p <- rbind(beta.p, x[[nam]]$beta2.p)
}
}
beta.pMed <- apply(beta.p, 2, median)
beta.pSd <- apply(beta.p, 2, sd)
beta.pUb <- apply(beta.p, 2, quantile, prob = (1-conf.level/2))
beta.pLb <- apply(beta.p, 2, quantile, prob = conf.level/2)
tbl2 <- matrix(NA, p2, 4)
rownames(tbl2) <- x$chain1$covNames2
tbl2[,1] <- beta.pMed
tbl2[,2] <- beta.pSd
tbl2[,3] <- exp(beta.pLb)
tbl2[,4] <- exp(beta.pUb)
tbl_beta <- rbind(tbl_beta, tbl2)
}
if(length(x$chain1$beta3.p) != 0){
p3 = dim(x$chain1$beta3.p)[2]
beta.p <- x$chain1$beta3.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
beta.p <- rbind(beta.p, x[[nam]]$beta3.p)
}
}
beta.pMed <- apply(beta.p, 2, median)
beta.pSd <- apply(beta.p, 2, sd)
beta.pUb <- apply(beta.p, 2, quantile, prob = (1-conf.level/2))
beta.pLb <- apply(beta.p, 2, quantile, prob = conf.level/2)
tbl3 <- matrix(NA, p3, 4)
rownames(tbl3) <- x$chain1$covNames3
tbl3[,1] <- beta.pMed
tbl3[,2] <- beta.pSd
tbl3[,3] <- exp(beta.pLb)
tbl3[,4] <- exp(beta.pUb)
tbl_beta <- rbind(tbl_beta, tbl3)
}
}
if(x$class[2] == "Surv")
{
tbl_beta <- NULL
if(length(x$chain1$beta.p) != 0){
p = dim(x$chain1$beta.p)[2]
beta.p <- x$chain1$beta.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
beta.p <- rbind(beta.p, x[[nam]]$beta.p)
}
}
beta.pMed <- apply(beta.p, 2, median)
beta.pSd <- apply(beta.p, 2, sd)
beta.pUb <- apply(beta.p, 2, quantile, prob = (1-conf.level/2))
beta.pLb <- apply(beta.p, 2, quantile, prob = conf.level/2)
tbl_beta <- matrix(NA, p, 4)
rownames(tbl_beta) <- x$chain1$covNames
tbl_beta[,1] <- beta.pMed
tbl_beta[,2] <- beta.pSd
tbl_beta[,3] <- exp(beta.pLb)
tbl_beta[,4] <- exp(beta.pUb)
}
}
if(!is.null(tbl_beta))
{
cat("\nRegression coefficients:\n")
colnames(tbl_beta) <- c( "Estimate", "SD", "LL", "UL")
print(round(tbl_beta, digits=digits))
if(x$class[2] == "ID")
{
cat("\nNote: Covariates are arranged in order of transition number, 1->3.\n")
}
}
invisible()
}
####
## SUMMARY METHOD
####
##
summary.Freq_HReg <- function(object, digits=3, alpha=0.05, ...)
{
conf.level = alpha
obj <- object
##
logEst <- obj$estimate
logSE <- sqrt(diag(obj$Finv))
results <- cbind(logEst, logEst - abs(qnorm(conf.level/2, 0, 1))*logSE, logEst + abs(qnorm(conf.level/2, 0, 1))*logSE)
##
if(obj$class[2] == "Surv")
{
nP <- length(obj$estimate)-2
##
#cat("\nRegression coefficients:\n")
if(nP > 0)
{
output.coef <- matrix(results[-c(1:2),], nrow=nP)
dimnames(output.coef) <- list(unique(obj$myLabels[-c(1:2)]), c("beta", "LL", "UL"))
}else
{
output.coef <- NULL
}
##
#cat("\nBaseline hazard function components:\n")
output.h0 <- results[c(1:2),]
dimnames(output.h0) <- list(c("Weibull: log-kappa", "Weibull: log-alpha"), c("h-PM", "LL", "UL"))
##
value <- list(coef=output.coef, h0=output.h0, code=obj$code, logLike=obj$logLike, nP=nrow(results), class=obj$class)
}
##
if(obj$class[2] == "ID")
{
##
nP.0 <- ifelse(obj$frailty, 7, 6)
nP.1 <- obj$nP[1]
nP.2 <- obj$nP[2]
nP.3 <- obj$nP[3]
##
beta.names <- unique(obj$myLabels[-c(1:nP.0)])
nP <- length(beta.names)
##
#cat("\nRegression coefficients:\n")
output <- matrix(NA, nrow=nP, ncol=9)
dimnames(output) <- list(beta.names, c("beta1", "LL", "UL", "beta2", "LL", "UL", "beta3", "LL", "UL"))
for(i in 1:nP)
{
if(nP.1 != 0)
{
for(j in 1:nP.1) if(obj$myLabels[nP.0+j] == beta.names[i]) output[i,1:3] <- results[nP.0+j,]
}
if(nP.2 != 0)
{
for(j in 1:nP.2) if(obj$myLabels[nP.0+nP.1+j] == beta.names[i]) output[i,4:6] <- results[nP.0+nP.1+j,]
}
if(nP.3 != 0)
{
for(j in 1:nP.3) if(obj$myLabels[nP.0+nP.1+nP.2+j] == beta.names[i]) output[i,7:9] <- results[nP.0+nP.1+nP.2+j,]
}
}
output.coef <- output
##
#cat("\nVariance of frailties:\n")
output <- matrix(NA, nrow=1, ncol=3)
dimnames(output) <- list(c("theta"), c("Estimate", "LL", "UL"))
if(obj$frailty == TRUE) output[1,] <- exp(results[7,])
if(obj$frailty == FALSE) output[1,] <- rep(NA, 3)
output.theta <- output
##
#cat("\nBaseline hazard function components:\n")
output <- matrix(NA, nrow=2, ncol=9)
dimnames(output) <- list(c("Weibull: log-kappa", "Weibull: log-alpha"), c("h1-PM", "LL", "UL", "h2-PM", "LL", "UL", "h3-PM", "LL", "UL"))
output[1,1:3] <- results[1,]
output[1,4:6] <- results[3,]
output[1,7:9] <- results[5,]
output[2,1:3] <- results[2,]
output[2,4:6] <- results[4,]
output[2,7:9] <- results[6,]
output.h0 <- output
##
value <- list(coef=output.coef, theta=output.theta, h0=output.h0, code=obj$code, logLike=obj$logLike, nP=nrow(results), class=obj$class, conf.level=conf.level)
}
class(value) <- "summ.Freq_HReg"
##
return(value)
}
summary.Bayes_HReg <- function(object, digits=3, alpha=0.05, ...)
{
conf.level = alpha
x <- object
nChain = x$setup$nChain
# convergence diagnostics
psrf <- NULL
if(nChain > 1){
if(x$class[2] == "ID")
{
theta <- x$chain1$theta.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
theta <- cbind(theta, x[[nam]]$theta.p)
}
psrftheta <- matrix(calcPSR(theta), 1, 1)
dimnames(psrftheta) <- list("", "")
beta.names <- unique(c(x$chain1$covNames1, x$chain1$covNames2, x$chain1$covNames3))
nP <- length(beta.names)
output <- matrix(NA, nrow=nP, ncol=3)
dimnames(output) <- list(beta.names, c("beta1", "beta2", "beta3"))
if(length(x$chain1$beta1.p) != 0){
#beta1
p1 = dim(x$chain1$beta1.p)[2]
psrfBeta1 <- rep(NA, p1)
for(j in 1:p1){
#namPara = paste("beta_", j, sep = "")
beta1 <- x$chain1$beta1[,j]
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
beta1 <- cbind(beta1, x[[nam]]$beta1[,j])
}
psrfBeta1[j] <- calcPSR(beta1)
}
for(i in 1:nP)
{
for(k in 1:p1) if(x$chain1$covNames1[k] == beta.names[i]) output[i,1] <- psrfBeta1[k]
}
}
if(length(x$chain1$beta2.p) != 0){
#beta2
p2 = dim(x$chain1$beta2.p)[2]
psrfBeta2 <- rep(NA, p2)
for(j in 1:p2){
#namPara = paste("beta_", j, sep = "")
beta2 <- x$chain1$beta2[,j]
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
beta2 <- cbind(beta2, x[[nam]]$beta2[,j])
}
psrfBeta2[j] <- calcPSR(beta2)
}
for(i in 1:nP)
{
for(k in 1:p2) if(x$chain1$covNames2[k] == beta.names[i]) output[i,2] <- psrfBeta2[k]
}
}
if(length(x$chain1$beta3.p) != 0){
#beta3
p3 = dim(x$chain1$beta3.p)[2]
psrfBeta3 <- rep(NA, p3)
for(j in 1:p3){
#namPara = paste("beta_", j, sep = "")
beta3 <- x$chain1$beta3[,j]
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
beta3 <- cbind(beta3, x[[nam]]$beta3[,j])
}
psrfBeta3[j] <- calcPSR(beta3)
}
for(i in 1:nP)
{
for(k in 1:p3) if(x$chain1$covNames3[k] == beta.names[i]) output[i,3] <- psrfBeta3[k]
}
}
psrfcoef <- NULL
if(nP > 0)
{
psrfcoef <- output
}
##
if(x$class[4] == "WB")
{
##
# alpha
alpha <- x$chain1$alpha1.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
alpha <- cbind(alpha, x[[nam]]$alpha1.p)
}
psrfAlpha1 <- calcPSR(alpha)
alpha <- x$chain1$alpha2.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
alpha <- cbind(alpha, x[[nam]]$alpha2.p)
}
psrfAlpha2 <- calcPSR(alpha)
alpha <- x$chain1$alpha3.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
alpha <- cbind(alpha, x[[nam]]$alpha3.p)
}
psrfAlpha3 <- calcPSR(alpha)
# kappa
kappa <- x$chain1$kappa1.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
kappa <- cbind(kappa, x[[nam]]$kappa1.p)
}
psrfKappa1 <- calcPSR(kappa)
kappa <- x$chain1$kappa2.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
kappa <- cbind(kappa, x[[nam]]$kappa2.p)
}
psrfKappa2 <- calcPSR(kappa)
kappa <- x$chain1$kappa3.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
kappa <- cbind(kappa, x[[nam]]$kappa3.p)
}
psrfKappa3 <- calcPSR(kappa)
bh <- matrix(c(psrfKappa1, psrfKappa2, psrfKappa3, psrfAlpha1, psrfAlpha2, psrfAlpha3), 2, 3, byrow = T)
dimnames(bh) <- list(c("kappa", "alpha"), c("h1", "h2", "h3"))
psrf <- list(theta=psrftheta, coef=psrfcoef, h0=bh)
}
if(x$class[4] == "PEM")
{
##
ntime1 = length(x$chain1$time_lambda1)
ntime2 = length(x$chain1$time_lambda2)
ntime3 = length(x$chain1$time_lambda3)
# lambda's
psrfLam1 <- rep(NA, ntime1)
for(j in 1:ntime1){
lambda1 <- x$chain1$lambda1.fin[,j]
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
lambda1 <- cbind(lambda1, x[[nam]]$lambda1.fin[,j])
}
psrfLam1[j] <- calcPSR(lambda1)
}
psrfLam2 <- rep(NA, ntime2)
for(j in 1:ntime2){
lambda2 <- x$chain1$lambda2.fin[,j]
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
lambda2 <- cbind(lambda2, x[[nam]]$lambda2.fin[,j])
}
psrfLam2[j] <- calcPSR(lambda2)
}
psrfLam3 <- rep(NA, ntime3)
for(j in 1:ntime3){
lambda3 <- x$chain1$lambda3.fin[,j]
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
lambda3 <- cbind(lambda3, x[[nam]]$lambda3.fin[,j])
}
psrfLam3[j] <- calcPSR(lambda3)
}
# mu_lam
mu <- x$chain1$mu_lam1.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
mu <- cbind(mu, x[[nam]]$mu_lam1.p)
}
psrfMu1 <- calcPSR(mu)
mu <- x$chain1$mu_lam2.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
mu <- cbind(mu, x[[nam]]$mu_lam2.p)
}
psrfMu2 <- calcPSR(mu)
mu <- x$chain1$mu_lam3.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
mu <- cbind(mu, x[[nam]]$mu_lam3.p)
}
psrfMu3 <- calcPSR(mu)
# sigSq_lam
sig <- x$chain1$sigSq_lam1.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
sig <- cbind(sig, x[[nam]]$sigSq_lam1.p)
}
psrfSig1 <- calcPSR(sig)
sig <- x$chain1$sigSq_lam2.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
sig <- cbind(sig, x[[nam]]$sigSq_lam2.p)
}
psrfSig2 <- calcPSR(sig)
sig <- x$chain1$sigSq_lam3.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
sig <- cbind(sig, x[[nam]]$sigSq_lam3.p)
}
psrfSig3 <- calcPSR(sig)
# J
J <- x$chain1$K1.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
J <- cbind(J, x[[nam]]$K1.p)
}
psrfJ1 <- calcPSR(J)
J <- x$chain1$K2.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
J <- cbind(J, x[[nam]]$K2.p)
}
psrfJ2 <- calcPSR(J)
J <- x$chain1$K3.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
J <- cbind(J, x[[nam]]$K3.p)
}
psrfJ3 <- calcPSR(J)
bh <- matrix(c(psrfMu1, psrfMu2, psrfMu3, psrfSig1, psrfSig2, psrfSig3, psrfJ1, psrfJ2, psrfJ3), 3, 3, byrow = T)
dimnames(bh) <- list(c("mu", "sigmaSq", "K"), c("h1", "h2", "h3"))
psrf <- list(theta=psrftheta, coef=psrfcoef, h0=bh, lambda1=psrfLam1, lambda2=psrfLam2, lambda3=psrfLam3)
}
}
if(x$class[2] == "Surv")
{
beta.names <- c(x$chain1$covNames)
nP <- length(beta.names)
output <- matrix(NA, nrow=nP, ncol=1)
dimnames(output) <- list(beta.names, c("beta"))
if(length(x$chain1$beta.p) != 0){
#beta
p = dim(x$chain1$beta.p)[2]
psrfBeta <- rep(NA, p)
for(j in 1:p){
beta <- x$chain1$beta.p[,j]
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
beta <- cbind(beta, x[[nam]]$beta.p[,j])
}
psrfBeta[j] <- calcPSR(beta)
}
for(i in 1:nP)
{
for(k in 1:p) if(x$chain1$covNames[k] == beta.names[i]) output[i,1] <- psrfBeta[k]
}
}
psrfcoef <- NULL
if(nP > 0)
{
psrfcoef <- output
}
if(x$class[4] == "WB")
{
##
# alpha
alpha <- x$chain1$alpha.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
alpha <- cbind(alpha, x[[nam]]$alpha.p)
}
psrfAlpha <- calcPSR(alpha)
# kappa
kappa <- x$chain1$kappa.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
kappa <- cbind(kappa, x[[nam]]$kappa.p)
}
psrfKappa <- calcPSR(kappa)
bh <- matrix(c(psrfKappa, psrfAlpha), 2, 1, byrow = T)
dimnames(bh) <- list(c("kappa", "alpha"), c("h"))
psrf <- list(coef=psrfcoef, h0=bh)
}
if(x$class[4] == "PEM")
{
##
ntime = length(x$chain1$time_lambda)
# lambda
psrfLam <- rep(NA, ntime)
for(j in 1:ntime){
namPara = paste("beta_", j, sep = "")
lambda <- x$chain1$lambda.fin[,j]
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
lambda <- cbind(lambda, x[[nam]]$lambda.fin[,j])
}
psrfLam[j] <- calcPSR(lambda)
}
# mu_lam
mu <- x$chain1$mu_lam.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
mu <- cbind(mu, x[[nam]]$mu_lam.p)
}
psrfMu <- calcPSR(mu)
# sigSq_lam
sig <- x$chain1$sigSq_lam.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
sig <- cbind(sig, x[[nam]]$sigSq_lam.p)
}
psrfSig <- calcPSR(sig)
# J
J <- x$chain1$K.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
J <- cbind(J, x[[nam]]$K.p)
}
psrfJ <- calcPSR(J)
bh <- matrix(c(psrfMu, psrfSig, psrfJ), 3, 1, byrow = T)
dimnames(bh) <- list(c("mu", "sigmaSq", "K"), c("h"))
psrf <- list(coef=psrfcoef, h0=bh, lambda=psrfLam)
}
}
}
# estimates
if(x$class[2] == "ID")
{
##
theta.p <- x$chain1$theta.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
theta.p <- rbind(theta.p, x[[nam]]$theta.p)
}
}
theta.pMed <- apply(theta.p, 2, median)
theta.pUb <- apply(theta.p, 2, quantile, prob = (1-conf.level/2))
theta.pLb <- apply(theta.p, 2, quantile, prob = conf.level/2)
tbl_theta <- matrix(NA, 1, 3)
dimnames(tbl_theta) <- list("", c( "theta", "LL", "UL"))
tbl_theta[,1] <- theta.pMed
tbl_theta[,2] <- theta.pLb
tbl_theta[,3] <- theta.pUb
##
beta.names <- unique(c(x$chain1$covNames1, x$chain1$covNames2, x$chain1$covNames3))
nP <- length(beta.names)
output <- matrix(NA, nrow=nP, ncol=9)
dimnames(output) <- list(beta.names, c("exp(beta1)", "LL", "UL", "exp(beta2)", "LL", "UL", "exp(beta3)", "LL", "UL"))
if(length(x$chain1$beta1.p) != 0){
#beta1
p1 = dim(x$chain1$beta1.p)[2]
beta.p <- x$chain1$beta1.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
beta.p <- rbind(beta.p, x[[nam]]$beta1.p)
}
}
beta.pMed <- apply(exp(beta.p), 2, median)
beta.pSd <- apply(exp(beta.p), 2, sd)
beta.pUb <- apply(exp(beta.p), 2, quantile, prob = (1-conf.level/2))
beta.pLb <- apply(exp(beta.p), 2, quantile, prob = conf.level/2)
tbl1 <- matrix(NA, p1, 3)
rownames(tbl1) <- x$chain1$covNames1
tbl1[,1] <- beta.pMed
tbl1[,2] <- beta.pLb
tbl1[,3] <- beta.pUb
for(i in 1:nP)
{
for(k in 1:p1) if(x$chain1$covNames1[k] == beta.names[i]) output[i,1:3] <- tbl1[k,]
}
}
if(length(x$chain1$beta2.p) != 0){
#beta2
p2 = dim(x$chain1$beta2.p)[2]
beta.p <- x$chain1$beta2.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
beta.p <- rbind(beta.p, x[[nam]]$beta2.p)
}
}
beta.pMed <- apply(exp(beta.p), 2, median)
beta.pSd <- apply(exp(beta.p), 2, sd)
beta.pUb <- apply(exp(beta.p), 2, quantile, prob = (1-conf.level/2))
beta.pLb <- apply(exp(beta.p), 2, quantile, prob = conf.level/2)
tbl2 <- matrix(NA, p2, 3)
rownames(tbl2) <- x$chain1$covNames2
tbl2[,1] <- beta.pMed
tbl2[,2] <- beta.pLb
tbl2[,3] <- beta.pUb
for(i in 1:nP)
{
for(k in 1:p2) if(x$chain1$covNames2[k] == beta.names[i]) output[i,4:6] <- tbl2[k,]
}
}
if(length(x$chain1$beta3.p) != 0){
#beta3
p3 = dim(x$chain1$beta3.p)[2]
beta.p <- x$chain1$beta3.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
beta.p <- rbind(beta.p, x[[nam]]$beta3.p)
}
}
beta.pMed <- apply(exp(beta.p), 2, median)
beta.pSd <- apply(exp(beta.p), 2, sd)
beta.pUb <- apply(exp(beta.p), 2, quantile, prob = (1-conf.level/2))
beta.pLb <- apply(exp(beta.p), 2, quantile, prob = conf.level/2)
tbl3 <- matrix(NA, p3, 3)
rownames(tbl3) <- x$chain1$covNames3
tbl3[,1] <- beta.pMed
tbl3[,2] <- beta.pLb
tbl3[,3] <- beta.pUb
for(i in 1:nP)
{
for(k in 1:p3) if(x$chain1$covNames3[k] == beta.names[i]) output[i,7:9] <- tbl3[k,]
}
}
output.coef <- NULL
if(nP > 0)
{
output.coef <- output
}
if(x$class[4] == "WB")
{
##
alpha.p <- x$chain1$alpha1.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
alpha.p <- rbind(alpha.p, x[[nam]]$alpha1.p)
}
}
alpha.pMed <- apply(log(alpha.p), 2, median)
alpha.pUb <- apply(log(alpha.p), 2, quantile, prob = (1-conf.level/2))
alpha.pLb <- apply(log(alpha.p), 2, quantile, prob = conf.level/2)
tbl_a1 <- c(alpha.pMed,alpha.pLb, alpha.pUb)
##
alpha.p <- x$chain1$alpha2.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
alpha.p <- rbind(alpha.p, x[[nam]]$alpha2.p)
}
}
alpha.pMed <- apply(log(alpha.p), 2, median)
alpha.pUb <- apply(log(alpha.p), 2, quantile, prob = (1-conf.level/2))
alpha.pLb <- apply(log(alpha.p), 2, quantile, prob = conf.level/2)
tbl_a2 <- c(alpha.pMed,alpha.pLb, alpha.pUb)
##
alpha.p <- x$chain1$alpha3.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
alpha.p <- rbind(alpha.p, x[[nam]]$alpha3.p)
}
}
alpha.pMed <- apply(log(alpha.p), 2, median)
alpha.pUb <- apply(log(alpha.p), 2, quantile, prob = (1-conf.level/2))
alpha.pLb <- apply(log(alpha.p), 2, quantile, prob = conf.level/2)
tbl_a3 <- c(alpha.pMed,alpha.pLb, alpha.pUb)
##
kappa.p <- x$chain1$kappa1.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
kappa.p <- rbind(kappa.p, x[[nam]]$kappa1.p)
}
}
kappa.pMed <- apply(log(kappa.p), 2, median)
kappa.pUb <- apply(log(kappa.p), 2, quantile, prob = (1-conf.level/2))
kappa.pLb <- apply(log(kappa.p), 2, quantile, prob = conf.level/2)
tbl_k1 <- c(kappa.pMed, kappa.pLb, kappa.pUb)
##
kappa.p <- x$chain1$kappa2.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
kappa.p <- rbind(kappa.p, x[[nam]]$kappa2.p)
}
}
kappa.pMed <- apply(log(kappa.p), 2, median)
kappa.pUb <- apply(log(kappa.p), 2, quantile, prob = (1-conf.level/2))
kappa.pLb <- apply(log(kappa.p), 2, quantile, prob = conf.level/2)
tbl_k2 <- c(kappa.pMed, kappa.pLb, kappa.pUb)
##
kappa.p <- x$chain1$kappa3.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
kappa.p <- rbind(kappa.p, x[[nam]]$kappa3.p)
}
}
kappa.pMed <- apply(log(kappa.p), 2, median)
kappa.pUb <- apply(log(kappa.p), 2, quantile, prob = (1-conf.level/2))
kappa.pLb <- apply(log(kappa.p), 2, quantile, prob = conf.level/2)
tbl_k3 <- c(kappa.pMed, kappa.pLb, kappa.pUb)
bh <- matrix(c(tbl_k1, tbl_k2, tbl_k3, tbl_a1, tbl_a2, tbl_a3), 2, 9, byrow = T)
dimnames(bh) <- list(c("Weibull: log-kappa", "Weibull: log-alpha"), c("h1-PM", "LL", "UL", "h2-PM", "LL", "UL", "h3-PM", "LL", "UL"))
value <- list(classFit=x$class, psrf=psrf, theta=tbl_theta, coef=output.coef, h0=bh)
}
if(x$class[4] == "PEM")
{
##
mu_lam.p <- x$chain1$mu_lam1.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
mu_lam.p <- rbind(mu_lam.p, x[[nam]]$mu_lam1.p)
}
}
mu_lam.pMed <- apply(mu_lam.p, 2, median)
mu_lam.pUb <- apply(mu_lam.p, 2, quantile, prob = (1-conf.level/2))
mu_lam.pLb <- apply(mu_lam.p, 2, quantile, prob = conf.level/2)
tbl_m1 <- c(mu_lam.pMed, mu_lam.pLb, mu_lam.pUb)
##
mu_lam.p <- x$chain1$mu_lam2.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
mu_lam.p <- rbind(mu_lam.p, x[[nam]]$mu_lam2.p)
}
}
mu_lam.pMed <- apply(mu_lam.p, 2, median)
mu_lam.pUb <- apply(mu_lam.p, 2, quantile, prob = (1-conf.level/2))
mu_lam.pLb <- apply(mu_lam.p, 2, quantile, prob = conf.level/2)
tbl_m2 <- c(mu_lam.pMed, mu_lam.pLb, mu_lam.pUb)
##
mu_lam.p <- x$chain1$mu_lam3.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
mu_lam.p <- rbind(mu_lam.p, x[[nam]]$mu_lam3.p)
}
}
mu_lam.pMed <- apply(mu_lam.p, 2, median)
mu_lam.pUb <- apply(mu_lam.p, 2, quantile, prob = (1-conf.level/2))
mu_lam.pLb <- apply(mu_lam.p, 2, quantile, prob = conf.level/2)
tbl_m3 <- c(mu_lam.pMed, mu_lam.pLb, mu_lam.pUb)
##
sigSq_lam.p <- x$chain1$sigSq_lam1.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
sigSq_lam.p <- rbind(sigSq_lam.p, x[[nam]]$sigSq_lam1.p)
}
}
sigSq_lam.pMed <- apply(sigSq_lam.p, 2, median)
sigSq_lam.pUb <- apply(sigSq_lam.p, 2, quantile, prob = (1-conf.level/2))
sigSq_lam.pLb <- apply(sigSq_lam.p, 2, quantile, prob = conf.level/2)
tbl_s1 <- c(sigSq_lam.pMed, sigSq_lam.pLb, sigSq_lam.pUb)
##
sigSq_lam.p <- x$chain1$sigSq_lam2.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
sigSq_lam.p <- rbind(sigSq_lam.p, x[[nam]]$sigSq_lam2.p)
}
}
sigSq_lam.pMed <- apply(sigSq_lam.p, 2, median)
sigSq_lam.pUb <- apply(sigSq_lam.p, 2, quantile, prob = (1-conf.level/2))
sigSq_lam.pLb <- apply(sigSq_lam.p, 2, quantile, prob = conf.level/2)
tbl_s2 <- c(sigSq_lam.pMed, sigSq_lam.pLb, sigSq_lam.pUb)
##
sigSq_lam.p <- x$chain1$sigSq_lam3.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
sigSq_lam.p <- rbind(sigSq_lam.p, x[[nam]]$sigSq_lam3.p)
}
}
sigSq_lam.pMed <- apply(sigSq_lam.p, 2, median)
sigSq_lam.pUb <- apply(sigSq_lam.p, 2, quantile, prob = (1-conf.level/2))
sigSq_lam.pLb <- apply(sigSq_lam.p, 2, quantile, prob = conf.level/2)
tbl_s3 <- c(sigSq_lam.pMed, sigSq_lam.pLb, sigSq_lam.pUb)
##
J.p <- x$chain1$K1.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
J.p <- rbind(J.p, x[[nam]]$K1.p)
}
}
J.pMed <- apply(J.p, 2, median)
J.pUb <- apply(J.p, 2, quantile, prob = (1-conf.level/2))
J.pLb <- apply(J.p, 2, quantile, prob = conf.level/2)
tbl_j1 <- c(J.pMed, J.pLb, J.pUb)
##
J.p <- x$chain1$K2.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
J.p <- rbind(J.p, x[[nam]]$K2.p)
}
}
J.pMed <- apply(J.p, 2, median)
J.pUb <- apply(J.p, 2, quantile, prob = (1-conf.level/2))
J.pLb <- apply(J.p, 2, quantile, prob = conf.level/2)
tbl_j2 <- c(J.pMed, J.pLb, J.pUb)
##
J.p <- x$chain1$K3.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
J.p <- rbind(J.p, x[[nam]]$K3.p)
}
}
J.pMed <- apply(J.p, 2, median)
J.pUb <- apply(J.p, 2, quantile, prob = (1-conf.level/2))
J.pLb <- apply(J.p, 2, quantile, prob = conf.level/2)
tbl_j3 <- c(J.pMed, J.pLb, J.pUb)
bh <- matrix(c(tbl_m1, tbl_m2, tbl_m3, tbl_s1, tbl_s2, tbl_s3, tbl_j1, tbl_j2, tbl_j3), 3, 9, byrow = T)
dimnames(bh) <- list(c("mu", "sigmaSq", "K"), c("h1-PM", "LL", "UL", "h2-PM", "LL", "UL", "h3-PM", "LL", "UL"))
##
lambda.p <- x$chain1$lambda1.fin
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
lambda.p <- rbind(lambda.p, x[[nam]]$lambda1.fin)
}
}
lambda.pMed <- apply(lambda.p, 2, median)
lambda.pUb <- apply(lambda.p, 2, quantile, prob = (1-conf.level/2))
lambda.pLb <- apply(lambda.p, 2, quantile, prob = conf.level/2)
lambda1 <- cbind(x$chain1$time_lambda1, lambda.pMed, lambda.pLb, lambda.pUb)
dimnames(lambda1) <- list(rep("", length(x$chain1$time_lambda1)), c("time", "lambda1-PM", "LL", "UL"))
##
lambda.p <- x$chain1$lambda2.fin
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
lambda.p <- rbind(lambda.p, x[[nam]]$lambda2.fin)
}
}
lambda.pMed <- apply(lambda.p, 2, median)
lambda.pUb <- apply(lambda.p, 2, quantile, prob = (1-conf.level/2))
lambda.pLb <- apply(lambda.p, 2, quantile, prob = conf.level/2)
lambda2 <- cbind(x$chain1$time_lambda2, lambda.pMed, lambda.pLb, lambda.pUb)
dimnames(lambda2) <- list(rep("", length(x$chain1$time_lambda2)), c("time", "lambda2-PM", "LL", "UL"))
##
lambda.p <- x$chain1$lambda3.fin
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
lambda.p <- rbind(lambda.p, x[[nam]]$lambda3.fin)
}
}
lambda.pMed <- apply(lambda.p, 2, median)
lambda.pUb <- apply(lambda.p, 2, quantile, prob = (1-conf.level/2))
lambda.pLb <- apply(lambda.p, 2, quantile, prob = conf.level/2)
lambda3 <- cbind(x$chain1$time_lambda3, lambda.pMed, lambda.pLb, lambda.pUb)
dimnames(lambda3) <- list(rep("", length(x$chain1$time_lambda3)), c("time", "lambda3-PM", "LL", "UL"))
value <- list(classFit=x$class, psrf=psrf, theta=tbl_theta, coef=output.coef, h0=bh, lambda1=lambda1, lambda2=lambda2, lambda3=lambda3)
}
if(x$class[3] == "Cor")
{
if(x$class[5] == "MVN")
{
nS <- dim(x$chain1$Sigma_V.p)[3]
Sigma <- array(NA, c(3,3, nS*nChain))
Sigma[,,1:nS] <- x$chain1$Sigma_V.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
Sigma[,,(nS*(i-1)+1):(nS*i)] <- x[[nam]]$Sigma_V.p
}
}
}
if(x$class[5] == "DPM")
{
##
tau.p <- x$chain1$tau.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
tau.p <- rbind(tau.p, x[[nam]]$tau.p)
}
}
tau.pMed <- apply(tau.p, 2, median)
tau.pUb <- apply(tau.p, 2, quantile, prob = (1-conf.level/2))
tau.pLb <- apply(tau.p, 2, quantile, prob = conf.level/2)
tbl_tau <- matrix(NA, 1, 3)
dimnames(tbl_tau) <- list("", c( "tau", "LL", "UL"))
tbl_tau[,1] <- tau.pMed
tbl_tau[,2] <- tau.pLb
tbl_tau[,3] <- tau.pUb
nS <- dim(x$chain1$Sigma.p)[3]
Sigma <- array(NA, c(3,3, nS*nChain))
Sigma[,,1:nS] <- calVar_DPM_MVN(x$chain1)
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
Sigma[,,(nS*(i-1)+1):(nS*i)] <- calVar_DPM_MVN(x[[nam]])
}
}
value$tau <- tbl_tau
}
Sigma.Med <- apply(Sigma, c(1,2), median)
Sigma.Sd <- apply(Sigma, c(1,2), sd)
Sigma.Ub <- apply(Sigma, c(1,2), quantile, prob = (1-conf.level/2))
Sigma.Lb <- apply(Sigma, c(1,2), quantile, prob = conf.level/2)
dimnames(Sigma.Med) <- list(c("", "", ""), c("Sigma_V-PM", "", ""))
dimnames(Sigma.Sd) <- list(c("", "", ""), c("Sigma_V-SD", "", ""))
dimnames(Sigma.Lb) <- list(c("", "", ""), c("Sigma_V-LL", "", ""))
dimnames(Sigma.Ub) <- list(c("", "", ""), c("Sigma_V-UL", "", ""))
value$Sigma.PM <- Sigma.Med
value$Sigma.SD <- Sigma.Sd
value$Sigma.UL <- Sigma.Ub
value$Sigma.LL <- Sigma.Lb
}
## DIC and LPML
value$DIC = x$DIC
value$LPML = x$LPML
}
if(x$class[2] == "Surv")
{
##
beta.names <- c(x$chain1$covNames)
nP <- length(beta.names)
output <- matrix(NA, nrow=nP, ncol=3)
dimnames(output) <- list(beta.names, c("exp(beta)", "LL", "UL"))
if(length(x$chain1$beta.p) != 0){
#beta
p = dim(x$chain1$beta.p)[2]
beta.p <- x$chain1$beta.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
beta.p <- rbind(beta.p, x[[nam]]$beta.p)
}
}
beta.pMed <- apply(exp(beta.p), 2, median)
beta.pSd <- apply(exp(beta.p), 2, sd)
beta.pUb <- apply(exp(beta.p), 2, quantile, prob = (1-conf.level/2))
beta.pLb <- apply(exp(beta.p), 2, quantile, prob = conf.level/2)
tbl <- matrix(NA, p, 3)
rownames(tbl) <- x$chain1$covNames
tbl[,1] <- beta.pMed
tbl[,2] <- beta.pLb
tbl[,3] <- beta.pUb
for(i in 1:nP)
{
for(k in 1:p) if(x$chain1$covNames[k] == beta.names[i]) output[i,1:3] <- tbl[k,]
}
}
output.coef <- NULL
if(nP > 0)
{
output.coef <- output
}
if(x$class[4] == "WB")
{
##
alpha.p <- x$chain1$alpha.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
alpha.p <- rbind(alpha.p, x[[nam]]$alpha.p)
}
}
alpha.pMed <- apply(log(alpha.p), 2, median)
alpha.pUb <- apply(log(alpha.p), 2, quantile, prob = (1-conf.level/2))
alpha.pLb <- apply(log(alpha.p), 2, quantile, prob = conf.level/2)
tbl_a <- c(alpha.pMed,alpha.pLb, alpha.pUb)
##
kappa.p <- x$chain1$kappa.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
kappa.p <- rbind(kappa.p, x[[nam]]$kappa.p)
}
}
kappa.pMed <- apply(log(kappa.p), 2, median)
kappa.pUb <- apply(log(kappa.p), 2, quantile, prob = (1-conf.level/2))
kappa.pLb <- apply(log(kappa.p), 2, quantile, prob = conf.level/2)
tbl_k <- c(kappa.pMed, kappa.pLb, kappa.pUb)
bh <- matrix(c(tbl_a, tbl_k), 2, 3, byrow = T)
dimnames(bh) <- list(c("Weibull: log-kappa", "Weibull: log-alpha"), c("h-PM", "LL", "UL"))
value <- list(coef=output.coef, h0=bh, psrf=psrf, classFit=x$class)
}
if(x$class[4] == "PEM")
{
##
mu_lam.p <- x$chain1$mu_lam.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
mu_lam.p <- rbind(mu_lam.p, x[[nam]]$mu_lam.p)
}
}
mu_lam.pMed <- apply(mu_lam.p, 2, median)
mu_lam.pUb <- apply(mu_lam.p, 2, quantile, prob = (1-conf.level/2))
mu_lam.pLb <- apply(mu_lam.p, 2, quantile, prob = conf.level/2)
tbl_m <- c(mu_lam.pMed, mu_lam.pLb, mu_lam.pUb)
##
sigSq_lam.p <- x$chain1$sigSq_lam.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
sigSq_lam.p <- rbind(sigSq_lam.p, x[[nam]]$sigSq_lam.p)
}
}
sigSq_lam.pMed <- apply(sigSq_lam.p, 2, median)
sigSq_lam.pUb <- apply(sigSq_lam.p, 2, quantile, prob = (1-conf.level/2))
sigSq_lam.pLb <- apply(sigSq_lam.p, 2, quantile, prob = conf.level/2)
tbl_s <- c(sigSq_lam.pMed, sigSq_lam.pLb, sigSq_lam.pUb)
##
J.p <- x$chain1$K.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
J.p <- rbind(J.p, x[[nam]]$K.p)
}
}
J.pMed <- apply(J.p, 2, median)
J.pUb <- apply(J.p, 2, quantile, prob = (1-conf.level/2))
J.pLb <- apply(J.p, 2, quantile, prob = conf.level/2)
tbl_j <- c(J.pMed, J.pLb, J.pUb)
bh <- matrix(c(tbl_m, tbl_s, tbl_j), 3, 3, byrow = T)
dimnames(bh) <- list(c("mu", "sigmaSq", "K"), c("h-PM", "LL", "UL"))
##
lambda.p <- x$chain1$lambda.fin
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
lambda.p <- rbind(lambda.p, x[[nam]]$lambda.fin)
}
}
lambda.pMed <- apply(lambda.p, 2, median)
lambda.pUb <- apply(lambda.p, 2, quantile, prob = (1-conf.level/2))
lambda.pLb <- apply(lambda.p, 2, quantile, prob = conf.level/2)
lambda <- cbind(x$chain1$time_lambda, lambda.pMed, lambda.pLb, lambda.pUb)
dimnames(lambda) <- list(rep("", length(x$chain1$time_lambda)), c("time", "lambda-PM", "LL", "UL"))
value <- list(coef=output.coef, h0=bh, psrf=psrf, lambda=lambda, classFit=x$class)
}
if(x$class[3] == "Cor")
{
if(x$class[5] == "Normal")
{
#sigmaV
sigV <- 1/x$chain1$zeta.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
sigV <- rbind(sigV, 1/x[[nam]]$zeta.p)
}
}
}
if(x$class[5] == "DPM")
{
##
tau.p <- x$chain1$tau.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
tau.p <- rbind(tau.p, x[[nam]]$tau.p)
}
}
tau.pMed <- apply(tau.p, 2, median)
tau.pUb <- apply(tau.p, 2, quantile, prob = (1-conf.level/2))
tau.pLb <- apply(tau.p, 2, quantile, prob = conf.level/2)
tbl_tau <- matrix(NA, 1, 3)
dimnames(tbl_tau) <- list("", c( "tau", "LL", "UL"))
tbl_tau[,1] <- tau.pMed
tbl_tau[,2] <- tau.pLb
tbl_tau[,3] <- tau.pUb
nS <- dim(x$chain1$zeta.p)[1]
sigV <- rep(NA, nS*nChain)
sigV[1:nS] <- calVar_DPM_Normal(x$chain1)
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
sigV[(nS*(i-1)+1):(nS*i)] <- calVar_DPM_Normal(x[[nam]])
}
}
value$tau <- tbl_tau
}
sigVMed <- median(sigV)
sigVSd <- sd(sigV)
sigVUb <- quantile(sigV, prob = (1-conf.level/2))
sigVLb <- quantile(sigV, prob = conf.level/2)
tbl_sigV <- matrix(NA, nrow=1, ncol=3)
tbl_sigV[,1] <- sigVMed
tbl_sigV[,2] <- sigVLb
tbl_sigV[,3] <- sigVUb
dimnames(tbl_sigV) <- list("", c("sigma_V-PM", "LL", "UL"))
value$sigma_V <- tbl_sigV
}
}
value$setup <- x$setup
value$conf.level <- conf.level
class(value) <- "summ.Bayes_HReg"
return(value)
}
summary.Bayes_AFT <- function(object, digits=3, alpha=0.05, ...)
{
conf.level = alpha
x <- object
nChain = x$setup$nChain
# convergence diagnostics
psrf <- NULL
if(nChain > 1){
if(x$class[2] == "ID")
{
theta <- x$chain1$theta.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
theta <- cbind(theta, x[[nam]]$theta.p)
}
psrftheta <- matrix(calcPSR(theta), 1, 1)
dimnames(psrftheta) <- list("", "")
beta.names <- unique(c(x$chain1$covNames1, x$chain1$covNames2, x$chain1$covNames3))
nP <- length(beta.names)
output <- matrix(NA, nrow=nP, ncol=3)
dimnames(output) <- list(beta.names, c("beta1", "beta2", "beta3"))
if(length(x$chain1$beta1.p) != 0){
#beta1
p1 = dim(x$chain1$beta1.p)[2]
psrfBeta1 <- rep(NA, p1)
for(j in 1:p1){
#namPara = paste("beta_", j, sep = "")
beta1 <- x$chain1$beta1[,j]
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
beta1 <- cbind(beta1, x[[nam]]$beta1[,j])
}
psrfBeta1[j] <- calcPSR(beta1)
}
for(i in 1:nP)
{
for(k in 1:p1) if(x$chain1$covNames1[k] == beta.names[i]) output[i,1] <- psrfBeta1[k]
}
}
if(length(x$chain1$beta2.p) != 0){
#beta2
p2 = dim(x$chain1$beta2.p)[2]
psrfBeta2 <- rep(NA, p2)
for(j in 1:p2){
#namPara = paste("beta_", j, sep = "")
beta2 <- x$chain1$beta2[,j]
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
beta2 <- cbind(beta2, x[[nam]]$beta2[,j])
}
psrfBeta2[j] <- calcPSR(beta2)
}
for(i in 1:nP)
{
for(k in 1:p2) if(x$chain1$covNames2[k] == beta.names[i]) output[i,2] <- psrfBeta2[k]
}
}
if(length(x$chain1$beta3.p) != 0){
#beta3
p3 = dim(x$chain1$beta3.p)[2]
psrfBeta3 <- rep(NA, p3)
for(j in 1:p3){
#namPara = paste("beta_", j, sep = "")
beta3 <- x$chain1$beta3[,j]
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
beta3 <- cbind(beta3, x[[nam]]$beta3[,j])
}
psrfBeta3[j] <- calcPSR(beta3)
}
for(i in 1:nP)
{
for(k in 1:p3) if(x$chain1$covNames3[k] == beta.names[i]) output[i,3] <- psrfBeta3[k]
}
}
psrfcoef <- NULL
if(nP > 0)
{
psrfcoef <- output
}
if(x$class[4] == "LN")
{
##
# mu
mu <- x$chain1$mu1.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
mu <- cbind(mu, x[[nam]]$mu1.p)
}
psrfmu1 <- calcPSR(mu)
mu <- x$chain1$mu2.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
mu <- cbind(mu, x[[nam]]$mu2.p)
}
psrfmu2 <- calcPSR(mu)
mu <- x$chain1$mu3.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
mu <- cbind(mu, x[[nam]]$mu3.p)
}
psrfmu3 <- calcPSR(mu)
# sigSq
sigSq <- x$chain1$sigSq1.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
sigSq <- cbind(sigSq, x[[nam]]$sigSq1.p)
}
psrfSigSq1 <- calcPSR(sigSq)
sigSq <- x$chain1$sigSq2.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
sigSq <- cbind(sigSq, x[[nam]]$sigSq2.p)
}
psrfSigSq2 <- calcPSR(sigSq)
sigSq <- x$chain1$sigSq3.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
sigSq <- cbind(sigSq, x[[nam]]$sigSq3.p)
}
psrfSigSq3 <- calcPSR(sigSq)
bh <- matrix(c(psrfmu1, psrfmu2, psrfmu3, psrfSigSq1, psrfSigSq2, psrfSigSq3), 2, 3, byrow = T)
dimnames(bh) <- list(c("mu", "sigmaSq"), c("g=1", "g=2", "g=3"))
psrf <- list(theta=psrftheta, coef=psrfcoef, h0=bh)
}
if(x$class[4] == "DPM")
{
##
# tau
tau <- x$chain1$tau1.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
tau <- cbind(tau, x[[nam]]$tau1.p)
}
psrfTau1 <- calcPSR(tau)
tau <- x$chain1$tau2.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
tau <- cbind(tau, x[[nam]]$tau2.p)
}
psrfTau2 <- calcPSR(tau)
tau <- x$chain1$tau3.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
tau <- cbind(tau, x[[nam]]$tau3.p)
}
psrfTau3 <- calcPSR(tau)
bh <- matrix(c(psrfTau1, psrfTau2, psrfTau3), 1, 3, byrow = T)
dimnames(bh) <- list(c("tau"), c("g=1", "g=2", "g=3"))
psrf <- list(theta=psrftheta, coef=psrfcoef, h0=bh)
}
}
if(x$class[2] == "Surv")
{
beta.names <- c(x$chain1$covNames)
nP <- length(beta.names)
output <- matrix(NA, nrow=nP, ncol=1)
dimnames(output) <- list(beta.names, c("beta"))
if(length(x$chain1$beta.p) != 0){
#beta
p = dim(x$chain1$beta.p)[2]
psrfBeta <- rep(NA, p)
for(j in 1:p){
beta <- x$chain1$beta.p[,j]
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
beta <- cbind(beta, x[[nam]]$beta.p[,j])
}
psrfBeta[j] <- calcPSR(beta)
}
for(i in 1:nP)
{
for(k in 1:p) if(x$chain1$covNames[k] == beta.names[i]) output[i,1] <- psrfBeta[k]
}
}
psrfcoef <- NULL
if(nP > 0)
{
psrfcoef <- output
}
if(x$class[4] == "LN")
{
##
# mu
mu <- x$chain1$mu.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
mu <- cbind(mu, x[[nam]]$mu.p)
}
psrfmu <- calcPSR(mu)
# sigSq
sigSq <- x$chain1$sigSq.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
sigSq <- cbind(sigSq, x[[nam]]$sigSq.p)
}
psrfSigSq <- calcPSR(sigSq)
bh <- matrix(c(psrfmu, psrfSigSq), 2, 1, byrow = T)
dimnames(bh) <- list(c("mu", "sigmaSq"), c("h"))
psrf <- list(coef=psrfcoef, h0=bh)
}
if(x$class[4] == "DPM")
{
##
# tau
tau <- x$chain1$tau.p
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
tau <- cbind(tau, x[[nam]]$tau.p)
}
psrfTau <- calcPSR(tau)
bh <- matrix(c(psrfTau), 1, 1, byrow = T)
dimnames(bh) <- list(c("tau"), c("h"))
psrf <- list(coef=psrfcoef, h0=bh)
}
}
}
# model comparison
if(nChain > 1){
if(x$class[2] == "ID")
{
# DIC and LPML
logLH_mean <- x$chain1$logLH_mean
LH_i_mean <- x$chain1$LH_i_mean
invLH_i_mean <- x$chain1$invLH_i_mean
for(i in 2:nChain){
nam <- paste("chain", i, sep = "")
logLH_mean <- cbind(logLH_mean, x[[nam]]$logLH_mean)
LH_i_mean <- rbind(LH_i_mean, x[[nam]]$LH_i_mean)
invLH_i_mean <- rbind(invLH_i_mean, x[[nam]]$invLH_i_mean)
}
dev <- -2*mean(logLH_mean)
DIC = 2*dev + 2*sum(log(apply(LH_i_mean, 2, mean)))
LPML = -sum(log(apply(invLH_i_mean, 2, mean)))
}
}
# estimates
if(x$class[2] == "ID")
{
##
theta.p <- x$chain1$theta.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
theta.p <- rbind(theta.p, x[[nam]]$theta.p)
}
}
theta.pMed <- apply(theta.p, 2, median)
theta.pUb <- apply(theta.p, 2, quantile, prob = (1-conf.level/2))
theta.pLb <- apply(theta.p, 2, quantile, prob = conf.level/2)
tbl_theta <- matrix(NA, 1, 3)
dimnames(tbl_theta) <- list("", c( "theta", "LL", "UL"))
tbl_theta[,1] <- theta.pMed
tbl_theta[,2] <- theta.pLb
tbl_theta[,3] <- theta.pUb
##
beta.names <- unique(c(x$chain1$covNames1, x$chain1$covNames2, x$chain1$covNames3))
nP <- length(beta.names)
output <- matrix(NA, nrow=nP, ncol=9)
dimnames(output) <- list(beta.names, c("exp(beta1)", "LL", "UL", "exp(beta2)", "LL", "UL", "exp(beta3)", "LL", "UL"))
if(length(x$chain1$beta1.p) != 0){
#beta1
p1 = dim(x$chain1$beta1.p)[2]
beta.p <- x$chain1$beta1.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
beta.p <- rbind(beta.p, x[[nam]]$beta1.p)
}
}
beta.pMed <- apply(exp(beta.p), 2, median)
beta.pSd <- apply(exp(beta.p), 2, sd)
beta.pUb <- apply(exp(beta.p), 2, quantile, prob = (1-conf.level/2))
beta.pLb <- apply(exp(beta.p), 2, quantile, prob = conf.level/2)
tbl1 <- matrix(NA, p1, 3)
rownames(tbl1) <- x$chain1$covNames1
tbl1[,1] <- beta.pMed
tbl1[,2] <- beta.pLb
tbl1[,3] <- beta.pUb
for(i in 1:nP)
{
for(k in 1:p1) if(x$chain1$covNames1[k] == beta.names[i]) output[i,1:3] <- tbl1[k,]
}
}
if(length(x$chain1$beta2.p) != 0){
#beta2
p2 = dim(x$chain1$beta2.p)[2]
beta.p <- x$chain1$beta2.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
beta.p <- rbind(beta.p, x[[nam]]$beta2.p)
}
}
beta.pMed <- apply(exp(beta.p), 2, median)
beta.pSd <- apply(exp(beta.p), 2, sd)
beta.pUb <- apply(exp(beta.p), 2, quantile, prob = (1-conf.level/2))
beta.pLb <- apply(exp(beta.p), 2, quantile, prob = conf.level/2)
tbl2 <- matrix(NA, p2, 3)
rownames(tbl2) <- x$chain1$covNames2
tbl2[,1] <- beta.pMed
tbl2[,2] <- beta.pLb
tbl2[,3] <- beta.pUb
for(i in 1:nP)
{
for(k in 1:p2) if(x$chain1$covNames2[k] == beta.names[i]) output[i,4:6] <- tbl2[k,]
}
}
if(length(x$chain1$beta3.p) != 0){
#beta3
p3 = dim(x$chain1$beta3.p)[2]
beta.p <- x$chain1$beta3.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
beta.p <- rbind(beta.p, x[[nam]]$beta3.p)
}
}
beta.pMed <- apply(exp(beta.p), 2, median)
beta.pSd <- apply(exp(beta.p), 2, sd)
beta.pUb <- apply(exp(beta.p), 2, quantile, prob = (1-conf.level/2))
beta.pLb <- apply(exp(beta.p), 2, quantile, prob = conf.level/2)
tbl3 <- matrix(NA, p3, 3)
rownames(tbl3) <- x$chain1$covNames3
tbl3[,1] <- beta.pMed
tbl3[,2] <- beta.pLb
tbl3[,3] <- beta.pUb
for(i in 1:nP)
{
for(k in 1:p3) if(x$chain1$covNames3[k] == beta.names[i]) output[i,7:9] <- tbl3[k,]
}
}
output.coef <- NULL
if(nP > 0)
{
output.coef <- output
}
if(x$class[4] == "LN")
{
##
sigSq.p <- x$chain1$sigSq1.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
sigSq.p <- rbind(sigSq.p, x[[nam]]$sigSq1.p)
}
}
sigSq.pMed <- apply(sigSq.p, 2, median)
sigSq.pUb <- apply(sigSq.p, 2, quantile, prob = (1-conf.level/2))
sigSq.pLb <- apply(sigSq.p, 2, quantile, prob = conf.level/2)
tbl_s1 <- c(sigSq.pMed,sigSq.pLb, sigSq.pUb)
##
sigSq.p <- x$chain1$sigSq2.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
sigSq.p <- rbind(sigSq.p, x[[nam]]$sigSq2.p)
}
}
sigSq.pMed <- apply(sigSq.p, 2, median)
sigSq.pUb <- apply(sigSq.p, 2, quantile, prob = (1-conf.level/2))
sigSq.pLb <- apply(sigSq.p, 2, quantile, prob = conf.level/2)
tbl_s2 <- c(sigSq.pMed,sigSq.pLb, sigSq.pUb)
##
sigSq.p <- x$chain1$sigSq3.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
sigSq.p <- rbind(sigSq.p, x[[nam]]$sigSq3.p)
}
}
sigSq.pMed <- apply(sigSq.p, 2, median)
sigSq.pUb <- apply(sigSq.p, 2, quantile, prob = (1-conf.level/2))
sigSq.pLb <- apply(sigSq.p, 2, quantile, prob = conf.level/2)
tbl_s3 <- c(sigSq.pMed,sigSq.pLb, sigSq.pUb)
##
mu.p <- x$chain1$mu1.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
mu.p <- rbind(mu.p, x[[nam]]$mu1.p)
}
}
mu.pMed <- apply(mu.p, 2, median)
mu.pUb <- apply(mu.p, 2, quantile, prob = (1-conf.level/2))
mu.pLb <- apply(mu.p, 2, quantile, prob = conf.level/2)
tbl_b1 <- c(mu.pMed,mu.pLb, mu.pUb)
##
mu.p <- x$chain1$mu2.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
mu.p <- rbind(mu.p, x[[nam]]$mu2.p)
}
}
mu.pMed <- apply(mu.p, 2, median)
mu.pUb <- apply(mu.p, 2, quantile, prob = (1-conf.level/2))
mu.pLb <- apply(mu.p, 2, quantile, prob = conf.level/2)
tbl_b2 <- c(mu.pMed,mu.pLb, mu.pUb)
##
mu.p <- x$chain1$mu3.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
mu.p <- rbind(mu.p, x[[nam]]$mu3.p)
}
}
mu.pMed <- apply(mu.p, 2, median)
mu.pUb <- apply(mu.p, 2, quantile, prob = (1-conf.level/2))
mu.pLb <- apply(mu.p, 2, quantile, prob = conf.level/2)
tbl_b3 <- c(mu.pMed,mu.pLb, mu.pUb)
bh <- matrix(c(tbl_b1, tbl_b2, tbl_b3, tbl_s1, tbl_s2, tbl_s3), 2, 9, byrow = T)
dimnames(bh) <- list(c("log-Normal: mu", "log-Normal: sigmaSq"), c("g=1: PM", "LL", "UL", "g=2: PM", "LL", "UL", "g=3: PM", "LL", "UL"))
value <- list(classFit=x$class, psrf=psrf, theta=tbl_theta, coef=output.coef, h0=bh, DIC=DIC, LPML=LPML)
}
if(x$class[4] == "DPM")
{
##
tau.p <- x$chain1$tau1.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
tau.p <- rbind(tau.p, x[[nam]]$tau1.p)
}
}
tau.pMed <- apply(tau.p, 2, median)
tau.pUb <- apply(tau.p, 2, quantile, prob = (1-conf.level/2))
tau.pLb <- apply(tau.p, 2, quantile, prob = conf.level/2)
tbl_t1 <- c(tau.pMed, tau.pLb, tau.pUb)
##
tau.p <- x$chain1$tau2.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
tau.p <- rbind(tau.p, x[[nam]]$tau2.p)
}
}
tau.pMed <- apply(tau.p, 2, median)
tau.pUb <- apply(tau.p, 2, quantile, prob = (1-conf.level/2))
tau.pLb <- apply(tau.p, 2, quantile, prob = conf.level/2)
tbl_t2 <- c(tau.pMed, tau.pLb, tau.pUb)
##
tau.p <- x$chain1$tau3.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
tau.p <- rbind(tau.p, x[[nam]]$tau3.p)
}
}
tau.pMed <- apply(tau.p, 2, median)
tau.pUb <- apply(tau.p, 2, quantile, prob = (1-conf.level/2))
tau.pLb <- apply(tau.p, 2, quantile, prob = conf.level/2)
tbl_t3 <- c(tau.pMed, tau.pLb, tau.pUb)
bh <- matrix(c(tbl_t1, tbl_t2, tbl_t3), 1, 9, byrow = T)
dimnames(bh) <- list(c("tau"), c("g=1: PM", "LL", "UL", "g=2: PM", "LL", "UL", "g=3: PM", "LL", "UL"))
value <- list(classFit=x$class, psrf=psrf, theta=tbl_theta, coef=output.coef, h0=bh, DIC=DIC, LPML=LPML)
}
}
if(x$class[2] == "Surv")
{
##
beta.names <- c(x$chain1$covNames)
nP <- length(beta.names)
output <- matrix(NA, nrow=nP, ncol=3)
dimnames(output) <- list(beta.names, c("exp(beta)", "LL", "UL"))
if(length(x$chain1$beta.p) != 0){
#beta
p = dim(x$chain1$beta.p)[2]
beta.p <- x$chain1$beta.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
beta.p <- rbind(beta.p, x[[nam]]$beta.p)
}
}
beta.pMed <- apply(exp(beta.p), 2, median)
beta.pSd <- apply(exp(beta.p), 2, sd)
beta.pUb <- apply(exp(beta.p), 2, quantile, prob = (1-conf.level/2))
beta.pLb <- apply(exp(beta.p), 2, quantile, prob = conf.level/2)
tbl <- matrix(NA, p, 3)
rownames(tbl) <- x$chain1$covNames
tbl[,1] <- beta.pMed
tbl[,2] <- beta.pLb
tbl[,3] <- beta.pUb
for(i in 1:nP)
{
for(k in 1:p) if(x$chain1$covNames[k] == beta.names[i]) output[i,1:3] <- tbl[k,]
}
}
output.coef <- NULL
if(nP > 0)
{
output.coef <- output
}
if(x$class[4] == "LN")
{
##
sigSq.p <- x$chain1$sigSq.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
sigSq.p <- rbind(sigSq.p, x[[nam]]$sigSq.p)
}
}
sigSq.pMed <- apply(sigSq.p, 2, median)
sigSq.pUb <- apply(sigSq.p, 2, quantile, prob = (1-conf.level/2))
sigSq.pLb <- apply(sigSq.p, 2, quantile, prob = conf.level/2)
tbl_s <- c(sigSq.pMed,sigSq.pLb, sigSq.pUb)
##
mu.p <- x$chain1$mu.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
mu.p <- rbind(mu.p, x[[nam]]$mu.p)
}
}
mu.pMed <- apply(mu.p, 2, median)
mu.pUb <- apply(mu.p, 2, quantile, prob = (1-conf.level/2))
mu.pLb <- apply(mu.p, 2, quantile, prob = conf.level/2)
tbl_b <- c(mu.pMed,mu.pLb, mu.pUb)
bh <- matrix(c(tbl_b, tbl_s), 2, 3, byrow = T)
dimnames(bh) <- list(c("log-Normal: mu", "log-Normal: sigmaSq"), c("PM", "LL", "UL"))
value <- list(coef=output.coef, h0=bh, psrf=psrf, classFit=x$class)
}
if(x$class[4] == "DPM")
{
##
tau.p <- x$chain1$tau.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
tau.p <- rbind(tau.p, x[[nam]]$tau.p)
}
}
tau.pMed <- apply(tau.p, 2, median)
tau.pUb <- apply(tau.p, 2, quantile, prob = (1-conf.level/2))
tau.pLb <- apply(tau.p, 2, quantile, prob = conf.level/2)
tbl_t <- c(tau.pMed, tau.pLb, tau.pUb)
bh <- matrix(c(tbl_t), 1, 3, byrow = T)
dimnames(bh) <- list(c("tau"), c("PM", "LL", "UL"))
value <- list(coef=output.coef, h0=bh, psrf=psrf, classFit=x$class)
}
}
value$setup <- x$setup
value$conf.level <- conf.level
class(value) <- "summ.Bayes_AFT"
return(value)
}
####
## PRINT.SUMMARY METHOD
####
##
print.summ.Freq_HReg <- function(x, digits=3, ...)
{
obj <- x
##
if(obj$class[2] == "Surv")
{
##
cat("\nAnalysis of independent univariate time-to-event data \n")
}
if(obj$class[2] == "ID")
{
##
cat("\nAnalysis of independent semi-competing risks data \n")
cat(obj$class[5], "assumption for h3\n")
}
cat("Confidence level: ", x$conf.level, "\n", sep = "")
##
if(!is.null(obj$coef))
{
##
cat("\nHazard ratios:\n")
print(round(exp(obj$coef), digits=digits))
}
##
if(obj$class[2] == "ID"){
cat("\nVariance of frailties:\n")
print(round(obj$theta, digits=digits))
}
##
cat("\nBaseline hazard function components:\n")
print(round(obj$h0, digits=digits))
##
invisible()
}
print.summ.Bayes_HReg <- function(x, digits=3, ...)
{
nChain = x$setup$nChain
if(x$classFit[2] == "ID")
{
if(x$classFit[3] == "Cor")
{
##
cat("\nAnalysis of cluster-correlated semi-competing risks data \n")
}
if(x$classFit[3] == "Ind")
{
##
cat("\nAnalysis of independent semi-competing risks data \n")
}
##
cat(x$setup$model, "assumption for h3\n")
##
cat("\n#####\n")
cat("\nDIC: ", round(x$DIC))
cat("\nLPML: ", round(x$LPML), "\n")
}
if(x$classFit[2] == "Surv")
{
if(x$classFit[3] == "Cor")
{
##
cat("\nAnalysis of cluster-correlated univariate time-to-event data \n")
}
if(x$classFit[3] == "Ind")
{
##
cat("\nAnalysis of independent univariate time-to-event data \n")
}
}
cat("Credibility level: ", x$conf.level, "\n", sep = "")
cat("\n#####\n")
if(!is.null(x$coef))
{
##
cat("\nHazard ratios:\n")
print(round(x$coef, digits=digits))
}
if(x$classFit[2] == "ID")
{
##
cat("\nVariance of frailties:\n")
print(round(x$theta, digits=digits))
}
##
cat("\nBaseline hazard function components:\n")
print(round(x$h0, digits=digits))
if(x$classFit[3] == "Cor")
{
if(x$classFit[5] == "DPM")
{
##
cat("\nPrecision parameter of DPM prior:\n")
print(round(x$tau, digits=digits))
}
if(x$classFit[2] == "ID")
{
##
cat("\nVariance-covariance matrix of cluster-specific random effects:\n")
print(round(x$Sigma.PM, digits=digits))
}
if(x$classFit[2] == "Surv")
{
##
cat("\nVariance of cluster-specific random effects:\n")
print(round(x$sigma_V, digits=digits))
}
}
invisible()
}
print.summ.Bayes_AFT <- function(x, digits=3, ...)
{
nChain = x$setup$nChain
if(x$classFit[2] == "ID")
{
if(x$classFit[3] == "Cor")
{
##
cat("\nAnalysis of cluster-correlated semi-competing risks data \n")
}
if(x$classFit[3] == "Ind")
{
##
cat("\nAnalysis of independent semi-competing risks data \n")
}
##
cat("\n#####\n")
cat("\nDIC: ", round(x$DIC))
cat("\nLPML: ", round(x$LPML), "\n")
}
if(x$classFit[2] == "Surv")
{
if(x$classFit[3] == "Cor")
{
##
cat("\nAnalysis of cluster-correlated univariate time-to-event data \n")
}
if(x$classFit[3] == "Ind")
{
##
cat("\nAnalysis of independent univariate time-to-event data \n")
}
}
cat("Credibility level: ", x$conf.level, "\n", sep = "")
cat("\n#####\n")
if(!is.null(x$coef))
{
##
cat("\nAcceleration factors:\n")
print(round(x$coef, digits=digits))
}
if(x$classFit[2] == "ID")
{
##
cat("\nVariance of frailties:\n")
print(round(x$theta, digits=digits))
}
##
cat("\nBaseline survival function components:\n")
print(round(x$h0, digits=digits))
invisible()
}
####
## COEF METHOD
####
##
coef.Bayes_HReg <- function(object, alpha=0.05, ...)
{
conf.level = alpha
x <- object
nChain = x$setup$nChain
# estimates
if(x$class[2] == "ID")
{
##
beta.names <- unique(c(x$chain1$covNames1, x$chain1$covNames2, x$chain1$covNames3))
nP <- length(beta.names)
output <- matrix(NA, nrow=nP, ncol=9)
dimnames(output) <- list(beta.names, c("beta1", "LL", "UL", "beta2", "LL", "UL", "beta3", "LL", "UL"))
if(length(x$chain1$beta1.p) != 0){
#beta1
p1 = dim(x$chain1$beta1.p)[2]
beta.p <- x$chain1$beta1.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
beta.p <- rbind(beta.p, x[[nam]]$beta1.p)
}
}
beta.pMed <- apply(beta.p, 2, median)
beta.pSd <- apply(beta.p, 2, sd)
beta.pUb <- apply(beta.p, 2, quantile, prob = (1-conf.level/2))
beta.pLb <- apply(beta.p, 2, quantile, prob = conf.level/2)
tbl1 <- matrix(NA, p1, 3)
rownames(tbl1) <- x$chain1$covNames1
tbl1[,1] <- beta.pMed
tbl1[,2] <- beta.pLb
tbl1[,3] <- beta.pUb
for(i in 1:nP)
{
for(k in 1:p1) if(x$chain1$covNames1[k] == beta.names[i]) output[i,1:3] <- tbl1[k,]
}
}
if(length(x$chain1$beta2.p) != 0){
#beta2
p2 = dim(x$chain1$beta2.p)[2]
beta.p <- x$chain1$beta2.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
beta.p <- rbind(beta.p, x[[nam]]$beta2.p)
}
}
beta.pMed <- apply(beta.p, 2, median)
beta.pSd <- apply(beta.p, 2, sd)
beta.pUb <- apply(beta.p, 2, quantile, prob = (1-conf.level/2))
beta.pLb <- apply(beta.p, 2, quantile, prob = conf.level/2)
tbl2 <- matrix(NA, p2, 3)
rownames(tbl2) <- x$chain1$covNames2
tbl2[,1] <- beta.pMed
tbl2[,2] <- beta.pLb
tbl2[,3] <- beta.pUb
for(i in 1:nP)
{
for(k in 1:p2) if(x$chain1$covNames2[k] == beta.names[i]) output[i,4:6] <- tbl2[k,]
}
}
if(length(x$chain1$beta3.p) != 0){
#beta3
p3 = dim(x$chain1$beta3.p)[2]
beta.p <- x$chain1$beta3.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
beta.p <- rbind(beta.p, x[[nam]]$beta3.p)
}
}
beta.pMed <- apply(beta.p, 2, median)
beta.pSd <- apply(beta.p, 2, sd)
beta.pUb <- apply(beta.p, 2, quantile, prob = (1-conf.level/2))
beta.pLb <- apply(beta.p, 2, quantile, prob = conf.level/2)
tbl3 <- matrix(NA, p3, 3)
rownames(tbl3) <- x$chain1$covNames3
tbl3[,1] <- beta.pMed
tbl3[,2] <- beta.pLb
tbl3[,3] <- beta.pUb
for(i in 1:nP)
{
for(k in 1:p3) if(x$chain1$covNames3[k] == beta.names[i]) output[i,7:9] <- tbl3[k,]
}
}
output.coef <- NULL
if(nP > 0)
{
output.coef <- output
}
}
if(x$class[2] == "Surv")
{
##
beta.names <- c(x$chain1$covNames)
nP <- length(beta.names)
output <- matrix(NA, nrow=nP, ncol=3)
dimnames(output) <- list(beta.names, c("exp(beta)", "LL", "UL"))
if(length(x$chain1$beta.p) != 0){
#beta
p = dim(x$chain1$beta.p)[2]
beta.p <- x$chain1$beta.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
beta.p <- rbind(beta.p, x[[nam]]$beta.p)
}
}
beta.pMed <- apply(beta.p, 2, median)
beta.pSd <- apply(beta.p, 2, sd)
beta.pUb <- apply(beta.p, 2, quantile, prob = (1-conf.level/2))
beta.pLb <- apply(beta.p, 2, quantile, prob = conf.level/2)
tbl <- matrix(NA, p, 3)
rownames(tbl) <- x$chain1$covNames
tbl[,1] <- beta.pMed
tbl[,2] <- beta.pLb
tbl[,3] <- beta.pUb
for(i in 1:nP)
{
for(k in 1:p) if(x$chain1$covNames[k] == beta.names[i]) output[i,1:3] <- tbl[k,]
}
}
output.coef <- NULL
if(nP > 0)
{
output.coef <- output
}
}
value <- output.coef
return(value)
}
coef.Bayes_AFT <- function(object, alpha=0.05, ...)
{
conf.level = alpha
x <- object
nChain = x$setup$nChain
# estimates
if(x$class[2] == "ID")
{
##
beta.names <- unique(c(x$chain1$covNames1, x$chain1$covNames2, x$chain1$covNames3))
nP <- length(beta.names)
output <- matrix(NA, nrow=nP, ncol=9)
dimnames(output) <- list(beta.names, c("beta1", "LL", "UL", "beta2", "LL", "UL", "beta3", "LL", "UL"))
if(length(x$chain1$beta1.p) != 0){
#beta1
p1 = dim(x$chain1$beta1.p)[2]
beta.p <- x$chain1$beta1.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
beta.p <- rbind(beta.p, x[[nam]]$beta1.p)
}
}
beta.pMed <- apply(beta.p, 2, median)
beta.pSd <- apply(beta.p, 2, sd)
beta.pUb <- apply(beta.p, 2, quantile, prob = (1-conf.level/2))
beta.pLb <- apply(beta.p, 2, quantile, prob = conf.level/2)
tbl1 <- matrix(NA, p1, 3)
rownames(tbl1) <- x$chain1$covNames1
tbl1[,1] <- beta.pMed
tbl1[,2] <- beta.pLb
tbl1[,3] <- beta.pUb
for(i in 1:nP)
{
for(k in 1:p1) if(x$chain1$covNames1[k] == beta.names[i]) output[i,1:3] <- tbl1[k,]
}
}
if(length(x$chain1$beta2.p) != 0){
#beta2
p2 = dim(x$chain1$beta2.p)[2]
beta.p <- x$chain1$beta2.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
beta.p <- rbind(beta.p, x[[nam]]$beta2.p)
}
}
beta.pMed <- apply(beta.p, 2, median)
beta.pSd <- apply(beta.p, 2, sd)
beta.pUb <- apply(beta.p, 2, quantile, prob = (1-conf.level/2))
beta.pLb <- apply(beta.p, 2, quantile, prob = conf.level/2)
tbl2 <- matrix(NA, p2, 3)
rownames(tbl2) <- x$chain1$covNames2
tbl2[,1] <- beta.pMed
tbl2[,2] <- beta.pLb
tbl2[,3] <- beta.pUb
for(i in 1:nP)
{
for(k in 1:p2) if(x$chain1$covNames2[k] == beta.names[i]) output[i,4:6] <- tbl2[k,]
}
}
if(length(x$chain1$beta3.p) != 0){
#beta3
p3 = dim(x$chain1$beta3.p)[2]
beta.p <- x$chain1$beta3.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
beta.p <- rbind(beta.p, x[[nam]]$beta3.p)
}
}
beta.pMed <- apply(beta.p, 2, median)
beta.pSd <- apply(beta.p, 2, sd)
beta.pUb <- apply(beta.p, 2, quantile, prob = (1-conf.level/2))
beta.pLb <- apply(beta.p, 2, quantile, prob = conf.level/2)
tbl3 <- matrix(NA, p3, 3)
rownames(tbl3) <- x$chain1$covNames3
tbl3[,1] <- beta.pMed
tbl3[,2] <- beta.pLb
tbl3[,3] <- beta.pUb
for(i in 1:nP)
{
for(k in 1:p3) if(x$chain1$covNames3[k] == beta.names[i]) output[i,7:9] <- tbl3[k,]
}
}
output.coef <- NULL
if(nP > 0)
{
output.coef <- output
}
}
if(x$class[2] == "Surv")
{
##
beta.names <- c(x$chain1$covNames)
nP <- length(beta.names)
output <- matrix(NA, nrow=nP, ncol=3)
dimnames(output) <- list(beta.names, c("beta", "LL", "UL"))
if(length(x$chain1$beta.p) != 0){
#beta
p = dim(x$chain1$beta.p)[2]
beta.p <- x$chain1$beta.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
beta.p <- rbind(beta.p, x[[nam]]$beta.p)
}
}
beta.pMed <- apply(beta.p, 2, median)
beta.pSd <- apply(beta.p, 2, sd)
beta.pUb <- apply(beta.p, 2, quantile, prob = (1-conf.level/2))
beta.pLb <- apply(beta.p, 2, quantile, prob = conf.level/2)
tbl <- matrix(NA, p, 3)
rownames(tbl) <- x$chain1$covNames
tbl[,1] <- beta.pMed
tbl[,2] <- beta.pLb
tbl[,3] <- beta.pUb
for(i in 1:nP)
{
for(k in 1:p) if(x$chain1$covNames[k] == beta.names[i]) output[i,1:3] <- tbl[k,]
}
}
output.coef <- NULL
if(nP > 0)
{
output.coef <- output
}
}
value <- output.coef
return(value)
}
coef.Freq_HReg <- function(object, alpha=0.05, ...)
{
conf.level = alpha
obj <- object
##
logEst <- obj$estimate
logSE <- sqrt(diag(obj$Finv))
results <- cbind(logEst, logEst - abs(qnorm(conf.level/2, 0, 1))*logSE, logEst + abs(qnorm(conf.level/2, 0, 1))*logSE)
##
if(obj$class[2] == "Surv")
{
##
#cat("\nRegression coefficients:\n")
if(length(obj$myLabels) >= 3)
{
output.coef <- results[-c(1:2),]
dimnames(output.coef) <- list(unique(obj$myLabels[-c(1:2)]), c("beta", "LL", "UL"))
}else
{
output.coef <- NULL
}
}
##
if(obj$class[2] == "ID")
{
##
nP.0 <- ifelse(obj$frailty, 7, 6)
nP.1 <- obj$nP[1]
nP.2 <- obj$nP[2]
nP.3 <- obj$nP[3]
##
beta.names <- unique(obj$myLabels[-c(1:nP.0)])
nP <- length(beta.names)
##
#cat("\nRegression coefficients:\n")
output <- matrix(NA, nrow=nP, ncol=9)
dimnames(output) <- list(beta.names, c("beta1", "LL", "UL", "beta2", "LL", "UL", "beta3", "LL", "UL"))
for(i in 1:nP)
{
if(nP.1 != 0)
{
for(j in 1:nP.1) if(obj$myLabels[nP.0+j] == beta.names[i]) output[i,1:3] <- results[nP.0+j,]
}
if(nP.2 != 0)
{
for(j in 1:nP.2) if(obj$myLabels[nP.0+nP.1+j] == beta.names[i]) output[i,4:6] <- results[nP.0+nP.1+j,]
}
if(nP.3 != 0)
{
for(j in 1:nP.3) if(obj$myLabels[nP.0+nP.1+nP.2+j] == beta.names[i]) output[i,7:9] <- results[nP.0+nP.1+nP.2+j,]
}
}
output.coef <- output
}
value <- output.coef
return(value)
}
####
## predict METHOD
####
##
predict.Bayes_HReg <- function(object, xnew=NULL, x1new=NULL, x2new=NULL, x3new=NULL, tseq=c(0, 5, 10), alpha = 0.05, ...)
{
conf.level = alpha
x <- object
nChain = x$setup$nChain
if(x$class[2] == "ID")
{
if(!is.null(xnew))
{
stop("'xnew' is for univariate models so it must be specified as NULL for semi-competing risks models")
}
if(is.null(x$chain1$beta1.p))
{
if(!is.null(x1new))
{
stop("x1new does not have the same length as the x1 specified in the Formula.")
}else
{
LP1 <- 0
}
}else if(length(x$chain1$beta1.p) != 0)
{
if(!is.null(x1new))
{
#beta1
p1 = dim(x$chain1$beta1.p)[2]
beta1.p <- x$chain1$beta1.p
if(nChain > 1)
{
for(i in 2:nChain)
{
nam <- paste("chain", i, sep="")
beta1.p <- rbind(beta1.p, x[[nam]]$beta1.p)
}
}
LP1 <- rowSums(beta1.p * matrix(x1new, nrow=dim(beta1.p)[1], ncol = dim(beta1.p)[2], byrow=T))
}else
{
LP1 <- 0
}
}
if(is.null(x$chain1$beta2.p))
{
if(!is.null(x2new))
{
stop("x2new does not have the same length as the x2 specified in the Formula.")
}else
{
LP2 <- 0
}
}else if(length(x$chain1$beta2.p) != 0){
if(!is.null(x2new))
{
#beta2
p2 = dim(x$chain1$beta2.p)[2]
beta2.p <- x$chain1$beta2.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
beta2.p <- rbind(beta2.p, x[[nam]]$beta2.p)
}
}
LP2 <- rowSums(beta2.p * matrix(x2new, nrow=dim(beta2.p)[1], ncol = dim(beta2.p)[2], byrow=T))
}else
{
LP2 <- 0
}
}
if(is.null(x$chain1$beta3.p))
{
if(!is.null(x3new))
{
stop("x3new does not have the same length as the x3 specified in the Formula.")
}else
{
LP3 <- 0
}
}else if(length(x$chain1$beta3.p) != 0){
if(!is.null(x3new))
{
#beta3
p3 = dim(x$chain1$beta3.p)[2]
beta3.p <- x$chain1$beta3.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
beta3.p <- rbind(beta3.p, x[[nam]]$beta3.p)
}
}
LP3 <- rowSums(beta3.p * matrix(x3new, nrow=dim(beta3.p)[1], ncol = dim(beta3.p)[2], byrow=T))
}else
{
LP3 <- 0
}
}
if(x$class[4] == "PEM")
{
time1 <- x$chain1$time_lambda1
time2 <- x$chain1$time_lambda2
time3 <- x$chain1$time_lambda3
time1hz <- time1
time2hz <- time2
time3hz <- time3
lambda1.fin <- x$chain1$lambda1.fin
lambda2.fin <- x$chain1$lambda2.fin
lambda3.fin <- x$chain1$lambda3.fin
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
lambda1.fin <- rbind(lambda1.fin, x[[nam]]$lambda1.fin)
lambda2.fin <- rbind(lambda2.fin, x[[nam]]$lambda2.fin)
lambda3.fin <- rbind(lambda3.fin, x[[nam]]$lambda3.fin)
}
}
BH1Med <- apply(exp(lambda1.fin)*exp(matrix(LP1, nrow=dim(lambda1.fin)[1], ncol=dim(lambda1.fin)[2], byrow=F)), 2, median)
BH1Ub <- apply(exp(lambda1.fin)*exp(matrix(LP1, nrow=dim(lambda1.fin)[1], ncol=dim(lambda1.fin)[2], byrow=F)), 2, quantile, prob = (1-conf.level/2))
BH1Lb <- apply(exp(lambda1.fin)*exp(matrix(LP1, nrow=dim(lambda1.fin)[1], ncol=dim(lambda1.fin)[2], byrow=F)), 2, quantile, prob = conf.level/2)
BH2Med <- apply(exp(lambda2.fin)*exp(matrix(LP2, nrow=dim(lambda2.fin)[1], ncol=dim(lambda2.fin)[2], byrow=F)), 2, median)
BH2Ub <- apply(exp(lambda2.fin)*exp(matrix(LP2, nrow=dim(lambda2.fin)[1], ncol=dim(lambda2.fin)[2], byrow=F)), 2, quantile, prob = (1-conf.level/2))
BH2Lb <- apply(exp(lambda2.fin)*exp(matrix(LP2, nrow=dim(lambda2.fin)[1], ncol=dim(lambda2.fin)[2], byrow=F)), 2, quantile, prob = conf.level/2)
BH3Med <- apply(exp(lambda3.fin)*exp(matrix(LP3, nrow=dim(lambda3.fin)[1], ncol=dim(lambda3.fin)[2], byrow=F)), 2, median)
BH3Ub <- apply(exp(lambda3.fin)*exp(matrix(LP3, nrow=dim(lambda3.fin)[1], ncol=dim(lambda3.fin)[2], byrow=F)), 2, quantile, prob = (1-conf.level/2))
BH3Lb <- apply(exp(lambda3.fin)*exp(matrix(LP3, nrow=dim(lambda3.fin)[1], ncol=dim(lambda3.fin)[2], byrow=F)), 2, quantile, prob = conf.level/2)
dif1 <- diff(c(0, time1hz))
dif2 <- diff(c(0, time2hz))
dif3 <- diff(c(0, time3hz))
BS1 <- matrix(NA, dim(lambda1.fin)[1], dim(lambda1.fin)[2])
for(i in 1:dim(lambda1.fin)[1])
{
BS1[i,] <- cumsum(exp(lambda1.fin[i,])* dif1)
}
BS2 <- matrix(NA, dim(lambda2.fin)[1], dim(lambda2.fin)[2])
for(i in 1:dim(lambda2.fin)[1])
{
BS2[i,] <- cumsum(exp(lambda2.fin[i,])* dif2)
}
BS3 <- matrix(NA, dim(lambda3.fin)[1], dim(lambda3.fin)[2])
for(i in 1:dim(lambda3.fin)[1])
{
BS3[i,] <- cumsum(exp(lambda3.fin[i,])* dif3)
}
BS1 <- exp(-BS1 * exp(matrix(LP1, nrow=dim(lambda1.fin)[1], ncol=dim(lambda1.fin)[2], byrow=F)))
BS2 <- exp(-BS2 * exp(matrix(LP2, nrow=dim(lambda2.fin)[1], ncol=dim(lambda2.fin)[2], byrow=F)))
BS3 <- exp(-BS3 * exp(matrix(LP3, nrow=dim(lambda3.fin)[1], ncol=dim(lambda3.fin)[2], byrow=F)))
BS1Med <- apply(BS1, 2, median)
BS1Ub <- apply(BS1, 2, quantile, prob = (1-conf.level/2))
BS1Lb <- apply(BS1, 2, quantile, prob = conf.level/2)
BS2Med <- apply(BS2, 2, median)
BS2Ub <- apply(BS2, 2, quantile, prob = (1-conf.level/2))
BS2Lb <- apply(BS2, 2, quantile, prob = conf.level/2)
BS3Med <- apply(BS3, 2, median)
BS3Ub <- apply(BS3, 2, quantile, prob = (1-conf.level/2))
BS3Lb <- apply(BS3, 2, quantile, prob = conf.level/2)
}
if(x$class[4] == "WB")
{
time1 <- time2 <- time3 <- seq(from=min(tseq), to=max(tseq), length=100)
nStore <- length(x$chain1$alpha1.p)
numSpl <- nStore * nChain
basehaz1 <- matrix(NA, numSpl, length(time1))
basehaz2 <- matrix(NA, numSpl, length(time2))
basehaz3 <- matrix(NA, numSpl, length(time3))
basesurv1 <- matrix(NA, numSpl, length(time1))
basesurv2 <- matrix(NA, numSpl, length(time2))
basesurv3 <- matrix(NA, numSpl, length(time3))
alpha1.p <- x$chain1$alpha1.p
alpha2.p <- x$chain1$alpha2.p
alpha3.p <- x$chain1$alpha3.p
kappa1.p <- x$chain1$kappa1.p
kappa2.p <- x$chain1$kappa2.p
kappa3.p <- x$chain1$kappa3.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
alpha1.p <- c(alpha1.p, x[[nam]]$alpha1.p)
alpha2.p <- c(alpha2.p, x[[nam]]$alpha2.p)
alpha3.p <- c(alpha3.p, x[[nam]]$alpha3.p)
kappa1.p <- c(kappa1.p, x[[nam]]$kappa1.p)
kappa2.p <- c(kappa2.p, x[[nam]]$kappa2.p)
kappa3.p <- c(kappa3.p, x[[nam]]$kappa3.p)
}
}
for(i in 1:numSpl){
basehaz1[i, ] <- alpha1.p[i] * kappa1.p[i] * time1^(alpha1.p[i] - 1)
basehaz2[i, ] <- alpha2.p[i] * kappa2.p[i] * time2^(alpha2.p[i] - 1)
basehaz3[i, ] <- alpha3.p[i] * kappa3.p[i] * time3^(alpha3.p[i] - 1)
basesurv1[i, ] <- kappa1.p[i] * time1^(alpha1.p[i])
basesurv2[i, ] <- kappa2.p[i] * time2^(alpha2.p[i])
basesurv3[i, ] <- kappa3.p[i] * time3^(alpha3.p[i])
}
basehaz1 <- basehaz1 * exp(matrix(LP1, nrow=dim(basehaz1)[1], ncol=dim(basehaz1)[2], byrow=F))
basehaz2 <- basehaz2 * exp(matrix(LP2, nrow=dim(basehaz2)[1], ncol=dim(basehaz2)[2], byrow=F))
basehaz3 <- basehaz3 * exp(matrix(LP3, nrow=dim(basehaz3)[1], ncol=dim(basehaz3)[2], byrow=F))
basesurv1 <- exp(-basesurv1 * exp(matrix(LP1, nrow=dim(basehaz1)[1], ncol=dim(basehaz1)[2], byrow=F)))
basesurv2 <- exp(-basesurv2 * exp(matrix(LP2, nrow=dim(basehaz2)[1], ncol=dim(basehaz2)[2], byrow=F)))
basesurv3 <- exp(-basesurv3 * exp(matrix(LP3, nrow=dim(basehaz3)[1], ncol=dim(basehaz3)[2], byrow=F)))
time1hz <- time1
time2hz <- time2
time3hz <- time3
if(tseq[1] == 0){
time1hz <- time1[-1]
time2hz <- time2[-1]
time3hz <- time3[-1]
basehaz1 <- basehaz1[,-1]
basehaz2 <- basehaz2[,-1]
basehaz3 <- basehaz3[,-1]
}
BH1Med <- apply(basehaz1, 2, median)
BH1Ub <- apply(basehaz1, 2, quantile, prob = (1-conf.level/2))
BH1Lb <- apply(basehaz1, 2, quantile, prob = conf.level/2)
BH2Med <- apply(basehaz2, 2, median)
BH2Ub <- apply(basehaz2, 2, quantile, prob = (1-conf.level/2))
BH2Lb <- apply(basehaz2, 2, quantile, prob = conf.level/2)
BH3Med <- apply(basehaz3, 2, median)
BH3Ub <- apply(basehaz3, 2, quantile, prob = (1-conf.level/2))
BH3Lb <- apply(basehaz3, 2, quantile, prob = conf.level/2)
BS1Med <- apply(basesurv1, 2, median)
BS1Ub <- apply(basesurv1, 2, quantile, prob = (1-conf.level/2))
BS1Lb <- apply(basesurv1, 2, quantile, prob = conf.level/2)
BS2Med <- apply(basesurv2, 2, median)
BS2Ub <- apply(basesurv2, 2, quantile, prob = (1-conf.level/2))
BS2Lb <- apply(basesurv2, 2, quantile, prob = conf.level/2)
BS3Med <- apply(basesurv3, 2, median)
BS3Ub <- apply(basesurv3, 2, quantile, prob = (1-conf.level/2))
BS3Lb <- apply(basesurv3, 2, quantile, prob = conf.level/2)
}
BH1_tbl <- data.frame(time1=time1hz, h.1=BH1Med, LL.1=BH1Lb, UL.1=BH1Ub)
BH2_tbl <- data.frame(time2=time2hz, h.2=BH2Med, LL.2=BH2Lb, UL.2=BH2Ub)
BH3_tbl <- data.frame(time3=time3hz, h.3=BH3Med, LL.3=BH3Lb, UL.3=BH3Ub)
BS1_tbl <- data.frame(time1=time1, S.1=BS1Med, LL.1=BS1Lb, UL.1=BS1Ub)
BS2_tbl <- data.frame(time2=time2, S.2=BS2Med, LL.2=BS2Lb, UL.2=BS2Ub)
BS3_tbl <- data.frame(time3=time3, S.3=BS3Med, LL.3=BS3Lb, UL.3=BS3Ub)
value <- list(h.1=BH1_tbl, h.2=BH2_tbl, h.3=BH3_tbl, S.1=BS1_tbl, S.2=BS2_tbl, S.3=BS3_tbl)
}
if(x$class[2] == "Surv")
{
if(!is.null(x1new))
{
stop("'x1new' is for semi-competing risks models so it must be specified as NULL for univariate models")
}
if(!is.null(x2new))
{
stop("'x2new' is for semi-competing risks models so it must be specified as NULL for univariate models")
}
if(!is.null(x3new))
{
stop("'x3new' is for semi-competing risks models so it must be specified as NULL for univariate models")
}
if(is.null(x$chain1$beta.p))
{
if(!is.null(xnew))
{
stop("xnew does not have the same length as the x specified in the Formula.")
}else
{
LP <- 0
}
}else if(length(x$chain1$beta.p) != 0)
{
if(!is.null(xnew))
{
#beta
p = dim(x$chain1$beta.p)[2]
beta.p <- x$chain1$beta.p
if(nChain > 1)
{
for(i in 2:nChain)
{
nam <- paste("chain", i, sep="")
beta.p <- rbind(beta.p, x[[nam]]$beta.p)
}
}
LP <- rowSums(beta.p * matrix(xnew, nrow=dim(beta.p)[1], ncol = dim(beta.p)[2], byrow=T))
}else
{
LP <- 0
}
}
if(x$class[4] == "PEM")
{
time <- x$chain1$time_lambda
timehz <- time
lambda.fin <- x$chain1$lambda.fin
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
lambda.fin <- rbind(lambda.fin, x[[nam]]$lambda.fin)
}
}
BHMed <- apply(exp(lambda.fin)*exp(matrix(LP, nrow=dim(lambda.fin)[1], ncol=dim(lambda.fin)[2], byrow=F)), 2, median)
BHUb <- apply(exp(lambda.fin)*exp(matrix(LP, nrow=dim(lambda.fin)[1], ncol=dim(lambda.fin)[2], byrow=F)), 2, quantile, prob = (1-conf.level/2))
BHLb <- apply(exp(lambda.fin)*exp(matrix(LP, nrow=dim(lambda.fin)[1], ncol=dim(lambda.fin)[2], byrow=F)), 2, quantile, prob = conf.level/2)
dif <- diff(c(0, timehz))
BS <- matrix(NA, dim(lambda.fin)[1], dim(lambda.fin)[2])
for(i in 1:dim(lambda.fin)[1])
{
BS[i,] <- cumsum(exp(lambda.fin[i,])* dif)
}
BS <- exp(-BS*exp(matrix(LP, nrow=dim(lambda.fin)[1], ncol=dim(lambda.fin)[2], byrow=F)))
BSMed <- apply(BS, 2, median)
BSUb <- apply(BS, 2, quantile, prob = (1-conf.level/2))
BSLb <- apply(BS, 2, quantile, prob = conf.level/2)
}
if(x$class[4] == "WB")
{
time <- seq(from=min(tseq), to=max(tseq), length=100)
nStore <- length(x$chain1$alpha.p)
numSpl <- nStore * nChain
basehaz <- matrix(NA, numSpl, length(time))
basesurv <- matrix(NA, numSpl, length(time))
alpha.p <- x$chain1$alpha.p
kappa.p <- x$chain1$kappa.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
alpha.p <- c(alpha.p, x[[nam]]$alpha.p)
kappa.p <- c(kappa.p, x[[nam]]$kappa.p)
}
}
for(i in 1:numSpl){
basehaz[i, ] <- alpha.p[i] * kappa.p[i] * time^(alpha.p[i] - 1)
basesurv[i, ] <- kappa.p[i] * time^(alpha.p[i])
}
basehaz <- basehaz * exp(matrix(LP, nrow=dim(basehaz)[1], ncol=dim(basehaz)[2], byrow=F))
basesurv <- exp(-basesurv * exp(matrix(LP, nrow=dim(basehaz)[1], ncol=dim(basehaz)[2], byrow=F)))
timehz <- time
if(tseq[1] == 0){
timehz <- time[-1]
basehaz <- basehaz[,-1]
}
BHMed <- apply(basehaz, 2, median)
BHUb <- apply(basehaz, 2, quantile, prob = (1-conf.level/2))
BHLb <- apply(basehaz, 2, quantile, prob = conf.level/2)
BSMed <- apply(basesurv, 2, median)
BSUb <- apply(basesurv, 2, quantile, prob = (1-conf.level/2))
BSLb <- apply(basesurv, 2, quantile, prob = conf.level/2)
}
BH_tbl <- data.frame(time=timehz, h=BHMed, LL=BHLb, UL=BHUb)
BS_tbl <- data.frame(time=time, S=BSMed, LL=BSLb, UL=BSUb)
value <- list(h=BH_tbl, S=BS_tbl)
}
value$xnew <- xnew
value$x1new <- x1new
value$x2new <- x2new
value$x3new <- x3new
value$tseq <- tseq
value$setup$model <- x$setup$model
value$class <- x$class
class(value) <- "pred.Bayes_HReg"
return(value)
}
plot.pred.Bayes_HReg <- function(x, plot.est="Haz", xlab=NULL, ylab=NULL, ...)
{
if(x$class[2] == "ID")
{
if(is.null(xlab))
{
xlab <- c("Time", "Time", "Time")
if(x$setup$model == "semi-Markov")
{
xlab[3] <- "Time since non-terminal event"
}
}
if(plot.est == "Haz")
{
if(is.null(ylab))
{
ylab <- "Hazard"
}
ygrid <- (max(x$h.1$UL.1, x$h.2$UL.2, x$h.3$UL.3) - 0)/5
ylim <- seq(from=0, to=max(x$h.1$UL.1, x$h.2$UL.2, x$h.3$UL.3), by=ygrid)
##
par(mfrow=c(1,3))
##
plot(c(0, max(x$h.1$time1)), range(ylim), xlab=xlab[1], ylab=ylab, type="n", main = expression(paste("Estimated ", h[1](t), "")), axes=FALSE)
if(x$class[4] == "PEM")
{
axis(1, at=c(0, max(x$h.1$time1)))
}
if(x$class[4] == "WB")
{
axis(1, at=x$tseq)
}
axis(2, at=round(ylim, 4))
lines(x$h.1$time1, x$h.1$h.1, col="blue", lwd=3)
lines(x$h.1$time1, x$h.1$UL.1, col="blue", lwd=3, lty=3)
lines(x$h.1$time1, x$h.1$LL.1, col="blue", lwd=3, lty=3)
##
plot(c(0, max(x$h.2$time2)), range(ylim), xlab=xlab[2], ylab=ylab, type="n", main = expression(paste("Estimated ", h[2](t), "")), axes=FALSE)
if(x$class[4] == "PEM")
{
axis(1, at=c(0, max(x$h.2$time2)))
}
if(x$class[4] == "WB")
{
axis(1, at=x$tseq)
}
axis(2, at=round(ylim, 4))
lines(x$h.2$time2, x$h.2$h.2, col="red", lwd=3)
lines(x$h.2$time2, x$h.2$UL.2, col="red", lwd=3, lty=3)
lines(x$h.2$time2, x$h.2$LL.2, col="red", lwd=3, lty=3)
##
plot(c(0, max(x$h.3$time3)), range(ylim), xlab=xlab[3], ylab=ylab, type="n", main = expression(paste("Estimated ", h[3](t), "")), axes=FALSE)
if(x$class[4] == "PEM")
{
axis(1, at=c(0, max(x$h.3$time3)))
}
if(x$class[4] == "WB")
{
axis(1, at=x$tseq)
}
axis(2, at=round(ylim, 4))
lines(x$h.3$time3, x$h.3$h.3, col="red", lwd=3)
lines(x$h.3$time3, x$h.3$UL.3, col="red", lwd=3, lty=3)
lines(x$h.3$time3, x$h.3$LL.3, col="red", lwd=3, lty=3)
}
if(plot.est == "Surv")
{
if(is.null(ylab))
{
ylab <- "Survival"
}
ylim <- seq(from=0, to=1, by=0.2)
##
par(mfrow=c(1,3))
##
plot(c(0, max(x$S.1$time1)), range(ylim), xlab=xlab[1], ylab=ylab, type="n", main = expression(paste("Estimated ", S[1](t), "")), axes=FALSE)
if(x$class[4] == "PEM")
{
axis(1, at=c(0, max(x$S.1$time1)))
}
if(x$class[4] == "WB")
{
axis(1, at=x$tseq)
}
axis(2, at=ylim)
if(x$S.1$time1[1] == 0)
{
lines(x$S.1$time1, x$S.1$S.1, col="blue", lwd=3)
lines(x$S.1$time1, x$S.1$UL.1, col="blue", lwd=3, lty=3)
lines(x$S.1$time1, x$S.1$LL.1, col="blue", lwd=3, lty=3)
}else
{
lines(unique(c(0, x$S.1$time1)), c(1, x$S.1$S.1), col="blue", lwd=3)
lines(unique(c(0, x$S.1$time1)), c(1, x$S.1$UL.1), col="blue", lwd=3, lty=3)
lines(unique(c(0, x$S.1$time1)), c(1, x$S.1$LL.1), col="blue", lwd=3, lty=3)
}
##
plot(c(0, max(x$S.2$time2)), range(ylim), xlab=xlab[2], ylab=ylab, type="n", main = expression(paste("Estimated ", S[2](t), "")), axes=FALSE)
if(x$class[4] == "PEM")
{
axis(1, at=c(0, max(x$S.2$time2)))
}
if(x$class[4] == "WB")
{
axis(1, at=x$tseq)
}
axis(2, at=ylim)
if(x$S.2$time2[1] == 0)
{
lines(x$S.2$time2, x$S.2$S.2, col="red", lwd=3)
lines(x$S.2$time2, x$S.2$UL.2, col="red", lwd=3, lty=3)
lines(x$S.2$time2, x$S.2$LL.2, col="red", lwd=3, lty=3)
}else
{
lines(unique(c(0, x$S.2$time2)), c(1, x$S.2$S.2), col="red", lwd=3)
lines(unique(c(0, x$S.2$time2)), c(1, x$S.2$UL.2), col="red", lwd=3, lty=3)
lines(unique(c(0, x$S.2$time2)), c(1, x$S.2$LL.2), col="red", lwd=3, lty=3)
}
##
plot(c(0, max(x$S.3$time3)), range(ylim), xlab=xlab[3], ylab=ylab, type="n", main = expression(paste("Estimated ", S[3](t), "")), axes=FALSE)
if(x$class[4] == "PEM")
{
axis(1, at=c(0, max(x$S.3$time3)))
}
if(x$class[4] == "WB")
{
axis(1, at=x$tseq)
}
axis(2, at=ylim)
if(x$S.3$time3[1] == 0)
{
lines(x$S.3$time3, x$S.3$S.3, col="red", lwd=3)
lines(x$S.3$time3, x$S.3$UL.3, col="red", lwd=3, lty=3)
lines(x$S.3$time3, x$S.3$LL.3, col="red", lwd=3, lty=3)
}else
{
lines(unique(c(0, x$S.3$time3)), c(1, x$S.3$S.3), col="red", lwd=3)
lines(unique(c(0, x$S.3$time3)), c(1, x$S.3$UL.3), col="red", lwd=3, lty=3)
lines(unique(c(0, x$S.3$time3)), c(1, x$S.3$LL.3), col="red", lwd=3, lty=3)
}
}
}
if(x$class[2] == "Surv")
{
if(is.null(xlab))
{
xlab <- "Time"
}
if(plot.est == "Haz")
{
if(is.null(ylab))
{
ylab <- "Hazard"
}
ygrid <- (max(x$h$UL) - 0)/5
ylim <- seq(from=0, to=max(x$h$UL), by=ygrid)
##
plot(c(0, max(x$h$time)), range(ylim), xlab=xlab, ylab=ylab, type="n", main = expression(paste("Estimated ", h(t), "")), axes=FALSE)
if(x$class[4] == "PEM")
{
axis(1, at=c(0, max(x$h$time)))
}
if(x$class[4] == "WB")
{
axis(1, at=x$tseq)
}
axis(2, at=round(ylim, 4))
lines(x$h$time, x$h$h, col="red", lwd=3)
lines(x$h$time, x$h$UL, col="red", lwd=3, lty=3)
lines(x$h$time, x$h$LL, col="red", lwd=3, lty=3)
}
if(plot.est == "Surv")
{
if(is.null(ylab))
{
ylab <- "Survival"
}
ylim <- seq(from=0, to=1, by=0.2)
##
plot(c(0, max(x$S$time)), range(ylim), xlab=xlab, ylab=ylab, type="n", main = expression(paste("Estimated ", S(t), "")), axes=FALSE)
if(x$class[4] == "PEM")
{
axis(1, at=c(0, max(x$S$time)))
}
if(x$class[4] == "WB")
{
axis(1, at=x$tseq)
}
axis(2, at=ylim)
if(x$S$time[1] == 0)
{
lines(x$S$time, x$S$S, col="red", lwd=3)
lines(x$S$time, x$S$UL, col="red", lwd=3, lty=3)
lines(x$S$time, x$S$LL, col="red", lwd=3, lty=3)
}else
{
lines(unique(c(0, x$S$time)), c(1, x$S$S), col="red", lwd=3)
lines(unique(c(0, x$S$time)), c(1, x$S$UL), col="red", lwd=3, lty=3)
lines(unique(c(0, x$S$time)), c(1, x$S$LL), col="red", lwd=3, lty=3)
}
}
}
}
predict.Bayes_AFT <- function(object, xnew=NULL, x1new=NULL, x2new=NULL, x3new=NULL, time, tseq=c(0, 5, 10), alpha=0.05, ...)
{
conf.level = alpha
x <- object
nChain = x$setup$nChain
if(x$class[2] == "ID")
{
if(!is.null(xnew))
{
stop("'xnew' is for univariate models so it must be specified as NULL for semi-competing risks models")
}
if(is.null(x$chain1$beta1.p))
{
if(!is.null(x1new))
{
stop("x1new does not have the same length as the x1 specified in the Formula.")
}else
{
LP1 <- 0
}
}else if(length(x$chain1$beta1.p) != 0)
{
if(!is.null(x1new))
{
#beta1
p1 = dim(x$chain1$beta1.p)[2]
beta1.p <- x$chain1$beta1.p
if(nChain > 1)
{
for(i in 2:nChain)
{
nam <- paste("chain", i, sep="")
beta1.p <- rbind(beta1.p, x[[nam]]$beta1.p)
}
}
LP1 <- rowSums(beta1.p * matrix(x1new, nrow=dim(beta1.p)[1], ncol = dim(beta1.p)[2], byrow=T))
}else
{
LP1 <- 0
}
}
if(is.null(x$chain1$beta2.p))
{
if(!is.null(x2new))
{
stop("x2new does not have the same length as the x2 specified in the Formula.")
}else
{
LP2 <- 0
}
}else if(length(x$chain1$beta2.p) != 0){
if(!is.null(x2new))
{
#beta2
p2 = dim(x$chain1$beta2.p)[2]
beta2.p <- x$chain1$beta2.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
beta2.p <- rbind(beta2.p, x[[nam]]$beta2.p)
}
}
LP2 <- rowSums(beta2.p * matrix(x2new, nrow=dim(beta2.p)[1], ncol = dim(beta2.p)[2], byrow=T))
}else
{
LP2 <- 0
}
}
if(is.null(x$chain1$beta3.p))
{
if(!is.null(x3new))
{
stop("x3new does not have the same length as the x3 specified in the Formula.")
}else
{
LP3 <- 0
}
}else if(length(x$chain1$beta3.p) != 0){
if(!is.null(x3new))
{
#beta3
p3 = dim(x$chain1$beta3.p)[2]
beta3.p <- x$chain1$beta3.p
if(nChain > 1){
for(i in 2:nChain){
nam <- paste("chain", i, sep="")
beta3.p <- rbind(beta3.p, x[[nam]]$beta3.p)
}
}
LP3 <- rowSums(beta3.p * matrix(x3new, nrow=dim(beta3.p)[1], ncol = dim(beta3.p)[2], byrow=T))
}else
{
LP3 <- 0
}
}
time1 <- time2 <- time3 <- time1hz <- time2hz<- time3hz <- time
if(time[1] == 0){
time1hz <- time1hz[-1]
time2hz <- time2hz[-1]
time3hz <- time3hz[-1]
}
if(x$class[4] == "LN")
{
bsbhN.1 <- BSBH.N(x, time=time1hz, 1, time.trunc=min(time), x1new)
bsbhN.2 <- BSBH.N(x, time=time2hz, 2, time.trunc=min(time), x2new)
bsbhN.3 <- BSBH.N(x, time=time3hz, 3, time.trunc=min(time), x3new)
basehaz1 <- bsbhN.1$BH
basehaz2 <- bsbhN.2$BH
basehaz3 <- bsbhN.3$BH
bsbhN.1 <- BSBH.N(x, time=time1, 1, time.trunc=min(time), x1new)
bsbhN.2 <- BSBH.N(x, time=time2, 2, time.trunc=min(time), x2new)
bsbhN.3 <- BSBH.N(x, time=time3, 3, time.trunc=min(time), x3new)
basesurv1 <- bsbhN.1$BS
basesurv2 <- bsbhN.2$BS
basesurv3 <- bsbhN.3$BS
}
if(x$class[4] == "DPM")
{
bsbhDP.1 <- BSBH.DP(x, time=time1hz, 1, time.trunc=min(time), x1new)
bsbhDP.2 <- BSBH.DP(x, time=time2hz, 1, time.trunc=min(time), x2new)
bsbhDP.3 <- BSBH.DP(x, time=time3hz, 1, time.trunc=min(time), x3new)
basehaz1 <- bsbhDP.1$BH
basehaz2 <- bsbhDP.2$BH
basehaz3 <- bsbhDP.3$BH
bsbhDP.1 <- BSBH.DP(x, time=time1, 1, time.trunc=min(time), x1new)
bsbhDP.2 <- BSBH.DP(x, time=time2, 1, time.trunc=min(time), x2new)
bsbhDP.3 <- BSBH.DP(x, time=time3, 1, time.trunc=min(time), x3new)
basesurv1 <- bsbhDP.1$BS
basesurv2 <- bsbhDP.2$BS
basesurv3 <- bsbhDP.3$BS
}
BH1Med <- apply(basehaz1, 2, median)
BH1Ub <- apply(basehaz1, 2, quantile, prob = (1-conf.level/2))
BH1Lb <- apply(basehaz1, 2, quantile, prob = conf.level/2)
BH2Med <- apply(basehaz2, 2, median)
BH2Ub <- apply(basehaz2, 2, quantile, prob = (1-conf.level/2))
BH2Lb <- apply(basehaz2, 2, quantile, prob = conf.level/2)
BH3Med <- apply(basehaz3, 2, median)
BH3Ub <- apply(basehaz3, 2, quantile, prob = (1-conf.level/2))
BH3Lb <- apply(basehaz3, 2, quantile, prob = conf.level/2)
BS1Med <- apply(basesurv1, 2, median)
BS1Ub <- apply(basesurv1, 2, quantile, prob = (1-conf.level/2))
BS1Lb <- apply(basesurv1, 2, quantile, prob = conf.level/2)
BS2Med <- apply(basesurv2, 2, median)
BS2Ub <- apply(basesurv2, 2, quantile, prob = (1-conf.level/2))
BS2Lb <- apply(basesurv2, 2, quantile, prob = conf.level/2)
BS3Med <- apply(basesurv3, 2, median)
BS3Ub <- apply(basesurv3, 2, quantile, prob = (1-conf.level/2))
BS3Lb <- apply(basesurv3, 2, quantile, prob = conf.level/2)
BH1_tbl <- data.frame(time=time1hz, h.1=BH1Med, LL.1=BH1Lb, UL.1=BH1Ub)
BH2_tbl <- data.frame(time=time2hz, h.2=BH2Med, LL.2=BH2Lb, UL.2=BH2Ub)
BH3_tbl <- data.frame(time=time3hz, h.3=BH3Med, LL.3=BH3Lb, UL.3=BH3Ub)
BS1_tbl <- data.frame(time=time1, S.1=BS1Med, LL.1=BS1Lb, UL.1=BS1Ub)
BS2_tbl <- data.frame(time=time2, S.2=BS2Med, LL.2=BS2Lb, UL.2=BS2Ub)
BS3_tbl <- data.frame(time=time3, S.3=BS3Med, LL.3=BS3Lb, UL.3=BS3Ub)
value <- list(h.1=BH1_tbl, h.2=BH2_tbl, h.3=BH3_tbl, S.1=BS1_tbl, S.2=BS2_tbl, S.3=BS3_tbl)
}
if(x$class[2] == "Surv")
{
if(!is.null(x1new))
{
stop("'x1new' is for semi-competing risks models so it must be specified as NULL for univariate models")
}
if(!is.null(x2new))
{
stop("'x2new' is for semi-competing risks models so it must be specified as NULL for univariate models")
}
if(!is.null(x3new))
{
stop("'x3new' is for semi-competing risks models so it must be specified as NULL for univariate models")
}
if(is.null(x$chain1$beta.p))
{
if(!is.null(xnew))
{
stop("xnew does not have the same length as the x specified in the Formula.")
}else
{
LP <- 0
}
}else if(length(x$chain1$beta.p) != 0)
{
if(!is.null(xnew))
{
#beta
p = dim(x$chain1$beta.p)[2]
beta.p <- x$chain1$beta.p
if(nChain > 1)
{
for(i in 2:nChain)
{
nam <- paste("chain", i, sep="")
beta.p <- rbind(beta.p, x[[nam]]$beta.p)
}
}
LP <- rowSums(beta.p * matrix(xnew, nrow=dim(beta.p)[1], ncol = dim(beta.p)[2], byrow=T))
}else
{
LP <- 0
}
}
timehz <- time
if(time[1] == 0){
timehz <- timehz[-1]
}
if(x$class[4] == "LN")
{
bsbhN.0 <- BSBH.N(x, time=timehz, 0, time.trunc=min(time), xnew)
basehaz <- bsbhN.0$BH
bsbhN.0 <- BSBH.N(x, time=time, 0, time.trunc=min(time), xnew)
basesurv <- bsbhN.0$BS
}
if(x$class[4] == "DPM")
{
bsbhDP.0 <- BSBH.DP(x, time=timehz, 0, time.trunc=min(time), xnew)
basehaz <- bsbhDP.0$BH
bsbhDP.0 <- BSBH.DP(x, time=time, 0, time.trunc=min(time), xnew)
basesurv <- bsbhDP.0$BS
}
BHMed <- apply(basehaz, 2, median)
BHUb <- apply(basehaz, 2, quantile, prob = (1-conf.level/2))
BHLb <- apply(basehaz, 2, quantile, prob = conf.level/2)
BSMed <- apply(basesurv, 2, median)
BSUb <- apply(basesurv, 2, quantile, prob = (1-conf.level/2))
BSLb <- apply(basesurv, 2, quantile, prob = conf.level/2)
BH_tbl <- data.frame(time=timehz, h=BHMed, LL=BHLb, UL=BHUb)
BS_tbl <- data.frame(time=time, S=BSMed, LL=BSLb, UL=BSUb)
value <- list(h=BH_tbl, S=BS_tbl)
}
value$xnew <- xnew
value$x1new <- x1new
value$x2new <- x2new
value$x3new <- x3new
value$tseq <- tseq
value$setup$model <- x$setup$model
value$class <- x$class
class(value) <- "pred.Bayes_AFT"
return(value)
}
plot.pred.Bayes_AFT <- function(x, plot.est="Haz", xlab=NULL, ylab=NULL, ...)
{
if(x$class[2] == "ID")
{
if(is.null(xlab))
{
xlab <- c("Time", "Time", "Time")
xlab[3] <- "Time since non-terminal event"
}
if(plot.est == "Haz")
{
if(is.null(ylab))
{
ylab <- "Hazard"
}
ygrid <- (max(x$h.1$UL.1, x$h.2$UL.2, x$h.3$UL.3) - 0)/5
ylim <- seq(from=0, to=max(x$h.1$UL.1, x$h.2$UL.2, x$h.3$UL.3), by=ygrid)
##
par(mfrow=c(1,3))
##
plot(c(0, max(x$h.1$time)), range(ylim), xlab=xlab[1], ylab=ylab, type="n", main = expression(paste("Estimated ", h[1](t), "")), axes=FALSE)
axis(1, at=x$tseq)
axis(2, at=round(ylim, 4))
lines(x$h.1$time, x$h.1$h.1, col="blue", lwd=3)
lines(x$h.1$time, x$h.1$UL.1, col="blue", lwd=3, lty=3)
lines(x$h.1$time, x$h.1$LL.1, col="blue", lwd=3, lty=3)
##
plot(c(0, max(x$h.2$time)), range(ylim), xlab=xlab[2], ylab=ylab, type="n", main = expression(paste("Estimated ", h[2](t), "")), axes=FALSE)
axis(1, at=x$tseq)
axis(2, at=round(ylim, 4))
lines(x$h.2$time, x$h.2$h.2, col="red", lwd=3)
lines(x$h.2$time, x$h.2$UL.2, col="red", lwd=3, lty=3)
lines(x$h.2$time, x$h.2$LL.2, col="red", lwd=3, lty=3)
##
plot(c(0, max(x$h.3$time)), range(ylim), xlab=xlab[3], ylab=ylab, type="n", main = expression(paste("Estimated ", h[3](t), "")), axes=FALSE)
axis(1, at=x$tseq)
axis(2, at=round(ylim, 4))
lines(x$h.3$time, x$h.3$h.3, col="red", lwd=3)
lines(x$h.3$time, x$h.3$UL.3, col="red", lwd=3, lty=3)
lines(x$h.3$time, x$h.3$LL.3, col="red", lwd=3, lty=3)
}
if(plot.est == "Surv")
{
if(is.null(ylab))
{
ylab <- "Survival"
}
ylim <- seq(from=0, to=1, by=0.2)
##
par(mfrow=c(1,3))
##
plot(c(0, max(x$S.1$time)), range(ylim), xlab=xlab[1], ylab=ylab, type="n", main = expression(paste("Estimated ", S[1](t), "")), axes=FALSE)
axis(1, at=x$tseq)
axis(2, at=ylim)
if(x$S.1$time[1] == 0)
{
lines(x$S.1$time, x$S.1$S.1, col="blue", lwd=3)
lines(x$S.1$time, x$S.1$UL.1, col="blue", lwd=3, lty=3)
lines(x$S.1$time, x$S.1$LL.1, col="blue", lwd=3, lty=3)
}else
{
lines(unique(c(0, x$S.1$time)), c(1, x$S.1$S.1), col="blue", lwd=3)
lines(unique(c(0, x$S.1$time)), c(1, x$S.1$UL.1), col="blue", lwd=3, lty=3)
lines(unique(c(0, x$S.1$time)), c(1, x$S.1$LL.1), col="blue", lwd=3, lty=3)
}
##
plot(c(0, max(x$S.2$time)), range(ylim), xlab=xlab[2], ylab=ylab, type="n", main = expression(paste("Estimated ", S[2](t), "")), axes=FALSE)
axis(1, at=x$tseq)
axis(2, at=ylim)
if(x$S.2$time[1] == 0)
{
lines(x$S.2$time, x$S.2$S.2, col="red", lwd=3)
lines(x$S.2$time, x$S.2$UL.2, col="red", lwd=3, lty=3)
lines(x$S.2$time, x$S.2$LL.2, col="red", lwd=3, lty=3)
}else
{
lines(unique(c(0, x$S.2$time)), c(1, x$S.2$S.2), col="red", lwd=3)
lines(unique(c(0, x$S.2$time)), c(1, x$S.2$UL.2), col="red", lwd=3, lty=3)
lines(unique(c(0, x$S.2$time)), c(1, x$S.2$LL.2), col="red", lwd=3, lty=3)
}
##
plot(c(0, max(x$S.3$time)), range(ylim), xlab=xlab[3], ylab=ylab, type="n", main = expression(paste("Estimated ", S[3](t), "")), axes=FALSE)
axis(1, at=x$tseq)
axis(2, at=ylim)
if(x$S.3$time[1] == 0)
{
lines(x$S.3$time, x$S.3$S.3, col="red", lwd=3)
lines(x$S.3$time, x$S.3$UL.3, col="red", lwd=3, lty=3)
lines(x$S.3$time, x$S.3$LL.3, col="red", lwd=3, lty=3)
}else
{
lines(unique(c(0, x$S.3$time)), c(1, x$S.3$S.3), col="red", lwd=3)
lines(unique(c(0, x$S.3$time)), c(1, x$S.3$UL.3), col="red", lwd=3, lty=3)
lines(unique(c(0, x$S.3$time)), c(1, x$S.3$LL.3), col="red", lwd=3, lty=3)
}
}
}
if(x$class[2] == "Surv")
{
if(is.null(xlab))
{
xlab <- "Time"
}
if(plot.est == "Haz")
{
if(is.null(ylab))
{
ylab <- "Hazard"
}
ygrid <- (max(x$h$UL) - 0)/5
ylim <- seq(from=0, to=max(x$h$UL), by=ygrid)
##
plot(c(0, max(x$h$time)), range(ylim), xlab=xlab, ylab=ylab, type="n", main = expression(paste("Estimated ", h(t), "")), axes=FALSE)
axis(1, at=x$tseq)
axis(2, at=round(ylim, 4))
lines(x$h$time, x$h$h, col="red", lwd=3)
lines(x$h$time, x$h$UL, col="red", lwd=3, lty=3)
lines(x$h$time, x$h$LL, col="red", lwd=3, lty=3)
}
if(plot.est == "Surv")
{
if(is.null(ylab))
{
ylab <- "Survival"
}
ylim <- seq(from=0, to=1, by=0.2)
##
plot(c(0, max(x$S$time)), range(ylim), xlab=xlab, ylab=ylab, type="n", main = expression(paste("Estimated ", S(t), "")), axes=FALSE)
axis(1, at=x$tseq)
axis(2, at=ylim)
if(x$S$time[1] == 0)
{
lines(x$S$time, x$S$S, col="red", lwd=3)
lines(x$S$time, x$S$UL, col="red", lwd=3, lty=3)
lines(x$S$time, x$S$LL, col="red", lwd=3, lty=3)
}else
{
lines(unique(c(0, x$S$time)), c(1, x$S$S), col="red", lwd=3)
lines(unique(c(0, x$S$time)), c(1, x$S$UL), col="red", lwd=3, lty=3)
lines(unique(c(0, x$S$time)), c(1, x$S$LL), col="red", lwd=3, lty=3)
}
}
}
}
predict.Freq_HReg <- function(object, xnew=NULL, x1new=NULL, x2new=NULL, x3new=NULL, tseq=c(0, 5, 10), alpha=0.05, ...)
{
conf.level = alpha
obj <- object
T2seq <- tseq
yLim <- NULL
##
## SEs based on the Delta method using log(-log(S0))
##
if(obj$class[2] == "Surv")
{
if(!is.null(x1new))
{
stop("'x1new' is for semi-competing risks models so it must be specified as NULL for univariate models")
}
if(!is.null(x2new))
{
stop("'x2new' is for semi-competing risks models so it must be specified as NULL for univariate models")
}
if(!is.null(x3new))
{
stop("'x3new' is for semi-competing risks models so it must be specified as NULL for univariate models")
}
T2 <- seq(from=min(T2seq), to=max(T2seq), length=100)
##
kappa <- exp(obj$estimate[1])
alpha <- exp(obj$estimate[2])
log_kappa <- obj$estimate[1]
log_alpha <- obj$estimate[2]
S0 <- exp(-(kappa*(T2)^alpha))
## Delta method based on log(-log(S0))
if(!is.null(xnew))
{
J <- cbind(1, exp(log_alpha)*log(T2), matrix(xnew, nrow=100, ncol=length(xnew), byrow=T))
Var.loglogS0 <- J %*% obj$Finv %*% t(J)
}else
{
J <- cbind(1, exp(log_alpha)*log(T2))
Var.loglogS0 <- J %*% obj$Finv[1:2, 1:2] %*% t(J)
}
se.loglogS0 <- sqrt(diag(Var.loglogS0))
se.loglogS0[is.na(se.loglogS0)] <- 0
LL <- S0^exp(-qnorm(conf.level/2)*se.loglogS0)
UL <- S0^exp(qnorm(conf.level/2)*se.loglogS0)
##
BS_tbl <- data.frame(time=T2, S=S0, LL=LL, UL=UL)
##
h0 <- alpha*kappa*(T2)^(alpha-1)
if(!is.null(xnew))
{
J <- cbind(h0, h0*(1+alpha*log(T2)), h0*matrix(xnew, nrow=100, ncol=length(xnew), byrow=T))
Var.h0 <- J %*% obj$Finv %*% t(J)
}else
{
J <- cbind(h0, h0*(1+alpha*log(T2)))
Var.h0 <- J %*% obj$Finv[1:2,1:2] %*% t(J)
}
se.h0 <- sqrt(diag(Var.h0))
se.h0[is.nan(se.h0)] <- 0
LLh0 <- h0 - qnorm(conf.level/2)*se.h0
ULh0 <- h0 + qnorm(conf.level/2)*se.h0
LLh0[LLh0 < 0] <- 0
T2h <- T2
if(T2[1] == 0)
{
T2h <- T2h[-1]
h0 <- h0[-1]
LLh0 <- LLh0[-1]
ULh0 <- ULh0[-1]
}
BH_tbl <- data.frame(time=T2h, h=h0, LL=LLh0, UL=ULh0)
value <- list(h=BH_tbl, S=BS_tbl)
}
##
if(obj$class[2] == "ID")
{
if(!is.null(xnew))
{
stop("'xnew' is for univariate models so it must be specified as NULL for semi-competing risks models")
}
nP = obj$nP
##
T2 <- seq(from=min(T2seq), to=max(T2seq), length=100)
##
kappa <- exp(obj$estimate[1])
alpha <- exp(obj$estimate[2])
log_alpha <- obj$estimate[2]
S0.1 <- exp(-kappa*(T2)^alpha)
if(!is.null(x1new) & nP[1] > 0)
{
J <- cbind(1, exp(log_alpha)*log(T2), matrix(x1new, nrow=100, ncol=length(x1new), byrow=T))
if(obj$frailty == TRUE)
{
Var.loglogS0 <- J %*% obj$Finv[c(1:2, 8:(8+nP[1]-1)), c(1:2, 8:(8+nP[1]-1))] %*% t(J)
}else if(obj$frailty == FALSE)
{
Var.loglogS0 <- J %*% obj$Finv[c(1:2, 7:(7+nP[1]-1)), c(1:2, 7:(7+nP[1]-1))] %*% t(J)
}
}else if(is.null(x1new) | nP[1] == 0)
{
J <- cbind(1, exp(log_alpha)*log(T2))
Var.loglogS0 <- J %*% obj$Finv[1:2,1:2] %*% t(J)
}
se.loglogS0 <- sqrt(diag(Var.loglogS0))
LL.1 <- S0.1^exp(-qnorm(conf.level/2)*se.loglogS0)
UL.1 <- S0.1^exp(qnorm(conf.level/2)*se.loglogS0)
##
h0.1 <- alpha*kappa*(T2)^(alpha-1)
if(!is.null(x1new) & nP[1] > 0)
{
J <- cbind(h0.1, h0.1*(1+alpha*log(T2)), h0.1*matrix(x1new, nrow=100, ncol=length(x1new), byrow=T))
if(obj$frailty == TRUE)
{
Var.h0.1 <- J %*% obj$Finv[c(1:2, 8:(8+nP[1]-1)), c(1:2, 8:(8+nP[1]-1))] %*% t(J)
}else if(obj$frailty == FALSE)
{
Var.h0.1 <- J %*% obj$Finv[c(1:2, 7:(7+nP[1]-1)), c(1:2, 7:(7+nP[1]-1))] %*% t(J)
}
}else if(is.null(x1new) | nP[1] == 0)
{
J <- cbind(h0.1, h0.1*(1+alpha*log(T2)))
Var.h0.1 <- J %*% obj$Finv[1:2,1:2] %*% t(J)
}
se.h0.1 <- sqrt(diag(Var.h0.1))
se.h0.1[is.nan(se.h0.1)] <- 0
LLh0.1 <- h0.1 - qnorm(conf.level/2)*se.h0.1
ULh0.1 <- h0.1 + qnorm(conf.level/2)*se.h0.1
LLh0.1[LLh0.1 < 0] <- 0
##
kappa <- exp(obj$estimate[3])
alpha <- exp(obj$estimate[4])
log_alpha <- obj$estimate[4]
S0.2 <- exp(-kappa*(T2)^alpha)
if(!is.null(x2new) & nP[2] > 0)
{
J <- cbind(1, exp(log_alpha)*log(T2), matrix(x2new, nrow=100, ncol=length(x2new), byrow=T))
if(obj$frailty == TRUE)
{
Var.loglogS0 <- J %*% obj$Finv[c(3:4, (8+nP[1]):(8+nP[1]+nP[2]-1)), c(3:4, (8+nP[1]):(8+nP[1]+nP[2]-1))] %*% t(J)
}else if(obj$frailty == FALSE)
{
Var.loglogS0 <- J %*% obj$Finv[c(3:4, (7+nP[1]):(7+nP[1]+nP[2]-1)), c(3:4, (7+nP[1]):(7+nP[1]+nP[2]-1))] %*% t(J)
}
}else if(is.null(x2new) | nP[2] == 0)
{
J <- cbind(1, exp(log_alpha)*log(T2))
Var.loglogS0 <- J %*% obj$Finv[3:4,3:4] %*% t(J)
}
se.loglogS0 <- sqrt(diag(Var.loglogS0))
LL.2 <- S0.2^exp(-qnorm(conf.level/2)*se.loglogS0)
UL.2 <- S0.2^exp(qnorm(conf.level/2)*se.loglogS0)
##
h0.2 <- alpha*kappa*(T2)^(alpha-1)
if(!is.null(x2new) & nP[2] > 0)
{
J <- cbind(h0.2, h0.2*(1+alpha*log(T2)), h0.2*matrix(x2new, nrow=100, ncol=length(x2new), byrow=T))
if(obj$frailty == TRUE)
{
Var.h0.2 <- J %*% obj$Finv[c(3:4, (8+nP[1]):(8+nP[1]+nP[2]-1)), c(3:4, (8+nP[1]):(8+nP[1]+nP[2]-1))] %*% t(J)
}else if(obj$frailty == FALSE)
{
Var.h0.2 <-J %*% obj$Finv[c(3:4, (7+nP[1]):(7+nP[1]+nP[2]-1)), c(3:4, (7+nP[1]):(7+nP[1]+nP[2]-1))] %*% t(J)
}
}else if(is.null(x2new) | nP[2] == 0)
{
J <- cbind(h0.2, h0.2*(1+alpha*log(T2)))
Var.h0.2 <- J %*% obj$Finv[3:4,3:4] %*% t(J)
}
se.h0.2 <- sqrt(diag(Var.h0.2))
se.h0.2[is.nan(se.h0.2)] <- 0
LLh0.2 <- h0.2 - qnorm(conf.level/2)*se.h0.2
ULh0.2 <- h0.2 + qnorm(conf.level/2)*se.h0.2
LLh0.2[LLh0.2 < 0] <- 0
##
kappa <- exp(obj$estimate[5])
alpha <- exp(obj$estimate[6])
log_alpha <- obj$estimate[6]
S0.3 <- exp(-kappa*(T2)^alpha)
if(!is.null(x3new) & nP[3] > 0)
{
J <- cbind(1, exp(log_alpha)*log(T2), matrix(x3new, nrow=100, ncol=length(x3new), byrow=T))
if(obj$frailty == TRUE)
{
Var.loglogS0 <- J %*% obj$Finv[c(5:6, (8+nP[1]+nP[2]):(8+nP[1]+nP[2]++nP[3]-1)), c(5:6, (8+nP[1]+nP[3]):(8+nP[1]+nP[2]+nP[3]-1))] %*% t(J)
}else if(obj$frailty == FALSE)
{
Var.loglogS0 <- J %*% obj$Finv[c(5:6, (7+nP[1]+nP[2]):(7+nP[1]+nP[2]++nP[3]-1)), c(5:6, (7+nP[1]+nP[3]):(7+nP[1]+nP[2]+nP[3]-1))] %*% t(J)
}
}else if(is.null(x3new) | nP[3] == 0)
{
J <- cbind(1, exp(log_alpha)*log(T2))
Var.loglogS0 <- J %*% obj$Finv[5:6,5:6] %*% t(J)
}
se.loglogS0 <- sqrt(diag(Var.loglogS0))
LL.3 <- S0.3^exp(-qnorm(conf.level/2)*se.loglogS0)
UL.3 <- S0.3^exp(qnorm(conf.level/2)*se.loglogS0)
##
h0.3 <- alpha*kappa*(T2)^(alpha-1)
if(!is.null(x3new) & nP[3] > 0)
{
J <- cbind(h0.3, h0.3*(1+alpha*log(T2)), h0.3*matrix(x3new, nrow=100, ncol=length(x3new), byrow=T))
if(obj$frailty == TRUE)
{
Var.h0.3 <- J %*% obj$Finv[c(5:6, (8+nP[1]+nP[2]):(8+nP[1]+nP[2]++nP[3]-1)), c(5:6, (8+nP[1]+nP[3]):(8+nP[1]+nP[2]+nP[3]-1))] %*% t(J)
}else if(obj$frailty == FALSE)
{
Var.h0.3 <- J %*% obj$Finv[c(5:6, (7+nP[1]+nP[2]):(7+nP[1]+nP[2]++nP[3]-1)), c(5:6, (7+nP[1]+nP[3]):(7+nP[1]+nP[2]+nP[3]-1))] %*% t(J)
}
}else if(is.null(x3new) | nP[3] == 0)
{
J <- cbind(h0.3, h0.3*(1+alpha*log(T2)))
Var.h0.3 <- J %*% obj$Finv[5:6,5:6] %*% t(J)
}
se.h0.3 <- sqrt(diag(Var.h0.3))
se.h0.3[is.nan(se.h0.3)] <- 0
LLh0.3 <- h0.3 - qnorm(conf.level/2)*se.h0.3
ULh0.3 <- h0.3 + qnorm(conf.level/2)*se.h0.3
LLh0.3[LLh0.3 < 0] <- 0
T2h <- T2
if(T2[1] == 0)
{
T2h <- T2h[-1]
h0.1 <- h0.1[-1]
LLh0.1 <- LLh0.1[-1]
ULh0.1 <- ULh0.1[-1]
h0.2 <- h0.2[-1]
LLh0.2 <- LLh0.2[-1]
ULh0.2 <- ULh0.2[-1]
h0.3 <- h0.3[-1]
LLh0.3 <- LLh0.3[-1]
ULh0.3 <- ULh0.3[-1]
}
BH1_tbl <- data.frame(time=T2h, h.1=h0.1, LL.1=LLh0.1, UL.1=ULh0.1)
BH2_tbl <- data.frame(time=T2h, h.2=h0.2, LL.2=LLh0.2, UL.2=ULh0.2)
BH3_tbl <- data.frame(time=T2h, h.3=h0.3, LL.3=LLh0.3, UL.3=ULh0.3)
BS1_tbl <- data.frame(time=T2, S.1=S0.1, LL.1=LL.1, UL.1=UL.1)
BS2_tbl <- data.frame(time=T2, S.2=S0.2, LL.2=LL.2, UL.2=UL.2)
BS3_tbl <- data.frame(time=T2, S.3=S0.3, LL.3=LL.3, UL.3=UL.3)
value <- list(h.1=BH1_tbl, h.2=BH2_tbl, h.3=BH3_tbl, S.1=BS1_tbl, S.2=BS2_tbl, S.3=BS3_tbl)
}
value$xnew <- xnew
value$x1new <- x1new
value$x2new <- x2new
value$x3new <- x3new
value$tseq <- tseq
value$setup$model <- obj$setup$model
value$class <- obj$class
class(value) <- "pred.Freq_HReg"
return(value)
}
plot.pred.Freq_HReg <- function(x, plot.est="Haz", xlab=NULL, ylab=NULL, ...)
{
obj <- x
T2seq <- x$tseq
yLim <- NULL
##
## SEs based on the Delta method using log(-log(S0))
##
if(obj$class[2] == "Surv")
{
##
if(is.null(yLim))
{
if(plot.est=="Surv")
{
yLim <- seq(from=0, to=1, by=0.2)
}
if(plot.est=="Haz")
{
grid <- (max(obj$h$UL) - min(obj$h$LL))/5
yLim <- seq(from=min(obj$h$LL), to=max(obj$h$UL), by=grid)
}
}
##
if(is.null(ylab))
{
if(plot.est=="Surv")
{
ylab <- "Survival"
}
if(plot.est=="Haz")
{
ylab <- "Hazard"
}
}
##
if(is.null(xlab)) xlab <- "Time"
##
if(plot.est == "Surv")
{
##
plot(range(T2seq), range(yLim), xlab=xlab, ylab=ylab, type="n", main = expression(paste("Estimated ", S(t), "")), axes=FALSE)
axis(1, at=T2seq)
axis(2, at=yLim)
lines(obj$S$time, obj$S$S, col="red", lwd=3)
lines(obj$S$time, obj$S$LL, col="red", lwd=3, lty=3)
lines(obj$S$time, obj$S$UL, col="red", lwd=3, lty=3)
}
if(plot.est == "Haz")
{
##
plot(range(T2seq), range(yLim), xlab=xlab, ylab=ylab, type="n", main = expression(paste("Estimated ", h(t), "")), axes=FALSE)
axis(1, at=T2seq)
axis(2, at=round(yLim, 4))
lines(obj$h$time, obj$h$h, col="red", lwd=3)
lines(obj$h$time, obj$h$LL, col="red", lwd=3, lty=3)
lines(obj$h$time, obj$h$UL, col="red", lwd=3, lty=3)
}
}
##
if(obj$class[2] == "ID")
{
##
if(is.null(ylab))
{
if(plot.est=="Surv")
{
ylab <- "Survival"
}
if(plot.est=="Haz")
{
ylab <- "Hazard"
}
}
##
if(is.null(xlab))
{
xlab <- c("Time", "Time", "Time")
if(obj$class[5] == "semi-Markov")
{
xlab[3] <- "Time since non-terminal event"
}
}
##
if(is.null(yLim))
{
if(plot.est=="Surv")
{
yLim <- seq(from=0, to=1, by=0.2)
}
if(plot.est=="Haz")
{
ygrid <- (max(x$h.1$UL.1, x$h.2$UL.2, x$h.3$UL.3) - 0)/5
yLim <- seq(from=0, to=max(x$h.1$UL.1, x$h.2$UL.2, x$h.3$UL.3), by=ygrid)
}
}
##
if(plot.est == "Surv")
{
##
par(mfrow=c(1,3))
##
plot(range(T2seq), range(yLim), xlab=xlab[1], ylab=ylab, type="n", main = expression(paste("Estimated ", S[1](t), "")), axes=FALSE)
axis(1, at=T2seq)
axis(2, at=yLim)
lines(obj$S.1$time, obj$S.1$S.1, col="blue", lwd=3)
lines(obj$S.1$time, obj$S.1$LL.1, col="blue", lwd=3, lty=3)
lines(obj$S.1$time, obj$S.1$UL.1, col="blue", lwd=3, lty=3)
##
plot(range(T2seq), range(yLim), xlab=xlab[2], ylab=ylab, type="n", main = expression(paste("Estimated ", S[2](t), "")), axes=FALSE)
axis(1, at=T2seq)
axis(2, at=yLim)
lines(obj$S.2$time, obj$S.2$S.2, col="red", lwd=3)
lines(obj$S.2$time, obj$S.2$LL.2, col="red", lwd=3, lty=3)
lines(obj$S.2$time, obj$S.2$UL.2, col="red", lwd=3, lty=3)
##
plot(range(T2seq), range(yLim), xlab=xlab[3], ylab=ylab, type="n", main = expression(paste("Estimated ", S[3](t), "")), axes=FALSE)
axis(1, at=T2seq)
axis(2, at=yLim)
lines(obj$S.3$time, obj$S.3$S.3, col="red", lwd=3)
lines(obj$S.3$time, obj$S.3$LL.3, col="red", lwd=3, lty=3)
lines(obj$S.3$time, obj$S.3$UL.3, col="red", lwd=3, lty=3)
}
if(plot.est == "Haz")
{
##
par(mfrow=c(1,3))
##
plot(range(T2seq), range(yLim), xlab=xlab[1], ylab=ylab, type="n", main = expression(paste("Estimated ", h[1](t), "")), axes=FALSE)
axis(1, at=T2seq)
axis(2, at=round(yLim, 4))
lines(obj$h.1$time, obj$h.1$h.1, col="blue", lwd=3)
lines(obj$h.1$time, obj$h.1$LL.1, col="blue", lwd=3, lty=3)
lines(obj$h.1$time, obj$h.1$UL.1, col="blue", lwd=3, lty=3)
##
plot(range(T2seq), range(yLim), xlab=xlab[2], ylab=ylab, type="n", main = expression(paste("Estimated ", h[2](t), "")), axes=FALSE)
axis(1, at=T2seq)
axis(2, at=round(yLim, 4))
lines(obj$h.2$time, obj$h.2$h.2, col="red", lwd=3)
lines(obj$h.2$time, obj$h.2$LL.2, col="red", lwd=3, lty=3)
lines(obj$h.2$time, obj$h.2$UL.2, col="red", lwd=3, lty=3)
##
plot(range(T2seq), range(yLim), xlab=xlab[3], ylab=ylab, type="n", main = expression(paste("Estimated ", h[3](t), "")), axes=FALSE)
axis(1, at=T2seq)
axis(2, at=round(yLim, 4))
lines(obj$h.3$time, obj$h.3$h.3, col="red", lwd=3)
lines(obj$h.3$time, obj$h.3$LL.3, col="red", lwd=3, lty=3)
lines(obj$h.3$time, obj$h.3$UL.3, col="red", lwd=3, lty=3)
}
}
##
invisible()
}
####
## VCOV METHOD
####
##
vcov.Freq_HReg <- function(object, ...)
{
obj <- object
value <- obj$Finv
dimnames(value) <- list(obj$myLabels, obj$myLabels)
if(obj$class[2] == "ID")
{
cat("\nNote: Covariates are arranged in order of transition number, 1->3.\n")
}
return(value)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.