fix bug in updateBH_cpp()

ocbe-uio · Mar 21, 2024 · 717092c · 717092c
1 parent 9fadb0e
commit 717092c
Show file tree

Hide file tree

Showing 10 changed files with 109 additions and 63 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: psbcSpeedUp
 Title: Penalized Semiparametric Bayesian Cox Models
-Version: 2.0.5
-Date: 2023-12-08
+Version: 2.0.6
+Date: 2024-03-21
 URL: https://github.com/ocbe-uio/psbcSpeedUp
 BugReports: https://github.com/ocbe-uio/psbcSpeedUp/issues
 Authors@R: c(person("Zhi", "Zhao", role=c("aut","cre"), email = "[email protected]"), 
@@ -23,7 +23,7 @@ Imports: Rcpp, xml2, ggplot2, GGally, MASS, survival, riskRegression,
 Suggests: knitr
 LazyData: true
 NeedsCompilation: yes
-Packaged: 2023-12-07 23:55:48 UTC; zhiz
+Packaged: 2024-03-21 17:28:53 UTC; zhiz
 Author: Zhi Zhao [aut, cre],
   Manuela Zucknick [aut],
   Maral Saadati [aut],

diff --git a/NEWS.md b/NEWS.md
@@ -6,6 +6,11 @@
 * add R functions for feature (stability) selection
 * extend C++ source code for implementing other shrinkage and group priors
 
+# psbcSpeedUp 2.0.6
+
+* Fixed a bug in function `updateBH_cpp()`
+* Other minor changes
+
 # psbcSpeedUp 2.0.5
 
 * Added R functions for survival predictions, e.g. time-dependent Brier scores, integrated Brier score

diff --git a/R/predict.psbcSpeedUp.R b/R/predict.psbcSpeedUp.R
@@ -19,8 +19,8 @@
 #' @param times time points at which to evaluate the risks. If \code{NULL}
 #' (default), the event/censoring times are used. If \code{type="brier"}, the
 #' largest one of the \code{times} is used
-#' @param type option to chose for predicting survival probabilities
-#' (\code{type="probability"}) and brier scores (\code{type="brier"})
+#' @param type option to chose for predicting survival probabilities (one of
+#' \code{c('hazard','cumhazard','survival')}) or brier scores (\code{type="brier"})
 #' @param method option to use the posterior mean (\code{"mean"}) of coefficients
 #' for prediction or Bayesian model averaging (\code{"BMA"}) for prediction
 #' @param \dots not used

diff --git a/README.md b/README.md
@@ -1,10 +1,7 @@
+# psbcSpeedUp
+
 [![CRAN](http://www.r-pkg.org/badges/version/psbcSpeedUp)](https://cran.r-project.org/package=psbcSpeedUp)
-[![r-universe](https://ocbe-uio.r-universe.dev/badges/psbcSpeedUp)](https://ocbe-uio.r-universe.dev/psbcSpeedUp)
-[![R-CMD-check](https://github.com/ocbe-uio/psbcSpeedUp/workflows/R-CMD-check/badge.svg)](https://github.com/ocbe-uio/psbcSpeedUp/actions)
-[![License](https://img.shields.io/badge/License-GPLv3-brightgreen.svg)](https://www.gnu.org/licenses/gpl-3.0)
-[![DOI](https://img.shields.io/badge/doi-10.1002%2Fbimj.201400160-brightgreen)](https://doi.org/10.1002/bimj.201400160)
 
-# psbcSpeedUp
 
 This is a C++ speed-up and extended version of the R-pakcage [psbcGroup](https://CRAN.R-project.org/package=psbcGroup).
 It implements the Bayesian Lasso Cox model ([Lee et al., 2011](https://doi.org/10.2202/1557-4679.1301)) and the Bayesian Lasso Cox with mandatory variables ([Zucknick et al., 2015](https://doi.org/10.1002/bimj.201400160)).
@@ -73,7 +70,7 @@ The function `psbcSpeedUp::plot()` can show the posterior mean and 95% credible
 plot(fitBayesCox)
 ```
 
-<img src="man/figures/README_plot_beta.png" width="70%" />
+<img src="man/figures/README_plot_beta.png" width="50%" />
 
 
 ### Plot time-dependent Brier scores
@@ -88,7 +85,7 @@ plotBrier(fitBayesCox, times = 80)
 IBS  0.2089742     0.109274
 ```
 
-<img src="man/figures/README_plot_brier.png" width="80%" />
+<img src="man/figures/README_plot_brier.png" width="50%" />
 
 
 ### Predict survival probabilities and cumulative hazards
@@ -113,7 +110,6 @@ predict(fitBayesCox, type = c("cumhazard", "survival"))
 ## 40000:         200 107.641  2.24e-01  7.99e-01
 ```
 
-
 ## References
 
 > Kyu Ha Lee, Sounak Chakraborty, Jianguo Sun (2011).

diff --git a/inst/doc/BayesCox.R b/inst/doc/BayesCox.R
@@ -1,5 +1,5 @@
 ## ----setup, include=FALSE-----------------------------------------------------
-knitr::opts_chunk$set(echo = TRUE, eval = TRUE)
+knitr::opts_chunk$set(echo = TRUE, eval = FALSE)
 options(rmarkdown.html_vignette.check_title = FALSE)
 
 ## ----install1, eval = FALSE---------------------------------------------------
@@ -10,28 +10,28 @@ options(rmarkdown.html_vignette.check_title = FALSE)
 #  remotes::install_github("ocbe-uio/psbcSpeedUp")
 
 ## ----results='hide', warning=FALSE--------------------------------------------
-# Load the example dataset
-data("exampleData", package = "psbcSpeedUp")
-p <- exampleData$p
-q <- exampleData$q
-survObj <- exampleData[1:3]
-
-# Set hyperparameters (see help file for specifying more hyperparameters)
-mypriorPara <- list('eta0'=0.02, 'kappa0'=1, 'c0'=2, 'r'=10/9, 'delta'=1e-05,
-'lambdaSq'=1, 'sigmaSq'= runif(1, 0.1, 10), 'beta.prop.var'=1, 'beta.clin.var'=1)
-
-# run Bayesian Lasso Cox
-library("psbcSpeedUp")
-set.seed(123)
-fitBayesCox <- psbcSpeedUp(survObj, p=p, q=q, hyperpar=mypriorPara,
-nIter=1000, burnin=500, outFilePath="/tmp")
+#  # Load the example dataset
+#  data("exampleData", package = "psbcSpeedUp")
+#  p <- exampleData$p
+#  q <- exampleData$q
+#  survObj <- exampleData[1:3]
+#  
+#  # Set hyperparameters (see help file for specifying more hyperparameters)
+#  mypriorPara <- list('eta0'=0.02, 'kappa0'=1, 'c0'=2, 'r'=10/9, 'delta'=1e-05,
+#  'lambdaSq'=1, 'sigmaSq'= runif(1, 0.1, 10), 'beta.prop.var'=1, 'beta.clin.var'=1)
+#  
+#  # run Bayesian Lasso Cox
+#  library("psbcSpeedUp")
+#  set.seed(123)
+#  fitBayesCox <- psbcSpeedUp(survObj, p=p, q=q, hyperpar=mypriorPara,
+#  nIter=1000, burnin=500, outFilePath="/tmp")
 
 ## ----fig.width=5, fig.height=8------------------------------------------------
-plot(fitBayesCox)
+#  plot(fitBayesCox)
 
 ## ----fig.width=6, fig.heigh=5-------------------------------------------------
-plotBrier(fitBayesCox, times = 80)
+#  plotBrier(fitBayesCox, times = 80)
 
 ## -----------------------------------------------------------------------------
-predict(fitBayesCox, type = c("cumhazard", "survival"))
+#  predict(fitBayesCox, type = c("cumhazard", "survival"))
 
diff --git a/inst/doc/BayesCox.Rmd b/inst/doc/BayesCox.Rmd
@@ -8,7 +8,7 @@ vignette: >
 ---
 
 ```{r setup, include=FALSE}
-knitr::opts_chunk$set(echo = TRUE, eval = TRUE)
+knitr::opts_chunk$set(echo = TRUE, eval = FALSE)
 options(rmarkdown.html_vignette.check_title = FALSE)
 ```
 
@@ -65,6 +65,7 @@ fitBayesCox <- psbcSpeedUp(survObj, p=p, q=q, hyperpar=mypriorPara,
 nIter=1000, burnin=500, outFilePath="/tmp")
 ```
 
+
 ### Plot posterior estimates of regression cofficients
 
 The function `psbcSpeedUp::plot()` can show the posterior mean and 95% credible intervals of regression coefficients.
@@ -73,13 +74,23 @@ The function `psbcSpeedUp::plot()` can show the posterior mean and 95% credible
 plot(fitBayesCox)
 ```
 
+<img src="../man/figures/README_plot_beta.png" width="70%" />
+
+
 ### Plot time-dependent Brier scores
 
 The function `psbcSpeedUp::plotBrier()` can show the time-dependent Brier scores based on posterior mean of coefficients or Bayesian model averaging.
 
 ```{r, fig.width=6, fig.heigh=5}
 plotBrier(fitBayesCox, times = 80)
 ```
+```
+##     Null.model Bayesian.Cox
+## IBS  0.2089742   0.08195375
+```
+
+<img src="../man/figures/README_plot_brier.png" width="90%" />
+
 
 ### Predict survival probabilities and cumulative hazards
 
@@ -88,7 +99,20 @@ The function `psbcSpeedUp::predict()` can estimate the survival probabilities an
 ```{r}
 predict(fitBayesCox, type = c("cumhazard", "survival"))
 ```
-
+```
+##        observation   times cumhazard  survival
+##     1:           1   0.264  1.14e-05  1.00e+00
+##     2:           2   0.264  4.66e-05  1.00e+00
+##     3:           3   0.264  6.07e-05  1.00e+00
+##     4:           4   0.264  1.71e-05  1.00e+00
+##     5:           5   0.264  7.55e-05  1.00e+00
+##    ---                                        
+## 39996:         196 107.641  2.76e+00  6.36e-02
+## 39997:         197 107.641  5.65e-01  5.68e-01
+## 39998:         198 107.641  5.37e+01  4.86e-24
+## 39999:         199 107.641  4.25e+02 2.25e-185
+## 40000:         200 107.641  1.89e-01  8.28e-01
+```
 
 ## References
 

diff --git a/inst/doc/BayesCox.html b/inst/doc/BayesCox.html
diff --git a/man/predict.psbcSpeedUp.Rd b/man/predict.psbcSpeedUp.Rd
diff --git a/src/psbc.cpp b/src/psbc.cpp
@@ -88,7 +88,7 @@ arma::vec PSBC::updateBH_cpp(arma::mat &ind_r_d_, arma::vec hPriorSh_, arma::vec
     for( unsigned int i=0; i<J_; ++i )
      //h_rate(i) += arma::accu( ind_r_d_.col(i) % arma::exp( xbeta_ ) );
         h_rate(i) += arma::dot( ind_r_d_.col(i), xbeta_ );*/
-    arma::vec h_rate = c0_ + sumMatProdVec(ind_r_d_, exp_xbeta);
+    arma::vec h_rate = c0_ + sumMatProdVec(ind_r_d_, arma::exp( xbeta_ ));
 
     // arma::vec shape = hPriorSh_ + d_;
     arma::vec h_ = arma::zeros<arma::vec>(J_);

diff --git a/vignettes/BayesCox.Rmd b/vignettes/BayesCox.Rmd
@@ -8,7 +8,7 @@ vignette: >
 ---
 
 ```{r setup, include=FALSE}
-knitr::opts_chunk$set(echo = TRUE, eval = TRUE)
+knitr::opts_chunk$set(echo = TRUE, eval = FALSE)
 options(rmarkdown.html_vignette.check_title = FALSE)
 ```
 
@@ -65,6 +65,7 @@ fitBayesCox <- psbcSpeedUp(survObj, p=p, q=q, hyperpar=mypriorPara,
 nIter=1000, burnin=500, outFilePath="/tmp")
 ```
 
+
 ### Plot posterior estimates of regression cofficients
 
 The function `psbcSpeedUp::plot()` can show the posterior mean and 95% credible intervals of regression coefficients.
@@ -73,13 +74,23 @@ The function `psbcSpeedUp::plot()` can show the posterior mean and 95% credible
 plot(fitBayesCox)
 ```
 
+<img src="../man/figures/README_plot_beta.png" width="70%" />
+
+
 ### Plot time-dependent Brier scores
 
 The function `psbcSpeedUp::plotBrier()` can show the time-dependent Brier scores based on posterior mean of coefficients or Bayesian model averaging.
 
 ```{r, fig.width=6, fig.heigh=5}
 plotBrier(fitBayesCox, times = 80)
 ```
+```
+##     Null.model Bayesian.Cox
+## IBS  0.2089742   0.08195375
+```
+
+<img src="../man/figures/README_plot_brier.png" width="90%" />
+
 
 ### Predict survival probabilities and cumulative hazards
 
@@ -88,7 +99,20 @@ The function `psbcSpeedUp::predict()` can estimate the survival probabilities an
 ```{r}
 predict(fitBayesCox, type = c("cumhazard", "survival"))
 ```
-
+```
+##        observation   times cumhazard  survival
+##     1:           1   0.264  1.14e-05  1.00e+00
+##     2:           2   0.264  4.66e-05  1.00e+00
+##     3:           3   0.264  6.07e-05  1.00e+00
+##     4:           4   0.264  1.71e-05  1.00e+00
+##     5:           5   0.264  7.55e-05  1.00e+00
+##    ---                                        
+## 39996:         196 107.641  2.76e+00  6.36e-02
+## 39997:         197 107.641  5.65e-01  5.68e-01
+## 39998:         198 107.641  5.37e+01  4.86e-24
+## 39999:         199 107.641  4.25e+02 2.25e-185
+## 40000:         200 107.641  1.89e-01  8.28e-01
+```
 
 ## References