From ce9da3fe593b2a70e599a01faa220b9fee2401d4 Mon Sep 17 00:00:00 2001 From: Waldir Leoncio Date: Thu, 15 Aug 2024 13:36:39 +0200 Subject: [PATCH] Translated `func_MCMC_graph()` to C++ (#16) * Added `src/.gitignore` * Added placeholder for `func_MCMC_graph_cpp()` (#11) * Added `cpp` flag (#11) * Updated docs * Added test file for `func_MCMC_graph_cpp()` (#11) * Translated a bit more of #11 * Added FIXME (#11) * Updated RoxygenNote version * Properly exporting `func_MCMC_graph_cpp()` (#11) * Fixed syntax (#11) * Updated unit test (#11) * Translated a bit more of #11 * Added commented R code to loop (#11) * Translated rest of `func_MCMC_graph()' (#11) * Added unit tests for `func_MCMC_graph()` (#11) * Put `func_MCMC()` progressbar under `verbose` (#14) Otherwise, the `pb` would be drawn even if the user sets `verbose = FALSE` on `BayesSurvive()` (which I assume is not the intended behaviour). * Updated NEWS.md (#11) --- DESCRIPTION | 2 +- NEWS.md | 1 + R/BayesSurvive.R | 7 +- R/RcppExports.R | 4 + R/func_MCMC.R | 15 +- R/func_MCMC_graph.R | 8 +- inst/doc/BayesCox.R | 4 - inst/doc/BayesCox.Rmd | 165 -------- inst/doc/BayesCox.html | 489 ---------------------- man/BayesSurvive.Rd | 5 +- man/func_MCMC.Rd | 5 +- man/func_MCMC_graph.Rd | 4 +- src/RcppExports.cpp | 16 + src/func_MCMC_graph_cpp.cpp | 183 ++++++++ src/init.c | 4 +- tests/testthat/test-func_MCMC_graph_cpp.R | 48 +++ 16 files changed, 287 insertions(+), 673 deletions(-) delete mode 100644 inst/doc/BayesCox.R delete mode 100755 inst/doc/BayesCox.Rmd delete mode 100644 inst/doc/BayesCox.html create mode 100644 src/func_MCMC_graph_cpp.cpp create mode 100644 tests/testthat/test-func_MCMC_graph_cpp.R diff --git a/DESCRIPTION b/DESCRIPTION index d0a4099..3e63c01 100755 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -15,7 +15,7 @@ License: GPL-3 VignetteBuilder: knitr Depends: R (>= 4.0) Encoding: UTF-8 -RoxygenNote: 7.3.1 +RoxygenNote: 7.3.2 LinkingTo: Rcpp, RcppArmadillo, testthat Imports: Rcpp, ggplot2, GGally, mvtnorm, survival, riskRegression, utils, stats diff --git a/NEWS.md b/NEWS.md index a080079..fd003f0 100644 --- a/NEWS.md +++ b/NEWS.md @@ -2,6 +2,7 @@ * Add units tests * Rename the output of MPM coefficients in function `coef.BayesSurvive()` +* Added `cpp` argument to `BayesSurvive()` to allow for faster computation using `Rcpp` # BayesSurvive 0.0.3 diff --git a/R/BayesSurvive.R b/R/BayesSurvive.R index 5476742..2a9edc1 100755 --- a/R/BayesSurvive.R +++ b/R/BayesSurvive.R @@ -39,6 +39,7 @@ #' output for parameters 'G', 'V', 'C' and 'Sig' in the graphical model #' if \code{MRF.G = FALSE} #' @param verbose logical value to display the progess of MCMC +#' @param cpp logical, whether to use C++ code for faster computation #' #' #' @return An object of class \code{BayesSurvive} is saved as @@ -108,7 +109,8 @@ BayesSurvive <- function(survObj, burnin = 0, thin = 1, output_graph_para = FALSE, - verbose = TRUE) { + verbose = TRUE, + cpp = FALSE) { # same number of covariates p in all subgroups p <- ifelse(is.list(survObj[[1]]), NCOL(survObj[[1]]$X), NCOL(survObj$X)) Beta.ini <- numeric(p) @@ -306,7 +308,8 @@ BayesSurvive <- function(survObj, MRF_2b = MRF2b, MRF_G = MRF.G, output_graph_para, - verbose + verbose, + cpp ) if (S == 1 && MRF.G) { diff --git a/R/RcppExports.R b/R/RcppExports.R index 8c70180..cbb6c63 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -1,6 +1,10 @@ # Generated by using Rcpp::compileAttributes() -> do not edit by hand # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 +func_MCMC_graph_cpp <- function(sobj, hyperpar, ini, S, method, MRF_2b) { + .Call(`_BayesSurvive_func_MCMC_graph_cpp`, sobj, hyperpar, ini, S, method, MRF_2b) +} + settingInterval_cpp <- function(y, delta_, s_, J_) { .Call(`_BayesSurvive_settingInterval_cpp`, y, delta_, s_, J_) } diff --git a/R/func_MCMC.R b/R/func_MCMC.R index 7a104bd..20e8933 100755 --- a/R/func_MCMC.R +++ b/R/func_MCMC.R @@ -27,6 +27,7 @@ #' output for parameters 'G', 'V', 'C' and 'Sig' in the graphical model #' if \code{MRF_G = FALSE} #' @param verbose logical value to display the progess of MCMC +#' @inheritParams BayesSurvive #' #' @return A list object saving the MCMC results with components including #' 'gamma.p', 'beta.p', 'h.p', 'gamma.margin', 'beta.margin', 's', 'eta0', @@ -38,7 +39,7 @@ func_MCMC <- function(survObj, hyperpar, initial, nIter, thin, burnin, S, method, MRF_2b, MRF_G, - output_graph_para, verbose) { + output_graph_para, verbose, cpp = FALSE) { # prior parameters for grouped data likelihood of Cox model if (method == "Pooled" && MRF_G) { # method = "Pooled" hyperpar$s <- sort(survObj$t[survObj$di == 1]) @@ -155,15 +156,17 @@ func_MCMC <- function(survObj, hyperpar, initial, # MCMC sampling # Initializes the progress bar - if (verbose) cat(" Running MCMC iterations ...\n") - pb <- txtProgressBar(min = 0, max = nIter, style = 3, width = 50, char = "=") + if (verbose) { + cat(" Running MCMC iterations ...\n") + pb <- txtProgressBar(min = 0, max = nIter, style = 3, width = 50, char = "=") + } for (M in 1:nIter) { # if (method %in% c("CoxBVSSL", "Sub-struct") || # (method == "Pooled" && !MRF_G)) { if (!MRF_G) { # update graph and precision matrix - network <- func_MCMC_graph(survObj, hyperpar, ini, S, method, MRF_2b) + network <- func_MCMC_graph(survObj, hyperpar, ini, S, method, MRF_2b, cpp) Sig.ini <- ini$Sig.ini <- network$Sig.ini # precision matrix? C.ini <- ini$C.ini <- network$C.ini @@ -283,9 +286,9 @@ func_MCMC <- function(survObj, hyperpar, initial, # } # Sets the progress bar to the current state - setTxtProgressBar(pb, M) + if (verbose) setTxtProgressBar(pb, M) } # the end of MCMC sampling - close(pb) # Close the connection of progress bar + if (verbose) close(pb) # Close the connection of progress bar if (S == 1 && MRF_G) { mcmcOutcome$gamma.margin <- mcmcOutcome$gamma.margin / (nIter - burnin) diff --git a/R/func_MCMC_graph.R b/R/func_MCMC_graph.R index ad345e5..7f3d045 100755 --- a/R/func_MCMC_graph.R +++ b/R/func_MCMC_graph.R @@ -15,6 +15,7 @@ #' @param method a method option from #' \code{c("Pooled", "CoxBVSSL", "Sub-struct")} #' @param MRF_2b two different b in MRF prior for subgraphs G_ss and G_rs +#' @inheritParams func_MCMC #' #' @return A list object with components "Sig" the updated covariance matrices, #' "G.ini" the updated graph, "V.ini" the updated variances for precision @@ -23,7 +24,11 @@ #' #' #' @export -func_MCMC_graph <- function(sobj, hyperpar, ini, S, method, MRF_2b) { +func_MCMC_graph <- function(sobj, hyperpar, ini, S, method, MRF_2b, cpp = FALSE) { + if (cpp) { + warning("This is not yet fully implemented. Please use cpp = FALSE for production") + return(func_MCMC_graph_cpp(sobj, hyperpar, ini, S, method, MRF_2b)) + } n <- sobj$n p <- sobj$p SSig <- sobj$SSig @@ -182,6 +187,5 @@ func_MCMC_graph <- function(sobj, hyperpar, ini, S, method, MRF_2b) { } } } - return(list(Sig.ini = Sig, G.ini = G, V.ini = V, C.ini = C)) } diff --git a/inst/doc/BayesCox.R b/inst/doc/BayesCox.R deleted file mode 100644 index 0402d91..0000000 --- a/inst/doc/BayesCox.R +++ /dev/null @@ -1,4 +0,0 @@ -## ----setup, include=FALSE----------------------------------------------------- -knitr::opts_chunk$set(echo = TRUE, eval = FALSE) -options(rmarkdown.html_vignette.check_title = FALSE) - diff --git a/inst/doc/BayesCox.Rmd b/inst/doc/BayesCox.Rmd deleted file mode 100755 index 0644c69..0000000 --- a/inst/doc/BayesCox.Rmd +++ /dev/null @@ -1,165 +0,0 @@ ---- -title: "Bayesian Cox Models with graph-structure priors" -output: rmarkdown::html_vignette -vignette: > - %\VignetteEngine{knitr::rmarkdown} - %\VignetteIndexEntry{Bayesian Cox Models with graph-structure priors} - \usepackage[utf8]{inputenc} ---- - -```{r setup, include=FALSE} -knitr::opts_chunk$set(echo = TRUE, eval = FALSE) -options(rmarkdown.html_vignette.check_title = FALSE) -``` - -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](https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-021-04483-z)) 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](https://CRAN.R-project.org/package=BayesSurvive) - -```r -install.packages("BayesSurvive") -``` - -Install the latest development version from [GitHub](https://github.com/ocbe-uio/BayesSurvive) - -```r -#install.packages("remotes") -remotes::install_github("ocbe-uio/BayesSurvive") -``` - - -## Examples - -### Simulate data - -```r -library("BayesSurvive") -# Load the example dataset -data("simData", package = "BayesSurvive") -dataset = list("X" = simData[[1]]$X, - "t" = simData[[1]]$time, - "di" = simData[[1]]$status) -``` - -### Run a Bayesian Cox model - -```r -## 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$. - -```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 -## -## 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 - -```r -hyperparPooled <- append(hyperparPooled, list("lambda" = 3, "nu0" = 0.05, "nu1" = 5)) -fit2 <- BayesSurvive(survObj = list(dataset), model.type = "Pooled", MRF.G = FALSE, - hyperpar = hyperparPooled, initial = initial, nIter = 10) -``` - -### Run a Bayesian Cox model with subgroups using fixed graph - -```r -# 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) -``` - -### Run a Bayesian Cox model with subgroups using graphical learning - -```r -fit4 <- BayesSurvive(survObj = dataset2, - hyperpar = hyperparPooled, initial = initial, - model.type="CoxBVSSL", MRF.G = FALSE, - nIter = 3, burnin = 0) -``` - -## 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](https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-021-04483-z). diff --git a/inst/doc/BayesCox.html b/inst/doc/BayesCox.html deleted file mode 100644 index d2f5bba..0000000 --- a/inst/doc/BayesCox.html +++ /dev/null @@ -1,489 +0,0 @@ - - - - - - - - - - - - - - -Bayesian Cox Models with graph-structure priors - - - - - - - - - - - - - - - - - - - - - - - - - - -

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.packages("BayesSurvive")
-

Install the latest development version from GitHub

-
#install.packages("remotes")
-remotes::install_github("ocbe-uio/BayesSurvive")
-
-
-

Examples

-
-

Simulate data

-
library("BayesSurvive")
-# Load the example dataset
-data("simData", package = "BayesSurvive")
-dataset = list("X" = simData[[1]]$X, 
-               "t" = simData[[1]]$time,
-               "di" = simData[[1]]$status)
-
-
-

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\).

