-
-The function `BayesSurvive::VS()` can show the index of selected variables by controlling Bayesian false discovery rate (FDR) at the level $\alpha = 0.05$.
-
-```r
-which( VS(fit, method = "FDR", threshold = 0.05) )
-```
-```{ .text .no-copy }
-#[1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 194
-```
-
-
-### Plot time-dependent Brier scores
-
-The function `BayesSurvive::plotBrier()` can show the time-dependent Brier scores based on posterior mean of coefficients or Bayesian model averaging.
-
-```r
-plotBrier(fit, survObj.new = dataset)
-```
-
-
-
-We can also use the function `BayesSurvive::predict()` to obtain the Brier score at time 8.5, the integrated Brier score (IBS) from time 0 to 8.5 and the index of prediction accuracy (IPA).
-
-```r
-predict(fit, survObj.new = dataset, times = 8.5)
-```
-```{ .text .no-copy }
-## Brier(t=8.5) IBS(t:0~8.5) IPA(t=8.5)
-## Null.model 0.2290318 0.08185316 0.0000000
-## Bayesian.Cox 0.1013692 0.02823275 0.5574011
-```
-
-### Predict survival probabilities and cumulative hazards
-
-The function `BayesSurvive::predict()` can estimate the survival probabilities and cumulative hazards.
-
-```r
-predict(fit, survObj.new = dataset, type = c("cumhazard", "survival"))
-```
-```{ .text .no-copy }
-# observation times cumhazard survival
-## Bayesian Cox Models with graph-structure -priors
- - - -This is a R/Rcpp package BayesSurvive for Bayesian -survival models with graph-structured selection priors for sparse -identification of high-dimensional features predictive of survival (Madjar -et al., 2021) and its extensions with the use of a fixed graph via a -Markov Random Field (MRF) prior for capturing known structure of -high-dimensional features, e.g. disease-specific pathways from the Kyoto -Encyclopedia of Genes and Genomes (KEGG) database.
-
-
-Installation
-Install the latest released version from CRAN
- -Install the latest development version from GitHub
- -
-
-Examples
-
-
-Simulate data
- -
-
-
-Run a Bayesian Cox model
-## Initial value: null model without covariates
-initial = list("gamma.ini" = rep(0, ncol(dataset$X)))
-# Prior parameters
-hyperparPooled = list(
- "c0" = 2, # prior of baseline hazard
- "tau" = 0.0375, # sd (spike) for coefficient prior
- "cb" = 20, # sd (slab) for coefficient prior
- "pi.ga" = 0.02, # prior variable selection probability for standard Cox models
- "a" = -4, # hyperparameter in MRF prior
- "b" = 0.1, # hyperparameter in MRF prior
- "G" = simData$G # hyperparameter in MRF prior
-)
-
-## run Bayesian Cox with graph-structured priors
-set.seed(123)
-fit <- BayesSurvive(survObj = dataset, model.type = "Pooled", MRF.G = TRUE,
- hyperpar = hyperparPooled, initial = initial,
- nIter = 200, burnin = 100)
-
-## show posterior mean of coefficients and 95% credible intervals
-library("GGally")
-plot(fit) +
- coord_flip() +
- theme(axis.text.x = element_text(angle = 90, size = 7))
The function BayesSurvive::VS() can show the index of
-selected variables by controlling Bayesian false discovery rate (FDR) at
-the level \(\alpha = 0.05\).
#[1] 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 194
-
-Plot time-dependent Brier scores
-The function BayesSurvive::plotBrier() can show the
-time-dependent Brier scores based on posterior mean of coefficients or
-Bayesian model averaging.
We can also use the function BayesSurvive::predict() to
-obtain the Brier score at time 8.5, the integrated Brier score (IBS)
-from time 0 to 8.5 and the index of prediction accuracy (IPA).
## Brier(t=8.5) IBS(t:0~8.5) IPA(t=8.5)
-## Null.model 0.2290318 0.08185316 0.0000000
-## Bayesian.Cox 0.1013692 0.02823275 0.5574011
-
-Predict survival probabilities and cumulative hazards
-The function BayesSurvive::predict() can estimate the
-survival probabilities and cumulative hazards.
# observation times cumhazard survival
-## <int> <num> <num> <num>
-## 1: 1 3.3 7.41e-05 1.00e+00
-## 2: 2 3.3 2.51e-01 7.78e-01
-## 3: 3 3.3 9.97e-07 1.00e+00
-## 4: 4 3.3 1.84e-03 9.98e-01
-## 5: 5 3.3 3.15e-04 1.00e+00
-## ---
-## 9996: 96 9.5 7.15e+00 7.88e-04
-## 9997: 97 9.5 3.92e+02 7.59e-171
-## 9998: 98 9.5 2.81e+00 6.02e-02
-## 9999: 99 9.5 3.12e+00 4.42e-02
-## 10000: 100 9.5 1.97e+01 2.79e-09
-
-Run a ‘Pooled’ Bayesian Cox model with graphical learning
- -
-
-
-
-Run a Bayesian Cox model with subgroups using fixed graph
-# specify a fixed joint graph between two subgroups
-hyperparPooled$G <- Matrix::bdiag(simData$G, simData$G)
-dataset2 <- simData[1:2]
-dataset2 <- lapply(dataset2, setNames, c("X", "t", "di", "X.unsc", "trueB"))
-fit3 <- BayesSurvive(survObj = dataset2,
- hyperpar = hyperparPooled, initial = initial,
- model.type="CoxBVSSL", MRF.G = TRUE,
- nIter = 10, burnin = 5)
-
-
-
-
-
-
-
-
-
-
-
-
diff --git a/src/func_MCMC_graph_cpp.cpp b/src/func_MCMC_graph_cpp.cpp
index 6b07d60..7d58509 100644
--- a/src/func_MCMC_graph_cpp.cpp
+++ b/src/func_MCMC_graph_cpp.cpp
@@ -84,7 +84,75 @@ Rcpp::List func_MCMC_graph_cpp(
arma::mat Sig_g = Sig[g];
// TODO: code i loop through genes
+ for (arma::uword i = 0; i < p; i++) {
+ // ind_noi <- setdiff(1:p, i)
+ // v_temp <- V_g[ind_noi, i]
+ // Sig11 <- Sig_g[ind_noi, ind_noi]
+ // Sig12 <- Sig_g[ind_noi, i]
+
+ // invC11 <- Sig11 - Sig12 %*% t(Sig12) / Sig_g[i, i] # Omega_11^(-1)
+
+ // Ci <- (S_g[i, i] + lambda) * invC11 + diag(1 / v_temp) # C^(-1)
+
+ // Ci <- (Ci + t(Ci)) / 2
+ // Ci_chol <- chol(Ci)
+
+ // mu_i <- -solve(Ci_chol, solve(t(Ci_chol), S_g[ind_noi, i]))
+ // beta <- mu_i + solve(Ci_chol, rnorm(p - 1))
+
+ // # Update of last column in Omega_gg
+ // C_g[ind_noi, i] <- C_g[i, ind_noi] <- beta # omega_12
+
+ // a_gam <- 0.5 * n_g + 1
+ // b_gam <- (S_g[i, i] + lambda) * 0.5
+ // gam <- rgamma(1, shape = a_gam, scale = 1 / b_gam)
+
+ // c <- t(beta) %*% invC11 %*% beta
+ // C_g[i, i] <- gam + c # omega_22
+
+ // # Below updating covariance matrix according to one-column change of precision matrix
+ // invC11beta <- invC11 %*% beta
+
+ // Sig_g[ind_noi, ind_noi] <- invC11 + invC11beta %*% t(invC11beta) / gam
+ // Sig_g[ind_noi, i] <- Sig_g[i, ind_noi] <- -invC11beta / gam
+ // Sig_g[i, i] <- 1 / gam
+
+ // # Update of variance matrix V_g and subgraph G_g
+
+ // if (i < p) {
+ // beta2 <- numeric(p)
+ // beta2[ind_noi] <- beta
+
+ // for (j in (i + 1):p) {
+ // # G where g_ss,ij=1 or 0:
+ // G.prop1 <- G.prop0 <- G.MRF
+ // G.prop1[(g - 1) * p + j, (g - 1) * p + i] <-
+ // G.prop1[(g - 1) * p + i, (g - 1) * p + j] <- b[1] # 1
+ // G.prop0[(g - 1) * p + j, (g - 1) * p + i] <-
+ // G.prop0[(g - 1) * p + i, (g - 1) * p + j] <- 0
+
+ // v0 <- V0[j, i]
+ // v1 <- V1[j, i]
+
+ // w1 <- -0.5 * log(v1) - 0.5 * beta2[j]^2 / v1 + log(pii)
+ // w2 <- -0.5 * log(v0) - 0.5 * beta2[j]^2 / v0 + log(1 - pii)
+
+ // wa <- w1 + (a * sum(gamma.ini) + t(gamma.ini) %*% G.prop1 %*% gamma.ini)
+ // wb <- w2 + (a * sum(gamma.ini) + t(gamma.ini) %*% G.prop0 %*% gamma.ini)
+
+ // w_max <- max(wa, wb)
+
+ // w <- exp(wa - w_max) / (exp(wa - w_max) + exp(wb - w_max))
+
+ // z <- (runif(1) < as.numeric(w)) # acceptance/ rejection of proposal
+ // v <- ifelse(z, v1, v0)
+ // V_g[j, i] <- V_g[i, j] <- v
+
+ // G_g[j, i] <- G_g[i, j] <- as.numeric(z)
+ // }
+ // }
+ }
V[g] = V_g;
C[g] = C_g;
G.submat(idx, idx) = G_g;
References
---Katrin Madjar, Manuela Zucknick, Katja Ickstadt, Jörg Rahnenführer -(2021). Combining heterogeneous subgroups with graph‐structured variable -selection priors for Cox regression. BMC Bioinformatics, -22(1):586. DOI: 10.1186/s12859-021-04483-z.
-