From 3e5e1a1ca65ff133e3b7e9491074af7857971345 Mon Sep 17 00:00:00 2001 From: Waldir Leoncio Date: Fri, 21 Jun 2024 12:41:41 +0200 Subject: [PATCH] Added unit tests from README (#15) * Added `src/.gitignore` * Added test skeleton Literally the output of `usethis::use_testthat()` * Added README as test script (#13) * Restyled file Output of `styler::style_file()`. * Added unit tests (solves #13) * Added gcno files to .gitignore * Added `inst/doc` to .gitignore Those are generated files that get removed whenever `devtools::check()` is ran. * Finished writing tests for #13 * Sped up #13 tests 3x * Resolved TODO, suppressed output on test (#13) * Reinstated test with `MRF.G = FALSE` (#13) This is necessary for `func_MCMC_graph()` to be called and work on #11 to continue. * Added test for `MRF.G = FALSE` (closes #13) * Relaxed unit test expectations So checks passes on macOS --- .gitignore | 4 +- DESCRIPTION | 5 +- inst/doc/BayesCox.R | 4 - inst/doc/BayesCox.Rmd | 165 ------------ inst/doc/BayesCox.html | 489 ---------------------------------- src/.gitignore | 3 + tests/testthat.R | 12 + tests/testthat/test-general.R | 98 +++++++ 8 files changed, 120 insertions(+), 660 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/.gitignore create mode 100644 tests/testthat.R create mode 100644 tests/testthat/test-general.R diff --git a/.gitignore b/.gitignore index dd52d9d..c0cf732 100644 --- a/.gitignore +++ b/.gitignore @@ -12,6 +12,8 @@ BayesSurv/src/symbols.rds BayesSurv/src/Makevars BayesSurv/BayesSurv.Rcheck/* BayesSurv.Rcheck/* +*.gcno +inst/doc BayesSurv/autom4te.cache BayesSurv/build @@ -19,4 +21,4 @@ BayesSurv/hyperpar.xml *.tar.gz *.log -*.status \ No newline at end of file +*.status diff --git a/DESCRIPTION b/DESCRIPTION index 27e091f..e75e802 100755 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -18,7 +18,9 @@ RoxygenNote: 7.3.1 LinkingTo: Rcpp, RcppArmadillo (>= 0.9.000) Imports: Rcpp, ggplot2, GGally, mvtnorm, survival, riskRegression, utils, stats -Suggests: knitr +Suggests: + knitr, + testthat (>= 3.0.0) LazyData: true NeedsCompilation: yes Packaged: 2024-06-05 07:00:29 UTC; zhiz @@ -28,3 +30,4 @@ Author: Zhi Zhao [aut, cre], Manuela Zucknick [ctb], Jörg Rahnenführer [ctb] Maintainer: Zhi Zhao +Config/testthat/edition: 3 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 1dd7402..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)) -``` - - - -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 143bfa0..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))
-

-

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/src/.gitignore b/src/.gitignore new file mode 100644 index 0000000..f60d2a7 --- /dev/null +++ b/src/.gitignore @@ -0,0 +1,3 @@ +*.o +*.so +Makevars diff --git a/tests/testthat.R b/tests/testthat.R new file mode 100644 index 0000000..5410aa3 --- /dev/null +++ b/tests/testthat.R @@ -0,0 +1,12 @@ +# This file is part of the standard setup for testthat. +# It is recommended that you do not modify it. +# +# Where should you do additional test configuration? +# Learn more about the roles of various files in: +# * https://r-pkgs.org/testing-design.html#sec-tests-files-overview +# * https://testthat.r-lib.org/articles/special-files.html + +library(testthat) +library(BayesSurvive) + +test_check("BayesSurvive") diff --git a/tests/testthat/test-general.R b/tests/testthat/test-general.R new file mode 100644 index 0000000..90bd95c --- /dev/null +++ b/tests/testthat/test-general.R @@ -0,0 +1,98 @@ +# 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 +) + +## 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 = 50, burnin = 0, + verbose = FALSE +) +pred_1 <- predict(fit, survObj.new = dataset, times = 8.5, verbose = FALSE) +pred_2 <- predict(fit, survObj.new = dataset, type = c("cumhazard", "survival")) + +test_that("fit has properly class and length", { + expect_s3_class(fit, "BayesSurvive") + expect_length(fit, 3L) + expect_length(fit$input, 9L) + expect_length(fit$input$hyperpar, 10L) + expect_length(fit$output, 17L) + expect_length(fit$output$survObj, 4L) +}) + +test_that("fit has expected values", { + tol <- 1e-3 + with(fit$output, { + expect_equal(eta0, c("(Intercept)" = 1.74e-5), tolerance = tol) + expect_equal(head(s, 4), c(3.2969, 3.3217, 4.0938, 4.4107), tolerance = tol) + expect_equal(head(survObj$t, 4), c(8.53, 4.09, 8.82, 6.09), tolerance = tol) + }) + expect_equal(which(VS(fit, method = "FDR", threshold = 0.6)), c(6, 12)) +}) + +test_that("predictions have expected values", { + tol <- 1e-1 + expect_equal( + head(pred_1$times), + c(0.00000000, 0.08585859, 0.17171717, 0.25757576, 0.34343434, 0.42929293), + tolerance = tol + ) + expect_equal( + head(pred_2$cumhazard[, 1], 5L), + c(1.782968e-04, 1.544469e-01, 9.193403e-07, 6.192432e-03, 4.701002e-04), + tolerance = tol + ) + expect_equal( + head(pred_2$survival[, 1]), + c(0.9998217, 0.8568890, 0.9999991, 0.9938267, 0.9995300, 0.9935771), + tolerance = tol + ) + expect_equal( + head(pred_2$times), + c(3.296921, 3.321787, 3.899565, 4.093849, 4.410782, 4.932070), + tolerance = tol + ) + expect_false(pred_2$se) + expect_false(pred_2$band) + expect_false(pred_2$diag) + expect_false(pred_2$baseline) +}) + +### Run a 'Pooled' Bayesian Cox model with graphical learning + +hyperparPooled <- append(hyperparPooled, list("lambda" = 3, "nu0" = 0.05, "nu1" = 5)) +set.seed(3346141) +fit2 <- BayesSurvive( + survObj = list(dataset), model.type = "Pooled", MRF.G = FALSE, + hyperpar = hyperparPooled, initial = initial, nIter = 3, verbose = FALSE +) +test_that("fit2 has expected values", { + tol <- 1e-3 + with(fit2$output, { + expect_equal(eta0[[1]], c("(Intercept)" = 1.74e-5), tolerance = tol) + expect_equal(head(s[[1]], 3), c(3.2969, 3.3217, 4.0938), tolerance = tol) + expect_equal(head(survObj[[1]]$t, 3), c(8.53, 4.09, 8.82), tolerance = tol) + }) + expect_equal(which(VS(fit2, method = "FDR", threshold = 0.8)), c(81, 182)) +})