R/BayesID_HReg.R

Defines functions BayesID_HReg

Documented in BayesID_HReg

BayesID_HReg <- function(Formula,
data,
id = NULL,
model = c("semi-Markov", "Weibull"),
hyperParams,
startValues,
mcmcParams,
na.action = "na.fail",
subset=NULL,
path = NULL)
{
    if(na.action != "na.fail" & na.action != "na.omit")
    {
        stop("na.action should be either na.fail or na.omit")
    }
    if(!is.null(id))
    {
        data$id <- id
        form2 <- as.Formula(paste(Formula[2], Formula[1], Formula[3], "| id", sep = ""))
    }else
    {
        form2 <- as.Formula(paste(Formula[2], Formula[1], Formula[3], sep = ""))
    }
    
    nChain <- length(startValues)
    for(i in 1:nChain)
    {
        nam <- paste("gam", i, sep = "")
        data[[nam]] <- startValues[[i]]$common$gamma.ji
        form2 <- as.Formula(paste(form2[2], form2[1], form2[3], "| ", nam, sep = ""))
    }
    
    data <- model.frame(form2, data=data, na.action = na.action, subset = subset)
    if(!is.null(id))
    {
        id <- data$id
    }
    
    
    for(i in 1:nChain)
    {
        nam <- paste("gam", i, sep = "")
        startValues[[i]]$common$gamma.ji <- data[[nam]]
    }
    
    mcmc        <- mcmcParams
    hyperP      <- hyperParams
    mcmcList    <- mcmc
    
    time1 <- model.part(Formula, data=data, lhs=1)
    time2 <- model.part(Formula, data=data, lhs=2)
    
    Y <- cbind(time1[1], time1[2], time2[1], time2[2])
    
    Xmat1 <- model.frame(formula(Formula, lhs=0, rhs=1), data=data)
    Xmat2 <- model.frame(formula(Formula, lhs=0, rhs=2), data=data)
    Xmat3 <- model.frame(formula(Formula, lhs=0, rhs=3), data=data)
    
    p1 <- ncol(Xmat1)
    p2 <- ncol(Xmat2)
    p3 <- ncol(Xmat3)
    
    nCov <- c(p1, p2, p3)
    
    if((mcmcList$run$numReps / mcmcList$run$thin * mcmcList$run$burninPerc) %% 1 == 0)
    {
        ## for independent semi-competing risks data
        
        if(is.null(id))
        {
            nChain <- length(startValues)
            
            model_h3    <- model[1]
            hz.type     <- model[2]
            
            ##
            
            if(p1 == 0 & p2 == 0 & p3 == 0){
                survData <- Y
            }
            if(p1 == 0 & p2 > 0 & p3 > 0){
                survData <- cbind(Y, Xmat2, Xmat3)
            }
            if(p1 > 0 & p2 == 0 & p3 > 0){
                survData <- cbind(Y, Xmat1, Xmat3)
            }
            if(p1 > 0 & p2 > 0 & p3 == 0){
                survData <- cbind(Y, Xmat1, Xmat2)
            }
            if(p1 > 0 & p2 == 0 & p3 == 0){
                survData <- cbind(Y, Xmat1)
            }
            if(p1 == 0 & p2 > 0 & p3 == 0){
                survData <- cbind(Y, Xmat2)
            }
            if(p1 == 0 & p2 == 0 & p3 > 0){
                survData <- cbind(Y, Xmat3)
            }
            if(p1 > 0 & p2 > 0 & p3 > 0){
                survData <- cbind(Y, Xmat1, Xmat2, Xmat3)
            }
            
            n    <- dim(survData)[1]
            
            #if(class(startValues) == "list" & length(startValues) == nChain){
            
            if(!is.null(path)){
                dir.create(paste(path), recursive = TRUE, showWarnings = FALSE)
            }
            
            ### setting hyperparameters
            
            if(hz.type == "Weibull")
            {
                hyperParams <- as.vector(c(hyperP$WB$WB.ab1, hyperP$WB$WB.ab2, hyperP$WB$WB.ab3, hyperP$WB$WB.cd1, hyperP$WB$WB.cd2, hyperP$WB$WB.cd3, hyperP$theta))
            }
            
            if(hz.type == "PEM")
            {
                hyperParams <- as.vector(c(hyperP$PEM$PEM.ab1, hyperP$PEM$PEM.ab2, hyperP$PEM$PEM.ab3, hyperP$PEM$PEM.alpha1, hyperP$PEM$PEM.alpha2, hyperP$PEM$PEM.alpha3, 1, 1, 1, hyperP$theta))
            }
            
            
            
            ### mcmc setting
            
            if(hz.type == "Weibull")
            {
                mcmc <- as.vector(c(mhProp_alpha1_var=mcmcList$tuning$mhProp_alphag_var[1], mhProp_alpha2_var=mcmcList$tuning$mhProp_alphag_var[2], mhProp_alpha3_var=mcmcList$tuning$mhProp_alphag_var[3], mhProp_theta_var=mcmcList$tuning$mhProp_theta_var, nGam_save=mcmcList$storage$nGam_save, numReps=mcmcList$run$numReps, thin=mcmcList$run$thin, burninPerc=mcmcList$run$burninPerc))
            }
            if(hz.type == "PEM")
            {
                mcmcList$tuning$sg_max <- c(1,1,1)
                mcmc <- as.vector(c(C1=mcmcList$tuning$Cg[1], C2=mcmcList$tuning$Cg[2], C3=mcmcList$tuning$Cg[3], delPert1=mcmcList$tuning$delPertg[1], delPert2=mcmcList$tuning$delPertg[2], delPert3=mcmcList$tuning$delPertg[3], rj.scheme = mcmcList$tuning$rj.scheme, K1_max=mcmcList$tuning$Kg_max[1], K2_max=mcmcList$tuning$Kg_max[2], K3_max=mcmcList$tuning$Kg_max[3], s1_max=mcmcList$tuning$sg_max[1], s2_max=mcmcList$tuning$sg_max[2], s3_max=mcmcList$tuning$sg_max[3], mhProp_theta_var=mcmcList$tuning$mhProp_theta_var, nGam_save=mcmcList$storage$nGam_save, numReps=mcmcList$run$numReps, thin=mcmcList$run$thin, burninPerc=mcmcList$run$burninPerc))
            }
            
            chain = 1
            ret <- list()
            
            while(chain <= nChain){
                
                cat("chain: ", chain, "\n")
                nam = paste("chain", chain, sep="")
                
                temp <- startValues[[chain]]
                
                
                ### setting starting values
                
                if(hz.type == "Weibull")
                {
                    startV <- as.vector(c(beta1=temp$common$beta1, beta2=temp$common$beta2, beta3=temp$common$beta3, temp$WB$WB.alpha, temp$WB$WB.kappa, theta=temp$common$theta, gamma=temp$common$gamma.ji))
                }
                if(hz.type == "PEM")
                {
                    startV <- as.vector(c(beta1=temp$common$beta1, beta2=temp$common$beta2, beta3=temp$common$beta3, theta=temp$common$theta, gamma=temp$common$gamma.ji))
                }
                
                
                # hz.type = "Weibull"
                
                if(hz.type == "Weibull"){
                    
                    nGam_save   <- mcmc[5]
                    numReps     <- mcmc[6]
                    thin        <- mcmc[7]
                    burninPerc  <- mcmc[8]
                    
                    nStore <- round(numReps/thin*(1-burninPerc))
                    
                    mcmcParams <- mcmc[1:4]
                    
                    # model_h3 = "Markov"
                    
                    if(model_h3 == "Markov"){
                        
                        ###
                        
                        mcmcRet <- .C("BweibScrmcmc",
                        survData         = as.double(as.matrix(survData)),
                        n                = as.integer(n),
                        p1                = as.integer(p1),
                        p2                = as.integer(p2),
                        p3                = as.integer(p3),
                        hyperParams     = as.double(hyperParams),
                        mcmcParams        = as.double(mcmcParams),
                        startValues     = as.double(startV),
                        numReps            = as.integer(numReps),
                        thin            = as.integer(thin),
                        burninPerc      = as.double(burninPerc),
                        nGam_save        = as.integer(nGam_save),
                        samples_beta1     = as.double(rep(0, nStore*p1)),
                        samples_beta2     = as.double(rep(0, nStore*p2)),
                        samples_beta3     = as.double(rep(0, nStore*p3)),
                        samples_alpha1     = as.double(rep(0, nStore*1)),
                        samples_alpha2     = as.double(rep(0, nStore*1)),
                        samples_alpha3     = as.double(rep(0, nStore*1)),
                        samples_kappa1     = as.double(rep(0, nStore*1)),
                        samples_kappa2     = as.double(rep(0, nStore*1)),
                        samples_kappa3     = as.double(rep(0, nStore*1)),
                        samples_theta     = as.double(rep(0, nStore*1)),
                        samples_gamma     = as.double(rep(0, nStore*nGam_save)),
                        samples_misc    = as.double(rep(0, (p1+p2+p3+4))),
                        dev                 = as.double(rep(0, nStore*1)),
                        invLH = as.double(rep(0, n)),
                        moveVec             = as.double(rep(0, numReps)))
                        
                        if(p1 > 0){
                            beta1.p         <- matrix(mcmcRet$samples_beta1, nrow = nStore, byrow = TRUE)
                        }
                        if(p1 == 0){
                            beta1.p         <- NULL
                        }
                        if(p2 > 0){
                            beta2.p         <- matrix(mcmcRet$samples_beta2, nrow = nStore, byrow = TRUE)
                        }
                        if(p2 == 0){
                            beta2.p         <- NULL
                        }
                        if(p3 > 0){
                            beta3.p         <- matrix(mcmcRet$samples_beta3, nrow = nStore, byrow = TRUE)
                        }
                        if(p3 == 0){
                            beta3.p         <- NULL
                        }
                        
                        alpha1.p         <- matrix(mcmcRet$samples_alpha1, nrow = nStore, byrow = TRUE)
                        alpha2.p         <- matrix(mcmcRet$samples_alpha2, nrow = nStore, byrow = TRUE)
                        alpha3.p         <- matrix(mcmcRet$samples_alpha3, nrow = nStore, byrow = TRUE)
                        kappa1.p         <- matrix(mcmcRet$samples_kappa1, nrow = nStore, byrow = TRUE)
                        kappa2.p         <- matrix(mcmcRet$samples_kappa2, nrow = nStore, byrow = TRUE)
                        kappa3.p         <- matrix(mcmcRet$samples_kappa3, nrow = nStore, byrow = TRUE)
                        theta.p         <- matrix(mcmcRet$samples_theta, nrow = nStore, byrow = TRUE)
                        gamma.p         <- matrix(mcmcRet$samples_gamma, nrow = nStore, byrow = TRUE)
                        
                        if(p1 > 0){
                            accept.beta1     <- as.vector(mcmcRet$samples_misc[1:(p1)])/sum(as.vector(mcmcRet$moveVec)==1)*p1
                        }
                        if(p1 == 0){
                            accept.beta1     <- NULL
                        }
                        if(p2 > 0){
                            accept.beta2     <- as.vector(mcmcRet$samples_misc[(p1+1):(p1+p2)])/sum(as.vector(mcmcRet$moveVec)==2)*p2
                        }
                        if(p2 == 0){
                            accept.beta2     <- NULL
                        }
                        if(p3 > 0){
                            accept.beta3     <- as.vector(mcmcRet$samples_misc[(p1+p2+1):(p1+p2+p3)])/sum(as.vector(mcmcRet$moveVec)==3)*p3
                        }
                        if(p3 == 0){
                            accept.beta3     <- NULL
                        }
                        
                        accept.alpha1    <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+1)])/sum(as.vector(mcmcRet$moveVec)==4)
                        accept.alpha2    <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+2)])/sum(as.vector(mcmcRet$moveVec)==5)
                        accept.alpha3    <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+3)])/sum(as.vector(mcmcRet$moveVec)==6)
                        accept.theta    <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+4)])/sum(as.vector(mcmcRet$moveVec)==10)
                        
                        if(nGam_save > 0 & !is.null(path)){
                            save(gamma.p, file = paste(path, "/gammaPch", chain, ".Rdata", sep = ""))
                        }
                        
                        if(p1 > 0){
                            covNames1 = colnames(Xmat1)
                        }
                        if(p1 == 0){
                            covNames1 = NULL
                        }
                        if(p2 > 0){
                            covNames2 = colnames(Xmat2)
                        }
                        if(p2 == 0){
                            covNames2 = NULL
                        }
                        if(p3 > 0){
                            covNames3 = colnames(Xmat3)
                        }
                        if(p3 == 0){
                            covNames3 = NULL
                        }
                        
                        ### posterior predictive checks ###
                        
                        ## 1. log pseudo-marginal likelihood
                        
                        invLH.p <- matrix(mcmcRet$invLH, nrow = n, byrow = TRUE)
                        
                        ## 2. deviance information criterion
                    
                        dev.p        <- matrix(mcmcRet$dev, nrow = nStore, byrow = TRUE)
                        
                        ret[[nam]] <- list(beta1.p = beta1.p, beta2.p = beta2.p, beta3.p = beta3.p, alpha1.p = alpha1.p, alpha2.p = alpha2.p, alpha3.p = alpha3.p, kappa1.p = kappa1.p, kappa2.p = kappa2.p, kappa3.p = kappa3.p, theta.p = theta.p, accept.beta1 = accept.beta1, accept.beta2 = accept.beta2, accept.beta3 = accept.beta3, accept.alpha1 = accept.alpha1, accept.alpha2 = accept.alpha2, accept.alpha3 = accept.alpha3, accept.theta = accept.theta, covNames1 = covNames1, covNames2 = covNames2, covNames3 = covNames3, dev.p=dev.p, invLH.p=invLH.p)
                    } # if(model_h3 == "Markov")
                    
                    # model_h3 = "semi-Markov"
                    
                    if(model_h3 == "semi-Markov"){
                        
                        ###
                        
                        mcmcRet <- .C("BweibScrSMmcmc",
                        survData         = as.double(as.matrix(survData)),
                        n                = as.integer(n),
                        p1                = as.integer(p1),
                        p2                = as.integer(p2),
                        p3                = as.integer(p3),
                        hyperParams     = as.double(hyperParams),
                        mcmcParams        = as.double(mcmcParams),
                        startValues     = as.double(startV),
                        numReps            = as.integer(numReps),
                        thin            = as.integer(thin),
                        burninPerc      = as.double(burninPerc),
                        nGam_save        = as.integer(nGam_save),
                        samples_beta1     = as.double(rep(0, nStore*p1)),
                        samples_beta2     = as.double(rep(0, nStore*p2)),
                        samples_beta3     = as.double(rep(0, nStore*p3)),
                        samples_alpha1     = as.double(rep(0, nStore*1)),
                        samples_alpha2     = as.double(rep(0, nStore*1)),
                        samples_alpha3     = as.double(rep(0, nStore*1)),
                        samples_kappa1     = as.double(rep(0, nStore*1)),
                        samples_kappa2     = as.double(rep(0, nStore*1)),
                        samples_kappa3     = as.double(rep(0, nStore*1)),
                        samples_theta     = as.double(rep(0, nStore*1)),
                        samples_gamma     = as.double(rep(0, nStore*nGam_save)),
                        samples_misc    = as.double(rep(0, (p1+p2+p3+4))),
                        dev                 = as.double(rep(0, nStore*1)),
                        invLH = as.double(rep(0, n)),
                        moveVec             = as.double(rep(0, numReps)))
                        
                        if(p1 > 0){
                            beta1.p         <- matrix(mcmcRet$samples_beta1, nrow = nStore, byrow = TRUE)
                        }
                        if(p1 == 0){
                            beta1.p         <- NULL
                        }
                        if(p2 > 0){
                            beta2.p         <- matrix(mcmcRet$samples_beta2, nrow = nStore, byrow = TRUE)
                        }
                        if(p2 == 0){
                            beta2.p         <- NULL
                        }
                        if(p3 > 0){
                            beta3.p         <- matrix(mcmcRet$samples_beta3, nrow = nStore, byrow = TRUE)
                        }
                        if(p3 == 0){
                            beta3.p         <- NULL
                        }
                        
                        alpha1.p         <- matrix(mcmcRet$samples_alpha1, nrow = nStore, byrow = TRUE)
                        alpha2.p         <- matrix(mcmcRet$samples_alpha2, nrow = nStore, byrow = TRUE)
                        alpha3.p         <- matrix(mcmcRet$samples_alpha3, nrow = nStore, byrow = TRUE)
                        kappa1.p         <- matrix(mcmcRet$samples_kappa1, nrow = nStore, byrow = TRUE)
                        kappa2.p         <- matrix(mcmcRet$samples_kappa2, nrow = nStore, byrow = TRUE)
                        kappa3.p         <- matrix(mcmcRet$samples_kappa3, nrow = nStore, byrow = TRUE)
                        theta.p         <- matrix(mcmcRet$samples_theta, nrow = nStore, byrow = TRUE)
                        gamma.p         <- matrix(mcmcRet$samples_gamma, nrow = nStore, byrow = TRUE)
                        
                        if(p1 > 0){
                            accept.beta1     <- as.vector(mcmcRet$samples_misc[1:(p1)])/sum(as.vector(mcmcRet$moveVec)==1)*p1
                        }
                        if(p1 == 0){
                            accept.beta1     <- NULL
                        }
                        if(p2 > 0){
                            accept.beta2     <- as.vector(mcmcRet$samples_misc[(p1+1):(p1+p2)])/sum(as.vector(mcmcRet$moveVec)==2)*p2
                        }
                        if(p2 == 0){
                            accept.beta2     <- NULL
                        }
                        if(p3 > 0){
                            accept.beta3     <- as.vector(mcmcRet$samples_misc[(p1+p2+1):(p1+p2+p3)])/sum(as.vector(mcmcRet$moveVec)==3)*p3
                        }
                        if(p3 == 0){
                            accept.beta3     <- NULL
                        }
                        
                        accept.alpha1    <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+1)])/sum(as.vector(mcmcRet$moveVec)==4)
                        accept.alpha2    <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+2)])/sum(as.vector(mcmcRet$moveVec)==5)
                        accept.alpha3    <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+3)])/sum(as.vector(mcmcRet$moveVec)==6)
                        accept.theta    <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+4)])/sum(as.vector(mcmcRet$moveVec)==10)
                        
                        if(nGam_save > 0 & !is.null(path)){
                            save(gamma.p, file = paste(path, "/gammaPch", chain, ".Rdata", sep = ""))
                        }
                        
                        if(p1 > 0){
                            covNames1 = colnames(Xmat1)
                        }
                        if(p1 == 0){
                            covNames1 = NULL
                        }
                        if(p2 > 0){
                            covNames2 = colnames(Xmat2)
                        }
                        if(p2 == 0){
                            covNames2 = NULL
                        }
                        if(p3 > 0){
                            covNames3 = colnames(Xmat3)
                        }
                        if(p3 == 0){
                            covNames3 = NULL
                        }
                        
                        ### posterior predictive checks ###
                        
                        ## 1. log pseudo-marginal likelihood
                        
                        invLH.p <- matrix(mcmcRet$invLH, nrow = n, byrow = TRUE)
                        
                        ## 2. deviance information criterion
                        
                        dev.p        <- matrix(mcmcRet$dev, nrow = nStore, byrow = TRUE)
                        
                        ret[[nam]] <- list(beta1.p = beta1.p, beta2.p = beta2.p, beta3.p = beta3.p, alpha1.p = alpha1.p, alpha2.p = alpha2.p, alpha3.p = alpha3.p, kappa1.p = kappa1.p, kappa2.p = kappa2.p, kappa3.p = kappa3.p, theta.p = theta.p, accept.beta1 = accept.beta1, accept.beta2 = accept.beta2, accept.beta3 = accept.beta3, accept.alpha1 = accept.alpha1, accept.alpha2 = accept.alpha2, accept.alpha3 = accept.alpha3, accept.theta = accept.theta, covNames1 = covNames1, covNames2 = covNames2, covNames3 = covNames3, dev.p=dev.p, invLH.p=invLH.p)
                    } # if(model_h3 == "semi-Markov")
                    
                    ## DIC and LPML
                    
                    be1.p <- ret$chain1$beta1.p
                    be2.p <- ret$chain1$beta2.p
                    be3.p <- ret$chain1$beta3.p
                    alp1.p <- ret$chain1$alpha1.p
                    alp2.p <- ret$chain1$alpha2.p
                    alp3.p <- ret$chain1$alpha3.p
                    kap1.p <- ret$chain1$kappa1.p
                    kap2.p <- ret$chain1$kappa2.p
                    kap3.p <- ret$chain1$kappa3.p
                    thet.p <- ret$chain1$theta.p
                    Dev.p <- ret$chain1$dev.p
                    InvLH.p <- ret$chain1$invLH.p
                    
                    if(nChain > 1){
                        for(i in 2:nChain){
                            nam <- paste("chain", i, sep="")
                            be1.p <- rbind(be1.p, ret[[nam]]$beta1.p)
                            be2.p <- rbind(be2.p, ret[[nam]]$beta2.p)
                            be3.p <- rbind(be3.p, ret[[nam]]$beta3.p)
                            alp1.p <- rbind(alp1.p, ret[[nam]]$alpha1.p)
                            alp2.p <- rbind(alp2.p, ret[[nam]]$alpha2.p)
                            alp3.p <- rbind(alp3.p, ret[[nam]]$alpha3.p)
                            kap1.p <- rbind(kap1.p, ret[[nam]]$kappa1.p)
                            kap2.p <- rbind(kap2.p, ret[[nam]]$kappa2.p)
                            kap3.p <- rbind(kap3.p, ret[[nam]]$kappa3.p)
                            thet.p <- rbind(thet.p, ret[[nam]]$theta.p)
                            Dev.p <- rbind(Dev.p, ret[[nam]]$dev.p)
                            InvLH.p <- cbind(InvLH.p, ret[[nam]]$invLH.p)
                        }
                    }
                    if(model_h3 == "Markov")
                    {
                        logLH_fin <- .C("BweibScr_logMLH_DIC",
                        survData         = as.double(as.matrix(survData)),
                        n                = as.integer(n),
                        p1                = as.integer(p1),
                        p2                = as.integer(p2),
                        p3                = as.integer(p3),
                        be1 = as.double(apply(be1.p, 2, mean)),
                        be2 = as.double(apply(be2.p, 2, mean)),
                        be3 = as.double(apply(be3.p, 2, mean)),
                        alp1 = as.double(apply(alp1.p, 2, mean)),
                        alp2 = as.double(apply(alp2.p, 2, mean)),
                        alp3 = as.double(apply(alp3.p, 2, mean)),
                        kap1 = as.double(apply(kap1.p, 2, mean)),
                        kap2 = as.double(apply(kap2.p, 2, mean)),
                        kap3 = as.double(apply(kap3.p, 2, mean)),
                        thet = as.double(apply(theta.p, 2, mean)),
                        val = as.double(0))$val
                    }else if(model_h3 == "semi-Markov")
                    {
                        logLH_fin <- .C("BweibScrSM_logMLH_DIC",
                        survData         = as.double(as.matrix(survData)),
                        n                = as.integer(n),
                        p1                = as.integer(p1),
                        p2                = as.integer(p2),
                        p3                = as.integer(p3),
                        be1 = as.double(apply(be1.p, 2, mean)),
                        be2 = as.double(apply(be2.p, 2, mean)),
                        be3 = as.double(apply(be3.p, 2, mean)),
                        alp1 = as.double(apply(alp1.p, 2, mean)),
                        alp2 = as.double(apply(alp2.p, 2, mean)),
                        alp3 = as.double(apply(alp3.p, 2, mean)),
                        kap1 = as.double(apply(kap1.p, 2, mean)),
                        kap2 = as.double(apply(kap2.p, 2, mean)),
                        kap3 = as.double(apply(kap3.p, 2, mean)),
                        thet = as.double(apply(theta.p, 2, mean)),
                        val = as.double(0))$val
                    }
                    
                    DIC = 2*mean(Dev.p) + 2 * logLH_fin
                    LPML = sum(log(1/apply(InvLH.p, 1, mean)))
                    
                    
                } # if(hz.type == "Weibull")
                
                # hz.type = "PEM"
                
                if(hz.type == "PEM"){
                    
                    C1 <- mcmc[1]
                    C2 <- mcmc[2]
                    C3 <- mcmc[3]
                    delPert1 <- mcmc[4]
                    delPert2 <- mcmc[5]
                    delPert3 <- mcmc[6]
                    rj.scheme <- mcmc[7]
                    K1_max <- mcmc[8]
                    K2_max <- mcmc[9]
                    K3_max <- mcmc[10]
                    s1_max <- max(temp$PEM$PEM.s1)
                    s2_max <- max(temp$PEM$PEM.s2)
                    s3_max <- max(temp$PEM$PEM.s3)
                    mhProp_theta_var <- mcmc[14]
                    nGam_save   <- mcmc[15]
                    numReps     <- mcmc[16]
                    thin        <- mcmc[17]
                    burninPerc  <- mcmc[18]
                    
                    nStore <- round(numReps/thin*(1-burninPerc))
                    
                    ## recommended when the unit of time is day
                    if(rj.scheme == 1){
                        s_propBI1 <- seq(1, s1_max, 1)
                        s_propBI1 <- s_propBI1[s_propBI1 < s1_max]
                        s_propBI2 <- seq(1, s2_max, 1)
                        s_propBI2 <- s_propBI2[s_propBI2 < s2_max]
                        s_propBI3 <- seq(1, s3_max, 1)
                        s_propBI3 <- s_propBI3[s_propBI3 < s3_max]
                    }
                    ## uniquely ordered failure time
                    if(rj.scheme == 2){
                        s_propBI1 <- sort(unique(survData[survData[,2]==1,1]))
                        s_propBI1 <- s_propBI1[s_propBI1 < s1_max]
                        s_propBI2 <- sort(unique(survData[(survData[,2]==0) & (survData[,4]==1),3]))
                        s_propBI2 <- s_propBI2[s_propBI2 < s2_max]
                        s_propBI3 <- sort(unique(survData[(survData[,2]==1) & (survData[,4]==1),3]))
                        s_propBI3 <- s_propBI3[s_propBI3 < s3_max]
                    }
                    
                    
                    num_s_propBI1=length(s_propBI1)
                    num_s_propBI2=length(s_propBI2)
                    num_s_propBI3=length(s_propBI3)
                    
                    time_lambda1 <- mcmcList$tuning$time_lambda1
                    time_lambda2 <- mcmcList$tuning$time_lambda2
                    time_lambda3 <- mcmcList$tuning$time_lambda3
                    
                    nTime_lambda1 <- length(time_lambda1)
                    nTime_lambda2 <- length(time_lambda2)
                    nTime_lambda3 <- length(time_lambda3)
                    
                    mcmcParams <- c(C1, C2, C3, delPert1, delPert2, delPert3, num_s_propBI1, num_s_propBI2, num_s_propBI3, K1_max, K2_max, K3_max, s1_max, s2_max, s3_max, nTime_lambda1, nTime_lambda2, nTime_lambda3, s_propBI1, s_propBI2, s_propBI3, time_lambda1, time_lambda2, time_lambda3, mhProp_theta_var)
                    
                    
                    s1          <- temp$PEM$PEM.s1
                    s2            <- temp$PEM$PEM.s2
                    s3            <- temp$PEM$PEM.s3
                    lambda1 <- temp$PEM$PEM.lambda1
                    lambda2 <- temp$PEM$PEM.lambda2
                    lambda3 <- temp$PEM$PEM.lambda3
                    
                    K1=temp$PEM$K1
                    K2=temp$PEM$K2
                    K3=temp$PEM$K3
                    mu_lam1=temp$PEM$PEM.mu_lam[1]
                    mu_lam2=temp$PEM$PEM.mu_lam[2]
                    mu_lam3=temp$PEM$PEM.mu_lam[3]
                    sigSq_lam1=temp$PEM$PEM.sigSq_lam[1]
                    sigSq_lam2=temp$PEM$PEM.sigSq_lam[2]
                    sigSq_lam3=temp$PEM$PEM.sigSq_lam[3]
                    
                    
                    # model_h3 = "Markov"
                    
                    if(model_h3 == "Markov"){
                        
                        ###
                        J_1 <- K1
                        J_2 <- K2
                        J_3 <- K3
                        
                        theta <- startV[(p1+p2+p3+1)]
                        gamma <- startV[(p1+p2+p3+2):(p1+p2+p3+n+1)]
                        
                        if(p1+p2+p3 > 0)
                        {
                            startV <- as.vector(c(startV[1:(p1+p2+p3)],
                            K1, K2, K3,
                            mu_lam1, mu_lam2, mu_lam3,
                            sigSq_lam1, sigSq_lam2, sigSq_lam3,
                            theta, gamma,
                            lambda1, lambda2, lambda3, s1, s2, s3))
                        }else
                        {
                            startV <- as.vector(c(K1, K2, K3,
                            mu_lam1, mu_lam2, mu_lam3,
                            sigSq_lam1, sigSq_lam2, sigSq_lam3,
                            theta, gamma,
                            lambda1, lambda2, lambda3, s1, s2, s3))
                        }
                        
                        mcmcRet <- .C("BpeScrmcmc",
                        survData            = as.double(as.matrix(survData)),
                        n                   = as.integer(n),
                        p1                    = as.integer(p1),
                        p2                    = as.integer(p2),
                        p3                    = as.integer(p3),
                        hyperParams         = as.double(hyperParams),
                        startValues         = as.double(startV),
                        mcmcParams          = as.double(mcmcParams),
                        numReps             = as.integer(numReps),
                        thin                = as.integer(thin),
                        burninPerc          = as.double(burninPerc),
                        nGam_save            = as.integer(nGam_save),
                        samples_beta1         = as.double(rep(0, nStore*p1)),
                        samples_beta2         = as.double(rep(0, nStore*p2)),
                        samples_beta3         = as.double(rep(0, nStore*p3)),
                        samples_mu_lam1     = as.double(rep(0, nStore*1)),
                        samples_mu_lam2     = as.double(rep(0, nStore*1)),
                        samples_mu_lam3     = as.double(rep(0, nStore*1)),
                        samples_sigSq_lam1    = as.double(rep(0, nStore*1)),
                        samples_sigSq_lam2    = as.double(rep(0, nStore*1)),
                        samples_sigSq_lam3    = as.double(rep(0, nStore*1)),
                        samples_J1          = as.double(rep(0, nStore*1)),
                        samples_J2          = as.double(rep(0, nStore*1)),
                        samples_J3          = as.double(rep(0, nStore*1)),
                        samples_s1          = as.double(rep(0, nStore*(K1_max + 1))),
                        samples_s2          = as.double(rep(0, nStore*(K2_max + 1))),
                        samples_s3          = as.double(rep(0, nStore*(K3_max + 1))),
                        samples_theta         = as.double(rep(0, nStore*1)),
                        samples_gamma         = as.double(rep(0, nStore*nGam_save)),
                        samples_misc        = as.double(rep(0, p1 + p2 + p3 + 7)),
                        lambda1_fin            = as.double(rep(0, nStore*nTime_lambda1)),
                        lambda2_fin            = as.double(rep(0, nStore*nTime_lambda2)),
                        lambda3_fin            = as.double(rep(0, nStore*nTime_lambda3)),
                        dev                 = as.double(rep(0, nStore*1)),
                        invLH = as.double(rep(0, n)),
                        moveVec             = as.double(rep(0, numReps)))
                        
                        if(p1 > 0){
                            beta1.p         <- matrix(mcmcRet$samples_beta1, nrow = nStore, byrow = TRUE)
                        }
                        if(p1 == 0){
                            beta1.p         <- NULL
                        }
                        if(p2 > 0){
                            beta2.p         <- matrix(mcmcRet$samples_beta2, nrow = nStore, byrow = TRUE)
                        }
                        if(p2 == 0){
                            beta2.p         <- NULL
                        }
                        if(p3 > 0){
                            beta3.p         <- matrix(mcmcRet$samples_beta3, nrow = nStore, byrow = TRUE)
                        }
                        if(p3 == 0){
                            beta3.p         <- NULL
                        }
                        
                        lambda1.fin     <- matrix(mcmcRet$lambda1_fin, nrow = nStore, byrow = TRUE)
                        lambda2.fin     <- matrix(mcmcRet$lambda2_fin, nrow = nStore, byrow = TRUE)
                        lambda3.fin     <- matrix(mcmcRet$lambda3_fin, nrow = nStore, byrow = TRUE)
                        gamma.p         <- matrix(mcmcRet$samples_gamma, nrow = nStore, byrow = TRUE)
                        theta.p         <- matrix(mcmcRet$samples_theta, nrow = nStore, byrow = TRUE)
                        mu_lam1.p         <- matrix(mcmcRet$samples_mu_lam1, nrow = nStore, byrow = TRUE)
                        mu_lam2.p         <- matrix(mcmcRet$samples_mu_lam2, nrow = nStore, byrow = TRUE)
                        mu_lam3.p         <- matrix(mcmcRet$samples_mu_lam3, nrow = nStore, byrow = TRUE)
                        sigSq_lam1.p     <- matrix(mcmcRet$samples_sigSq_lam1, nrow = nStore, byrow = TRUE)
                        sigSq_lam2.p     <- matrix(mcmcRet$samples_sigSq_lam2, nrow = nStore, byrow = TRUE)
                        sigSq_lam3.p     <- matrix(mcmcRet$samples_sigSq_lam3, nrow = nStore, byrow = TRUE)
                        K1.p             <- matrix(mcmcRet$samples_J1, nrow = nStore, byrow = TRUE)
                        K2.p             <- matrix(mcmcRet$samples_J2, nrow = nStore, byrow = TRUE)
                        K3.p             <- matrix(mcmcRet$samples_J3, nrow = nStore, byrow = TRUE)
                        s1.p             <- matrix(mcmcRet$samples_s1, nrow = nStore, byrow = TRUE)
                        s2.p             <- matrix(mcmcRet$samples_s2, nrow = nStore, byrow = TRUE)
                        s3.p             <- matrix(mcmcRet$samples_s3, nrow = nStore, byrow = TRUE)
                        
                        if(p1 > 0){
                            accept.beta1     <- as.vector(mcmcRet$samples_misc[1:(p1)])/sum(as.vector(mcmcRet$moveVec)==1)*p1
                        }
                        if(p1 == 0){
                            accept.beta1     <- NULL
                        }
                        if(p2 > 0){
                            accept.beta2     <- as.vector(mcmcRet$samples_misc[(p1+1):(p1+p2)])/sum(as.vector(mcmcRet$moveVec)==2)*p2
                        }
                        if(p2 == 0){
                            accept.beta2     <- NULL
                        }
                        if(p3 > 0){
                            accept.beta3     <- as.vector(mcmcRet$samples_misc[(p1+p2+1):(p1+p2+p3)])/sum(as.vector(mcmcRet$moveVec)==3)*p3
                        }
                        if(p3 == 0){
                            accept.beta3     <- NULL
                        }
                        accept.BI1        <- as.vector(mcmcRet$samples_misc[(p1+p2+p3)+1])/sum(as.vector(mcmcRet$moveVec)==12)
                        accept.DI1        <- as.vector(mcmcRet$samples_misc[(p1+p2+p3)+2])/sum(as.vector(mcmcRet$moveVec)==15)
                        accept.BI2        <- as.vector(mcmcRet$samples_misc[(p1+p2+p3)+3])/sum(as.vector(mcmcRet$moveVec)==13)
                        accept.DI2        <- as.vector(mcmcRet$samples_misc[(p1+p2+p3)+4])/sum(as.vector(mcmcRet$moveVec)==16)
                        accept.BI3        <- as.vector(mcmcRet$samples_misc[(p1+p2+p3)+5])/sum(as.vector(mcmcRet$moveVec)==14)
                        accept.DI3        <- as.vector(mcmcRet$samples_misc[(p1+p2+p3)+6])/sum(as.vector(mcmcRet$moveVec)==17)
                        accept.theta    <- as.vector(mcmcRet$samples_misc[(p1+p2+p3)+7])/sum(as.vector(mcmcRet$moveVec)==11)
                        
                        if(nGam_save > 0 & !is.null(path)){
                            save(gamma.p, file = paste(path, "/gammaPch", chain, ".Rdata", sep = ""))
                        }
                        
                        if(p1 > 0){
                            covNames1 = colnames(Xmat1)
                        }
                        if(p1 == 0){
                            covNames1 = NULL
                        }
                        
                        if(p2 > 0){
                            covNames2 = colnames(Xmat2)
                        }
                        if(p2 == 0){
                            covNames2 = NULL
                        }
                        if(p3 > 0){
                            covNames3 = colnames(Xmat3)
                        }
                        if(p3 == 0){
                            covNames3 = NULL
                        }
                        
                        ### posterior predictive checks ###
                        
                        ## 1. log pseudo-marginal likelihood
                        
                        invLH.p <- matrix(mcmcRet$invLH, nrow = n, byrow = TRUE)
                        
                        ## 2. deviance information criterion
                        
                        dev.p        <- matrix(mcmcRet$dev, nrow = nStore, byrow = TRUE)
                        
                        ret[[nam]] <- list(beta1.p = beta1.p, beta2.p = beta2.p, beta3.p = beta3.p, lambda1.fin = lambda1.fin, lambda2.fin = lambda2.fin, lambda3.fin = lambda3.fin, mu_lam1.p = mu_lam1.p, mu_lam2.p = mu_lam2.p, mu_lam3.p = mu_lam3.p, sigSq_lam1.p = sigSq_lam1.p, sigSq_lam2.p = sigSq_lam2.p, sigSq_lam3.p = sigSq_lam3.p, theta.p = theta.p, K1.p = K1.p, K2.p = K2.p, K3.p = K3.p, s1.p = s1.p, s2.p = s2.p, s3.p = s3.p, accept.beta1 = accept.beta1, accept.beta2 = accept.beta2, accept.beta3 = accept.beta3, accept.BI1 = accept.BI1, accept.BI2 = accept.BI2, accept.BI3 = accept.BI3, accept.DI1 = accept.DI1, accept.DI2 = accept.DI2, accept.DI3 = accept.DI3, accept.theta = accept.theta, time_lambda1 = time_lambda1, time_lambda2 = time_lambda2, time_lambda3 = time_lambda3, covNames1 = covNames1, covNames2 = covNames2, covNames3 = covNames3, dev.p=dev.p, invLH.p=invLH.p)
                        
                    } #if(model_h3 == "Markov")
                    
                    # model_h3 = "semi-Markov"
                    
                    if(model_h3 == "semi-Markov"){
                        
                        ###
                        
                        J_1 <- K1
                        J_2 <- K2
                        J_3 <- K3
                        
                        theta <- startV[(p1+p2+p3+1)]
                        gamma <- startV[(p1+p2+p3+2):(p1+p2+p3+n+1)]
                        
                        if(p1+p2+p3 > 0)
                        {
                            startV <- as.vector(c(startV[1:(p1+p2+p3)],
                            K1, K2, K3,
                            mu_lam1, mu_lam2, mu_lam3,
                            sigSq_lam1, sigSq_lam2, sigSq_lam3,
                            theta, gamma,
                            lambda1, lambda2, lambda3, s1, s2, s3))
                        }else
                        {
                            startV <- as.vector(c(K1, K2, K3,
                            mu_lam1, mu_lam2, mu_lam3,
                            sigSq_lam1, sigSq_lam2, sigSq_lam3,
                            theta, gamma,
                            lambda1, lambda2, lambda3, s1, s2, s3))
                        }
                        
                        mcmcRet <- .C("BpeScrSMmcmc",
                        survData            = as.double(as.matrix(survData)),
                        n                   = as.integer(n),
                        p1                    = as.integer(p1),
                        p2                    = as.integer(p2),
                        p3                    = as.integer(p3),
                        hyperParams         = as.double(hyperParams),
                        startValues         = as.double(startV),
                        mcmcParams          = as.double(mcmcParams),
                        numReps             = as.integer(numReps),
                        thin                = as.integer(thin),
                        burninPerc          = as.double(burninPerc),
                        nGam_save            = as.integer(nGam_save),
                        samples_beta1         = as.double(rep(0, nStore*p1)),
                        samples_beta2         = as.double(rep(0, nStore*p2)),
                        samples_beta3         = as.double(rep(0, nStore*p3)),
                        samples_mu_lam1     = as.double(rep(0, nStore*1)),
                        samples_mu_lam2     = as.double(rep(0, nStore*1)),
                        samples_mu_lam3     = as.double(rep(0, nStore*1)),
                        samples_sigSq_lam1    = as.double(rep(0, nStore*1)),
                        samples_sigSq_lam2    = as.double(rep(0, nStore*1)),
                        samples_sigSq_lam3    = as.double(rep(0, nStore*1)),
                        samples_J1          = as.double(rep(0, nStore*1)),
                        samples_J2          = as.double(rep(0, nStore*1)),
                        samples_J3          = as.double(rep(0, nStore*1)),
                        samples_s1          = as.double(rep(0, nStore*(K1_max + 1))),
                        samples_s2          = as.double(rep(0, nStore*(K2_max + 1))),
                        samples_s3          = as.double(rep(0, nStore*(K3_max + 1))),
                        samples_theta         = as.double(rep(0, nStore*1)),
                        samples_gamma         = as.double(rep(0, nStore*nGam_save)),
                        samples_misc        = as.double(rep(0, p1 + p2 + p3 + 7)),
                        lambda1_fin            = as.double(rep(0, nStore*nTime_lambda1)),
                        lambda2_fin            = as.double(rep(0, nStore*nTime_lambda2)),
                        lambda3_fin            = as.double(rep(0, nStore*nTime_lambda3)),
                        dev                 = as.double(rep(0, nStore*1)),
                        invLH = as.double(rep(0, n)),
                        moveVec             = as.double(rep(0, numReps)))
                        
                        if(p1 > 0){
                            beta1.p         <- matrix(mcmcRet$samples_beta1, nrow = nStore, byrow = TRUE)
                        }
                        if(p1 == 0){
                            beta1.p         <- NULL
                        }
                        if(p2 > 0){
                            beta2.p         <- matrix(mcmcRet$samples_beta2, nrow = nStore, byrow = TRUE)
                        }
                        if(p2 == 0){
                            beta2.p         <- NULL
                        }
                        if(p3 > 0){
                            beta3.p         <- matrix(mcmcRet$samples_beta3, nrow = nStore, byrow = TRUE)
                        }
                        if(p3 == 0){
                            beta3.p         <- NULL
                        }
                        
                        lambda1.fin     <- matrix(mcmcRet$lambda1_fin, nrow = nStore, byrow = TRUE)
                        lambda2.fin     <- matrix(mcmcRet$lambda2_fin, nrow = nStore, byrow = TRUE)
                        lambda3.fin     <- matrix(mcmcRet$lambda3_fin, nrow = nStore, byrow = TRUE)
                        gamma.p         <- matrix(mcmcRet$samples_gamma, nrow = nStore, byrow = TRUE)
                        theta.p         <- matrix(mcmcRet$samples_theta, nrow = nStore, byrow = TRUE)
                        mu_lam1.p         <- matrix(mcmcRet$samples_mu_lam1, nrow = nStore, byrow = TRUE)
                        mu_lam2.p         <- matrix(mcmcRet$samples_mu_lam2, nrow = nStore, byrow = TRUE)
                        mu_lam3.p         <- matrix(mcmcRet$samples_mu_lam3, nrow = nStore, byrow = TRUE)
                        sigSq_lam1.p     <- matrix(mcmcRet$samples_sigSq_lam1, nrow = nStore, byrow = TRUE)
                        sigSq_lam2.p     <- matrix(mcmcRet$samples_sigSq_lam2, nrow = nStore, byrow = TRUE)
                        sigSq_lam3.p     <- matrix(mcmcRet$samples_sigSq_lam3, nrow = nStore, byrow = TRUE)
                        K1.p             <- matrix(mcmcRet$samples_J1, nrow = nStore, byrow = TRUE)
                        K2.p             <- matrix(mcmcRet$samples_J2, nrow = nStore, byrow = TRUE)
                        K3.p             <- matrix(mcmcRet$samples_J3, nrow = nStore, byrow = TRUE)
                        s1.p             <- matrix(mcmcRet$samples_s1, nrow = nStore, byrow = TRUE)
                        s2.p             <- matrix(mcmcRet$samples_s2, nrow = nStore, byrow = TRUE)
                        s3.p             <- matrix(mcmcRet$samples_s3, nrow = nStore, byrow = TRUE)
                        
                        if(p1 > 0){
                            accept.beta1     <- as.vector(mcmcRet$samples_misc[1:(p1)])/sum(as.vector(mcmcRet$moveVec)==1)*p1
                        }
                        if(p1 == 0){
                            accept.beta1     <- NULL
                        }
                        if(p2 > 0){
                            accept.beta2     <- as.vector(mcmcRet$samples_misc[(p1+1):(p1+p2)])/sum(as.vector(mcmcRet$moveVec)==2)*p2
                        }
                        if(p2 == 0){
                            accept.beta2     <- NULL
                        }
                        if(p3 > 0){
                            accept.beta3     <- as.vector(mcmcRet$samples_misc[(p1+p2+1):(p1+p2+p3)])/sum(as.vector(mcmcRet$moveVec)==3)*p3
                        }
                        if(p3 == 0){
                            accept.beta3     <- NULL
                        }
                        accept.BI1        <- as.vector(mcmcRet$samples_misc[(p1+p2+p3)+1])/sum(as.vector(mcmcRet$moveVec)==12)
                        accept.DI1        <- as.vector(mcmcRet$samples_misc[(p1+p2+p3)+2])/sum(as.vector(mcmcRet$moveVec)==15)
                        accept.BI2        <- as.vector(mcmcRet$samples_misc[(p1+p2+p3)+3])/sum(as.vector(mcmcRet$moveVec)==13)
                        accept.DI2        <- as.vector(mcmcRet$samples_misc[(p1+p2+p3)+4])/sum(as.vector(mcmcRet$moveVec)==16)
                        accept.BI3        <- as.vector(mcmcRet$samples_misc[(p1+p2+p3)+5])/sum(as.vector(mcmcRet$moveVec)==14)
                        accept.DI3        <- as.vector(mcmcRet$samples_misc[(p1+p2+p3)+6])/sum(as.vector(mcmcRet$moveVec)==17)
                        accept.theta    <- as.vector(mcmcRet$samples_misc[(p1+p2+p3)+7])/sum(as.vector(mcmcRet$moveVec)==11)
                        
                        
                        if(nGam_save > 0 & !is.null(path)){
                            save(gamma.p, file = paste(path, "/gammaPch", chain, ".Rdata", sep = ""))
                        }
                        
                        if(p1 > 0){
                            covNames1 = colnames(Xmat1)
                        }
                        if(p1 == 0){
                            covNames1 = NULL
                        }
                        
                        if(p2 > 0){
                            covNames2 = colnames(Xmat2)
                        }
                        if(p2 == 0){
                            covNames2 = NULL
                        }
                        if(p3 > 0){
                            covNames3 = colnames(Xmat3)
                        }
                        if(p3 == 0){
                            covNames3 = NULL
                        }
                        
                        ### posterior predictive checks ###
                        
                        ## 1. log pseudo-marginal likelihood
                        
                        invLH.p <- matrix(mcmcRet$invLH, nrow = n, byrow = TRUE)
                        
                        ## 2. deviance information criterion
                        
                        dev.p        <- matrix(mcmcRet$dev, nrow = nStore, byrow = TRUE)
                        
                        
                        ret[[nam]] <- list(beta1.p = beta1.p, beta2.p = beta2.p, beta3.p = beta3.p, lambda1.fin = lambda1.fin, lambda2.fin = lambda2.fin, lambda3.fin = lambda3.fin, mu_lam1.p = mu_lam1.p, mu_lam2.p = mu_lam2.p, mu_lam3.p = mu_lam3.p, sigSq_lam1.p = sigSq_lam1.p, sigSq_lam2.p = sigSq_lam2.p, sigSq_lam3.p = sigSq_lam3.p, theta.p = theta.p, K1.p = K1.p, K2.p = K2.p, K3.p = K3.p, s1.p = s1.p, s2.p = s2.p, s3.p = s3.p, accept.beta1 = accept.beta1, accept.beta2 = accept.beta2, accept.beta3 = accept.beta3, accept.BI1 = accept.BI1, accept.BI2 = accept.BI2, accept.BI3 = accept.BI3, accept.DI1 = accept.DI1, accept.DI2 = accept.DI2, accept.DI3 = accept.DI3, accept.theta = accept.theta, time_lambda1 = time_lambda1, time_lambda2 = time_lambda2, time_lambda3 = time_lambda3, covNames1 = covNames1, covNames2 = covNames2, covNames3 = covNames3, dev.p=dev.p, invLH.p=invLH.p)
                        
                    }    # if(model_h3 == "semi-Markov")
                    
                    ## DIC and LPML
                    
                    be1.p <- ret$chain1$beta1.p
                    be2.p <- ret$chain1$beta2.p
                    be3.p <- ret$chain1$beta3.p
                    thet.p <- ret$chain1$theta.p
                    lam1.p <- ret$chain1$lambda1.fin
                    lam2.p <- ret$chain1$lambda2.fin
                    lam3.p <- ret$chain1$lambda3.fin
                    Dev.p <- ret$chain1$dev.p
                    InvLH.p <- ret$chain1$invLH.p
                    
                    if(nChain > 1){
                        for(i in 2:nChain){
                            nam <- paste("chain", i, sep="")
                            be1.p <- rbind(be1.p, ret[[nam]]$beta1.p)
                            be2.p <- rbind(be2.p, ret[[nam]]$beta2.p)
                            be3.p <- rbind(be3.p, ret[[nam]]$beta3.p)
                            thet.p <- rbind(thet.p, ret[[nam]]$theta.p)
                            lam1.p <- rbind(lam1.p, ret[[nam]]$lambda1.fin)
                            lam2.p <- rbind(lam2.p, ret[[nam]]$lambda2.fin)
                            lam3.p <- rbind(lam3.p, ret[[nam]]$lambda3.fin)
                            Dev.p <- rbind(Dev.p, ret[[nam]]$dev.p)
                            InvLH.p <- cbind(InvLH.p, ret[[nam]]$invLH.p)
                        }
                    }
                    if(model_h3 == "Markov")
                    {
                        logLH_fin <- .C("BpeScr_logMLH_DIC",
                        survData         = as.double(as.matrix(survData)),
                        n                = as.integer(n),
                        p1                = as.integer(p1),
                        p2                = as.integer(p2),
                        p3                = as.integer(p3),
                        be1 = as.double(apply(be1.p, 2, mean)),
                        be2 = as.double(apply(be2.p, 2, mean)),
                        be3 = as.double(apply(be3.p, 2, mean)),
                        thet = as.double(apply(theta.p, 2, mean)),
                        lam1 = as.double(apply(lam1.p, 2, mean)),
                        lam2 = as.double(apply(lam2.p, 2, mean)),
                        lam3 = as.double(apply(lam3.p, 2, mean)),
                        s_1 = as.double(time_lambda1),
                        s_2 = as.double(time_lambda2),
                        s_3 = as.double(time_lambda3),
                        K_1                = as.integer(nTime_lambda1),
                        K_2                = as.integer(nTime_lambda2),
                        K_3                = as.integer(nTime_lambda3),
                        val = as.double(0))$val
                    }else if(model_h3 == "semi-Markov")
                    {
                        logLH_fin <- .C("BpeScrSM_logMLH_DIC",
                        survData         = as.double(as.matrix(survData)),
                        n                = as.integer(n),
                        p1                = as.integer(p1),
                        p2                = as.integer(p2),
                        p3                = as.integer(p3),
                        be1 = as.double(apply(be1.p, 2, mean)),
                        be2 = as.double(apply(be2.p, 2, mean)),
                        be3 = as.double(apply(be3.p, 2, mean)),
                        thet = as.double(apply(theta.p, 2, mean)),
                        lam1 = as.double(apply(lam1.p, 2, mean)),
                        lam2 = as.double(apply(lam2.p, 2, mean)),
                        lam3 = as.double(apply(lam3.p, 2, mean)),
                        s_1 = as.double(time_lambda1),
                        s_2 = as.double(time_lambda2),
                        s_3 = as.double(time_lambda3),
                        K_1                = as.integer(nTime_lambda1),
                        K_2                = as.integer(nTime_lambda2),
                        K_3                = as.integer(nTime_lambda3),
                        val = as.double(0))$val
                    }
                    
                    DIC = 2*mean(Dev.p) + 2 * logLH_fin
                    LPML = sum(log(1/apply(InvLH.p, 1, mean)))
                    
                } # if(hz.type == "PEM")
                
                chain = chain + 1
                
            }    # while(chain <= nChain)
            
            ret[["setup"]]    <- list(Formula=Formula, nCov = nCov, hyperParams = hyperParams, startValues = startValues, mcmcParams = mcmcParams, nGam_save = nGam_save, numReps = numReps, thin = thin, path = path, burninPerc = burninPerc, hz.type = hz.type, model = model_h3, nChain = nChain)
            
            ret[["DIC"]] <- DIC
            ret[["LPML"]] <- LPML            
            
            if(hz.type == "Weibull")
            {
                ret$class <- c("Bayes_HReg", "ID", "Ind", "WB")
            }
            if(hz.type == "PEM")
            {
                ret$class <- c("Bayes_HReg", "ID", "Ind", "PEM")
            }
            
            class(ret) <- "Bayes_HReg"
            
            return(ret)
            
        }
        
        
        
        
        ## for cluster-correlated semi-competing risks data
        
        if(!is.null(id))
        {
            nChain <- length(startValues)
            
            model_h3     <- model[1]
            hz.type     <- model[2]
            re.type     <- model[3]
            
            ##
            survData <- cbind(Y, id, Xmat1, Xmat2, Xmat3)
            
            n    <- dim(survData)[1]
            
            J    <- length(unique(id))
            
            nj    <- rep(NA, J)
            
            for(i in 1:J){
                nj[i]    <- length(which(id == i))
            }
            
            
            #if(class(startValues) == "list" & length(startValues) == nChain){
            
            if(!is.null(path)){
                dir.create(paste(path), recursive = TRUE, showWarnings = FALSE)
            }
            
            
            
            ### setting hyperparameters
            
            if(hz.type == "Weibull" & re.type == "MVN")
            {
                hyperParams <- as.vector(c(hyperP$WB$WB.ab1, hyperP$WB$WB.ab2, hyperP$WB$WB.ab3, hyperP$WB$WB.cd1, hyperP$WB$WB.cd2, hyperP$WB$WB.cd3, hyperP$theta, hyperP$MVN$rho_v, hyperP$MVN$Psi_v))
            }
            if(hz.type == "Weibull" & re.type == "DPM")
            {
                hyperParams <- as.vector(c(hyperP$WB$WB.ab1, hyperP$WB$WB.ab2, hyperP$WB$WB.ab3, hyperP$WB$WB.cd1, hyperP$WB$WB.cd2, hyperP$WB$WB.cd3, hyperP$theta, rep(0, 3), 1, hyperP$DPM$Psi0, hyperP$DPM$rho0, hyperP$DPM$aTau, hyperP$DPM$bTau))
            }
            if(hz.type == "PEM" & re.type == "MVN")
            {
                hyperParams <- as.vector(c(hyperP$PEM$PEM.ab1, hyperP$PEM$PEM.ab2, hyperP$PEM$PEM.ab3, hyperP$PEM$PEM.alpha1, hyperP$PEM$PEM.alpha2, hyperP$PEM$PEM.alpha3, 1, 1, 1, hyperP$theta, hyperP$MVN$rho_v, hyperP$MVN$Psi_v))
            }
            if(hz.type == "PEM" & re.type == "DPM")
            {
                hyperParams <- as.vector(c(hyperP$PEM$PEM.ab1, hyperP$PEM$PEM.ab2, hyperP$PEM$PEM.ab3, hyperP$PEM$PEM.alpha1, hyperP$PEM$PEM.alpha2, hyperP$PEM$PEM.alpha3, 1, 1, 1, hyperP$theta, rep(0, 3), 1, hyperP$DPM$Psi0, hyperP$DPM$rho0, hyperP$DPM$aTau, hyperP$DPM$bTau))
            }
            
            ### mcmc setting
            
            if(hz.type == "Weibull")
            {
                mcmc <- as.vector(c(mhProp_alpha1_var=mcmcList$tuning$mhProp_alphag_var[1], mhProp_alpha2_var=mcmcList$tuning$mhProp_alphag_var[2], mhProp_alpha3_var=mcmcList$tuning$mhProp_alphag_var[3], mhProp_theta_var=mcmcList$tuning$mhProp_theta_var, mhProp_V1_var=mcmcList$tuning$mhProp_Vg_var[1], mhProp_V2_var=mcmcList$tuning$mhProp_Vg_var[2], mhProp_V3_var=mcmcList$tuning$mhProp_Vg_var[3], nGam_save=mcmcList$storage$nGam_save, numReps=mcmcList$run$numReps, thin=mcmcList$run$thin, burninPerc=mcmcList$run$burninPerc, storeV=mcmcList$storage$storeV))
            }
            if(hz.type == "PEM")
            {
                mcmcList$tuning$sg_max <- c(1,1,1)
                mcmc <- as.vector(c(C1=mcmcList$tuning$Cg[1], C2=mcmcList$tuning$Cg[2], C3=mcmcList$tuning$Cg[3], delPert1=mcmcList$tuning$delPertg[1], delPert2=mcmcList$tuning$delPertg[2], delPert3=mcmcList$tuning$delPertg[3], rj.scheme = mcmcList$tuning$rj.scheme, K1_max=mcmcList$tuning$Kg_max[1], K2_max=mcmcList$tuning$Kg_max[2], K3_max=mcmcList$tuning$Kg_max[3], s1_max=mcmcList$tuning$sg_max[1], s2_max=mcmcList$tuning$sg_max[2], s3_max=mcmcList$tuning$sg_max[3], mhProp_theta_var=mcmcList$tuning$mhProp_theta_var, mhProp_V1_var=mcmcList$tuning$mhProp_Vg_var[1], mhProp_V2_var=mcmcList$tuning$mhProp_Vg_var[2], mhProp_V3_var=mcmcList$tuning$mhProp_Vg_var[3], nGam_save=mcmcList$storage$nGam_save, numReps=mcmcList$run$numReps, thin=mcmcList$run$thin, burninPerc=mcmcList$run$burninPerc, storeV=mcmcList$storage$storeV))
            }
            
            chain = 1
            ret <- list()
            
            while(chain <= nChain){
                
                cat("chain: ", chain, "\n")
                nam = paste("chain", chain, sep="")
                
                temp <- startValues[[chain]]
                
                ### setting starting values
                
                if(hz.type == "Weibull" & re.type == "MVN")
                {
                    startV <- as.vector(c(beta1=temp$common$beta1, beta2=temp$common$beta2, beta3=temp$common$beta3, temp$WB$WB.alpha, temp$WB$WB.kappa, theta=temp$common$theta, gamma=temp$common$gamma.ji, V1=temp$common$V.j1, V2=temp$common$V.j2, V3=temp$common$V.j3, Sigma_V=temp$MVN$MVN.SigmaV))
                }
                if(hz.type == "Weibull" & re.type == "DPM")
                {
                    startV <- as.vector(c(beta1=temp$common$beta1, beta2=temp$common$beta2, beta3=temp$common$beta3, temp$WB$WB.alpha, temp$WB$WB.kappa, temp$common$theta, gamma=temp$common$gamma.ji, V1=temp$common$V.j1, V2=temp$common$V.j2, V3=temp$common$V.j3, class=temp$DPM$DPM.class, tau=temp$DPM$DPM.tau))
                }
                if(hz.type == "PEM" & re.type == "MVN")
                {
                    startV <- as.vector(c(beta1=temp$common$beta1, beta2=temp$common$beta2, beta3=temp$common$beta3, theta=temp$common$theta, gamma=temp$common$gamma.ji, V1=temp$common$V.j1, V2=temp$common$V.j2, V3=temp$common$V.j3, Sigma_V=temp$MVN$MVN.SigmaV))
                }
                if(hz.type == "PEM" & re.type == "DPM")
                {
                    startV <- as.vector(c(beta1=temp$common$beta1, beta2=temp$common$beta2, beta3=temp$common$beta3, theta=temp$common$theta, gamma=temp$common$gamma.ji, V1=temp$common$V.j1, V2=temp$common$V.j2, V3=temp$common$V.j3, class=temp$DPM$DPM.class, tau=temp$DPM$DPM.tau))
                }
                
                
                # hz.type = "Weibull"
                
                if(hz.type == "Weibull"){
                    
                    nGam_save   <- mcmc[8]
                    numReps     <- mcmc[9]
                    thin        <- mcmc[10]
                    burninPerc  <- mcmc[11]
                    storeV      <- mcmc[12:14]
                    
                    nStore <- round(numReps/thin*(1-burninPerc))
                    
                    # re.type = "MVN"
                    
                    if(re.type == "MVN"){
                        
                        # model = "Markov"
                        
                        #######################################
                        ############ Weibull-MVN-M ############
                        #######################################
                        
                        if(model_h3 == "Markov"){
                            
                            ###
                            
                            startValues1 <- startV[1:(p1+p2+p3+n+7)]
                            startValues2 <- startV[(p1+p2+p3+n+8):length(startV)]
                            startV <- c(startValues1, 1, 1, startValues2)
                            
                            mcmc1 <- mcmc[1:4]
                            mcmc2 <- mcmc[5:7]
                            
                            mcmcNew <- c(mcmc1, 1, mcmc2, n)
                            
                            mcmcRet <- .C("BweibMvnCorScrmcmc",
                            survData         = as.double(as.matrix(survData)),
                            n                = as.integer(n),
                            p1                = as.integer(p1),
                            p2                = as.integer(p2),
                            p3                = as.integer(p3),
                            J                = as.integer(J),
                            nj                = as.double(nj),
                            hyperParams     = as.double(hyperParams),
                            mcmc        = as.double(mcmcNew),
                            startValues     = as.double(startV),
                            numReps            = as.integer(numReps),
                            thin            = as.integer(thin),
                            burninPerc      = as.double(burninPerc),
                            nGam_save        = as.integer(nGam_save),
                            samples_beta1     = as.double(rep(0, nStore*p1)),
                            samples_beta2     = as.double(rep(0, nStore*p2)),
                            samples_beta3     = as.double(rep(0, nStore*p3)),
                            samples_alpha1     = as.double(rep(0, nStore*1)),
                            samples_alpha2     = as.double(rep(0, nStore*1)),
                            samples_alpha3     = as.double(rep(0, nStore*1)),
                            samples_kappa1     = as.double(rep(0, nStore*1)),
                            samples_kappa2     = as.double(rep(0, nStore*1)),
                            samples_kappa3     = as.double(rep(0, nStore*1)),
                            samples_nu2     = as.double(rep(0, nStore*1)),
                            samples_nu3     = as.double(rep(0, nStore*1)),
                            samples_theta     = as.double(rep(0, nStore*1)),
                            samples_V1        = as.double(rep(0, nStore*J)),
                            samples_V2        = as.double(rep(0, nStore*J)),
                            samples_V3        = as.double(rep(0, nStore*J)),
                            samples_Sigma_V    = as.double(rep(0, nStore*3*3)),
                            samples_gamma     = as.double(rep(0, nStore*nGam_save)),
                            samples_misc    = as.double(rep(0, (p1+p2+p3+6+n+1+n+J+J+J))),
                            gammaP            = as.double(rep(0, n)),
                            dev                 = as.double(rep(0, nStore*1)),
                            invLH = as.double(rep(0, n)),
                            logLH_fin = as.double(0),
                            lpml                 = as.double(rep(0, nStore*1)),
                            lpml2                 = as.double(rep(0, nStore*1)),
                            moveVec             = as.double(rep(0, numReps)))
                            
                            if(p1 > 0){
                                beta1.p         <- matrix(mcmcRet$samples_beta1, nrow = nStore, byrow = TRUE)
                            }
                            if(p1 == 0){
                                beta1.p         <- NULL
                            }
                            if(p2 > 0){
                                beta2.p         <- matrix(mcmcRet$samples_beta2, nrow = nStore, byrow = TRUE)
                            }
                            if(p2 == 0){
                                beta2.p         <- NULL
                            }
                            if(p3 > 0){
                                beta3.p         <- matrix(mcmcRet$samples_beta3, nrow = nStore, byrow = TRUE)
                            }
                            if(p3 == 0){
                                beta3.p         <- NULL
                            }
                            
                            alpha1.p         <- matrix(mcmcRet$samples_alpha1, nrow = nStore, byrow = TRUE)
                            alpha2.p         <- matrix(mcmcRet$samples_alpha2, nrow = nStore, byrow = TRUE)
                            alpha3.p         <- matrix(mcmcRet$samples_alpha3, nrow = nStore, byrow = TRUE)
                            kappa1.p         <- matrix(mcmcRet$samples_kappa1, nrow = nStore, byrow = TRUE)
                            kappa2.p         <- matrix(mcmcRet$samples_kappa2, nrow = nStore, byrow = TRUE)
                            kappa3.p         <- matrix(mcmcRet$samples_kappa3, nrow = nStore, byrow = TRUE)
                            theta.p         <- matrix(mcmcRet$samples_theta, nrow = nStore, byrow = TRUE)
                            gamma.p         <- matrix(mcmcRet$samples_gamma, nrow = nStore, byrow = TRUE)
                            
                            V1.p            <- matrix(mcmcRet$samples_V1, nrow = nStore, byrow = TRUE)
                            V2.p            <- matrix(mcmcRet$samples_V2, nrow = nStore, byrow = TRUE)
                            V3.p            <- matrix(mcmcRet$samples_V3, nrow = nStore, byrow = TRUE)
                            Sigma_V.p        <- array(as.vector(mcmcRet$samples_Sigma_V), c(3, 3, nStore))
                            
                            if(p1 > 0){
                                accept.beta1     <- as.vector(mcmcRet$samples_misc[1:(p1)])/sum(as.vector(mcmcRet$moveVec)==1)*p1
                            }
                            if(p1 == 0){
                                accept.beta1     <- NULL
                            }
                            if(p2 > 0){
                                accept.beta2     <- as.vector(mcmcRet$samples_misc[(p1+1):(p1+p2)])/sum(as.vector(mcmcRet$moveVec)==2)*p2
                            }
                            if(p2 == 0){
                                accept.beta2     <- NULL
                            }
                            if(p3 > 0){
                                accept.beta3     <- as.vector(mcmcRet$samples_misc[(p1+p2+1):(p1+p2+p3)])/sum(as.vector(mcmcRet$moveVec)==3)*p3
                            }
                            if(p3 == 0){
                                accept.beta3     <- NULL
                            }
                            
                            accept.alpha1    <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+1)])/sum(as.vector(mcmcRet$moveVec)==4)
                            accept.alpha2    <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+2)])/sum(as.vector(mcmcRet$moveVec)==5)
                            accept.alpha3    <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+3)])/sum(as.vector(mcmcRet$moveVec)==6)
                            accept.theta    <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+4)])/sum(as.vector(mcmcRet$moveVec)==10)
                            accept.V1       <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+7+n+n+1):(p1+p2+p3+7+n+n+J)])/sum(as.vector(mcmcRet$moveVec)==13)
                            accept.V2       <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+7+n+n+J+1):(p1+p2+p3+7+n+n+J+J)])/sum(as.vector(mcmcRet$moveVec)==14)
                            accept.V3       <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+7+n+n+J+J+1):(p1+p2+p3+7+n+n+J+J+J)])/sum(as.vector(mcmcRet$moveVec)==15)
                            
                            V1summary <- as.matrix(apply(V1.p, 2, summary))
                            V1summary <- rbind(V1summary, apply(V1.p, 2, quantile, prob = 0.975))
                            V1summary <- rbind(V1summary, apply(V1.p, 2, quantile, prob = 0.025))
                            V1summary <- rbind(V1summary, apply(V1.p, 2, sd))
                            rownames(V1summary)[7:9] <- c("0.975", "0.025", "sd")
                            
                            V2summary <- as.matrix(apply(V2.p, 2, summary))
                            V2summary <- rbind(V2summary, apply(V2.p, 2, quantile, prob = 0.975))
                            V2summary <- rbind(V2summary, apply(V2.p, 2, quantile, prob = 0.025))
                            V2summary <- rbind(V2summary, apply(V2.p, 2, sd))
                            rownames(V2summary)[7:9] <- c("0.975", "0.025", "sd")
                            
                            V3summary <- as.matrix(apply(V3.p, 2, summary))
                            V3summary <- rbind(V3summary, apply(V3.p, 2, quantile, prob = 0.975))
                            V3summary <- rbind(V3summary, apply(V3.p, 2, quantile, prob = 0.025))
                            V3summary <- rbind(V3summary, apply(V3.p, 2, sd))
                            rownames(V3summary)[7:9] <- c("0.975", "0.025", "sd")
                            
                            if(storeV[1] == TRUE & !is.null(path))
                            {
                                save(V1.p, file = paste(path, "V1Pch", chain, ".RData", sep = ""))
                            }
                            if(storeV[2] == TRUE & !is.null(path))
                            {
                                save(V2.p, file = paste(path, "V2Pch", chain, ".RData", sep = ""))
                            }
                            if(storeV[3] == TRUE & !is.null(path))
                            {
                                save(V3.p, file = paste(path, "V3Pch", chain, ".RData", sep = ""))
                            }
                            
                            if(nGam_save > 0 & !is.null(path)){
                                save(gamma.p, file = paste(path, "/gammaPch", chain, ".Rdata", sep = ""))
                            }
                            
                            if(p1 > 0){
                                covNames1 = colnames(Xmat1)
                            }
                            if(p1 == 0){
                                covNames1 = NULL
                            }
                            
                            if(p2 > 0){
                                covNames2 = colnames(Xmat2)
                            }
                            if(p2 == 0){
                                covNames2 = NULL
                            }
                            if(p3 > 0){
                                covNames3 = colnames(Xmat3)
                            }
                            if(p3 == 0){
                                covNames3 = NULL
                            }
                            
                            ### posterior predictive checks ###
                            
                            ## 1. log pseudo-marginal likelihood
                            
                            invLH.p <- matrix(mcmcRet$invLH, nrow = n, byrow = TRUE)

                            ## 2. deviance information criterion
                            
                            gamma_mean <- matrix(mcmcRet$gammaP, nrow = n, byrow = TRUE)
                            dev.p        <- matrix(mcmcRet$dev, nrow = nStore, byrow = TRUE)
                            
                            ret[[nam]] <- list(beta1.p = beta1.p, beta2.p = beta2.p, beta3.p = beta3.p, alpha1.p = alpha1.p, alpha2.p = alpha2.p, alpha3.p = alpha3.p, kappa1.p = kappa1.p, kappa2.p = kappa2.p, kappa3.p = kappa3.p, theta.p = theta.p, Sigma_V.p = Sigma_V.p, accept.beta1 = accept.beta1, accept.beta2 = accept.beta2, accept.beta3 = accept.beta3, accept.alpha1 = accept.alpha1, accept.alpha2 = accept.alpha2, accept.alpha3 = accept.alpha3, accept.theta = accept.theta, accept.V1 = accept.V1, accept.V2 = accept.V2, accept.V3 = accept.V3, covNames1 = covNames1, covNames2 = covNames2, covNames3 = covNames3, V1sum = V1summary, V2sum = V2summary, V3sum = V3summary, gamma_mean = gamma_mean, dev.p=dev.p, invLH.p=invLH.p)
                            
                        } ## end: if Weibull-MVN-M
                        
                        
                        # model = "semi-Markov"
                        
                        #######################################
                        ############ 2 Weibull-MVN-SM #########
                        #######################################
                        
                        if(model_h3 == "semi-Markov"){
                            
                            ###
                            
                            startValues1 <- startV[1:(p1+p2+p3+n+7)]
                            startValues2 <- startV[(p1+p2+p3+n+8):length(startV)]
                            startV <- c(startValues1, 1, 1, startValues2)
                            
                            mcmc1 <- mcmc[1:4]
                            mcmc2 <- mcmc[5:7]
                            
                            mcmcNew <- c(mcmc1, 1, mcmc2, n)
                            
                            
                            mcmcRet <- .C("BweibMvnCorScrSMmcmc",
                            survData         = as.double(as.matrix(survData)),
                            n                = as.integer(n),
                            p1                = as.integer(p1),
                            p2                = as.integer(p2),
                            p3                = as.integer(p3),
                            J                = as.integer(J),
                            nj                = as.double(nj),
                            hyperParams     = as.double(hyperParams),
                            mcmc        = as.double(mcmcNew),
                            startValues     = as.double(startV),
                            numReps            = as.integer(numReps),
                            thin            = as.integer(thin),
                            burninPerc      = as.double(burninPerc),
                            nGam_save        = as.integer(nGam_save),
                            samples_beta1     = as.double(rep(0, nStore*p1)),
                            samples_beta2     = as.double(rep(0, nStore*p2)),
                            samples_beta3     = as.double(rep(0, nStore*p3)),
                            samples_alpha1     = as.double(rep(0, nStore*1)),
                            samples_alpha2     = as.double(rep(0, nStore*1)),
                            samples_alpha3     = as.double(rep(0, nStore*1)),
                            samples_kappa1     = as.double(rep(0, nStore*1)),
                            samples_kappa2     = as.double(rep(0, nStore*1)),
                            samples_kappa3     = as.double(rep(0, nStore*1)),
                            samples_nu2     = as.double(rep(0, nStore*1)),
                            samples_nu3     = as.double(rep(0, nStore*1)),
                            samples_theta     = as.double(rep(0, nStore*1)),
                            samples_V1        = as.double(rep(0, nStore*J)),
                            samples_V2        = as.double(rep(0, nStore*J)),
                            samples_V3        = as.double(rep(0, nStore*J)),
                            samples_Sigma_V    = as.double(rep(0, nStore*3*3)),
                            samples_gamma     = as.double(rep(0, nStore*nGam_save)),
                            samples_misc    = as.double(rep(0, (p1+p2+p3+6+n+1+n+J+J+J))),
                            gammaP            = as.double(rep(0, n)),
                            dev                 = as.double(rep(0, nStore*1)),
                            invLH = as.double(rep(0, n)),
                            logLH_fin = as.double(0),
                            lpml                 = as.double(rep(0, nStore*1)),
                            lpml2                 = as.double(rep(0, nStore*1)),
                            moveVec             = as.double(rep(0, numReps)))
                            
                            
                            
                            if(p1 > 0){
                                beta1.p         <- matrix(mcmcRet$samples_beta1, nrow = nStore, byrow = TRUE)
                            }
                            if(p1 == 0){
                                beta1.p         <- NULL
                            }
                            if(p2 > 0){
                                beta2.p         <- matrix(mcmcRet$samples_beta2, nrow = nStore, byrow = TRUE)
                            }
                            if(p2 == 0){
                                beta2.p         <- NULL
                            }
                            if(p3 > 0){
                                beta3.p         <- matrix(mcmcRet$samples_beta3, nrow = nStore, byrow = TRUE)
                            }
                            if(p3 == 0){
                                beta3.p         <- NULL
                            }
                            
                            alpha1.p         <- matrix(mcmcRet$samples_alpha1, nrow = nStore, byrow = TRUE)
                            alpha2.p         <- matrix(mcmcRet$samples_alpha2, nrow = nStore, byrow = TRUE)
                            alpha3.p         <- matrix(mcmcRet$samples_alpha3, nrow = nStore, byrow = TRUE)
                            kappa1.p         <- matrix(mcmcRet$samples_kappa1, nrow = nStore, byrow = TRUE)
                            kappa2.p         <- matrix(mcmcRet$samples_kappa2, nrow = nStore, byrow = TRUE)
                            kappa3.p         <- matrix(mcmcRet$samples_kappa3, nrow = nStore, byrow = TRUE)
                            theta.p         <- matrix(mcmcRet$samples_theta, nrow = nStore, byrow = TRUE)
                            gamma.p         <- matrix(mcmcRet$samples_gamma, nrow = nStore, byrow = TRUE)
                            V1.p            <- matrix(mcmcRet$samples_V1, nrow = nStore, byrow = TRUE)
                            V2.p            <- matrix(mcmcRet$samples_V2, nrow = nStore, byrow = TRUE)
                            V3.p            <- matrix(mcmcRet$samples_V3, nrow = nStore, byrow = TRUE)
                            Sigma_V.p        <- array(as.vector(mcmcRet$samples_Sigma_V), c(3, 3, nStore))
                            
                            if(p1 > 0){
                                accept.beta1     <- as.vector(mcmcRet$samples_misc[1:(p1)])/sum(as.vector(mcmcRet$moveVec)==1)*p1
                            }
                            if(p1 == 0){
                                accept.beta1     <- NULL
                            }
                            if(p2 > 0){
                                accept.beta2     <- as.vector(mcmcRet$samples_misc[(p1+1):(p1+p2)])/sum(as.vector(mcmcRet$moveVec)==2)*p2
                            }
                            if(p2 == 0){
                                accept.beta2     <- NULL
                            }
                            if(p3 > 0){
                                accept.beta3     <- as.vector(mcmcRet$samples_misc[(p1+p2+1):(p1+p2+p3)])/sum(as.vector(mcmcRet$moveVec)==3)*p3
                            }
                            if(p3 == 0){
                                accept.beta3     <- NULL
                            }
                            
                            accept.alpha1    <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+1)])/sum(as.vector(mcmcRet$moveVec)==4)
                            accept.alpha2    <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+2)])/sum(as.vector(mcmcRet$moveVec)==5)
                            accept.alpha3    <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+3)])/sum(as.vector(mcmcRet$moveVec)==6)
                            accept.theta    <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+4)])/sum(as.vector(mcmcRet$moveVec)==10)
                            accept.V1       <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+7+n+n+1):(p1+p2+p3+7+n+n+J)])/sum(as.vector(mcmcRet$moveVec)==13)
                            accept.V2       <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+7+n+n+J+1):(p1+p2+p3+7+n+n+J+J)])/sum(as.vector(mcmcRet$moveVec)==14)
                            accept.V3       <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+7+n+n+J+J+1):(p1+p2+p3+7+n+n+J+J+J)])/sum(as.vector(mcmcRet$moveVec)==15)
                            
                            V1summary <- as.matrix(apply(V1.p, 2, summary))
                            V1summary <- rbind(V1summary, apply(V1.p, 2, quantile, prob = 0.975))
                            V1summary <- rbind(V1summary, apply(V1.p, 2, quantile, prob = 0.025))
                            V1summary <- rbind(V1summary, apply(V1.p, 2, sd))
                            rownames(V1summary)[7:9] <- c("0.975", "0.025", "sd")
                            
                            V2summary <- as.matrix(apply(V2.p, 2, summary))
                            V2summary <- rbind(V2summary, apply(V2.p, 2, quantile, prob = 0.975))
                            V2summary <- rbind(V2summary, apply(V2.p, 2, quantile, prob = 0.025))
                            V2summary <- rbind(V2summary, apply(V2.p, 2, sd))
                            rownames(V2summary)[7:9] <- c("0.975", "0.025", "sd")
                            
                            V3summary <- as.matrix(apply(V3.p, 2, summary))
                            V3summary <- rbind(V3summary, apply(V3.p, 2, quantile, prob = 0.975))
                            V3summary <- rbind(V3summary, apply(V3.p, 2, quantile, prob = 0.025))
                            V3summary <- rbind(V3summary, apply(V3.p, 2, sd))
                            rownames(V3summary)[7:9] <- c("0.975", "0.025", "sd")
                            
                            if(storeV[1] == TRUE & !is.null(path))
                            {
                                save(V1.p, file = paste(path, "V1Pch", chain, ".RData", sep = ""))
                            }
                            if(storeV[2] == TRUE & !is.null(path))
                            {
                                save(V2.p, file = paste(path, "V2Pch", chain, ".RData", sep = ""))
                            }
                            if(storeV[3] == TRUE & !is.null(path))
                            {
                                save(V3.p, file = paste(path, "V3Pch", chain, ".RData", sep = ""))
                            }
                            
                            if(nGam_save > 0 & !is.null(path)){
                                save(gamma.p, file = paste(path, "/gammaPch", chain, ".Rdata", sep = ""))
                            }
                            
                            if(p1 > 0){
                                covNames1 = colnames(Xmat1)
                            }
                            if(p1 == 0){
                                covNames1 = NULL
                            }
                            
                            if(p2 > 0){
                                covNames2 = colnames(Xmat2)
                            }
                            if(p2 == 0){
                                covNames2 = NULL
                            }
                            if(p3 > 0){
                                covNames3 = colnames(Xmat3)
                            }
                            if(p3 == 0){
                                covNames3 = NULL
                            }
                            
                            
                            
                            ### posterior predictive checks ###
                            
                            ## 1. log pseudo-marginal likelihood
                            
                            invLH.p <- matrix(mcmcRet$invLH, nrow = n, byrow = TRUE)

                            ## 2. deviance information criterion
                            
                            gamma_mean <- matrix(mcmcRet$gammaP, nrow = n, byrow = TRUE)
                            dev.p        <- matrix(mcmcRet$dev, nrow = nStore, byrow = TRUE)

                            ret[[nam]] <- list(beta1.p = beta1.p, beta2.p = beta2.p, beta3.p = beta3.p, alpha1.p = alpha1.p, alpha2.p = alpha2.p, alpha3.p = alpha3.p, kappa1.p = kappa1.p, kappa2.p = kappa2.p, kappa3.p = kappa3.p, theta.p = theta.p, Sigma_V.p = Sigma_V.p, accept.beta1 = accept.beta1, accept.beta2 = accept.beta2, accept.beta3 = accept.beta3, accept.alpha1 = accept.alpha1, accept.alpha2 = accept.alpha2, accept.alpha3 = accept.alpha3, accept.theta = accept.theta, accept.V1 = accept.V1, accept.V2 = accept.V2, accept.V3 = accept.V3, covNames1 = covNames1, covNames2 = covNames2, covNames3 = covNames3, V1sum = V1summary, V2sum = V2summary, V3sum = V3summary, gamma_mean = gamma_mean, dev.p=dev.p, invLH.p=invLH.p)
                            
                            
                        }  ## end: if Weibull-MVN-SM
                        
                        
                        ## DIC and LPML
                        
                        be1.p <- ret$chain1$beta1.p
                        be2.p <- ret$chain1$beta2.p
                        be3.p <- ret$chain1$beta3.p
                        alp1.p <- ret$chain1$alpha1.p
                        alp2.p <- ret$chain1$alpha2.p
                        alp3.p <- ret$chain1$alpha3.p
                        kap1.p <- ret$chain1$kappa1.p
                        kap2.p <- ret$chain1$kappa2.p
                        kap3.p <- ret$chain1$kappa3.p
                        thet.p <- ret$chain1$theta.p
                        V1me <- ret$chain1$V1sum[4,]
                        V2me <- ret$chain1$V2sum[4,]
                        V3me <- ret$chain1$V3sum[4,]
                        Dev.p <- ret$chain1$dev.p
                        InvLH.p <- ret$chain1$invLH.p
                        
                        if(nChain > 1){
                            for(i in 2:nChain){
                                nam <- paste("chain", i, sep="")
                                be1.p <- rbind(be1.p, ret[[nam]]$beta1.p)
                                be2.p <- rbind(be2.p, ret[[nam]]$beta2.p)
                                be3.p <- rbind(be3.p, ret[[nam]]$beta3.p)
                                alp1.p <- rbind(alp1.p, ret[[nam]]$alpha1.p)
                                alp2.p <- rbind(alp2.p, ret[[nam]]$alpha2.p)
                                alp3.p <- rbind(alp3.p, ret[[nam]]$alpha3.p)
                                kap1.p <- rbind(kap1.p, ret[[nam]]$kappa1.p)
                                kap2.p <- rbind(kap2.p, ret[[nam]]$kappa2.p)
                                kap3.p <- rbind(kap3.p, ret[[nam]]$kappa3.p)
                                thet.p <- rbind(thet.p, ret[[nam]]$theta.p)
                                V1me <- rbind(V1me, ret[[nam]]$V1sum[4,])
                                V2me <- rbind(V2me, ret[[nam]]$V2sum[4,])
                                V3me <- rbind(V3me, ret[[nam]]$V3sum[4,])
                                Dev.p <- rbind(Dev.p, ret[[nam]]$dev.p)
                                InvLH.p <- cbind(InvLH.p, ret[[nam]]$invLH.p)
                            }
                        }
                        if(model_h3 == "Markov")
                        {
                            logLH_fin <- .C("BweibMvnCorScr_logMLH_DIC",
                            survData         = as.double(as.matrix(survData)),
                            n                = as.integer(n),
                            p1                = as.integer(p1),
                            p2                = as.integer(p2),
                            p3                = as.integer(p3),
                            J                = as.integer(J),
                            be1 = as.double(apply(be1.p, 2, mean)),
                            be2 = as.double(apply(be2.p, 2, mean)),
                            be3 = as.double(apply(be3.p, 2, mean)),
                            alp1 = as.double(apply(alp1.p, 2, mean)),
                            alp2 = as.double(apply(alp2.p, 2, mean)),
                            alp3 = as.double(apply(alp3.p, 2, mean)),
                            kap1 = as.double(apply(kap1.p, 2, mean)),
                            kap2 = as.double(apply(kap2.p, 2, mean)),
                            kap3 = as.double(apply(kap3.p, 2, mean)),
                            thet = as.double(apply(theta.p, 2, mean)),
                            V_1 = as.double(apply(V1me, 2, mean)),
                            V_2 = as.double(apply(V2me, 2, mean)),
                            V_3 = as.double(apply(V3me, 2, mean)),
                            val = as.double(0))$val
                        }else if(model_h3 == "semi-Markov")
                        {
                            logLH_fin <- .C("BweibMvnCorScrSM_logMLH_DIC",
                            survData         = as.double(as.matrix(survData)),
                            n                = as.integer(n),
                            p1                = as.integer(p1),
                            p2                = as.integer(p2),
                            p3                = as.integer(p3),
                            J                = as.integer(J),
                            be1 = as.double(apply(be1.p, 2, mean)),
                            be2 = as.double(apply(be2.p, 2, mean)),
                            be3 = as.double(apply(be3.p, 2, mean)),
                            alp1 = as.double(apply(alp1.p, 2, mean)),
                            alp2 = as.double(apply(alp2.p, 2, mean)),
                            alp3 = as.double(apply(alp3.p, 2, mean)),
                            kap1 = as.double(apply(kap1.p, 2, mean)),
                            kap2 = as.double(apply(kap2.p, 2, mean)),
                            kap3 = as.double(apply(kap3.p, 2, mean)),
                            thet = as.double(apply(theta.p, 2, mean)),
                            V_1 = as.double(apply(V1me, 2, mean)),
                            V_2 = as.double(apply(V2me, 2, mean)),
                            V_3 = as.double(apply(V3me, 2, mean)),
                            val = as.double(0))$val
                        }
                        
                        DIC = 2*mean(Dev.p) + 2 * logLH_fin
                        LPML = sum(log(1/apply(InvLH.p, 1, mean)))
                        
                    } ## end: if Weibull-MVN
                    
                    
                    # re.type = "DPM"
                    
                    if(re.type == "DPM"){
                        
                        # model = "Markov"
                        
                        #######################################
                        ############ Weibull-DPM-M ############
                        #######################################
                        
                        if(model_h3 == "Markov"){
                            
                            ##
                            startValues1 <- startV[1:(p1+p2+p3+n+7)]
                            startValues2 <- startV[(p1+p2+p3+n+8):length(startV)]
                            startV <- c(startValues1, 1, 1, startValues2)
                            
                            mcmc1 <- mcmc[1:4]
                            mcmc2 <- mcmc[5:7]
                            
                            mcmcNew <- c(mcmc1, 1, mcmc2, n)
                            
                            
                            mcmcRet <- .C("BweibDpCorScrmcmc",
                            survData         = as.double(as.matrix(survData)),
                            n                = as.integer(n),
                            p1                = as.integer(p1),
                            p2                = as.integer(p2),
                            p3                = as.integer(p3),
                            J                = as.integer(J),
                            nj                = as.double(nj),
                            hyperParams     = as.double(hyperParams),
                            mcmc        = as.double(mcmcNew),
                            startValues     = as.double(startV),
                            numReps            = as.integer(numReps),
                            thin            = as.integer(thin),
                            burninPerc      = as.double(burninPerc),
                            nGam_save        = as.integer(nGam_save),
                            samples_beta1     = as.double(rep(0, nStore*p1)),
                            samples_beta2     = as.double(rep(0, nStore*p2)),
                            samples_beta3     = as.double(rep(0, nStore*p3)),
                            samples_alpha1     = as.double(rep(0, nStore*1)),
                            samples_alpha2     = as.double(rep(0, nStore*1)),
                            samples_alpha3     = as.double(rep(0, nStore*1)),
                            samples_kappa1     = as.double(rep(0, nStore*1)),
                            samples_kappa2     = as.double(rep(0, nStore*1)),
                            samples_kappa3     = as.double(rep(0, nStore*1)),
                            samples_nu2     = as.double(rep(0, nStore*1)),
                            samples_nu3     = as.double(rep(0, nStore*1)),
                            samples_theta     = as.double(rep(0, nStore*1)),
                            samples_V1        = as.double(rep(0, nStore*J)),
                            samples_V2        = as.double(rep(0, nStore*J)),
                            samples_V3        = as.double(rep(0, nStore*J)),
                            samples_c        = as.double(rep(0, nStore*J)),
                            samples_mu        = as.double(rep(0, nStore*3*J)),
                            samples_Sigma    = as.double(rep(0, nStore*3*3*J)),
                            samples_tau     = as.double(rep(0, nStore*1)),
                            samples_gamma     = as.double(rep(0, nStore*nGam_save)),
                            samples_misc    = as.double(rep(0, (p1+p2+p3+6+n+1+n+J))),
                            gammaP            = as.double(rep(0, n)),
                            dev                 = as.double(rep(0, nStore*1)),
                            invLH = as.double(rep(0, n)),
                            logLH_fin = as.double(0),
                            lpml                 = as.double(rep(0, nStore*1)),
                            lpml2                 = as.double(rep(0, nStore*1)),
                            moveVec             = as.double(rep(0, numReps)))
                            
                            
                            
                            if(p1 > 0){
                                beta1.p         <- matrix(mcmcRet$samples_beta1, nrow = nStore, byrow = TRUE)
                            }
                            if(p1 == 0){
                                beta1.p         <- NULL
                            }
                            if(p2 > 0){
                                beta2.p         <- matrix(mcmcRet$samples_beta2, nrow = nStore, byrow = TRUE)
                            }
                            if(p2 == 0){
                                beta2.p         <- NULL
                            }
                            if(p3 > 0){
                                beta3.p         <- matrix(mcmcRet$samples_beta3, nrow = nStore, byrow = TRUE)
                            }
                            if(p3 == 0){
                                beta3.p         <- NULL
                            }
                            
                            alpha1.p         <- matrix(mcmcRet$samples_alpha1, nrow = nStore, byrow = TRUE)
                            alpha2.p         <- matrix(mcmcRet$samples_alpha2, nrow = nStore, byrow = TRUE)
                            alpha3.p         <- matrix(mcmcRet$samples_alpha3, nrow = nStore, byrow = TRUE)
                            kappa1.p         <- matrix(mcmcRet$samples_kappa1, nrow = nStore, byrow = TRUE)
                            kappa2.p         <- matrix(mcmcRet$samples_kappa2, nrow = nStore, byrow = TRUE)
                            kappa3.p         <- matrix(mcmcRet$samples_kappa3, nrow = nStore, byrow = TRUE)
                            nu2.p           <- matrix(mcmcRet$samples_nu2, nrow = nStore, byrow = TRUE)
                            nu3.p           <- matrix(mcmcRet$samples_nu3, nrow = nStore, byrow = TRUE)
                            theta.p         <- matrix(mcmcRet$samples_theta, nrow = nStore, byrow = TRUE)
                            gamma.p         <- matrix(mcmcRet$samples_gamma, nrow = nStore, byrow = TRUE)
                            
                            V1.p            <- matrix(mcmcRet$samples_V1, nrow = nStore, byrow = TRUE)
                            V2.p            <- matrix(mcmcRet$samples_V2, nrow = nStore, byrow = TRUE)
                            V3.p            <- matrix(mcmcRet$samples_V3, nrow = nStore, byrow = TRUE)
                            c.p            <- matrix(mcmcRet$samples_c, nrow = nStore, byrow = TRUE)
                            mu.p            <- matrix(mcmcRet$samples_mu, nrow = nStore, byrow = TRUE)
                            Sigma.p         <- array(as.vector(mcmcRet$samples_Sigma), c(3, 3 * J, nStore))
                            tau.p           <- matrix(mcmcRet$samples_tau, nrow = nStore, byrow = TRUE)
                            
                            if(p1 > 0){
                                accept.beta1     <- as.vector(mcmcRet$samples_misc[1:(p1)])/sum(as.vector(mcmcRet$moveVec)==1)*p1
                            }
                            if(p1 == 0){
                                accept.beta1     <- NULL
                            }
                            if(p2 > 0){
                                accept.beta2     <- as.vector(mcmcRet$samples_misc[(p1+1):(p1+p2)])/sum(as.vector(mcmcRet$moveVec)==2)*p2
                            }
                            if(p2 == 0){
                                accept.beta2     <- NULL
                            }
                            if(p3 > 0){
                                accept.beta3     <- as.vector(mcmcRet$samples_misc[(p1+p2+1):(p1+p2+p3)])/sum(as.vector(mcmcRet$moveVec)==3)*p3
                            }
                            if(p3 == 0){
                                accept.beta3     <- NULL
                            }
                            
                            accept.alpha1    <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+1)])/sum(as.vector(mcmcRet$moveVec)==4)
                            accept.alpha2    <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+2)])/sum(as.vector(mcmcRet$moveVec)==5)
                            accept.alpha3    <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+3)])/sum(as.vector(mcmcRet$moveVec)==6)
                            accept.theta    <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+4)])/sum(as.vector(mcmcRet$moveVec)==10)
                            accept.nu2      <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+5)])
                            accept.nu3      <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+6)])
                            accept.gamma    <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+6+1):(p1+p2+p3+6+n)])
                            lastChgProp     <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+6+n+1)])
                            mhGam_chk       <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+6+n+1+1):(p1+p2+p3+6+n+1+n)])
                            accept.V        <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+7+n+n+1):(p1+p2+p3+7+n+n+J)])/sum(as.vector(mcmcRet$moveVec)==13)/3
                            
                            
                            V1summary <- as.matrix(apply(V1.p, 2, summary))
                            V1summary <- rbind(V1summary, apply(V1.p, 2, quantile, prob = 0.975))
                            V1summary <- rbind(V1summary, apply(V1.p, 2, quantile, prob = 0.025))
                            V1summary <- rbind(V1summary, apply(V1.p, 2, sd))
                            rownames(V1summary)[7:9] <- c("0.975", "0.025", "sd")
                            
                            V2summary <- as.matrix(apply(V2.p, 2, summary))
                            V2summary <- rbind(V2summary, apply(V2.p, 2, quantile, prob = 0.975))
                            V2summary <- rbind(V2summary, apply(V2.p, 2, quantile, prob = 0.025))
                            V2summary <- rbind(V2summary, apply(V2.p, 2, sd))
                            rownames(V2summary)[7:9] <- c("0.975", "0.025", "sd")
                            
                            V3summary <- as.matrix(apply(V3.p, 2, summary))
                            V3summary <- rbind(V3summary, apply(V3.p, 2, quantile, prob = 0.975))
                            V3summary <- rbind(V3summary, apply(V3.p, 2, quantile, prob = 0.025))
                            V3summary <- rbind(V3summary, apply(V3.p, 2, sd))
                            rownames(V3summary)[7:9] <- c("0.975", "0.025", "sd")
                            
                            
                            if(storeV[1] == TRUE & !is.null(path))
                            {
                                save(V1.p, file = paste(path, "V1Pch", chain, ".RData", sep = ""))
                            }
                            if(storeV[2] == TRUE & !is.null(path))
                            {
                                save(V2.p, file = paste(path, "V2Pch", chain, ".RData", sep = ""))
                            }
                            if(storeV[3] == TRUE & !is.null(path))
                            {
                                save(V3.p, file = paste(path, "V3Pch", chain, ".RData", sep = ""))
                            }
                            
                            if(nGam_save > 0 & !is.null(path)){
                                save(gamma.p, file = paste(path, "/gammaPch", chain, ".Rdata", sep = ""))
                            }
                            
                            
                            if(p1 > 0){
                                covNames1 = colnames(Xmat1)
                            }
                            if(p1 == 0){
                                covNames1 = NULL
                            }
                            
                            if(p2 > 0){
                                covNames2 = colnames(Xmat2)
                            }
                            if(p2 == 0){
                                covNames2 = NULL
                            }
                            if(p3 > 0){
                                covNames3 = colnames(Xmat3)
                            }
                            if(p3 == 0){
                                covNames3 = NULL
                            }
                            
                            
                            ### posterior predictive checks ###
                            
                            ## 1. log pseudo-marginal likelihood
                            
                            invLH.p <- matrix(mcmcRet$invLH, nrow = n, byrow = TRUE)

                            ## 2. deviance information criterion
                            
                            gamma_mean <- matrix(mcmcRet$gammaP, nrow = n, byrow = TRUE)
                            dev.p        <- matrix(mcmcRet$dev, nrow = nStore, byrow = TRUE)
                            
                            ret[[nam]] <- list(beta1.p = beta1.p, beta2.p = beta2.p, beta3.p = beta3.p, alpha1.p = alpha1.p, alpha2.p = alpha2.p, alpha3.p = alpha3.p, kappa1.p = kappa1.p, kappa2.p = kappa2.p, kappa3.p = kappa3.p, theta.p = theta.p, V1.p = V1.p, V2.p = V2.p, V3.p = V3.p, class.p = c.p, mu.p = mu.p, Sigma.p = Sigma.p, tau.p = tau.p, accept.beta1 = accept.beta1, accept.beta2 = accept.beta2, accept.beta3 = accept.beta3, accept.alpha1 = accept.alpha1, accept.alpha2 = accept.alpha2, accept.alpha3 = accept.alpha3, accept.theta = accept.theta, accept.V = accept.V, covNames1 = covNames1, covNames2 = covNames2, covNames3 = covNames3, V1sum = V1summary, V2sum = V2summary, V3sum = V3summary, gamma_mean = gamma_mean, dev.p=dev.p, invLH.p=invLH.p)
                            
                            
                        }  ## end: if Weibull-DPM-M
                        
                        
                        # model = "semi-Markov"
                        
                        #######################################
                        ############ Weibull-DPM-SM ############
                        #######################################
                        
                        if(model_h3 == "semi-Markov"){
                            
                            ###
                            startValues1 <- startV[1:(p1+p2+p3+n+7)]
                            startValues2 <- startV[(p1+p2+p3+n+8):length(startV)]
                            startV <- c(startValues1, 1, 1, startValues2)
                            
                            mcmc1 <- mcmc[1:4]
                            mcmc2 <- mcmc[5:7]
                            
                            mcmcNew <- c(mcmc1, 1, mcmc2, n)
                            
                            mcmcRet <- .C("BweibDpCorScrSMmcmc",
                            survData         = as.double(as.matrix(survData)),
                            n                = as.integer(n),
                            p1                = as.integer(p1),
                            p2                = as.integer(p2),
                            p3                = as.integer(p3),
                            J                = as.integer(J),
                            nj                = as.double(nj),
                            hyperParams     = as.double(hyperParams),
                            mcmc        = as.double(mcmcNew),
                            startValues     = as.double(startV),
                            numReps            = as.integer(numReps),
                            thin            = as.integer(thin),
                            burninPerc      = as.double(burninPerc),
                            nGam_save        = as.integer(nGam_save),
                            samples_beta1     = as.double(rep(0, nStore*p1)),
                            samples_beta2     = as.double(rep(0, nStore*p2)),
                            samples_beta3     = as.double(rep(0, nStore*p3)),
                            samples_alpha1     = as.double(rep(0, nStore*1)),
                            samples_alpha2     = as.double(rep(0, nStore*1)),
                            samples_alpha3     = as.double(rep(0, nStore*1)),
                            samples_kappa1     = as.double(rep(0, nStore*1)),
                            samples_kappa2     = as.double(rep(0, nStore*1)),
                            samples_kappa3     = as.double(rep(0, nStore*1)),
                            samples_nu2     = as.double(rep(0, nStore*1)),
                            samples_nu3     = as.double(rep(0, nStore*1)),
                            samples_theta     = as.double(rep(0, nStore*1)),
                            samples_V1        = as.double(rep(0, nStore*J)),
                            samples_V2        = as.double(rep(0, nStore*J)),
                            samples_V3        = as.double(rep(0, nStore*J)),
                            samples_c        = as.double(rep(0, nStore*J)),
                            samples_mu        = as.double(rep(0, nStore*3*J)),
                            samples_Sigma    = as.double(rep(0, nStore*3*3*J)),
                            samples_tau     = as.double(rep(0, nStore*1)),
                            samples_gamma     = as.double(rep(0, nStore*nGam_save)),
                            samples_misc    = as.double(rep(0, (p1+p2+p3+6+n+1+n+J))),
                            gammaP            = as.double(rep(0, n)),
                            dev                 = as.double(rep(0, nStore*1)),
                            invLH = as.double(rep(0, n)),
                            logLH_fin = as.double(0),
                            lpml                 = as.double(rep(0, nStore*1)),
                            lpml2                 = as.double(rep(0, nStore*1)),
                            moveVec             = as.double(rep(0, numReps)))
                            
                            
                            if(p1 > 0){
                                beta1.p         <- matrix(mcmcRet$samples_beta1, nrow = nStore, byrow = TRUE)
                            }
                            if(p1 == 0){
                                beta1.p         <- NULL
                            }
                            if(p2 > 0){
                                beta2.p         <- matrix(mcmcRet$samples_beta2, nrow = nStore, byrow = TRUE)
                            }
                            if(p2 == 0){
                                beta2.p         <- NULL
                            }
                            if(p3 > 0){
                                beta3.p         <- matrix(mcmcRet$samples_beta3, nrow = nStore, byrow = TRUE)
                            }
                            if(p3 == 0){
                                beta3.p         <- NULL
                            }
                            
                            alpha1.p         <- matrix(mcmcRet$samples_alpha1, nrow = nStore, byrow = TRUE)
                            alpha2.p         <- matrix(mcmcRet$samples_alpha2, nrow = nStore, byrow = TRUE)
                            alpha3.p         <- matrix(mcmcRet$samples_alpha3, nrow = nStore, byrow = TRUE)
                            kappa1.p         <- matrix(mcmcRet$samples_kappa1, nrow = nStore, byrow = TRUE)
                            kappa2.p         <- matrix(mcmcRet$samples_kappa2, nrow = nStore, byrow = TRUE)
                            kappa3.p         <- matrix(mcmcRet$samples_kappa3, nrow = nStore, byrow = TRUE)
                            nu2.p           <- matrix(mcmcRet$samples_nu2, nrow = nStore, byrow = TRUE)
                            nu3.p           <- matrix(mcmcRet$samples_nu3, nrow = nStore, byrow = TRUE)
                            theta.p         <- matrix(mcmcRet$samples_theta, nrow = nStore, byrow = TRUE)
                            gamma.p         <- matrix(mcmcRet$samples_gamma, nrow = nStore, byrow = TRUE)
                            
                            V1.p            <- matrix(mcmcRet$samples_V1, nrow = nStore, byrow = TRUE)
                            V2.p            <- matrix(mcmcRet$samples_V2, nrow = nStore, byrow = TRUE)
                            V3.p            <- matrix(mcmcRet$samples_V3, nrow = nStore, byrow = TRUE)
                            c.p            <- matrix(mcmcRet$samples_c, nrow = nStore, byrow = TRUE)
                            mu.p            <- matrix(mcmcRet$samples_mu, nrow = nStore, byrow = TRUE)
                            Sigma.p         <- array(as.vector(mcmcRet$samples_Sigma), c(3, 3 * J, nStore))
                            tau.p           <- matrix(mcmcRet$samples_tau, nrow = nStore, byrow = TRUE)
                            
                            if(p1 > 0){
                                accept.beta1     <- as.vector(mcmcRet$samples_misc[1:(p1)])/sum(as.vector(mcmcRet$moveVec)==1)*p1
                            }
                            if(p1 == 0){
                                accept.beta1     <- NULL
                            }
                            if(p2 > 0){
                                accept.beta2     <- as.vector(mcmcRet$samples_misc[(p1+1):(p1+p2)])/sum(as.vector(mcmcRet$moveVec)==2)*p2
                            }
                            if(p2 == 0){
                                accept.beta2     <- NULL
                            }
                            if(p3 > 0){
                                accept.beta3     <- as.vector(mcmcRet$samples_misc[(p1+p2+1):(p1+p2+p3)])/sum(as.vector(mcmcRet$moveVec)==3)*p3
                            }
                            if(p3 == 0){
                                accept.beta3     <- NULL
                            }
                            
                            accept.alpha1    <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+1)])/sum(as.vector(mcmcRet$moveVec)==4)
                            accept.alpha2    <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+2)])/sum(as.vector(mcmcRet$moveVec)==5)
                            accept.alpha3    <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+3)])/sum(as.vector(mcmcRet$moveVec)==6)
                            accept.theta    <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+4)])/sum(as.vector(mcmcRet$moveVec)==10)
                            accept.nu2      <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+5)])
                            accept.nu3      <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+6)])
                            accept.gamma    <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+6+1):(p1+p2+p3+6+n)])
                            lastChgProp     <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+6+n+1)])
                            mhGam_chk       <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+6+n+1+1):(p1+p2+p3+6+n+1+n)])
                            accept.V        <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+7+n+n+1):(p1+p2+p3+7+n+n+J)])/sum(as.vector(mcmcRet$moveVec)==13)/3
                            
                            V1summary <- as.matrix(apply(V1.p, 2, summary))
                            V1summary <- rbind(V1summary, apply(V1.p, 2, quantile, prob = 0.975))
                            V1summary <- rbind(V1summary, apply(V1.p, 2, quantile, prob = 0.025))
                            V1summary <- rbind(V1summary, apply(V1.p, 2, sd))
                            rownames(V1summary)[7:9] <- c("0.975", "0.025", "sd")
                            
                            V2summary <- as.matrix(apply(V2.p, 2, summary))
                            V2summary <- rbind(V2summary, apply(V2.p, 2, quantile, prob = 0.975))
                            V2summary <- rbind(V2summary, apply(V2.p, 2, quantile, prob = 0.025))
                            V2summary <- rbind(V2summary, apply(V2.p, 2, sd))
                            rownames(V2summary)[7:9] <- c("0.975", "0.025", "sd")
                            
                            V3summary <- as.matrix(apply(V3.p, 2, summary))
                            V3summary <- rbind(V3summary, apply(V3.p, 2, quantile, prob = 0.975))
                            V3summary <- rbind(V3summary, apply(V3.p, 2, quantile, prob = 0.025))
                            V3summary <- rbind(V3summary, apply(V3.p, 2, sd))
                            rownames(V3summary)[7:9] <- c("0.975", "0.025", "sd")
                            
                            if(storeV[1] == TRUE & !is.null(path))
                            {
                                save(V1.p, file = paste(path, "V1Pch", chain, ".RData", sep = ""))
                            }
                            if(storeV[2] == TRUE & !is.null(path))
                            {
                                save(V2.p, file = paste(path, "V2Pch", chain, ".RData", sep = ""))
                            }
                            if(storeV[3] == TRUE & !is.null(path))
                            {
                                save(V3.p, file = paste(path, "V3Pch", chain, ".RData", sep = ""))
                            }
                            
                            if(nGam_save > 0 & !is.null(path)){
                                save(gamma.p, file = paste(path, "/gammaPch", chain, ".Rdata", sep = ""))
                            }
                            
                            if(p1 > 0){
                                covNames1 = colnames(Xmat1)
                            }
                            if(p1 == 0){
                                covNames1 = NULL
                            }
                            
                            if(p2 > 0){
                                covNames2 = colnames(Xmat2)
                            }
                            if(p2 == 0){
                                covNames2 = NULL
                            }
                            if(p3 > 0){
                                covNames3 = colnames(Xmat3)
                            }
                            if(p3 == 0){
                                covNames3 = NULL
                            }
                            
                            
                            ### posterior predictive checks ###
                            
                            ## 1. log pseudo-marginal likelihood
                            
                            invLH.p <- matrix(mcmcRet$invLH, nrow = n, byrow = TRUE)
                            
                            ## 2. deviance information criterion
                            
                            gamma_mean <- matrix(mcmcRet$gammaP, nrow = n, byrow = TRUE)
                            dev.p        <- matrix(mcmcRet$dev, nrow = nStore, byrow = TRUE)
                            
                            ret[[nam]] <- list(beta1.p = beta1.p, beta2.p = beta2.p, beta3.p = beta3.p, alpha1.p = alpha1.p, alpha2.p = alpha2.p, alpha3.p = alpha3.p, kappa1.p = kappa1.p, kappa2.p = kappa2.p, kappa3.p = kappa3.p, theta.p = theta.p, V1.p = V1.p, V2.p = V2.p, V3.p = V3.p, class.p = c.p, mu.p = mu.p, Sigma.p = Sigma.p, tau.p = tau.p, accept.beta1 = accept.beta1, accept.beta2 = accept.beta2, accept.beta3 = accept.beta3, accept.alpha1 = accept.alpha1, accept.alpha2 = accept.alpha2, accept.alpha3 = accept.alpha3, accept.theta = accept.theta, accept.V = accept.V, covNames1 = covNames1, covNames2 = covNames2, covNames3 = covNames3, V1sum = V1summary, V2sum = V2summary, V3sum = V3summary, gamma_mean = gamma_mean, dev.p=dev.p, invLH.p=invLH.p)
                            
                            
                        }  ## end: if Weibull-DPM-SM
                        
                        ## DIC and LPML
                        
                        be1.p <- ret$chain1$beta1.p
                        be2.p <- ret$chain1$beta2.p
                        be3.p <- ret$chain1$beta3.p
                        alp1.p <- ret$chain1$alpha1.p
                        alp2.p <- ret$chain1$alpha2.p
                        alp3.p <- ret$chain1$alpha3.p
                        kap1.p <- ret$chain1$kappa1.p
                        kap2.p <- ret$chain1$kappa2.p
                        kap3.p <- ret$chain1$kappa3.p
                        thet.p <- ret$chain1$theta.p
                        V1me <- ret$chain1$V1sum[4,]
                        V2me <- ret$chain1$V2sum[4,]
                        V3me <- ret$chain1$V3sum[4,]
                        Dev.p <- ret$chain1$dev.p
                        InvLH.p <- ret$chain1$invLH.p
                        
                        if(nChain > 1){
                            for(i in 2:nChain){
                                nam <- paste("chain", i, sep="")
                                be1.p <- rbind(be1.p, ret[[nam]]$beta1.p)
                                be2.p <- rbind(be2.p, ret[[nam]]$beta2.p)
                                be3.p <- rbind(be3.p, ret[[nam]]$beta3.p)
                                alp1.p <- rbind(alp1.p, ret[[nam]]$alpha1.p)
                                alp2.p <- rbind(alp2.p, ret[[nam]]$alpha2.p)
                                alp3.p <- rbind(alp3.p, ret[[nam]]$alpha3.p)
                                kap1.p <- rbind(kap1.p, ret[[nam]]$kappa1.p)
                                kap2.p <- rbind(kap2.p, ret[[nam]]$kappa2.p)
                                kap3.p <- rbind(kap3.p, ret[[nam]]$kappa3.p)
                                thet.p <- rbind(thet.p, ret[[nam]]$theta.p)
                                V1me <- rbind(V1me, ret[[nam]]$V1sum[4,])
                                V2me <- rbind(V2me, ret[[nam]]$V2sum[4,])
                                V3me <- rbind(V3me, ret[[nam]]$V3sum[4,])
                                Dev.p <- rbind(Dev.p, ret[[nam]]$dev.p)
                                InvLH.p <- cbind(InvLH.p, ret[[nam]]$invLH.p)
                            }
                        }
                        if(model_h3 == "Markov")
                        {
                            logLH_fin <- .C("BweibDpCorScr_logMLH_DIC",
                            survData         = as.double(as.matrix(survData)),
                            n                = as.integer(n),
                            p1                = as.integer(p1),
                            p2                = as.integer(p2),
                            p3                = as.integer(p3),
                            J                = as.integer(J),
                            be1 = as.double(apply(be1.p, 2, mean)),
                            be2 = as.double(apply(be2.p, 2, mean)),
                            be3 = as.double(apply(be3.p, 2, mean)),
                            alp1 = as.double(apply(alp1.p, 2, mean)),
                            alp2 = as.double(apply(alp2.p, 2, mean)),
                            alp3 = as.double(apply(alp3.p, 2, mean)),
                            kap1 = as.double(apply(kap1.p, 2, mean)),
                            kap2 = as.double(apply(kap2.p, 2, mean)),
                            kap3 = as.double(apply(kap3.p, 2, mean)),
                            thet = as.double(apply(theta.p, 2, mean)),
                            V_1 = as.double(apply(V1me, 2, mean)),
                            V_2 = as.double(apply(V2me, 2, mean)),
                            V_3 = as.double(apply(V3me, 2, mean)),
                            val = as.double(0))$val
                        }else if(model_h3 == "semi-Markov")
                        {
                            logLH_fin <- .C("BweibDpCorScrSM_logMLH_DIC",
                            survData         = as.double(as.matrix(survData)),
                            n                = as.integer(n),
                            p1                = as.integer(p1),
                            p2                = as.integer(p2),
                            p3                = as.integer(p3),
                            J                = as.integer(J),
                            be1 = as.double(apply(be1.p, 2, mean)),
                            be2 = as.double(apply(be2.p, 2, mean)),
                            be3 = as.double(apply(be3.p, 2, mean)),
                            alp1 = as.double(apply(alp1.p, 2, mean)),
                            alp2 = as.double(apply(alp2.p, 2, mean)),
                            alp3 = as.double(apply(alp3.p, 2, mean)),
                            kap1 = as.double(apply(kap1.p, 2, mean)),
                            kap2 = as.double(apply(kap2.p, 2, mean)),
                            kap3 = as.double(apply(kap3.p, 2, mean)),
                            thet = as.double(apply(theta.p, 2, mean)),
                            V_1 = as.double(apply(V1me, 2, mean)),
                            V_2 = as.double(apply(V2me, 2, mean)),
                            V_3 = as.double(apply(V3me, 2, mean)),
                            val = as.double(0))$val
                        }
                        
                        DIC = 2*mean(Dev.p) + 2 * logLH_fin
                        LPML = sum(log(1/apply(InvLH.p, 1, mean)))
                        
                    } ## end: if Weibull-DPM
                    
                } ## end: if Weibull
                
                
                # hz.type = "PEM"
                
                if(hz.type == "PEM"){

                    C1 <- mcmc[1]
                    C2 <- mcmc[2]
                    C3 <- mcmc[3]
                    delPert1 <- mcmc[4]
                    delPert2 <- mcmc[5]
                    delPert3 <- mcmc[6]
                    rj.scheme <- mcmc[7]
                    K1_max <- mcmc[8]
                    K2_max <- mcmc[9]
                    K3_max <- mcmc[10]
                    s1_max <- max(temp$PEM$PEM.s1)
                    s2_max <- max(temp$PEM$PEM.s2)
                    s3_max <- max(temp$PEM$PEM.s3)
                    mhProp_theta_var <- mcmc[14]
                    mhProp_V1_var   <- mcmc[15]
                    mhProp_V2_var   <- mcmc[16]
                    mhProp_V3_var   <- mcmc[17]
                    nGam_save   <- mcmc[18]
                    numReps     <- mcmc[19]
                    thin        <- mcmc[20]
                    burninPerc  <- mcmc[21]
                    storeV      <- mcmc[22:24]
                    
                    nStore <- round(numReps/thin*(1-burninPerc))
                    
                    
                    ## recommended when the unit of time is day
                    if(rj.scheme == 1){
                        s_propBI1 <- seq(1, s1_max, 1)
                        s_propBI1 <- s_propBI1[s_propBI1 < s1_max]
                        s_propBI2 <- seq(1, s2_max, 1)
                        s_propBI2 <- s_propBI2[s_propBI2 < s2_max]
                        s_propBI3 <- seq(1, s3_max, 1)
                        s_propBI3 <- s_propBI3[s_propBI3 < s3_max]
                    }
                    ## uniquely ordered failure time
                    if(rj.scheme == 2){
                        s_propBI1 <- sort(unique(survData[survData[,2]==1,1]))
                        s_propBI1 <- s_propBI1[s_propBI1 < s1_max]
                        s_propBI2 <- sort(unique(survData[(survData[,2]==0) & (survData[,4]==1),3]))
                        s_propBI2 <- s_propBI2[s_propBI2 < s2_max]
                        s_propBI3 <- sort(unique(survData[(survData[,2]==1) & (survData[,4]==1),3]))
                        s_propBI3 <- s_propBI3[s_propBI3 < s3_max]
                    }
                    
                    num_s_propBI1=length(s_propBI1)
                    num_s_propBI2=length(s_propBI2)
                    num_s_propBI3=length(s_propBI3)
                    
                    time_lambda1 <- mcmcList$tuning$time_lambda1
                    time_lambda2 <- mcmcList$tuning$time_lambda2
                    time_lambda3 <- mcmcList$tuning$time_lambda3
                    
                    nTime_lambda1 <- length(time_lambda1)
                    nTime_lambda2 <- length(time_lambda2)
                    nTime_lambda3 <- length(time_lambda3)
                    
                    mcmcParams <- c(C1, C2, C3, delPert1, delPert2, delPert3, num_s_propBI1, num_s_propBI2, num_s_propBI3, K1_max, K2_max, K3_max, s1_max, s2_max, s3_max, nTime_lambda1, nTime_lambda2, nTime_lambda3, s_propBI1, s_propBI2, s_propBI3, time_lambda1, time_lambda2, time_lambda3, mhProp_theta_var, mhProp_V1_var, mhProp_V2_var, mhProp_V3_var)
                    
                    s1          <- temp$PEM$PEM.s1
                    s2            <- temp$PEM$PEM.s2
                    s3            <- temp$PEM$PEM.s3
                    
                    lambda1 <- temp$PEM$PEM.lambda1
                    lambda2 <- temp$PEM$PEM.lambda1
                    lambda3 <- temp$PEM$PEM.lambda3
                    
                    K1=temp$PEM$K1
                    K2=temp$PEM$K2
                    K3=temp$PEM$K3
                    mu_lam1=temp$PEM$PEM.mu_lam[1]
                    mu_lam2=temp$PEM$PEM.mu_lam[2]
                    mu_lam3=temp$PEM$PEM.mu_lam[3]
                    sigSq_lam1=temp$PEM$PEM.sigSq_lam[1]
                    sigSq_lam2=temp$PEM$PEM.sigSq_lam[2]
                    sigSq_lam3=temp$PEM$PEM.sigSq_lam[3]
                    
                    
                    # re.type = "MVN"
                    
                    if(re.type == "MVN"){
                        
                        # model = "Markov"
                        
                        #######################################
                        ############ PEM-MVN-M ############
                        #######################################
                        
                        if(model_h3 == "Markov"){
                            
                            ###
                            K_1 <- K1
                            K_2 <- K2
                            K_3 <- K3
                            
                            theta <- startV[(p1+p2+p3+1)]
                            gamma <- startV[(p1+p2+p3+2):(p1+p2+p3+n+1)]
                            
                            if(p1+p2+p3 > 0)
                            {
                                startV <- as.vector(c(startV[1:(p1+p2+p3)],
                                K1, K2, K3,
                                mu_lam1, mu_lam2, mu_lam3,
                                sigSq_lam1, sigSq_lam2, sigSq_lam3,
                                theta, gamma,
                                lambda1, lambda2, lambda3, s1, s2, s3,
                                startV[(p1+p2+p3+n+1+1):length(startV)]))
                            }else
                            {
                                startV <- as.vector(c(K1, K2, K3,
                                mu_lam1, mu_lam2, mu_lam3,
                                sigSq_lam1, sigSq_lam2, sigSq_lam3,
                                theta, gamma,
                                lambda1, lambda2, lambda3, s1, s2, s3,
                                startV[(p1+p2+p3+n+1+1):length(startV)]))
                            }
                            
                            startValues1 <- startV[1:(p1+p2+p3+10+n+2*(K_1+1)+2*(K_2+1)+2*(K_3+1))]
                            startValues2 <- startV[(p1+p2+p3+10+n+2*(K_1+1)+2*(K_2+1)+2*(K_3+1)+1):length(startV)]
                            startV <- c(startValues1, 1, 1, startValues2)
                            
                            mcmcParams1 <- mcmcParams[1:(18+num_s_propBI1+num_s_propBI2+num_s_propBI3+nTime_lambda1+nTime_lambda2+nTime_lambda3+1)]
                            mcmcParams2 <- mcmcParams[(18+num_s_propBI1+num_s_propBI2+num_s_propBI3+nTime_lambda1+nTime_lambda2+nTime_lambda3+2):(length(mcmcParams)-1)]
                            
                            mcmcParams <- c(mcmcParams1, 1, mcmcParams2, n)
                        
                            mcmcRet <- .C("BpeMvnCorScrmcmc",
                            survData         = as.double(as.matrix(survData)),
                            n                = as.integer(n),
                            p1                = as.integer(p1),
                            p2                = as.integer(p2),
                            p3                = as.integer(p3),
                            J                = as.integer(J),
                            nj                = as.double(nj),
                            hyperParams     = as.double(hyperParams),
                            startValues     = as.double(startV),
                            mcmcParams        = as.double(mcmcParams),
                            numReps            = as.integer(numReps),
                            thin            = as.integer(thin),
                            burninPerc      = as.double(burninPerc),
                            nGam_save        = as.integer(nGam_save),
                            samples_beta1     = as.double(rep(0, nStore*p1)),
                            samples_beta2     = as.double(rep(0, nStore*p2)),
                            samples_beta3     = as.double(rep(0, nStore*p3)),
                            samples_mu_lam1     = as.double(rep(0, nStore*1)),
                            samples_mu_lam2     = as.double(rep(0, nStore*1)),
                            samples_mu_lam3     = as.double(rep(0, nStore*1)),
                            samples_sigSq_lam1    = as.double(rep(0, nStore*1)),
                            samples_sigSq_lam2    = as.double(rep(0, nStore*1)),
                            samples_sigSq_lam3    = as.double(rep(0, nStore*1)),
                            samples_K1          = as.double(rep(0, nStore*1)),
                            samples_K2          = as.double(rep(0, nStore*1)),
                            samples_K3          = as.double(rep(0, nStore*1)),
                            samples_s1          = as.double(rep(0, nStore*(K1_max + 1))),
                            samples_s2          = as.double(rep(0, nStore*(K2_max + 1))),
                            samples_s3          = as.double(rep(0, nStore*(K3_max + 1))),
                            samples_nu2     = as.double(rep(0, nStore*1)),
                            samples_nu3     = as.double(rep(0, nStore*1)),
                            samples_theta     = as.double(rep(0, nStore*1)),
                            samples_V1        = as.double(rep(0, nStore*J)),
                            samples_V2        = as.double(rep(0, nStore*J)),
                            samples_V3        = as.double(rep(0, nStore*J)),
                            samples_Sigma_V    = as.double(rep(0, nStore*3*3)),
                            samples_gamma     = as.double(rep(0, nStore*nGam_save)),
                            samples_gamma_last = as.double(rep(0, n)),
                            samples_misc    = as.double(rep(0, (p1+p2+p3+9+n+1+n+J+J+J))),
                            lambda1_fin            = as.double(rep(0, nStore*nTime_lambda1)),
                            lambda2_fin            = as.double(rep(0, nStore*nTime_lambda2)),
                            lambda3_fin            = as.double(rep(0, nStore*nTime_lambda3)),
                            gammaP            = as.double(rep(0, n)),
                            dev                 = as.double(rep(0, nStore*1)),
                            invLH = as.double(rep(0, n)),
                            logLH_fin = as.double(0),
                            lpml                 = as.double(rep(0, nStore*1)),
                            lpml2                 = as.double(rep(0, nStore*1)),
                            moveVec             = as.double(rep(0, numReps)))

                            if(p1 > 0){
                                beta1.p         <- matrix(mcmcRet$samples_beta1, nrow = nStore, byrow = TRUE)
                            }
                            if(p1 == 0){
                                beta1.p         <- NULL
                            }
                            if(p2 > 0){
                                beta2.p         <- matrix(mcmcRet$samples_beta2, nrow = nStore, byrow = TRUE)
                            }
                            if(p2 == 0){
                                beta2.p         <- NULL
                            }
                            if(p3 > 0){
                                beta3.p         <- matrix(mcmcRet$samples_beta3, nrow = nStore, byrow = TRUE)
                            }
                            if(p3 == 0){
                                beta3.p         <- NULL
                            }
                            
                            lambda1.fin     <- matrix(mcmcRet$lambda1_fin, nrow = nStore, byrow = TRUE)
                            lambda2.fin     <- matrix(mcmcRet$lambda2_fin, nrow = nStore, byrow = TRUE)
                            lambda3.fin     <- matrix(mcmcRet$lambda3_fin, nrow = nStore, byrow = TRUE)
                            
                            mu_lam1.p         <- matrix(mcmcRet$samples_mu_lam1, nrow = nStore, byrow = TRUE)
                            mu_lam2.p         <- matrix(mcmcRet$samples_mu_lam2, nrow = nStore, byrow = TRUE)
                            mu_lam3.p         <- matrix(mcmcRet$samples_mu_lam3, nrow = nStore, byrow = TRUE)
                            sigSq_lam1.p     <- matrix(mcmcRet$samples_sigSq_lam1, nrow = nStore, byrow = TRUE)
                            sigSq_lam2.p     <- matrix(mcmcRet$samples_sigSq_lam2, nrow = nStore, byrow = TRUE)
                            sigSq_lam3.p     <- matrix(mcmcRet$samples_sigSq_lam3, nrow = nStore, byrow = TRUE)
                            
                            K1.p             <- matrix(mcmcRet$samples_K1, nrow = nStore, byrow = TRUE)
                            K2.p             <- matrix(mcmcRet$samples_K2, nrow = nStore, byrow = TRUE)
                            K3.p             <- matrix(mcmcRet$samples_K3, nrow = nStore, byrow = TRUE)
                            s1.p             <- matrix(mcmcRet$samples_s1, nrow = nStore, byrow = TRUE)
                            s2.p             <- matrix(mcmcRet$samples_s2, nrow = nStore, byrow = TRUE)
                            s3.p             <- matrix(mcmcRet$samples_s3, nrow = nStore, byrow = TRUE)
                            
                            theta.p         <- matrix(mcmcRet$samples_theta, nrow = nStore, byrow = TRUE)
                            gamma.p         <- matrix(mcmcRet$samples_gamma, nrow = nStore, byrow = TRUE)
                            
                            V1.p            <- matrix(mcmcRet$samples_V1, nrow = nStore, byrow = TRUE)
                            V2.p            <- matrix(mcmcRet$samples_V2, nrow = nStore, byrow = TRUE)
                            V3.p            <- matrix(mcmcRet$samples_V3, nrow = nStore, byrow = TRUE)
                            Sigma_V.p        <- array(as.vector(mcmcRet$samples_Sigma_V), c(3, 3, nStore))
                            
                            if(p1 > 0){
                                accept.beta1     <- as.vector(mcmcRet$samples_misc[1:(p1)])/sum(as.vector(mcmcRet$moveVec)==1)*p1
                            }
                            if(p1 == 0){
                                accept.beta1     <- NULL
                            }
                            if(p2 > 0){
                                accept.beta2     <- as.vector(mcmcRet$samples_misc[(p1+1):(p1+p2)])/sum(as.vector(mcmcRet$moveVec)==2)*p2
                            }
                            if(p2 == 0){
                                accept.beta2     <- NULL
                            }
                            if(p3 > 0){
                                accept.beta3     <- as.vector(mcmcRet$samples_misc[(p1+p2+1):(p1+p2+p3)])/sum(as.vector(mcmcRet$moveVec)==3)*p3
                            }
                            if(p3 == 0){
                                accept.beta3     <- NULL
                            }
                            
                            accept.BI1        <- as.vector(mcmcRet$samples_misc[(p1+p2+p3)+1])/sum(as.vector(mcmcRet$moveVec)==12)
                            accept.DI1        <- as.vector(mcmcRet$samples_misc[(p1+p2+p3)+2])/sum(as.vector(mcmcRet$moveVec)==15)
                            accept.BI2        <- as.vector(mcmcRet$samples_misc[(p1+p2+p3)+3])/sum(as.vector(mcmcRet$moveVec)==13)
                            accept.DI2        <- as.vector(mcmcRet$samples_misc[(p1+p2+p3)+4])/sum(as.vector(mcmcRet$moveVec)==16)
                            accept.BI3        <- as.vector(mcmcRet$samples_misc[(p1+p2+p3)+5])/sum(as.vector(mcmcRet$moveVec)==14)
                            accept.DI3        <- as.vector(mcmcRet$samples_misc[(p1+p2+p3)+6])/sum(as.vector(mcmcRet$moveVec)==17)
                            accept.theta    <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+7)])/sum(as.vector(mcmcRet$moveVec)==11)
                            accept.V1       <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+10+n+n+1):(p1+p2+p3+10+n+n+J)])/sum(as.vector(mcmcRet$moveVec)==18)
                            accept.V2       <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+10+n+n+J+1):(p1+p2+p3+10+n+n+J+J)])/sum(as.vector(mcmcRet$moveVec)==19)
                            accept.V3       <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+10+n+n+J+J+1):(p1+p2+p3+10+n+n+J+J+J)])/sum(as.vector(mcmcRet$moveVec)==20)
                            

                            V1summary <- as.matrix(apply(V1.p, 2, summary))
                            V1summary <- rbind(V1summary, apply(V1.p, 2, quantile, prob = 0.975))
                            V1summary <- rbind(V1summary, apply(V1.p, 2, quantile, prob = 0.025))
                            V1summary <- rbind(V1summary, apply(V1.p, 2, sd))
                            rownames(V1summary)[7:9] <- c("0.975", "0.025", "sd")
                            
                            V2summary <- as.matrix(apply(V2.p, 2, summary))
                            V2summary <- rbind(V2summary, apply(V2.p, 2, quantile, prob = 0.975))
                            V2summary <- rbind(V2summary, apply(V2.p, 2, quantile, prob = 0.025))
                            V2summary <- rbind(V2summary, apply(V2.p, 2, sd))
                            rownames(V2summary)[7:9] <- c("0.975", "0.025", "sd")
                            
                            V3summary <- as.matrix(apply(V3.p, 2, summary))
                            V3summary <- rbind(V3summary, apply(V3.p, 2, quantile, prob = 0.975))
                            V3summary <- rbind(V3summary, apply(V3.p, 2, quantile, prob = 0.025))
                            V3summary <- rbind(V3summary, apply(V3.p, 2, sd))
                            rownames(V3summary)[7:9] <- c("0.975", "0.025", "sd")
                            
                            if(storeV[1] == TRUE & !is.null(path))
                            {
                                save(V1.p, file = paste(path, "V1Pch", chain, ".RData", sep = ""))
                            }
                            if(storeV[2] == TRUE & !is.null(path))
                            {
                                save(V2.p, file = paste(path, "V2Pch", chain, ".RData", sep = ""))
                            }
                            if(storeV[3] == TRUE & !is.null(path))
                            {
                                save(V3.p, file = paste(path, "V3Pch", chain, ".RData", sep = ""))
                            }
                            
                            if(nGam_save > 0 & !is.null(path)){
                                save(gamma.p, file = paste(path, "/gammaPch", chain, ".Rdata", sep = ""))
                            }
                            
                            if(p1 > 0){
                                covNames1 = colnames(Xmat1)
                            }
                            if(p1 == 0){
                                covNames1 = NULL
                            }
                            
                            if(p2 > 0){
                                covNames2 = colnames(Xmat2)
                            }
                            if(p2 == 0){
                                covNames2 = NULL
                            }
                            if(p3 > 0){
                                covNames3 = colnames(Xmat3)
                            }
                            if(p3 == 0){
                                covNames3 = NULL
                            }
                            
                            
                            
                            ### posterior predictive checks ###
                            
                            ## 1. log pseudo-marginal likelihood
                            
                            invLH.p <- matrix(mcmcRet$invLH, nrow = n, byrow = TRUE)
                            
                            ## 2. deviance information criterion
                            
                            gamma_mean <- matrix(mcmcRet$gammaP, nrow = n, byrow = TRUE)
                            dev.p        <- matrix(mcmcRet$dev, nrow = nStore, byrow = TRUE)
       
                            ret[[nam]] <- list(beta1.p = beta1.p, beta2.p = beta2.p, beta3.p = beta3.p, lambda1.fin = lambda1.fin, lambda2.fin = lambda2.fin, lambda3.fin = lambda3.fin, mu_lam1.p = mu_lam1.p, mu_lam2.p = mu_lam2.p, mu_lam3.p = mu_lam3.p, sigSq_lam1.p = sigSq_lam1.p, sigSq_lam2.p = sigSq_lam2.p, sigSq_lam3.p = sigSq_lam3.p, K1.p = K1.p, K2.p = K2.p, K3.p = K3.p, s1.p = s1.p, s2.p = s2.p, s3.p = s3.p, theta.p = theta.p, Sigma_V.p = Sigma_V.p, accept.beta1 = accept.beta1, accept.beta2 = accept.beta2, accept.beta3 = accept.beta3, accept.BI1 = accept.BI1, accept.BI2 = accept.BI2, accept.BI3 = accept.BI3, accept.DI1 = accept.DI1, accept.DI2 = accept.DI2, accept.DI3 = accept.DI3, accept.theta = accept.theta, time_lambda1 = time_lambda1, time_lambda2 = time_lambda2, time_lambda3 = time_lambda3, accept.V1 = accept.V1, accept.V2 = accept.V2, accept.V3 = accept.V3, covNames1 = covNames1, covNames2 = covNames2, covNames3 = covNames3, V1sum = V1summary, V2sum = V2summary, V3sum = V3summary, gamma_mean = gamma_mean, dev.p=dev.p, invLH.p=invLH.p)
                            
                            
                        }   ## end: if PEM-MVN-M
                        
                        # model = "semi-Markov"
                        
                        #######################################
                        ############ PEM-MVN-SM ############
                        #######################################
                        
                        if(model_h3 == "semi-Markov"){
                            
                            ###
                            
                            K_1 <- K1
                            K_2 <- K2
                            K_3 <- K3
                            
                            
                            theta <- startV[(p1+p2+p3+1)]
                            gamma <- startV[(p1+p2+p3+2):(p1+p2+p3+n+1)]
                            
                            if(p1+p2+p3 > 0)
                            {
                                startV <- as.vector(c(startV[1:(p1+p2+p3)],
                                K1, K2, K3,
                                mu_lam1, mu_lam2, mu_lam3,
                                sigSq_lam1, sigSq_lam2, sigSq_lam3,
                                theta, gamma,
                                lambda1, lambda2, lambda3, s1, s2, s3,
                                startV[(p1+p2+p3+n+1+1):length(startV)]))
                            }else
                            {
                                startV <- as.vector(c(K1, K2, K3,
                                mu_lam1, mu_lam2, mu_lam3,
                                sigSq_lam1, sigSq_lam2, sigSq_lam3,
                                theta, gamma,
                                lambda1, lambda2, lambda3, s1, s2, s3,
                                startV[(p1+p2+p3+n+1+1):length(startV)]))
                            }
                            
                            
                            startValues1 <- startV[1:(p1+p2+p3+10+n+2*(K_1+1)+2*(K_2+1)+2*(K_3+1))]
                            startValues2 <- startV[(p1+p2+p3+10+n+2*(K_1+1)+2*(K_2+1)+2*(K_3+1)+1):length(startV)]
                            startV <- c(startValues1, 1, 1, startValues2)
                            
                            mcmcParams1 <- mcmcParams[1:(18+num_s_propBI1+num_s_propBI2+num_s_propBI3+nTime_lambda1+nTime_lambda2+nTime_lambda3+1)]
                            mcmcParams2 <- mcmcParams[(18+num_s_propBI1+num_s_propBI2+num_s_propBI3+nTime_lambda1+nTime_lambda2+nTime_lambda3+2):(length(mcmcParams)-1)]
                            
                            mcmcParams <- c(mcmcParams1, 1, mcmcParams2, n)
                            
                            
                            
                            mcmcRet <- .C("BpeMvnCorScrSMmcmc",
                            survData         = as.double(as.matrix(survData)),
                            n                = as.integer(n),
                            p1                = as.integer(p1),
                            p2                = as.integer(p2),
                            p3                = as.integer(p3),
                            J                = as.integer(J),
                            nj                = as.double(nj),
                            hyperParams     = as.double(hyperParams),
                            startValues     = as.double(startV),
                            mcmcParams        = as.double(mcmcParams),
                            numReps            = as.integer(numReps),
                            thin            = as.integer(thin),
                            burninPerc      = as.double(burninPerc),
                            nGam_save        = as.integer(nGam_save),
                            samples_beta1     = as.double(rep(0, nStore*p1)),
                            samples_beta2     = as.double(rep(0, nStore*p2)),
                            samples_beta3     = as.double(rep(0, nStore*p3)),
                            samples_mu_lam1     = as.double(rep(0, nStore*1)),
                            samples_mu_lam2     = as.double(rep(0, nStore*1)),
                            samples_mu_lam3     = as.double(rep(0, nStore*1)),
                            samples_sigSq_lam1    = as.double(rep(0, nStore*1)),
                            samples_sigSq_lam2    = as.double(rep(0, nStore*1)),
                            samples_sigSq_lam3    = as.double(rep(0, nStore*1)),
                            samples_K1          = as.double(rep(0, nStore*1)),
                            samples_K2          = as.double(rep(0, nStore*1)),
                            samples_K3          = as.double(rep(0, nStore*1)),
                            samples_s1          = as.double(rep(0, nStore*(K1_max + 1))),
                            samples_s2          = as.double(rep(0, nStore*(K2_max + 1))),
                            samples_s3          = as.double(rep(0, nStore*(K3_max + 1))),
                            samples_nu2     = as.double(rep(0, nStore*1)),
                            samples_nu3     = as.double(rep(0, nStore*1)),
                            samples_theta     = as.double(rep(0, nStore*1)),
                            samples_V1        = as.double(rep(0, nStore*J)),
                            samples_V2        = as.double(rep(0, nStore*J)),
                            samples_V3        = as.double(rep(0, nStore*J)),
                            samples_Sigma_V    = as.double(rep(0, nStore*3*3)),
                            samples_gamma     = as.double(rep(0, nStore*nGam_save)),
                            samples_gamma_last = as.double(rep(0, n)),
                            samples_misc    = as.double(rep(0, (p1+p2+p3+9+n+1+n+J+J+J))),
                            lambda1_fin            = as.double(rep(0, nStore*nTime_lambda1)),
                            lambda2_fin            = as.double(rep(0, nStore*nTime_lambda2)),
                            lambda3_fin            = as.double(rep(0, nStore*nTime_lambda3)),
                            gammaP            = as.double(rep(0, n)),
                            dev                 = as.double(rep(0, nStore*1)),
                            invLH = as.double(rep(0, n)),
                            logLH_fin = as.double(0),
                            lpml                 = as.double(rep(0, nStore*1)),
                            moveVec             = as.double(rep(0, numReps)))
                            
                            
                            
                            if(p1 > 0){
                                beta1.p         <- matrix(mcmcRet$samples_beta1, nrow = nStore, byrow = TRUE)
                            }
                            if(p1 == 0){
                                beta1.p         <- NULL
                            }
                            if(p2 > 0){
                                beta2.p         <- matrix(mcmcRet$samples_beta2, nrow = nStore, byrow = TRUE)
                            }
                            if(p2 == 0){
                                beta2.p         <- NULL
                            }
                            if(p3 > 0){
                                beta3.p         <- matrix(mcmcRet$samples_beta3, nrow = nStore, byrow = TRUE)
                            }
                            if(p3 == 0){
                                beta3.p         <- NULL
                            }
                            
                            lambda1.fin     <- matrix(mcmcRet$lambda1_fin, nrow = nStore, byrow = TRUE)
                            lambda2.fin     <- matrix(mcmcRet$lambda2_fin, nrow = nStore, byrow = TRUE)
                            lambda3.fin     <- matrix(mcmcRet$lambda3_fin, nrow = nStore, byrow = TRUE)
                            
                            mu_lam1.p         <- matrix(mcmcRet$samples_mu_lam1, nrow = nStore, byrow = TRUE)
                            mu_lam2.p         <- matrix(mcmcRet$samples_mu_lam2, nrow = nStore, byrow = TRUE)
                            mu_lam3.p         <- matrix(mcmcRet$samples_mu_lam3, nrow = nStore, byrow = TRUE)
                            sigSq_lam1.p     <- matrix(mcmcRet$samples_sigSq_lam1, nrow = nStore, byrow = TRUE)
                            sigSq_lam2.p     <- matrix(mcmcRet$samples_sigSq_lam2, nrow = nStore, byrow = TRUE)
                            sigSq_lam3.p     <- matrix(mcmcRet$samples_sigSq_lam3, nrow = nStore, byrow = TRUE)
                            
                            K1.p             <- matrix(mcmcRet$samples_K1, nrow = nStore, byrow = TRUE)
                            K2.p             <- matrix(mcmcRet$samples_K2, nrow = nStore, byrow = TRUE)
                            K3.p             <- matrix(mcmcRet$samples_K3, nrow = nStore, byrow = TRUE)
                            s1.p             <- matrix(mcmcRet$samples_s1, nrow = nStore, byrow = TRUE)
                            s2.p             <- matrix(mcmcRet$samples_s2, nrow = nStore, byrow = TRUE)
                            s3.p             <- matrix(mcmcRet$samples_s3, nrow = nStore, byrow = TRUE)
                            
                            theta.p         <- matrix(mcmcRet$samples_theta, nrow = nStore, byrow = TRUE)
                            gamma.p         <- matrix(mcmcRet$samples_gamma, nrow = nStore, byrow = TRUE)
                            
                            V1.p            <- matrix(mcmcRet$samples_V1, nrow = nStore, byrow = TRUE)
                            V2.p            <- matrix(mcmcRet$samples_V2, nrow = nStore, byrow = TRUE)
                            V3.p            <- matrix(mcmcRet$samples_V3, nrow = nStore, byrow = TRUE)
                            Sigma_V.p        <- array(as.vector(mcmcRet$samples_Sigma_V), c(3, 3, nStore))
                            
                            if(p1 > 0){
                                accept.beta1     <- as.vector(mcmcRet$samples_misc[1:(p1)])/sum(as.vector(mcmcRet$moveVec)==1)*p1
                            }
                            if(p1 == 0){
                                accept.beta1     <- NULL
                            }
                            if(p2 > 0){
                                accept.beta2     <- as.vector(mcmcRet$samples_misc[(p1+1):(p1+p2)])/sum(as.vector(mcmcRet$moveVec)==2)*p2
                            }
                            if(p2 == 0){
                                accept.beta2     <- NULL
                            }
                            if(p3 > 0){
                                accept.beta3     <- as.vector(mcmcRet$samples_misc[(p1+p2+1):(p1+p2+p3)])/sum(as.vector(mcmcRet$moveVec)==3)*p3
                            }
                            if(p3 == 0){
                                accept.beta3     <- NULL
                            }
                            
                            accept.BI1        <- as.vector(mcmcRet$samples_misc[(p1+p2+p3)+1])/sum(as.vector(mcmcRet$moveVec)==12)
                            accept.DI1        <- as.vector(mcmcRet$samples_misc[(p1+p2+p3)+2])/sum(as.vector(mcmcRet$moveVec)==15)
                            accept.BI2        <- as.vector(mcmcRet$samples_misc[(p1+p2+p3)+3])/sum(as.vector(mcmcRet$moveVec)==13)
                            accept.DI2        <- as.vector(mcmcRet$samples_misc[(p1+p2+p3)+4])/sum(as.vector(mcmcRet$moveVec)==16)
                            accept.BI3        <- as.vector(mcmcRet$samples_misc[(p1+p2+p3)+5])/sum(as.vector(mcmcRet$moveVec)==14)
                            accept.DI3        <- as.vector(mcmcRet$samples_misc[(p1+p2+p3)+6])/sum(as.vector(mcmcRet$moveVec)==17)
                            accept.theta    <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+7)])/sum(as.vector(mcmcRet$moveVec)==11)
                            accept.V1       <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+10+n+n+1):(p1+p2+p3+10+n+n+J)])/sum(as.vector(mcmcRet$moveVec)==18)
                            accept.V2       <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+10+n+n+J+1):(p1+p2+p3+10+n+n+J+J)])/sum(as.vector(mcmcRet$moveVec)==19)
                            accept.V3       <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+10+n+n+J+J+1):(p1+p2+p3+10+n+n+J+J+J)])/sum(as.vector(mcmcRet$moveVec)==20)
                            
                            
                            V1summary <- as.matrix(apply(V1.p, 2, summary))
                            V1summary <- rbind(V1summary, apply(V1.p, 2, quantile, prob = 0.975))
                            V1summary <- rbind(V1summary, apply(V1.p, 2, quantile, prob = 0.025))
                            V1summary <- rbind(V1summary, apply(V1.p, 2, sd))
                            rownames(V1summary)[7:9] <- c("0.975", "0.025", "sd")
                            
                            V2summary <- as.matrix(apply(V2.p, 2, summary))
                            V2summary <- rbind(V2summary, apply(V2.p, 2, quantile, prob = 0.975))
                            V2summary <- rbind(V2summary, apply(V2.p, 2, quantile, prob = 0.025))
                            V2summary <- rbind(V2summary, apply(V2.p, 2, sd))
                            rownames(V2summary)[7:9] <- c("0.975", "0.025", "sd")
                            
                            V3summary <- as.matrix(apply(V3.p, 2, summary))
                            V3summary <- rbind(V3summary, apply(V3.p, 2, quantile, prob = 0.975))
                            V3summary <- rbind(V3summary, apply(V3.p, 2, quantile, prob = 0.025))
                            V3summary <- rbind(V3summary, apply(V3.p, 2, sd))
                            rownames(V3summary)[7:9] <- c("0.975", "0.025", "sd")
                            
                            if(storeV[1] == TRUE & !is.null(path))
                            {
                                save(V1.p, file = paste(path, "V1Pch", chain, ".RData", sep = ""))
                            }
                            if(storeV[2] == TRUE & !is.null(path))
                            {
                                save(V2.p, file = paste(path, "V2Pch", chain, ".RData", sep = ""))
                            }
                            if(storeV[3] == TRUE & !is.null(path))
                            {
                                save(V3.p, file = paste(path, "V3Pch", chain, ".RData", sep = ""))
                            }
                            
                            if(nGam_save > 0 & !is.null(path)){
                                save(gamma.p, file = paste(path, "/gammaPch", chain, ".Rdata", sep = ""))
                            }
                            
                            if(p1 > 0){
                                covNames1 = colnames(Xmat1)
                            }
                            if(p1 == 0){
                                covNames1 = NULL
                            }
                            
                            if(p2 > 0){
                                covNames2 = colnames(Xmat2)
                            }
                            if(p2 == 0){
                                covNames2 = NULL
                            }
                            if(p3 > 0){
                                covNames3 = colnames(Xmat3)
                            }
                            if(p3 == 0){
                                covNames3 = NULL
                            }
                            
                            
                            
                            ### posterior predictive checks ###
                            
                            ## 1. log pseudo-marginal likelihood
                            
                            invLH.p <- matrix(mcmcRet$invLH, nrow = n, byrow = TRUE)
                            
                            ## 2. deviance information criterion
                            
                            gamma_mean <- matrix(mcmcRet$gammaP, nrow = n, byrow = TRUE)
                            dev.p        <- matrix(mcmcRet$dev, nrow = nStore, byrow = TRUE)
                            
                            ret[[nam]] <- list(beta1.p = beta1.p, beta2.p = beta2.p, beta3.p = beta3.p, lambda1.fin = lambda1.fin, lambda2.fin = lambda2.fin, lambda3.fin = lambda3.fin, mu_lam1.p = mu_lam1.p, mu_lam2.p = mu_lam2.p, mu_lam3.p = mu_lam3.p, sigSq_lam1.p = sigSq_lam1.p, sigSq_lam2.p = sigSq_lam2.p, sigSq_lam3.p = sigSq_lam3.p, K1.p = K1.p, K2.p = K2.p, K3.p = K3.p, s1.p = s1.p, s2.p = s2.p, s3.p = s3.p, theta.p = theta.p, Sigma_V.p = Sigma_V.p, accept.beta1 = accept.beta1, accept.beta2 = accept.beta2, accept.beta3 = accept.beta3, accept.BI1 = accept.BI1, accept.BI2 = accept.BI2, accept.BI3 = accept.BI3, accept.DI1 = accept.DI1, accept.DI2 = accept.DI2, accept.DI3 = accept.DI3, accept.theta = accept.theta, time_lambda1 = time_lambda1, time_lambda2 = time_lambda2, time_lambda3 = time_lambda3, accept.V1 = accept.V1, accept.V2 = accept.V2, accept.V3 = accept.V3, covNames1 = covNames1, covNames2 = covNames2, covNames3 = covNames3, V1sum = V1summary, V2sum = V2summary, V3sum = V3summary, gamma_mean = gamma_mean, dev.p=dev.p, invLH.p=invLH.p)
                            
                            
                        }   ## end: if PEM-MVN-SM
                        
                        ## DIC and LPML
                        
                        be1.p <- ret$chain1$beta1.p
                        be2.p <- ret$chain1$beta2.p
                        be3.p <- ret$chain1$beta3.p
                        thet.p <- ret$chain1$theta.p
                        lam1.p <- ret$chain1$lambda1.fin
                        lam2.p <- ret$chain1$lambda2.fin
                        lam3.p <- ret$chain1$lambda3.fin
                        V1me <- ret$chain1$V1sum[4,]
                        V2me <- ret$chain1$V2sum[4,]
                        V3me <- ret$chain1$V3sum[4,]
                        Dev.p <- ret$chain1$dev.p
                        InvLH.p <- ret$chain1$invLH.p
                        
                        if(nChain > 1){
                            for(i in 2:nChain){
                                nam <- paste("chain", i, sep="")
                                be1.p <- rbind(be1.p, ret[[nam]]$beta1.p)
                                be2.p <- rbind(be2.p, ret[[nam]]$beta2.p)
                                be3.p <- rbind(be3.p, ret[[nam]]$beta3.p)
                                thet.p <- rbind(thet.p, ret[[nam]]$theta.p)
                                lam1.p <- rbind(lam1.p, ret[[nam]]$lambda1.fin)
                                lam2.p <- rbind(lam2.p, ret[[nam]]$lambda2.fin)
                                lam3.p <- rbind(lam3.p, ret[[nam]]$lambda3.fin)
                                V1me <- rbind(V1me, ret[[nam]]$V1sum[4,])
                                V2me <- rbind(V2me, ret[[nam]]$V2sum[4,])
                                V3me <- rbind(V3me, ret[[nam]]$V3sum[4,])
                                Dev.p <- rbind(Dev.p, ret[[nam]]$dev.p)
                                InvLH.p <- cbind(InvLH.p, ret[[nam]]$invLH.p)
                            }
                        }
                        if(model_h3 == "Markov")
                        {
                            logLH_fin <- .C("BpeMvnCorScr_logMLH_DIC",
                            survData         = as.double(as.matrix(survData)),
                            n                = as.integer(n),
                            p1                = as.integer(p1),
                            p2                = as.integer(p2),
                            p3                = as.integer(p3),
                            J                = as.integer(J),
                            be1 = as.double(apply(be1.p, 2, mean)),
                            be2 = as.double(apply(be2.p, 2, mean)),
                            be3 = as.double(apply(be3.p, 2, mean)),
                            thet = as.double(apply(theta.p, 2, mean)),
                            lam1 = as.double(apply(lam1.p, 2, mean)),
                            lam2 = as.double(apply(lam2.p, 2, mean)),
                            lam3 = as.double(apply(lam3.p, 2, mean)),
                            s_1 = as.double(time_lambda1),
                            s_2 = as.double(time_lambda2),
                            s_3 = as.double(time_lambda3),
                            V_1 = as.double(apply(V1me, 2, mean)),
                            V_2 = as.double(apply(V2me, 2, mean)),
                            V_3 = as.double(apply(V3me, 2, mean)),
                            K_1                = as.integer(nTime_lambda1),
                            K_2                = as.integer(nTime_lambda2),
                            K_3                = as.integer(nTime_lambda3),
                            val = as.double(0))$val
                        }else if(model_h3 == "semi-Markov")
                        {
                            logLH_fin <- .C("BpeMvnCorScrSM_logMLH_DIC",
                            survData         = as.double(as.matrix(survData)),
                            n                = as.integer(n),
                            p1                = as.integer(p1),
                            p2                = as.integer(p2),
                            p3                = as.integer(p3),
                            J                = as.integer(J),
                            be1 = as.double(apply(be1.p, 2, mean)),
                            be2 = as.double(apply(be2.p, 2, mean)),
                            be3 = as.double(apply(be3.p, 2, mean)),
                            thet = as.double(apply(theta.p, 2, mean)),
                            lam1 = as.double(apply(lam1.p, 2, mean)),
                            lam2 = as.double(apply(lam2.p, 2, mean)),
                            lam3 = as.double(apply(lam3.p, 2, mean)),
                            s_1 = as.double(time_lambda1),
                            s_2 = as.double(time_lambda2),
                            s_3 = as.double(time_lambda3),
                            V_1 = as.double(apply(V1me, 2, mean)),
                            V_2 = as.double(apply(V2me, 2, mean)),
                            V_3 = as.double(apply(V3me, 2, mean)),
                            K_1                = as.integer(nTime_lambda1),
                            K_2                = as.integer(nTime_lambda2),
                            K_3                = as.integer(nTime_lambda3),
                            val = as.double(0))$val
                        }
                        
                        DIC = 2*mean(Dev.p) + 2 * logLH_fin
                        LPML = sum(log(1/apply(InvLH.p, 1, mean)))
                        
                    }   ## end: if PEM-MVN
                    
                    
                    # re.type = "DPM"
                    
                    if(re.type == "DPM"){
                        
                        # model = "Markov"
                        
                        #######################################
                        ############ PEM-DPM-M ############
                        #######################################
                        
                        if(model_h3 == "Markov"){
                            
                            ###
                            
                            K_1 <- K1
                            K_2 <- K2
                            K_3 <- K3
                            
                            theta <- startV[(p1+p2+p3+1)]
                            gamma <- startV[(p1+p2+p3+2):(p1+p2+p3+n+1)]
                            
                            if(p1+p2+p3 > 0)
                            {
                                startV <- as.vector(c(startV[1:(p1+p2+p3)],
                                K1, K2, K3,
                                mu_lam1, mu_lam2, mu_lam3,
                                sigSq_lam1, sigSq_lam2, sigSq_lam3,
                                theta, gamma,
                                lambda1, lambda2, lambda3, s1, s2, s3,
                                startV[(p1+p2+p3+n+1+1):length(startV)]))
                            }else
                            {
                                startV <- as.vector(c(K1, K2, K3,
                                mu_lam1, mu_lam2, mu_lam3,
                                sigSq_lam1, sigSq_lam2, sigSq_lam3,
                                theta, gamma,
                                lambda1, lambda2, lambda3, s1, s2, s3,
                                startV[(p1+p2+p3+n+1+1):length(startV)]))
                            }
                            
                            
                            startValues1 <- startV[1:(p1+p2+p3+10+n+2*(K_1+1)+2*(K_2+1)+2*(K_3+1))]
                            startValues2 <- startV[(p1+p2+p3+10+n+2*(K_1+1)+2*(K_2+1)+2*(K_3+1)+1):length(startV)]
                            startV <- c(startValues1, 1, 1, startValues2)
                            
                            mcmcParams1 <- mcmcParams[1:(18+num_s_propBI1+num_s_propBI2+num_s_propBI3+nTime_lambda1+nTime_lambda2+nTime_lambda3+1)]
                            mcmcParams2 <- mcmcParams[(18+num_s_propBI1+num_s_propBI2+num_s_propBI3+nTime_lambda1+nTime_lambda2+nTime_lambda3+2):(length(mcmcParams)-1)]
                            
                            mcmcParams <- c(mcmcParams1, 1, mcmcParams2, n)
                            
                            
                            mcmcRet <- .C("BpeDpCorScrmcmc",
                            survData         = as.double(as.matrix(survData)),
                            n                = as.integer(n),
                            p1                = as.integer(p1),
                            p2                = as.integer(p2),
                            p3                = as.integer(p3),
                            J                = as.integer(J),
                            nj                = as.double(nj),
                            hyperParams     = as.double(hyperParams),
                            startValues     = as.double(startV),
                            mcmcParams        = as.double(mcmcParams),
                            numReps            = as.integer(numReps),
                            thin            = as.integer(thin),
                            burninPerc      = as.double(burninPerc),
                            nGam_save        = as.integer(nGam_save),
                            samples_beta1     = as.double(rep(0, nStore*p1)),
                            samples_beta2     = as.double(rep(0, nStore*p2)),
                            samples_beta3     = as.double(rep(0, nStore*p3)),
                            samples_mu_lam1     = as.double(rep(0, nStore*1)),
                            samples_mu_lam2     = as.double(rep(0, nStore*1)),
                            samples_mu_lam3     = as.double(rep(0, nStore*1)),
                            samples_sigSq_lam1    = as.double(rep(0, nStore*1)),
                            samples_sigSq_lam2    = as.double(rep(0, nStore*1)),
                            samples_sigSq_lam3    = as.double(rep(0, nStore*1)),
                            samples_K1          = as.double(rep(0, nStore*1)),
                            samples_K2          = as.double(rep(0, nStore*1)),
                            samples_K3          = as.double(rep(0, nStore*1)),
                            samples_s1          = as.double(rep(0, nStore*(K1_max + 1))),
                            samples_s2          = as.double(rep(0, nStore*(K2_max + 1))),
                            samples_s3          = as.double(rep(0, nStore*(K3_max + 1))),
                            samples_nu2     = as.double(rep(0, nStore*1)),
                            samples_nu3     = as.double(rep(0, nStore*1)),
                            samples_theta     = as.double(rep(0, nStore*1)),
                            samples_V1        = as.double(rep(0, nStore*J)),
                            samples_V2        = as.double(rep(0, nStore*J)),
                            samples_V3        = as.double(rep(0, nStore*J)),
                            samples_c        = as.double(rep(0, nStore*J)),
                            samples_mu        = as.double(rep(0, nStore*3*J)),
                            samples_Sigma    = as.double(rep(0, nStore*3*3*J)),
                            samples_tau     = as.double(rep(0, nStore*1)),
                            samples_gamma     = as.double(rep(0, nStore*nGam_save)),
                            samples_gamma_last = as.double(rep(0, n)),
                            samples_misc    = as.double(rep(0, (p1+p2+p3+9+n+1+n+J))),
                            lambda1_fin            = as.double(rep(0, nStore*nTime_lambda1)),
                            lambda2_fin            = as.double(rep(0, nStore*nTime_lambda2)),
                            lambda3_fin            = as.double(rep(0, nStore*nTime_lambda3)),
                            gammaP            = as.double(rep(0, n)),
                            dev                 = as.double(rep(0, nStore*1)),
                            invLH = as.double(rep(0, n)),
                            logLH_fin = as.double(0),
                            lpml                 = as.double(rep(0, nStore*1)),
                            moveVec             = as.double(rep(0, numReps)))
                            
                            
                            
                            if(p1 > 0){
                                beta1.p         <- matrix(mcmcRet$samples_beta1, nrow = nStore, byrow = TRUE)
                            }
                            if(p1 == 0){
                                beta1.p         <- NULL
                            }
                            if(p2 > 0){
                                beta2.p         <- matrix(mcmcRet$samples_beta2, nrow = nStore, byrow = TRUE)
                            }
                            if(p2 == 0){
                                beta2.p         <- NULL
                            }
                            if(p3 > 0){
                                beta3.p         <- matrix(mcmcRet$samples_beta3, nrow = nStore, byrow = TRUE)
                            }
                            if(p3 == 0){
                                beta3.p         <- NULL
                            }
                            
                            lambda1.fin     <- matrix(mcmcRet$lambda1_fin, nrow = nStore, byrow = TRUE)
                            lambda2.fin     <- matrix(mcmcRet$lambda2_fin, nrow = nStore, byrow = TRUE)
                            lambda3.fin     <- matrix(mcmcRet$lambda3_fin, nrow = nStore, byrow = TRUE)
                            
                            mu_lam1.p         <- matrix(mcmcRet$samples_mu_lam1, nrow = nStore, byrow = TRUE)
                            mu_lam2.p         <- matrix(mcmcRet$samples_mu_lam2, nrow = nStore, byrow = TRUE)
                            mu_lam3.p         <- matrix(mcmcRet$samples_mu_lam3, nrow = nStore, byrow = TRUE)
                            sigSq_lam1.p     <- matrix(mcmcRet$samples_sigSq_lam1, nrow = nStore, byrow = TRUE)
                            sigSq_lam2.p     <- matrix(mcmcRet$samples_sigSq_lam2, nrow = nStore, byrow = TRUE)
                            sigSq_lam3.p     <- matrix(mcmcRet$samples_sigSq_lam3, nrow = nStore, byrow = TRUE)
                            
                            K1.p             <- matrix(mcmcRet$samples_K1, nrow = nStore, byrow = TRUE)
                            K2.p             <- matrix(mcmcRet$samples_K2, nrow = nStore, byrow = TRUE)
                            K3.p             <- matrix(mcmcRet$samples_K3, nrow = nStore, byrow = TRUE)
                            s1.p             <- matrix(mcmcRet$samples_s1, nrow = nStore, byrow = TRUE)
                            s2.p             <- matrix(mcmcRet$samples_s2, nrow = nStore, byrow = TRUE)
                            s3.p             <- matrix(mcmcRet$samples_s3, nrow = nStore, byrow = TRUE)
                            
                            theta.p         <- matrix(mcmcRet$samples_theta, nrow = nStore, byrow = TRUE)
                            gamma.p         <- matrix(mcmcRet$samples_gamma, nrow = nStore, byrow = TRUE)
                            
                            V1.p            <- matrix(mcmcRet$samples_V1, nrow = nStore, byrow = TRUE)
                            V2.p            <- matrix(mcmcRet$samples_V2, nrow = nStore, byrow = TRUE)
                            V3.p            <- matrix(mcmcRet$samples_V3, nrow = nStore, byrow = TRUE)
                            c.p            <- matrix(mcmcRet$samples_c, nrow = nStore, byrow = TRUE)
                            mu.p            <- matrix(mcmcRet$samples_mu, nrow = nStore, byrow = TRUE)
                            Sigma.p         <- array(as.vector(mcmcRet$samples_Sigma), c(3, 3 * J, nStore))
                            tau.p           <- matrix(mcmcRet$samples_tau, nrow = nStore, byrow = TRUE)
                            
                            if(p1 > 0){
                                accept.beta1     <- as.vector(mcmcRet$samples_misc[1:(p1)])/sum(as.vector(mcmcRet$moveVec)==1)*p1
                            }
                            if(p1 == 0){
                                accept.beta1     <- NULL
                            }
                            if(p2 > 0){
                                accept.beta2     <- as.vector(mcmcRet$samples_misc[(p1+1):(p1+p2)])/sum(as.vector(mcmcRet$moveVec)==2)*p2
                            }
                            if(p2 == 0){
                                accept.beta2     <- NULL
                            }
                            if(p3 > 0){
                                accept.beta3     <- as.vector(mcmcRet$samples_misc[(p1+p2+1):(p1+p2+p3)])/sum(as.vector(mcmcRet$moveVec)==3)*p3
                            }
                            if(p3 == 0){
                                accept.beta3     <- NULL
                            }
                            
                            accept.BI1        <- as.vector(mcmcRet$samples_misc[(p1+p2+p3)+1])/sum(as.vector(mcmcRet$moveVec)==12)
                            accept.DI1        <- as.vector(mcmcRet$samples_misc[(p1+p2+p3)+2])/sum(as.vector(mcmcRet$moveVec)==15)
                            accept.BI2        <- as.vector(mcmcRet$samples_misc[(p1+p2+p3)+3])/sum(as.vector(mcmcRet$moveVec)==13)
                            accept.DI2        <- as.vector(mcmcRet$samples_misc[(p1+p2+p3)+4])/sum(as.vector(mcmcRet$moveVec)==16)
                            accept.BI3        <- as.vector(mcmcRet$samples_misc[(p1+p2+p3)+5])/sum(as.vector(mcmcRet$moveVec)==14)
                            accept.DI3        <- as.vector(mcmcRet$samples_misc[(p1+p2+p3)+6])/sum(as.vector(mcmcRet$moveVec)==17)
                            accept.theta    <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+7)])/sum(as.vector(mcmcRet$moveVec)==11)
                            accept.V       <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+10+n+n+1):(p1+p2+p3+10+n+n+J)])/sum(as.vector(mcmcRet$moveVec)==18)/3
                            
                            
                            V1summary <- as.matrix(apply(V1.p, 2, summary))
                            V1summary <- rbind(V1summary, apply(V1.p, 2, quantile, prob = 0.975))
                            V1summary <- rbind(V1summary, apply(V1.p, 2, quantile, prob = 0.025))
                            V1summary <- rbind(V1summary, apply(V1.p, 2, sd))
                            rownames(V1summary)[7:9] <- c("0.975", "0.025", "sd")
                            
                            V2summary <- as.matrix(apply(V2.p, 2, summary))
                            V2summary <- rbind(V2summary, apply(V2.p, 2, quantile, prob = 0.975))
                            V2summary <- rbind(V2summary, apply(V2.p, 2, quantile, prob = 0.025))
                            V2summary <- rbind(V2summary, apply(V2.p, 2, sd))
                            rownames(V2summary)[7:9] <- c("0.975", "0.025", "sd")
                            
                            V3summary <- as.matrix(apply(V3.p, 2, summary))
                            V3summary <- rbind(V3summary, apply(V3.p, 2, quantile, prob = 0.975))
                            V3summary <- rbind(V3summary, apply(V3.p, 2, quantile, prob = 0.025))
                            V3summary <- rbind(V3summary, apply(V3.p, 2, sd))
                            rownames(V3summary)[7:9] <- c("0.975", "0.025", "sd")
                            
                            if(storeV[1] == TRUE & !is.null(path))
                            {
                                save(V1.p, file = paste(path, "V1Pch", chain, ".RData", sep = ""))
                            }
                            if(storeV[2] == TRUE & !is.null(path))
                            {
                                save(V2.p, file = paste(path, "V2Pch", chain, ".RData", sep = ""))
                            }
                            if(storeV[3] == TRUE & !is.null(path))
                            {
                                save(V3.p, file = paste(path, "V3Pch", chain, ".RData", sep = ""))
                            }
                            
                            if(nGam_save > 0 & !is.null(path)){
                                save(gamma.p, file = paste(path, "/gammaPch", chain, ".Rdata", sep = ""))
                            }
                            
                            
                            if(p1 > 0){
                                covNames1 = colnames(Xmat1)
                            }
                            if(p1 == 0){
                                covNames1 = NULL
                            }
                            
                            if(p2 > 0){
                                covNames2 = colnames(Xmat2)
                            }
                            if(p2 == 0){
                                covNames2 = NULL
                            }
                            if(p3 > 0){
                                covNames3 = colnames(Xmat3)
                            }
                            if(p3 == 0){
                                covNames3 = NULL
                            }
                            
                            
                            
                            ### posterior predictive checks ###
                            
                            ## 1. log pseudo-marginal likelihood
                            
                            invLH.p <- matrix(mcmcRet$invLH, nrow = n, byrow = TRUE)
                            
                            ## 2. deviance information criterion
                            
                            gamma_mean <- matrix(mcmcRet$gammaP, nrow = n, byrow = TRUE)
                            dev.p        <- matrix(mcmcRet$dev, nrow = nStore, byrow = TRUE)

                            
                            ret[[nam]] <- list(beta1.p = beta1.p, beta2.p = beta2.p, beta3.p = beta3.p, lambda1.fin = lambda1.fin, lambda2.fin = lambda2.fin, lambda3.fin = lambda3.fin, mu_lam1.p = mu_lam1.p, mu_lam2.p = mu_lam2.p, mu_lam3.p = mu_lam3.p, sigSq_lam1.p = sigSq_lam1.p, sigSq_lam2.p = sigSq_lam2.p, sigSq_lam3.p = sigSq_lam3.p, K1.p = K1.p, K2.p = K2.p, K3.p = K3.p, s1.p = s1.p, s2.p = s2.p, s3.p = s3.p, theta.p = theta.p, class.p = c.p, mu.p = mu.p, Sigma.p = Sigma.p, tau.p = tau.p, accept.beta1 = accept.beta1, accept.beta2 = accept.beta2, accept.beta3 = accept.beta3, accept.BI1 = accept.BI1, accept.BI2 = accept.BI2, accept.BI3 = accept.BI3, accept.DI1 = accept.DI1, accept.DI2 = accept.DI2, accept.DI3 = accept.DI3, accept.theta = accept.theta, time_lambda1 = time_lambda1, time_lambda2 = time_lambda2, time_lambda3 = time_lambda3, accept.V = accept.V, covNames1 = covNames1, covNames2 = covNames2, covNames3 = covNames3, V1sum = V1summary, V2sum = V2summary, V3sum = V3summary, gamma_mean = gamma_mean, dev.p=dev.p, invLH.p=invLH.p)
                            
                            
                        }   ## end: if PEM-DPM-M
                        
                        
                        
                        # model = "semi-Markov"
                        
                        #######################################
                        ############ PEM-DPM-SM ############
                        #######################################
                        
                        if(model_h3 == "semi-Markov"){
                            
                            ###
                            
                            K_1 <- K1
                            K_2 <- K2
                            K_3 <- K3
                            
                            theta <- startV[(p1+p2+p3+1)]
                            gamma <- startV[(p1+p2+p3+2):(p1+p2+p3+n+1)]
                            
                            if(p1+p2+p3 > 0)
                            {
                                startV <- as.vector(c(startV[1:(p1+p2+p3)],
                                K1, K2, K3,
                                mu_lam1, mu_lam2, mu_lam3,
                                sigSq_lam1, sigSq_lam2, sigSq_lam3,
                                theta, gamma,
                                lambda1, lambda2, lambda3, s1, s2, s3,
                                startV[(p1+p2+p3+n+1+1):length(startV)]))
                            }else
                            {
                                startV <- as.vector(c(K1, K2, K3,
                                mu_lam1, mu_lam2, mu_lam3,
                                sigSq_lam1, sigSq_lam2, sigSq_lam3,
                                theta, gamma,
                                lambda1, lambda2, lambda3, s1, s2, s3,
                                startV[(p1+p2+p3+n+1+1):length(startV)]))
                            }
                            
                            startValues1 <- startV[1:(p1+p2+p3+10+n+2*(K_1+1)+2*(K_2+1)+2*(K_3+1))]
                            startValues2 <- startV[(p1+p2+p3+10+n+2*(K_1+1)+2*(K_2+1)+2*(K_3+1)+1):length(startV)]
                            startV <- c(startValues1, 1, 1, startValues2)
                            
                            mcmcParams1 <- mcmcParams[1:(18+num_s_propBI1+num_s_propBI2+num_s_propBI3+nTime_lambda1+nTime_lambda2+nTime_lambda3+1)]
                            mcmcParams2 <- mcmcParams[(18+num_s_propBI1+num_s_propBI2+num_s_propBI3+nTime_lambda1+nTime_lambda2+nTime_lambda3+2):(length(mcmcParams)-1)]
                            
                            mcmcParams <- c(mcmcParams1, 1, mcmcParams2, n)
                            
                            mcmcRet <- .C("BpeDpCorScrSMmcmc",
                            survData         = as.double(as.matrix(survData)),
                            n                = as.integer(n),
                            p1                = as.integer(p1),
                            p2                = as.integer(p2),
                            p3                = as.integer(p3),
                            J                = as.integer(J),
                            nj                = as.double(nj),
                            hyperParams     = as.double(hyperParams),
                            startValues     = as.double(startV),
                            mcmcParams        = as.double(mcmcParams),
                            numReps            = as.integer(numReps),
                            thin            = as.integer(thin),
                            burninPerc      = as.double(burninPerc),
                            nGam_save        = as.integer(nGam_save),
                            samples_beta1     = as.double(rep(0, nStore*p1)),
                            samples_beta2     = as.double(rep(0, nStore*p2)),
                            samples_beta3     = as.double(rep(0, nStore*p3)),
                            samples_mu_lam1     = as.double(rep(0, nStore*1)),
                            samples_mu_lam2     = as.double(rep(0, nStore*1)),
                            samples_mu_lam3     = as.double(rep(0, nStore*1)),
                            samples_sigSq_lam1    = as.double(rep(0, nStore*1)),
                            samples_sigSq_lam2    = as.double(rep(0, nStore*1)),
                            samples_sigSq_lam3    = as.double(rep(0, nStore*1)),
                            samples_K1          = as.double(rep(0, nStore*1)),
                            samples_K2          = as.double(rep(0, nStore*1)),
                            samples_K3          = as.double(rep(0, nStore*1)),
                            samples_s1          = as.double(rep(0, nStore*(K1_max + 1))),
                            samples_s2          = as.double(rep(0, nStore*(K2_max + 1))),
                            samples_s3          = as.double(rep(0, nStore*(K3_max + 1))),
                            samples_nu2     = as.double(rep(0, nStore*1)),
                            samples_nu3     = as.double(rep(0, nStore*1)),
                            samples_theta     = as.double(rep(0, nStore*1)),
                            samples_V1        = as.double(rep(0, nStore*J)),
                            samples_V2        = as.double(rep(0, nStore*J)),
                            samples_V3        = as.double(rep(0, nStore*J)),
                            samples_c        = as.double(rep(0, nStore*J)),
                            samples_mu        = as.double(rep(0, nStore*3*J)),
                            samples_Sigma    = as.double(rep(0, nStore*3*3*J)),
                            samples_tau     = as.double(rep(0, nStore*1)),
                            samples_gamma     = as.double(rep(0, nStore*nGam_save)),
                            samples_gamma_last = as.double(rep(0, n)),
                            samples_misc    = as.double(rep(0, (p1+p2+p3+9+n+1+n+J))),
                            lambda1_fin            = as.double(rep(0, nStore*nTime_lambda1)),
                            lambda2_fin            = as.double(rep(0, nStore*nTime_lambda2)),
                            lambda3_fin            = as.double(rep(0, nStore*nTime_lambda3)),
                            gammaP            = as.double(rep(0, n)),
                            dev                 = as.double(rep(0, nStore*1)),
                            invLH = as.double(rep(0, n)),
                            logLH_fin = as.double(0),
                            lpml                 = as.double(rep(0, nStore*1)),
                            moveVec             = as.double(rep(0, numReps)))
                            
                            
                            if(p1 > 0){
                                beta1.p         <- matrix(mcmcRet$samples_beta1, nrow = nStore, byrow = TRUE)
                            }
                            if(p1 == 0){
                                beta1.p         <- NULL
                            }
                            if(p2 > 0){
                                beta2.p         <- matrix(mcmcRet$samples_beta2, nrow = nStore, byrow = TRUE)
                            }
                            if(p2 == 0){
                                beta2.p         <- NULL
                            }
                            if(p3 > 0){
                                beta3.p         <- matrix(mcmcRet$samples_beta3, nrow = nStore, byrow = TRUE)
                            }
                            if(p3 == 0){
                                beta3.p         <- NULL
                            }
                            
                            lambda1.fin     <- matrix(mcmcRet$lambda1_fin, nrow = nStore, byrow = TRUE)
                            lambda2.fin     <- matrix(mcmcRet$lambda2_fin, nrow = nStore, byrow = TRUE)
                            lambda3.fin     <- matrix(mcmcRet$lambda3_fin, nrow = nStore, byrow = TRUE)
                            
                            mu_lam1.p         <- matrix(mcmcRet$samples_mu_lam1, nrow = nStore, byrow = TRUE)
                            mu_lam2.p         <- matrix(mcmcRet$samples_mu_lam2, nrow = nStore, byrow = TRUE)
                            mu_lam3.p         <- matrix(mcmcRet$samples_mu_lam3, nrow = nStore, byrow = TRUE)
                            sigSq_lam1.p     <- matrix(mcmcRet$samples_sigSq_lam1, nrow = nStore, byrow = TRUE)
                            sigSq_lam2.p     <- matrix(mcmcRet$samples_sigSq_lam2, nrow = nStore, byrow = TRUE)
                            sigSq_lam3.p     <- matrix(mcmcRet$samples_sigSq_lam3, nrow = nStore, byrow = TRUE)
                            
                            K1.p             <- matrix(mcmcRet$samples_K1, nrow = nStore, byrow = TRUE)
                            K2.p             <- matrix(mcmcRet$samples_K2, nrow = nStore, byrow = TRUE)
                            K3.p             <- matrix(mcmcRet$samples_K3, nrow = nStore, byrow = TRUE)
                            s1.p             <- matrix(mcmcRet$samples_s1, nrow = nStore, byrow = TRUE)
                            s2.p             <- matrix(mcmcRet$samples_s2, nrow = nStore, byrow = TRUE)
                            s3.p             <- matrix(mcmcRet$samples_s3, nrow = nStore, byrow = TRUE)
                            
                            theta.p         <- matrix(mcmcRet$samples_theta, nrow = nStore, byrow = TRUE)
                            gamma.p         <- matrix(mcmcRet$samples_gamma, nrow = nStore, byrow = TRUE)
                            
                            V1.p            <- matrix(mcmcRet$samples_V1, nrow = nStore, byrow = TRUE)
                            V2.p            <- matrix(mcmcRet$samples_V2, nrow = nStore, byrow = TRUE)
                            V3.p            <- matrix(mcmcRet$samples_V3, nrow = nStore, byrow = TRUE)
                            c.p            <- matrix(mcmcRet$samples_c, nrow = nStore, byrow = TRUE)
                            mu.p            <- matrix(mcmcRet$samples_mu, nrow = nStore, byrow = TRUE)
                            Sigma.p         <- array(as.vector(mcmcRet$samples_Sigma), c(3, 3 * J, nStore))
                            tau.p           <- matrix(mcmcRet$samples_tau, nrow = nStore, byrow = TRUE)
                            
                            if(p1 > 0){
                                accept.beta1     <- as.vector(mcmcRet$samples_misc[1:(p1)])/sum(as.vector(mcmcRet$moveVec)==1)*p1
                            }
                            if(p1 == 0){
                                accept.beta1     <- NULL
                            }
                            if(p2 > 0){
                                accept.beta2     <- as.vector(mcmcRet$samples_misc[(p1+1):(p1+p2)])/sum(as.vector(mcmcRet$moveVec)==2)*p2
                            }
                            if(p2 == 0){
                                accept.beta2     <- NULL
                            }
                            if(p3 > 0){
                                accept.beta3     <- as.vector(mcmcRet$samples_misc[(p1+p2+1):(p1+p2+p3)])/sum(as.vector(mcmcRet$moveVec)==3)*p3
                            }
                            if(p3 == 0){
                                accept.beta3     <- NULL
                            }
                            
                            accept.BI1        <- as.vector(mcmcRet$samples_misc[(p1+p2+p3)+1])/sum(as.vector(mcmcRet$moveVec)==12)
                            accept.DI1        <- as.vector(mcmcRet$samples_misc[(p1+p2+p3)+2])/sum(as.vector(mcmcRet$moveVec)==15)
                            accept.BI2        <- as.vector(mcmcRet$samples_misc[(p1+p2+p3)+3])/sum(as.vector(mcmcRet$moveVec)==13)
                            accept.DI2        <- as.vector(mcmcRet$samples_misc[(p1+p2+p3)+4])/sum(as.vector(mcmcRet$moveVec)==16)
                            accept.BI3        <- as.vector(mcmcRet$samples_misc[(p1+p2+p3)+5])/sum(as.vector(mcmcRet$moveVec)==14)
                            accept.DI3        <- as.vector(mcmcRet$samples_misc[(p1+p2+p3)+6])/sum(as.vector(mcmcRet$moveVec)==17)
                            accept.theta    <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+7)])/sum(as.vector(mcmcRet$moveVec)==11)
                            accept.V       <- as.vector(mcmcRet$samples_misc[(p1+p2+p3+10+n+n+1):(p1+p2+p3+10+n+n+J)])/sum(as.vector(mcmcRet$moveVec)==18)/3
                            
                            
                            V1summary <- as.matrix(apply(V1.p, 2, summary))
                            V1summary <- rbind(V1summary, apply(V1.p, 2, quantile, prob = 0.975))
                            V1summary <- rbind(V1summary, apply(V1.p, 2, quantile, prob = 0.025))
                            V1summary <- rbind(V1summary, apply(V1.p, 2, sd))
                            rownames(V1summary)[7:9] <- c("0.975", "0.025", "sd")
                            
                            V2summary <- as.matrix(apply(V2.p, 2, summary))
                            V2summary <- rbind(V2summary, apply(V2.p, 2, quantile, prob = 0.975))
                            V2summary <- rbind(V2summary, apply(V2.p, 2, quantile, prob = 0.025))
                            V2summary <- rbind(V2summary, apply(V2.p, 2, sd))
                            rownames(V2summary)[7:9] <- c("0.975", "0.025", "sd")
                            
                            V3summary <- as.matrix(apply(V3.p, 2, summary))
                            V3summary <- rbind(V3summary, apply(V3.p, 2, quantile, prob = 0.975))
                            V3summary <- rbind(V3summary, apply(V3.p, 2, quantile, prob = 0.025))
                            V3summary <- rbind(V3summary, apply(V3.p, 2, sd))
                            rownames(V3summary)[7:9] <- c("0.975", "0.025", "sd")
                            
                            if(storeV[1] == TRUE & !is.null(path))
                            {
                                save(V1.p, file = paste(path, "V1Pch", chain, ".RData", sep = ""))
                            }
                            if(storeV[2] == TRUE & !is.null(path))
                            {
                                save(V2.p, file = paste(path, "V2Pch", chain, ".RData", sep = ""))
                            }
                            if(storeV[3] == TRUE & !is.null(path))
                            {
                                save(V3.p, file = paste(path, "V3Pch", chain, ".RData", sep = ""))
                            }
                            
                            if(nGam_save > 0 & !is.null(path)){
                                save(gamma.p, file = paste(path, "/gammaPch", chain, ".Rdata", sep = ""))
                            }
                            
                            if(p1 > 0){
                                covNames1 = colnames(Xmat1)
                            }
                            if(p1 == 0){
                                covNames1 = NULL
                            }
                            
                            if(p2 > 0){
                                covNames2 = colnames(Xmat2)
                            }
                            if(p2 == 0){
                                covNames2 = NULL
                            }
                            if(p3 > 0){
                                covNames3 = colnames(Xmat3)
                            }
                            if(p3 == 0){
                                covNames3 = NULL
                            }
                            
                            
                            
                            
                            ### posterior predictive checks ###
                            
                            ## 1. log pseudo-marginal likelihood
                            
                            invLH.p <- matrix(mcmcRet$invLH, nrow = n, byrow = TRUE)
                            
                            ## 2. deviance information criterion
                            
                            gamma_mean <- matrix(mcmcRet$gammaP, nrow = n, byrow = TRUE)
                            dev.p        <- matrix(mcmcRet$dev, nrow = nStore, byrow = TRUE)
                            
                            ret[[nam]] <- list(beta1.p = beta1.p, beta2.p = beta2.p, beta3.p = beta3.p, lambda1.fin = lambda1.fin, lambda2.fin = lambda2.fin, lambda3.fin = lambda3.fin, mu_lam1.p = mu_lam1.p, mu_lam2.p = mu_lam2.p, mu_lam3.p = mu_lam3.p, sigSq_lam1.p = sigSq_lam1.p, sigSq_lam2.p = sigSq_lam2.p, sigSq_lam3.p = sigSq_lam3.p, K1.p = K1.p, K2.p = K2.p, K3.p = K3.p, s1.p = s1.p, s2.p = s2.p, s3.p = s3.p, theta.p = theta.p, class.p = c.p, mu.p = mu.p, Sigma.p = Sigma.p, tau.p = tau.p, accept.beta1 = accept.beta1, accept.beta2 = accept.beta2, accept.beta3 = accept.beta3, accept.BI1 = accept.BI1, accept.BI2 = accept.BI2, accept.BI3 = accept.BI3, accept.DI1 = accept.DI1, accept.DI2 = accept.DI2, accept.DI3 = accept.DI3, accept.theta = accept.theta, time_lambda1 = time_lambda1, time_lambda2 = time_lambda2, time_lambda3 = time_lambda3, accept.V = accept.V, covNames1 = covNames1, covNames2 = covNames2, covNames3 = covNames3, V1sum = V1summary, V2sum = V2summary, V3sum = V3summary, gamma_mean = gamma_mean, dev.p=dev.p, invLH.p=invLH.p)
                            
                            
                        }   ## end: if PEM-DPM-SM
                        
                        ## DIC and LPML
                        
                        be1.p <- ret$chain1$beta1.p
                        be2.p <- ret$chain1$beta2.p
                        be3.p <- ret$chain1$beta3.p
                        thet.p <- ret$chain1$theta.p
                        lam1.p <- ret$chain1$lambda1.fin
                        lam2.p <- ret$chain1$lambda2.fin
                        lam3.p <- ret$chain1$lambda3.fin
                        V1me <- ret$chain1$V1sum[4,]
                        V2me <- ret$chain1$V2sum[4,]
                        V3me <- ret$chain1$V3sum[4,]
                        Dev.p <- ret$chain1$dev.p
                        InvLH.p <- ret$chain1$invLH.p
                        
                        if(nChain > 1){
                            for(i in 2:nChain){
                                nam <- paste("chain", i, sep="")
                                be1.p <- rbind(be1.p, ret[[nam]]$beta1.p)
                                be2.p <- rbind(be2.p, ret[[nam]]$beta2.p)
                                be3.p <- rbind(be3.p, ret[[nam]]$beta3.p)
                                thet.p <- rbind(thet.p, ret[[nam]]$theta.p)
                                lam1.p <- rbind(lam1.p, ret[[nam]]$lambda1.fin)
                                lam2.p <- rbind(lam2.p, ret[[nam]]$lambda2.fin)
                                lam3.p <- rbind(lam3.p, ret[[nam]]$lambda3.fin)
                                V1me <- rbind(V1me, ret[[nam]]$V1sum[4,])
                                V2me <- rbind(V2me, ret[[nam]]$V2sum[4,])
                                V3me <- rbind(V3me, ret[[nam]]$V3sum[4,])
                                Dev.p <- rbind(Dev.p, ret[[nam]]$dev.p)
                                InvLH.p <- cbind(InvLH.p, ret[[nam]]$invLH.p)
                            }
                        }
                        if(model_h3 == "Markov")
                        {
                            logLH_fin <- .C("BpeDpCorScr_logMLH_DIC",
                            survData         = as.double(as.matrix(survData)),
                            n                = as.integer(n),
                            p1                = as.integer(p1),
                            p2                = as.integer(p2),
                            p3                = as.integer(p3),
                            J                = as.integer(J),
                            be1 = as.double(apply(be1.p, 2, mean)),
                            be2 = as.double(apply(be2.p, 2, mean)),
                            be3 = as.double(apply(be3.p, 2, mean)),
                            thet = as.double(apply(theta.p, 2, mean)),
                            lam1 = as.double(apply(lam1.p, 2, mean)),
                            lam2 = as.double(apply(lam2.p, 2, mean)),
                            lam3 = as.double(apply(lam3.p, 2, mean)),
                            s_1 = as.double(time_lambda1),
                            s_2 = as.double(time_lambda2),
                            s_3 = as.double(time_lambda3),
                            V_1 = as.double(apply(V1me, 2, mean)),
                            V_2 = as.double(apply(V2me, 2, mean)),
                            V_3 = as.double(apply(V3me, 2, mean)),
                            K_1                = as.integer(nTime_lambda1),
                            K_2                = as.integer(nTime_lambda2),
                            K_3                = as.integer(nTime_lambda3),
                            val = as.double(0))$val
                        }else if(model_h3 == "semi-Markov")
                        {
                            logLH_fin <- .C("BpeDpCorScrSM_logMLH_DIC",
                            survData         = as.double(as.matrix(survData)),
                            n                = as.integer(n),
                            p1                = as.integer(p1),
                            p2                = as.integer(p2),
                            p3                = as.integer(p3),
                            J                = as.integer(J),
                            be1 = as.double(apply(be1.p, 2, mean)),
                            be2 = as.double(apply(be2.p, 2, mean)),
                            be3 = as.double(apply(be3.p, 2, mean)),
                            thet = as.double(apply(theta.p, 2, mean)),
                            lam1 = as.double(apply(lam1.p, 2, mean)),
                            lam2 = as.double(apply(lam2.p, 2, mean)),
                            lam3 = as.double(apply(lam3.p, 2, mean)),
                            s_1 = as.double(time_lambda1),
                            s_2 = as.double(time_lambda2),
                            s_3 = as.double(time_lambda3),
                            V_1 = as.double(apply(V1me, 2, mean)),
                            V_2 = as.double(apply(V2me, 2, mean)),
                            V_3 = as.double(apply(V3me, 2, mean)),
                            K_1                = as.integer(nTime_lambda1),
                            K_2                = as.integer(nTime_lambda2),
                            K_3                = as.integer(nTime_lambda3),
                            val = as.double(0))$val
                        }
                        
                        DIC = 2*mean(Dev.p) + 2 * logLH_fin
                        LPML = sum(log(1/apply(InvLH.p, 1, mean)))
                        
                        
                    }   ## end: if PEM-DPM
                    
                    
                }   ## end: if PEM
                
                
                chain = chain + 1
                
                
            }    ## end: while(chain <= nChain)
            
            
            ret[["setup"]]    <- list(Formula=Formula, nCov = nCov, hyperParams = hyperParams, startValues = startValues, mcmc = mcmc, nGam_save = nGam_save, numReps = numReps, thin = thin, path = path, burninPerc = burninPerc, hz.type = hz.type, re.type = re.type, model = model_h3, nChain = nChain, survData=survData, p1=p1, p2=p2, p3=p3, J=J)
            ret[["DIC"]] <- DIC
            ret[["LPML"]] <- LPML
            
            if(hz.type == "Weibull")
            {
                if(re.type == "MVN")
                {
                    ret$class <- c("Bayes_HReg", "ID", "Cor", "WB", "MVN")
                }
                if(re.type == "DPM")
                {
                    ret$class <- c("Bayes_HReg", "ID", "Cor", "WB", "DPM")
                }
            }
            if(hz.type == "PEM")
            {
                if(re.type == "MVN")
                {
                    ret$class <- c("Bayes_HReg", "ID", "Cor", "PEM", "MVN")
                }
                if(re.type == "DPM")
                {
                    ret$class <- c("Bayes_HReg", "ID", "Cor", "PEM", "DPM")
                }
            }
            
            class(ret) <- "Bayes_HReg"
            
            return(ret)
            
        }
        
        
    }
    else{
        warning(" (numReps * burninPerc) must be divisible by (thin)")
    }
    
    
    
}# end of function "BayesID"

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.