# mcmc function MCMC routine for the strucatural equation model ```r mcmc(ct, u1, u2, X, nburnin = 1000, nmc = 2000, nthin = 1) ``` ## Arguments - `ct`: survival response, a $n*2$ matrix with first column as response and second column as right censored indicator, 1 is event time and 0 is right censored. - `u1`: Matix of predictors from the first platform, dimension $n*q_1$ - `u2`: Matix of predictors from the first platform, dimension $n*q_2$ - `X`: Matrix of covariates, dimension $n*p$. - `nburnin`: number of burnin samples - `nmc`: number of markov chain samples - `nthin`: thinning parameter. Default is 1 (no thinning) ## Returns - **pMean.beta.t**: - **pMean.beta.t**: - **pMean.alpha.t**: - **pMean.alpha.t**: - **pMean.phi.t**: - **pMean.phi.t**: - **pMean.alpha.u1**: - **pMean.alpha.u2**: - **pMean.alpha.u2**: - **pMean.phi.u1**: - **pMean.eta1**: - **pMean.eta2**: - **pMean.sigma.t.square**: - **pMean.sigma.u1.square**: - **pMean.sigma.u2.square**: - **alpha.t.samples**: - **phi.t.samples**: - **beta1.t.samples**: - **beta2.t.samples**: - **beta.t.samples**: - **alpha.u1.samples**: - **alpha.u2.samples**: - **phi.u1.samples**: - **phi.u2.samples**: - **eta1.samples**: - **eta2.samples**: - **sigma.t.square.samples**: - **sigma.u1.square.samples**: - **sigma.u2.square.samples**: - **pMean.logt.hat**: - **DIC**: - **WAIC**: ## Examples ```r require(MASS) # for random number from multivariate normal distribution n <- 100 # number of individuals p <- 5 # number of variables q1 <- 20 # dimension of the response q2 <- 20 # dimension of the response ngrid <- 1000 nburnin <- 100 nmc <- 200 nthin <- 5 niter <- nburnin + nmc effsamp <- (niter - nburnin)/nthin alpha.tt <- runif(n = 1, min = -1, max = 1) # intercept term alpha.u1t <- runif(n = 1, min = -1, max = 1) # intercept term alpha.u2t <- runif(n = 1, min = -1, max = 1) # intercept term beta.tt <- runif(n = p, min = -1, max = 1) # regression parameter gamma1.t <- runif(n = q1, min = -1, max = 1) gamma2.t <- runif(n = q2, min = -1, max = 1) phi.tt <- 1 phi.u1t <- 1 phi.u2t <- 1 sigma.tt <- 1 sigma.u1t <- 1 sigma.u2t <- 1 sigma.etat1 <- 1 sigma.etat2 <- 1 x <- mvrnorm(n = n, mu = rep(0, p), Sigma = diag(p)) eta2 <- rnorm(n = 1, mean = 0, sd = sigma.etat2) eta1 <- rnorm(n = 1, mean = eta2, sd = sigma.etat1) logt <- rnorm(n = n, mean = alpha.tt + x %*% beta.tt + eta1 * phi.tt, sd = sigma.tt) u1 <- matrix(rnorm(n = n * q1, mean = alpha.u1t + eta1 * phi.u1t, sd = sigma.u1t), nrow = n, ncol = q1) u2 <- matrix(rnorm(n = n * q2, mean = alpha.u2t + eta2 * phi.u2t, sd = sigma.u2t), nrow = n, ncol = q2) logt <- rnorm(n = n, mean = alpha.tt + x %*% beta.tt + u1 %*% gamma1.t + u2 %*% gamma2.t, sd = sigma.tt) # Survival time generation T <- exp(logt) # AFT model C <- rgamma(n, shape = 1, rate = 1) # 50% censor time <- pmin(T, C) # observed time is min of censored and true status = time == T # set to 1 if event is observed 1 - sum(status)/length(T) # censoring rate censor.rate <- 1 - sum(status)/length(T) # censoring rate censor.rate summary(C) summary(T) ct <- as.matrix(cbind(time = time, status = status)) # censored time logt.grid <- seq(from = min(logt) - 1, to = max(logt) + 1, length.out = ngrid) index1 <- which(ct[, 2] == 1) # which are NOT censored ct1 <- ct[index1, ] posterior.fit.sem <- mcmc(ct, u1, u2, x, nburnin = nburnin, nmc = nmc, nthin = nthin) pMean.beta.t <- posterior.fit.sem$pMean.beta.t pMean.alpha.t <- posterior.fit.sem$pMean.alpha.t pMean.phi.t <- posterior.fit.sem$pMean.phi.t pMean.alpha.u1 <- posterior.fit.sem$pMean.alpha.u1 pMean.alpha.u2 <- posterior.fit.sem$pMean.alpha.u2 pMean.phi.u1 <- posterior.fit.sem$pMean.phi.u1 pMean.phi.u2 <- posterior.fit.sem$pMean.phi.u2 pMean.eta1 <- posterior.fit.sem$pMean.eta1 pMean.eta2 <- posterior.fit.sem$pMean.eta2 pMean.logt.hat <- posterior.fit.sem$posterior.fit.sem pMean.sigma.t.square <- posterior.fit.sem$pMean.sigma.t.square pMean.sigma.u1.square <- posterior.fit.sem$pMean.sigma.u1.square pMean.sigma.u2.square <- posterior.fit.sem$pMean.sigma.u2.square pMean.logt.hat <- posterior.fit.sem$pMean.logt.hat DIC.sem <- posterior.fit.sem$DIC WAIC.sem <- posterior.fit.sem$WAIC mse.sem <- mean(pMean.logt.hat[index1] - log(ct1[, 1]))^2 alpha.t.samples <- posterior.fit.sem$alpha.t.samples beta1.t.samples <- posterior.fit.sem$beta1.t.samples beta2.t.samples <- posterior.fit.sem$beta2.t.samples beta.t.samples <- posterior.fit.sem$beta.t.samples phi.t.samples <- posterior.fit.sem$phi.t.samples alpha.u1.samples <- posterior.fit.sem$alpha.u1.samples alpha.u2.samples <- posterior.fit.sem$alpha.u2.samples phi.u1.samples <- posterior.fit.sem$phi.u1.samples phi.u2.samples <- posterior.fit.sem$phi.u2.samples sigma.t.square.samples <- posterior.fit.sem$sigma.t.square.samples sigma.u1.square.samples <- posterior.fit.sem$sigma.u1.square.samples sigma.u2.square.samples <- posterior.fit.sem$sigma.u2.square.samples eta1.samples <- posterior.fit.sem$eta1.samples eta2.samples <- posterior.fit.sem$eta2.samples inv.cpo <- matrix(0, nrow = effsamp, ncol = n) # this will store inverse cpo values log.cpo <- rep(0, n) # this will store log cpo for(iter in 1:effsamp) # Post burn in { inv.cpo[iter, ] <- 1/(dnorm(ct[, 1], mean = alpha.t.samples[iter] + x %*% beta.t.samples[, iter] + + eta1.samples[iter] * phi.t.samples[iter], sd = sqrt(sigma.t.square.samples[iter]))^ct[, 2] * pnorm(ct[, 1], mean = alpha.t.samples[iter] + x %*% beta.t.samples[, iter] + + eta1.samples[iter] * phi.t.samples[iter], sd = sqrt(sigma.t.square.samples[iter]), lower.tail = FALSE)^(1 - ct[, 2])) } # End of iter loop for (i in 1:n){ log.cpo[i] <- -log(mean(inv.cpo[, i])) # You average invcpo[i] over the iterations, # then take 1/average and then take log. # Hence the negative sign in the log } lpml.sem <- sum(log.cpo) DIC.sem WAIC.sem mse.sem ``` ## References Maity, A. K., Lee, S. C., Mallick, B. K., & Sarkar, T. R. (2020). Bayesian structural equation modeling in multiple omics data integration with application to circadian genes. Bioinformatics, 36(13), 3951-3958.