Nothing
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"
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.