-
which( VS(fit, method = "FDR", threshold = 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.

-
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).

-
predict(fit, survObj.new = dataset, times = 8.5)
-
##               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.

-
predict(fit, survObj.new = dataset, type = c("cumhazard", "survival"))
-
#        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

-
hyperparPooled <- append(hyperparPooled, list("lambda" = 3, "nu0" = 0.05, "nu1" = 5))
-fit2 <- BayesSurvive(survObj = list(dataset), model.type = "Pooled", MRF.G = FALSE,
-                     hyperpar = hyperparPooled, initial = initial, nIter = 10)
-
-
-

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)
-
-
-

Run a Bayesian Cox model with subgroups using graphical -learning

-
fit4 <- BayesSurvive(survObj = dataset2, 
-                     hyperpar = hyperparPooled, initial = initial, 
-                     model.type="CoxBVSSL", MRF.G = FALSE, 
-                     nIter = 3, burnin = 0)
-
-
-
-

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.

-
-
- - - - - - - - - - - diff --git a/man/BayesSurvive.Rd b/man/BayesSurvive.Rd index f31e9db..38122c1 100644 --- a/man/BayesSurvive.Rd +++ b/man/BayesSurvive.Rd @@ -16,7 +16,8 @@ BayesSurvive( burnin = 0, thin = 1, output_graph_para = FALSE, - verbose = TRUE + verbose = TRUE, + cpp = FALSE ) } \arguments{ @@ -60,6 +61,8 @@ output for parameters 'G', 'V', 'C' and 'Sig' in the graphical model if \code{MRF.G = FALSE}} \item{verbose}{logical value to display the progess of MCMC} + +\item{cpp}{logical, whether to use C++ code for faster computation} } \value{ An object of class \code{BayesSurvive} is saved as diff --git a/man/func_MCMC.Rd b/man/func_MCMC.Rd index acf3153..b47a411 100644 --- a/man/func_MCMC.Rd +++ b/man/func_MCMC.Rd @@ -16,7 +16,8 @@ func_MCMC( MRF_2b, MRF_G, output_graph_para, - verbose + verbose, + cpp = FALSE ) } \arguments{ @@ -50,6 +51,8 @@ output for parameters 'G', 'V', 'C' and 'Sig' in the graphical model if \code{MRF_G = FALSE}} \item{verbose}{logical value to display the progess of MCMC} + +\item{cpp}{logical, whether to use C++ code for faster computation} } \value{ A list object saving the MCMC results with components including diff --git a/man/func_MCMC_graph.Rd b/man/func_MCMC_graph.Rd index 6910971..e610f3f 100644 --- a/man/func_MCMC_graph.Rd +++ b/man/func_MCMC_graph.Rd @@ -4,7 +4,7 @@ \alias{func_MCMC_graph} \title{Function to learn MRF graph} \usage{ -func_MCMC_graph(sobj, hyperpar, ini, S, method, MRF_2b) +func_MCMC_graph(sobj, hyperpar, ini, S, method, MRF_2b, cpp = FALSE) } \arguments{ \item{sobj}{a list containing observed data from \code{n} subjects; @@ -20,6 +20,8 @@ func_MCMC_graph(sobj, hyperpar, ini, S, method, MRF_2b) \code{c("Pooled", "CoxBVSSL", "Sub-struct")}} \item{MRF_2b}{two different b in MRF prior for subgraphs G_ss and G_rs} + +\item{cpp}{logical, whether to use C++ code for faster computation} } \value{ A list object with components "Sig" the updated covariance matrices, diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index f13b2ea..a93b284 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -11,6 +11,22 @@ Rcpp::Rostream& Rcpp::Rcout = Rcpp::Rcpp_cout_get(); Rcpp::Rostream& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get(); #endif +// func_MCMC_graph_cpp +Rcpp::List func_MCMC_graph_cpp(const Rcpp::List sobj, const Rcpp::List hyperpar, const Rcpp::List ini, const uint S, const std::string method, const bool MRF_2b); +RcppExport SEXP _BayesSurvive_func_MCMC_graph_cpp(SEXP sobjSEXP, SEXP hyperparSEXP, SEXP iniSEXP, SEXP SSEXP, SEXP methodSEXP, SEXP MRF_2bSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< const Rcpp::List >::type sobj(sobjSEXP); + Rcpp::traits::input_parameter< const Rcpp::List >::type hyperpar(hyperparSEXP); + Rcpp::traits::input_parameter< const Rcpp::List >::type ini(iniSEXP); + Rcpp::traits::input_parameter< const uint >::type S(SSEXP); + Rcpp::traits::input_parameter< const std::string >::type method(methodSEXP); + Rcpp::traits::input_parameter< const bool >::type MRF_2b(MRF_2bSEXP); + rcpp_result_gen = Rcpp::wrap(func_MCMC_graph_cpp(sobj, hyperpar, ini, S, method, MRF_2b)); + return rcpp_result_gen; +END_RCPP +} // settingInterval_cpp Rcpp::List settingInterval_cpp(const arma::vec y, const arma::vec delta_, const arma::vec s_, const unsigned int J_); RcppExport SEXP _BayesSurvive_settingInterval_cpp(SEXP ySEXP, SEXP delta_SEXP, SEXP s_SEXP, SEXP J_SEXP) { diff --git a/src/func_MCMC_graph_cpp.cpp b/src/func_MCMC_graph_cpp.cpp new file mode 100644 index 0000000..c80e527 --- /dev/null +++ b/src/func_MCMC_graph_cpp.cpp @@ -0,0 +1,183 @@ +#include +#include +// [[Rcpp::depends(RcppArmadillo)]] + +arma::mat construct_G_MRF(arma::mat G, arma::vec b, uint S, int p, bool MRF_2b) { + arma::mat G_MRF(G.n_rows, G.n_cols); + if (MRF_2b) { + // TODO: develop test for this case. Off-by-one minefield! + // two different values for b in MRF prior for subgraphs G_ss and G_rs + arma::uvec p_seq = arma::regspace(0, p - 1); + for (uint g = 0; g < S; g++) { + // b1 * G_ss + arma::uvec g_seq(p, arma::fill::value(g * p)); + arma::uvec idx = g_seq + p_seq; + G_MRF.submat(idx, idx) *= b(0); + } + for (uint g = 0; g < S - 1; g++) { + // b2 * G_rs + for (uint r = g; r < S - 1; r++) { + arma::uvec g_seq(p, arma::fill::value(g * p)); + arma::uvec r_seq(p, arma::fill::value(r * p)); + arma::uvec idx = g_seq + p_seq; + arma::uvec idx2 = r_seq + p_seq; + G_MRF.submat(idx, idx2) *= b(1); + G_MRF.submat(idx2, idx) *= b(1); + } + } + } else { + // one value for b in MRF prior for all subgraphs + G_MRF = b[0] * G; + } + return G_MRF; +} + +// [[Rcpp::export]] +Rcpp::List func_MCMC_graph_cpp( + const Rcpp::List sobj, + const Rcpp::List hyperpar, + const Rcpp::List ini, + const uint S, + const std::string method, + const bool MRF_2b +) { + // Extracting data + Rcpp::List n = sobj["n"]; + uint p = Rcpp::as(sobj["p"]); + Rcpp::List SSig = sobj["SSig"]; + + double pii = Rcpp::as(hyperpar["pi.G"]); + double a = Rcpp::as(hyperpar["a"]); + arma::vec b(2); + double lambda = Rcpp::as(hyperpar["lambda"]); + arma::mat V0 = Rcpp::as(hyperpar["V0"]); + arma::mat V1 = Rcpp::as(hyperpar["V1"]); + + arma::mat G = Rcpp::as(ini["G.ini"]); + Rcpp::List V = Rcpp::as(ini["V.ini"]); + Rcpp::List Sig = Rcpp::as(ini["Sig.ini"]); + Rcpp::List C = Rcpp::as(ini["C.ini"]); + + Rcpp::List gamma_ini_list = Rcpp::as(ini["gamma.ini"]); + arma::vec gamma_ini = Rcpp::as(gamma_ini_list[0]); + + if (MRF_2b) { + // two different values for b in MRF prior for subgraphs G_ss and G_rs + b = Rcpp::as(hyperpar["b"]); + } else { + b = {hyperpar["b"], hyperpar["b"]}; + } + arma::mat G_MRF = construct_G_MRF(G, b, S, p, MRF_2b); + + // Update of precision matrix and graph within each subgroup + // (analogous to SSSL algorithm for concentration (precision) graph models in + // function 'BayesGGM_SSVS_FixedV0V1' (Wang, 2015)) + for (uint g = 0; g < S; g++) { // loop through subgroups + arma::mat V_g = V[g]; + arma::mat C_g = C[g]; + + arma::uvec p_seq = arma::regspace(0, p - 1); + arma::uvec g_seq(p, arma::fill::value(g * p)); + arma::uvec idx = g_seq + p_seq; + arma::mat G_g = G.submat(idx, idx); + uint n_g = Rcpp::as(n[g]); + arma::mat S_g = SSig[g]; + arma::mat Sig_g = Sig[g]; + + // TODO: code i loop through genes + for (arma::uword i = 0; i < p; i++) { + arma::uvec ind_noi = arma::regspace(0, p - 1); + ind_noi.shed_row(i); + arma::vec v_temp = V_g.submat(ind_noi, arma::uvec({i})); + + arma::mat Sig11 = Sig_g.submat(ind_noi, ind_noi); + arma::vec Sig12 = Sig_g.submat(ind_noi, arma::uvec({i})); + + // Omega_11^(-1) + arma::mat invC11 = Sig11 - Sig12 * Sig12.t() / Sig_g(i, i); + + // C^(-1) + arma::mat Ci = (S_g(i, i) + lambda) * invC11 + arma::diagmat(1 / v_temp); + + Ci = 0.5 * (Ci + Ci.t()); + arma::mat Ci_chol = arma::chol(Ci); + + arma::mat mu_i = -arma::solve(Ci_chol, arma::solve(Ci_chol.t(), S_g.submat(ind_noi, arma::uvec({i})))); + arma::mat beta = mu_i + arma::solve(Ci_chol, arma::randn(p - 1)); + + // Update of last column in Omega_gg + C_g.submat(ind_noi, arma::uvec({i})) = beta; + C_g.submat(arma::uvec({i}), ind_noi) = beta.t(); + + double a_gam = 0.5 * n_g + 1; + double b_gam = (S_g(i, i) + lambda) * 0.5; + double gam = arma::as_scalar(arma::randg(1, arma::distr_param(a_gam, 1 / b_gam))); + + double c = arma::as_scalar(beta.t() * invC11 * beta); + C_g(i, i) = gam + c; + + // Below updating covariance matrix according to one-column change of precision matrix + // invC11beta <- invC11 %*% beta + + arma::mat invC11beta = invC11 * beta; + + Sig_g.submat(ind_noi, ind_noi) = invC11 + invC11beta * invC11beta.t() / gam; + Sig_g.submat(ind_noi, arma::uvec({i})) = -invC11beta / gam; + Sig_g.submat(arma::uvec({i}), ind_noi) = arma::trans(-invC11beta / gam); + Sig_g(i, i) = 1 / gam; + + // Update of variance matrix V_g and subgraph G_g + + if (i < p) { + arma::vec beta2(p); + beta2.elem(ind_noi) = beta; + + for (arma::uword j = i + 1; j < p; j++) { + // G where g_ss,ij=1 or 0: + arma::mat G_prop1 = G_MRF; + arma::mat G_prop0 = G_MRF; + arma::uword gpj = g * p + j; + arma::uword gpi = g * p + i; + G_prop1(gpj, gpi) = b(0); + G_prop1(gpi, gpj) = b(0); + G_prop0(gpj, gpi) = 0; + G_prop0(gpi, gpj) = 0; + + double v0 = V0(j, i); + double v1 = V1(j, i); + + double w1 = -0.5 * log(v1) - 0.5 * std::pow(beta2(j), 2) / v1 + log(pii); + double w2 = -0.5 * log(v0) - 0.5 * std::pow(beta2(j), 2) / v0 + log(1 - pii); + + double wa = arma::as_scalar(w1 + (a * arma::sum(gamma_ini) + gamma_ini.t() * G_prop1 * gamma_ini)); + double wb = arma::as_scalar(w2 + (a * arma::sum(gamma_ini) + gamma_ini.t() * G_prop0 * gamma_ini)); + + double w_max = std::max(wa, wb); + + double w = exp(wa - w_max) / (exp(wa - w_max) + exp(wb - w_max)); + + bool z = Rcpp::runif(1)[0] < w; + double v = z ? v1 : v0; + V_g(j, i) = v; + V_g(i, j) = v; + + G_g(j, i) = static_cast(z); + G_g(i, j) = static_cast(z); + } + } + } + V[g] = V_g; + C[g] = C_g; + G.submat(idx, idx) = G_g; + Sig[g] = Sig_g; + } + + // Assembling output + Rcpp::List out = Rcpp::List::create( + Rcpp::Named("Sig.ini") = Sig, + Rcpp::Named("G.ini") = G, + Rcpp::Named("V.ini") = V, + Rcpp::Named("C.ini") = C + ); + return out; +} diff --git a/src/init.c b/src/init.c index ac8aef4..d03c0ac 100644 --- a/src/init.c +++ b/src/init.c @@ -3,7 +3,7 @@ #include // for NULL #include -/* FIXME: +/* FIXME: Check these declarations against the C/Fortran source code. */ @@ -15,6 +15,7 @@ extern SEXP _BayesSurvive_sumMatProdVec(SEXP, SEXP); extern SEXP _BayesSurvive_updateBH_cpp(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); extern SEXP _BayesSurvive_updateBH_list_cpp(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); extern SEXP _BayesSurvive_updateRP_genomic_cpp(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); +extern SEXP _BayesSurvive_func_MCMC_graph_cpp(SEXP, SEXP, SEXP, SEXP, SEXP, SEXP); static const R_CallMethodDef CallEntries[] = { {"_BayesSurvive_calJpost_helper_cpp", (DL_FUNC) &_BayesSurvive_calJpost_helper_cpp, 9}, @@ -24,6 +25,7 @@ static const R_CallMethodDef CallEntries[] = { {"_BayesSurvive_updateBH_cpp", (DL_FUNC) &_BayesSurvive_updateBH_cpp, 7}, {"_BayesSurvive_updateBH_list_cpp", (DL_FUNC) &_BayesSurvive_updateBH_list_cpp, 7}, {"_BayesSurvive_updateRP_genomic_cpp", (DL_FUNC) &_BayesSurvive_updateRP_genomic_cpp, 12}, + {"_BayesSurvive_func_MCMC_graph_cpp", (DL_FUNC) &_BayesSurvive_func_MCMC_graph_cpp, 6}, {NULL, NULL, 0} }; diff --git a/tests/testthat/test-func_MCMC_graph_cpp.R b/tests/testthat/test-func_MCMC_graph_cpp.R new file mode 100644 index 0000000..509f82a --- /dev/null +++ b/tests/testthat/test-func_MCMC_graph_cpp.R @@ -0,0 +1,48 @@ +# Load the example dataset +dataset <- list( + "X" = simData[[1]]$X, + "t" = simData[[1]]$time, + "di" = simData[[1]]$status +) + +# 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 + "lambda" = 3, + "nu0" = 0.05, + "nu1" = 5 +) + +# Run a 'Pooled' Bayesian Cox model with graphical learning + +set.seed(3346141) +BayesSurvive_wrap <- function(use_cpp = FALSE) { + suppressWarnings( + BayesSurvive( + survObj = list(dataset), model.type = "Pooled", MRF.G = FALSE, + hyperpar = hyperparPooled, initial = initial, nIter = 3, + verbose = FALSE, cpp = use_cpp + ) + ) +} +fit_R <- BayesSurvive_wrap(use_cpp = FALSE) +fit_C <- BayesSurvive_wrap(use_cpp = TRUE) + +test_that("R and C++ objects are similar", { + expect_equal(fit_R$call, fit_C$call) + expect_equal(fit_R$input, fit_C$input) + for (obj in names(fit_R$output)[2]) { + expect_equal(fit_R$output[[obj]], fit_C$output[[obj]], tolerance = 1) + } +})