Description Usage Arguments Value References Examples
MCMC routine for the strucatural equation model
1 | mcmc(ct, u1, u2, X, nburnin = 1000, nmc = 2000, nthin = 1)
|
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) |
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 |
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.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 | 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
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.