R/methods.R

Defines functions vcov.Freq_HReg plot.pred.Freq_HReg predict.Freq_HReg plot.pred.Bayes_AFT predict.Bayes_AFT plot.pred.Bayes_HReg predict.Bayes_HReg coef.Freq_HReg coef.Bayes_AFT coef.Bayes_HReg print.summ.Bayes_AFT print.summ.Bayes_HReg print.summ.Freq_HReg summary.Bayes_AFT summary.Bayes_HReg summary.Freq_HReg print.Bayes_AFT print.Bayes_HReg print.Freq_HReg

Documented in coef.Bayes_AFT coef.Bayes_HReg coef.Freq_HReg plot.pred.Bayes_AFT plot.pred.Bayes_HReg plot.pred.Freq_HReg predict.Bayes_AFT predict.Bayes_HReg predict.Freq_HReg print.Bayes_AFT print.Bayes_HReg print.Freq_HReg print.summ.Bayes_AFT print.summ.Bayes_HReg print.summ.Freq_HReg summary.Bayes_AFT summary.Bayes_HReg summary.Freq_HReg vcov.Freq_HReg

####
## 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)
}

Try the SemiCompRisks package in your browser

Any scripts or data that you put into this service are public.

SemiCompRisks documentation built on Feb. 3, 2021, 5:06 p.m.