From a680ff00607cb9bd07eeb10662f6f8a8f911a407 Mon Sep 17 00:00:00 2001 From: zhiz Date: Tue, 11 Jun 2024 10:07:16 +0200 Subject: [PATCH] change folder structure --- {BayesSUR/.github => .github}/.gitignore | 0 BayesSUR/.github/workflows/R-CMD-check.yaml | 54 --- BayesSUR/NEWS.md | 50 --- BayesSUR/README.md | 424 ------------------ BayesSUR/DESCRIPTION => DESCRIPTION | 0 BayesSUR/LICENSE => LICENSE | 0 Makefile | 51 --- BayesSUR/NAMESPACE => NAMESPACE | 0 NEWS.md | 16 +- {BayesSUR/R => R}/BayesSUR.R | 0 {BayesSUR/R => R}/RcppExports.R | 0 {BayesSUR/R => R}/coef.BayesSUR.R | 0 {BayesSUR/R => R}/datasets.R | 0 {BayesSUR/R => R}/elpd.R | 0 {BayesSUR/R => R}/exampleEQTL.R | 0 {BayesSUR/R => R}/exampleGDSC.R | 0 {BayesSUR/R => R}/fitted.BayesSUR.R | 0 {BayesSUR/R => R}/getEstimator.R | 0 {BayesSUR/R => R}/plot.BayesSUR.R | 0 {BayesSUR/R => R}/plotCPO.R | 0 {BayesSUR/R => R}/plotEstimator.R | 0 {BayesSUR/R => R}/plotGraph.R | 0 {BayesSUR/R => R}/plotMCMCdiag.R | 0 {BayesSUR/R => R}/plotManhattan.R | 0 {BayesSUR/R => R}/plotNetwork.R | 0 {BayesSUR/R => R}/predict.BayesSUR.R | 0 {BayesSUR/R => R}/print.BayesSUR.R | 0 {BayesSUR/R => R}/summary.BayesSUR.R | 0 {BayesSUR/R => R}/targetGene.R | 0 README.md | 7 +- {BayesSUR/build => build}/vignette.rds | Bin check_and_build.R | 2 - BayesSUR/cleanup => cleanup | 0 BayesSUR/configure => configure | 0 BayesSUR/configure.ac => configure.ac | 0 {BayesSUR/data => data}/datalist | 0 {BayesSUR/data => data}/exampleEQTL.rda | Bin {BayesSUR/data => data}/exampleGDSC.rda | Bin {BayesSUR/data => data}/targetGene.rda | Bin {BayesSUR/inst => inst}/CITATION | 0 {BayesSUR/inst => inst}/doc/BayesSUR-RE.R | 0 {BayesSUR/inst => inst}/doc/BayesSUR-RE.Rmd | 0 {BayesSUR/inst => inst}/doc/BayesSUR-RE.html | 0 {BayesSUR/inst => inst}/doc/BayesSUR.pdf | Bin {BayesSUR/inst => inst}/doc/BayesSUR.pdf.asis | 0 main.cpp | 335 -------------- {BayesSUR/man => man}/BayesSUR.Rd | 0 {BayesSUR/man => man}/BayesSUR_internal.Rd | 0 {BayesSUR/man => man}/coef.BayesSUR.Rd | 0 {BayesSUR/man => man}/elpd.Rd | 0 {BayesSUR/man => man}/exampleEQTL.Rd | 0 {BayesSUR/man => man}/exampleGDSC.Rd | 0 {BayesSUR/man => man}/figures/figure2.png | Bin {BayesSUR/man => man}/fitted.BayesSUR.Rd | 0 {BayesSUR/man => man}/getEstimator.Rd | 0 {BayesSUR/man => man}/plot.BayesSUR.Rd | 0 {BayesSUR/man => man}/plotCPO.Rd | 0 {BayesSUR/man => man}/plotEstimator.Rd | 0 {BayesSUR/man => man}/plotGraph.Rd | 0 {BayesSUR/man => man}/plotMCMCdiag.Rd | 0 {BayesSUR/man => man}/plotManhattan.Rd | 0 {BayesSUR/man => man}/plotNetwork.Rd | 0 {BayesSUR/man => man}/predict.BayesSUR.Rd | 0 {BayesSUR/man => man}/print.BayesSUR.Rd | 0 {BayesSUR/man => man}/summary.BayesSUR.Rd | 0 {BayesSUR/man => man}/targetGene.Rd | 0 {BayesSUR/src => src}/BayesSUR.cpp | 0 {BayesSUR/src => src}/ESS_Atom.h | 0 {BayesSUR/src => src}/ESS_Sampler.h | 0 {BayesSUR/src => src}/HRR_Chain.cpp | 0 {BayesSUR/src => src}/HRR_Chain.h | 0 {BayesSUR/src => src}/Makevars.in | 0 {BayesSUR/src => src}/Makevars.win | 0 {BayesSUR/src => src}/Parameter_types.h | 0 {BayesSUR/src => src}/RcppExports.cpp | 0 {BayesSUR/src => src}/SUR_Chain.cpp | 0 {BayesSUR/src => src}/SUR_Chain.h | 0 {BayesSUR/src => src}/distr.cpp | 0 {BayesSUR/src => src}/distr.h | 0 {BayesSUR/src => src}/drive.cpp | 0 {BayesSUR/src => src}/drive.h | 0 {BayesSUR/src => src}/global.cpp | 0 {BayesSUR/src => src}/global.h | 0 {BayesSUR/src => src}/junction_tree.cpp | 0 {BayesSUR/src => src}/junction_tree.h | 0 {BayesSUR/src => src}/pugiconfig.hpp | 0 {BayesSUR/src => src}/pugixml.cpp | 0 {BayesSUR/src => src}/pugixml.hpp | 0 {BayesSUR/src => src}/utils.cpp | 0 {BayesSUR/src => src}/utils.h | 0 test.R | 52 --- .../vignettes => vignettes}/BayesSUR-RE.Rmd | 0 .../vignettes => vignettes}/BayesSUR.pdf | Bin .../vignettes => vignettes}/BayesSUR.pdf.asis | 0 94 files changed, 20 insertions(+), 971 deletions(-) rename {BayesSUR/.github => .github}/.gitignore (100%) delete mode 100644 BayesSUR/.github/workflows/R-CMD-check.yaml delete mode 100644 BayesSUR/NEWS.md delete mode 100644 BayesSUR/README.md rename BayesSUR/DESCRIPTION => DESCRIPTION (100%) rename BayesSUR/LICENSE => LICENSE (100%) delete mode 100644 Makefile rename BayesSUR/NAMESPACE => NAMESPACE (100%) rename {BayesSUR/R => R}/BayesSUR.R (100%) rename {BayesSUR/R => R}/RcppExports.R (100%) rename {BayesSUR/R => R}/coef.BayesSUR.R (100%) rename {BayesSUR/R => R}/datasets.R (100%) rename {BayesSUR/R => R}/elpd.R (100%) rename {BayesSUR/R => R}/exampleEQTL.R (100%) rename {BayesSUR/R => R}/exampleGDSC.R (100%) rename {BayesSUR/R => R}/fitted.BayesSUR.R (100%) rename {BayesSUR/R => R}/getEstimator.R (100%) rename {BayesSUR/R => R}/plot.BayesSUR.R (100%) rename {BayesSUR/R => R}/plotCPO.R (100%) rename {BayesSUR/R => R}/plotEstimator.R (100%) rename {BayesSUR/R => R}/plotGraph.R (100%) rename {BayesSUR/R => R}/plotMCMCdiag.R (100%) rename {BayesSUR/R => R}/plotManhattan.R (100%) rename {BayesSUR/R => R}/plotNetwork.R (100%) rename {BayesSUR/R => R}/predict.BayesSUR.R (100%) rename {BayesSUR/R => R}/print.BayesSUR.R (100%) rename {BayesSUR/R => R}/summary.BayesSUR.R (100%) rename {BayesSUR/R => R}/targetGene.R (100%) rename {BayesSUR/build => build}/vignette.rds (100%) delete mode 100644 check_and_build.R rename BayesSUR/cleanup => cleanup (100%) rename BayesSUR/configure => configure (100%) rename BayesSUR/configure.ac => configure.ac (100%) rename {BayesSUR/data => data}/datalist (100%) rename {BayesSUR/data => data}/exampleEQTL.rda (100%) rename {BayesSUR/data => data}/exampleGDSC.rda (100%) rename {BayesSUR/data => data}/targetGene.rda (100%) rename {BayesSUR/inst => inst}/CITATION (100%) rename {BayesSUR/inst => inst}/doc/BayesSUR-RE.R (100%) rename {BayesSUR/inst => inst}/doc/BayesSUR-RE.Rmd (100%) rename {BayesSUR/inst => inst}/doc/BayesSUR-RE.html (100%) rename {BayesSUR/inst => inst}/doc/BayesSUR.pdf (100%) rename {BayesSUR/inst => inst}/doc/BayesSUR.pdf.asis (100%) delete mode 100644 main.cpp rename {BayesSUR/man => man}/BayesSUR.Rd (100%) rename {BayesSUR/man => man}/BayesSUR_internal.Rd (100%) rename {BayesSUR/man => man}/coef.BayesSUR.Rd (100%) rename {BayesSUR/man => man}/elpd.Rd (100%) rename {BayesSUR/man => man}/exampleEQTL.Rd (100%) rename {BayesSUR/man => man}/exampleGDSC.Rd (100%) rename {BayesSUR/man => man}/figures/figure2.png (100%) rename {BayesSUR/man => man}/fitted.BayesSUR.Rd (100%) rename {BayesSUR/man => man}/getEstimator.Rd (100%) rename {BayesSUR/man => man}/plot.BayesSUR.Rd (100%) rename {BayesSUR/man => man}/plotCPO.Rd (100%) rename {BayesSUR/man => man}/plotEstimator.Rd (100%) rename {BayesSUR/man => man}/plotGraph.Rd (100%) rename {BayesSUR/man => man}/plotMCMCdiag.Rd (100%) rename {BayesSUR/man => man}/plotManhattan.Rd (100%) rename {BayesSUR/man => man}/plotNetwork.Rd (100%) rename {BayesSUR/man => man}/predict.BayesSUR.Rd (100%) rename {BayesSUR/man => man}/print.BayesSUR.Rd (100%) rename {BayesSUR/man => man}/summary.BayesSUR.Rd (100%) rename {BayesSUR/man => man}/targetGene.Rd (100%) rename {BayesSUR/src => src}/BayesSUR.cpp (100%) rename {BayesSUR/src => src}/ESS_Atom.h (100%) rename {BayesSUR/src => src}/ESS_Sampler.h (100%) rename {BayesSUR/src => src}/HRR_Chain.cpp (100%) rename {BayesSUR/src => src}/HRR_Chain.h (100%) rename {BayesSUR/src => src}/Makevars.in (100%) rename {BayesSUR/src => src}/Makevars.win (100%) rename {BayesSUR/src => src}/Parameter_types.h (100%) rename {BayesSUR/src => src}/RcppExports.cpp (100%) rename {BayesSUR/src => src}/SUR_Chain.cpp (100%) rename {BayesSUR/src => src}/SUR_Chain.h (100%) rename {BayesSUR/src => src}/distr.cpp (100%) rename {BayesSUR/src => src}/distr.h (100%) rename {BayesSUR/src => src}/drive.cpp (100%) rename {BayesSUR/src => src}/drive.h (100%) rename {BayesSUR/src => src}/global.cpp (100%) rename {BayesSUR/src => src}/global.h (100%) rename {BayesSUR/src => src}/junction_tree.cpp (100%) rename {BayesSUR/src => src}/junction_tree.h (100%) rename {BayesSUR/src => src}/pugiconfig.hpp (100%) rename {BayesSUR/src => src}/pugixml.cpp (100%) rename {BayesSUR/src => src}/pugixml.hpp (100%) rename {BayesSUR/src => src}/utils.cpp (100%) rename {BayesSUR/src => src}/utils.h (100%) delete mode 100644 test.R rename {BayesSUR/vignettes => vignettes}/BayesSUR-RE.Rmd (100%) rename {BayesSUR/vignettes => vignettes}/BayesSUR.pdf (100%) rename {BayesSUR/vignettes => vignettes}/BayesSUR.pdf.asis (100%) diff --git a/BayesSUR/.github/.gitignore b/.github/.gitignore similarity index 100% rename from BayesSUR/.github/.gitignore rename to .github/.gitignore diff --git a/BayesSUR/.github/workflows/R-CMD-check.yaml b/BayesSUR/.github/workflows/R-CMD-check.yaml deleted file mode 100644 index c4008370..00000000 --- a/BayesSUR/.github/workflows/R-CMD-check.yaml +++ /dev/null @@ -1,54 +0,0 @@ -# Workflow derived from https://github.com/r-lib/actions/tree/v2/examples -# Need help debugging build failures? Start at https://github.com/r-lib/actions#where-to-find-help -# -# NOTE: This workflow is overkill for most R packages and -# check-standard.yaml is likely a better choice. -# usethis::use_github_action("check-standard") will install it. -on: - push: - branches: [main, master] - pull_request: - branches: [main, master] - -name: R-CMD-check - -jobs: - R-CMD-check: - runs-on: ${{ matrix.config.os }} - - name: ${{ matrix.config.os }} (${{ matrix.config.r }}) - - strategy: - fail-fast: false - matrix: - config: - - { os: macos-latest, r: "release" } - - { os: ubuntu-latest, r: "devel", http-user-agent: "release" } - - { os: ubuntu-latest, r: "release" } - - { os: ubuntu-latest, r: "oldrel-1" } - # - { os: ubuntu-latest, r: "oldrel-2" } - # - { os: ubuntu-latest, r: "oldrel-3" } - - env: - GITHUB_PAT: ${{ secrets.GITHUB_TOKEN }} - R_KEEP_PKG_SOURCE: yes - - steps: - - uses: actions/checkout@v3 - - - uses: r-lib/actions/setup-pandoc@v2 - - - uses: r-lib/actions/setup-r@v2 - with: - r-version: ${{ matrix.config.r }} - http-user-agent: ${{ matrix.config.http-user-agent }} - use-public-rspm: true - - - uses: r-lib/actions/setup-r-dependencies@v2 - with: - extra-packages: any::rcmdcheck - needs: check - - - uses: r-lib/actions/check-r-package@v2 - with: - upload-snapshots: true diff --git a/BayesSUR/NEWS.md b/BayesSUR/NEWS.md deleted file mode 100644 index 6ef5b356..00000000 --- a/BayesSUR/NEWS.md +++ /dev/null @@ -1,50 +0,0 @@ -# BayesSUR 2.1-9 (or BayesSUR 2.2-0 for CRAN) - -* Fixed a bug about the use of temperature parameter in `HRR_Chain.cpp`, also minor update in `SUR_Chain.cpp` - -# BayesSUR 2.1-8 - -* Added argument `beta.type` in function `plotEstimator()` to plot MPM coefficient estimates -* Fixed argument `mrfG` index issue in function `BayesSUR()` - -# BayesSUR 2.1-7 - -* Increased the threshold of predictor dimension to 100000 (previously 5000) when pre-computing XtX -* Minor changes of the parameter updates of the HRR models, especially in function `step()` - -# BayesSUR 2.1-6 - -* Fixed issues of knitr VignetteBuilder and README - -# BayesSUR 2.1-5 - -* Updated citation, added README and one more vignette - -# BayesSUR 2.1-4 - -* Fixed a documentation issue by adding `@alias BayesSUR-package` - -# BayesSUR 2.1-3 - -* Updated C++11 specification to current default of C++17 -* Cleaned R scripts - -# BayesSUR 2.1-2 - -* Fixed gcc-UBSAN error `reference binding to null pointer` in SUR_Chain.cpp:3818 - -# BayesSUR 2.1-1 - -* Fixed R session aborted fatal error - -# BayesSUR 2.1-0 - -* Fixed an issue with `omp_set_num_threads` that did not work - -# BayesSUR 2.0-1 - -* Code update - -# BayesSUR 2.0-0 - -* Major update, introducing random effects for mandatory variables without variable selection. diff --git a/BayesSUR/README.md b/BayesSUR/README.md deleted file mode 100644 index 67db7dd0..00000000 --- a/BayesSUR/README.md +++ /dev/null @@ -1,424 +0,0 @@ - -# BayesSUR - - - -[![CRAN](http://www.r-pkg.org/badges/version/BayesSUR)](https://cran.r-project.org/package=BayesSUR) -[![r-universe](https://mbant.r-universe.dev/badges/BayesSUR)](https://mbant.r-universe.dev/BayesSUR) -[![R-CMD-check](https://github.com/zhizuio/BayesSUR/BayesSUR/workflows/R-CMD-check/badge.svg)](https://github.com/zhizuio/BayesSUR/actions) -[![License](https://img.shields.io/badge/License-MIT-yellow.svg)](https://opensource.org/licenses/MIT) - - - -This repository contains a [new and improved](inst/doc/BayesSUR.pdf) R package for high-dimensional multivariate Bayesian variable and covariance selection in linear regression, started as an interface to the [Bayesian SSUR](https://github.com/mbant/Bayesian_SSUR) C++-only, UNIX-specific, code. -This R package includes high-dimensional Bayesian SUR models from [Bottolo et al. (2021)](https://doi.org/10.1111/rssc.12490), [Zhao et al. (2021)](https://doi.org/10.18637/jss.v100.i11) and [Zhao et al. (2024)](https://doi.org/10.1093/jrsssc/qlad102). -See the package vignettes [`BayesSUR.pdf`](inst/doc/BayesSUR.pdf) for more information and an additional example below for the BayesSUR model with random effects. - -## Installation - -Install the latest released version from [CRAN](https://CRAN.R-project.org/package=BayesSUR) - -```r -install.packages("BayesSUR") -``` - -Install the latest development version from GitHub - -```r -#install.packages("remotes") -remotes::install_github("mbant/BayesSUR/BayesSUR") -``` - -## Additional example - -The BayesSUR model has been extended to include mandatory variables by assigning Gaussian priors as random effects rather than spike-and-slab priors, named as **SSUR-MRF with random effects** in [Zhao et al. 2023](https://doi.org/10.1093/jrsssc/qlad102). -The R code for the simulated data and real data analyses in [Zhao et al. 2023](https://doi.org/10.1093/jrsssc/qlad102) can be found at the GitHub repository [BayesSUR-RE](https://github.com/zhizuio/BayesSUR-RE). - -Here, we show a simulation example to run the BayesSUR mdoel with random effects. - -### Simulate data - -We design a network as the following figure (a) to construct a complex structure between $20$ response variables and $300$ predictors. -It assumes that the responses are divided into six groups, and the first $120$ predictors are divided into nine groups. - - - -_**Figure**: True relationships between response variables and predictors. -(a) Network structure between $\mathbf Y$ and $\mathbf X$. -(b) Spare latent indicator variable $\Gamma$ for the associations between $\mathbf Y$ and $\mathbf X$ in the SUR model. -Black blocks indicate nonzero coefficients and white blocks indicate zero coefficients. -(c) Additional structure in the residual covariance matrix between response variables not explained by $\mathbf X\mathbf B$. -Black blocks indicate correlated residuals of the corresponding response variables and white blocks indicate uncorrelated residuals of the corresponding response variables._ - -
- -Load the simulation function `sim.ssur()` as follows. - -```{r} -sim.ssur <- function(n, s, p, t0 = 0, seed = 123, mv = TRUE, - t.df = Inf, random.intercept = 0, intercept = TRUE) { - # set seed to fix coefficients - set.seed(7193) - sd_b <- 1 - mu_b <- 1 - b <- matrix(rnorm((p + ifelse(t0 == 0, 1, 0)) * s, mu_b, sd_b), p + ifelse(t0 == 0, 1, 0), s) - - # design groups and pathways of Gamma matrix - gamma <- matrix(FALSE, p + ifelse(t0 == 0, 1, 0), s) - if (t0 == 0) gamma[1, ] <- TRUE - gamma[2:6 - ifelse(t0 == 0, 0, 1), 1:5] <- TRUE - gamma[11:21 - ifelse(t0 == 0, 0, 1), 6:12] <- TRUE - gamma[31:51 - ifelse(t0 == 0, 0, 1), 1:5] <- TRUE - gamma[31:51 - ifelse(t0 == 0, 0, 1), 13:15] <- TRUE - gamma[52:61 - ifelse(t0 == 0, 0, 1), 1:12] <- TRUE - gamma[71:91 - ifelse(t0 == 0, 0, 1), 6:15] <- TRUE - gamma[111:121 - ifelse(t0 == 0, 0, 1), 1:15] <- TRUE - gamma[122 - ifelse(t0 == 0, 0, 1), 16:18] <- TRUE - gamma[123 - ifelse(t0 == 0, 0, 1), 19] <- TRUE - gamma[124 - ifelse(t0 == 0, 0, 1), 20] <- TRUE - - G_kron <- matrix(0, s * p, s * p) - G_m <- bdiag(matrix(1, ncol = 5, nrow = 5), - matrix(1, ncol = 7, nrow = 7), - matrix(1, ncol = 8, nrow = 8)) - G_p <- bdiag(matrix(1, ncol = 5, nrow = 5), diag(3), - matrix(1, ncol = 11, nrow = 11), diag(9), - matrix(1, ncol = 21, nrow = 21), - matrix(1, ncol = 10, nrow = 10), diag(9), - matrix(1, ncol = 21, nrow = 21), diag(19), - matrix(1, ncol = 11, nrow = 11), diag(181)) - G_kron <- kronecker(G_m, G_p) - - combn11 <- combn(rep((1:5 - 1) * p, each = length(1:5)) + - rep(1:5, times = length(1:5)), 2) - combn12 <- combn(rep((1:5 - 1) * p, each = length(30:60)) + - rep(30:60, times = length(1:5)), 2) - combn13 <- combn(rep((1:5 - 1) * p, each = length(110:120)) + - rep(110:120, times = length(1:5)), 2) - combn21 <- combn(rep((6:12 - 1) * p, each = length(10:20)) + - rep(10:20, times = length(6:12)), 2) - combn22 <- combn(rep((6:12 - 1) * p, each = length(51:60)) + - rep(51:60, times = length(6:12)), 2) - combn23 <- combn(rep((6:12 - 1) * p, each = length(70:90)) + - rep(70:90, times = length(6:12)), 2) - combn24 <- combn(rep((6:12 - 1) * p, each = length(110:120)) + - rep(110:120, times = length(6:12)), 2) - combn31 <- combn(rep((13:15 - 1) * p, each = length(30:50)) + - rep(30:50, times = length(13:15)), 2) - combn32 <- combn(rep((13:15 - 1) * p, each = length(70:90)) + - rep(70:90, times = length(13:15)), 2) - combn33 <- combn(rep((13:15 - 1) * p, each = length(110:120)) + - rep(110:120, times = length(13:15)), 2) - combn4 <- combn(rep((16:18 - 1) * p, each = length(121)) + - rep(121, times = length(16:18)), 2) - combn5 <- matrix(rep((19 - 1) * p, each = length(122)) + - rep(122, times = length(19)), nrow = 1, ncol = 2) - combn6 <- matrix(rep((20 - 1) * p, each = length(123)) + - rep(123, times = length(20)), nrow = 1, ncol = 2) - - combnAll <- rbind(t(combn11), t(combn12), t(combn13), - t(combn21), t(combn22), t(combn23), t(combn24), - t(combn31), t(combn32), t(combn33), - t(combn4), combn5, combn6) - - set.seed(seed + 7284) - sd_x <- 1 - x <- matrix(rnorm(n * p, 0, sd_x), n, p) - - if (t0 == 0 & intercept) x <- cbind(rep(1, n), x) - if (!intercept) { - gamma <- gamma[-1, ] - b <- b[-1, ] - } - xb <- matrix(NA, n, s) - if (mv) { - for (i in 1:s) { - if (sum(gamma[, i]) >= 1) { - if (sum(gamma[, i]) == 1) { - xb[, i] <- x[, gamma[, i]] * b[gamma[, i], i] - } else { - xb[, i] <- x[, gamma[, i]] %*% b[gamma[, i], i] - } - } else { - xb[, i] <- sapply(1:s, function(i) rep(1, n) * b[1, i]) - } - } - } else { - if (sum(gamma) >= 1) { - xb <- x[, gamma] %*% b[gamma, ] - } else { - xb <- sapply(1:s, function(i) rep(1, n) * b[1, i]) - } - } - - corr_param <- 0.9 - M <- matrix(corr_param, s, s) - diag(M) <- rep(1, s) - - ## wanna make it decomposable - Prime <- list(c(1:(s * .4), (s * .8):s), - c((s * .4):(s * .6)), - c((s * .65):(s * .75)), - c((s * .8):s)) - G <- matrix(0, s, s) - for (i in 1:length(Prime)) { - G[Prime[[i]], Prime[[i]]] <- 1 - } - - # check - dimnames(G) <- list(1:s, 1:s) - length(gRbase::mcsMAT(G - diag(s))) > 0 - - var <- solve(BDgraph::rgwish(n = 1, adj = G, b = 3, D = M)) - - # change seeds to add randomness on error - set.seed(seed + 8493) - sd_err <- 0.5 - if (is.infinite(t.df)) { - err <- matrix(rnorm(n * s, 0, sd_err), n, s) %*% chol(as.matrix(var)) - } else { - err <- matrix(rt(n * s, t.df), n, s) %*% chol(as.matrix(var)) - } - - if (t0 == 0) { - b.re <- NA - z <- NA - y <- xb + err - if (random.intercept != 0) { - y <- y + matrix(rnorm(n * s, 0, sqrt(random.intercept)), n, s) - } - - z <- sample(1:4, n, replace = T, prob = rep(1 / 4, 4)) - - return(list(y = y, x = x, b = b, gamma = gamma, z = model.matrix(~ factor(z) + 0)[, ], - b.re = b.re, Gy = G, mrfG = combnAll)) - } else { - # add random effects - z <- t(rmultinom(n, size = 1, prob = c(.1, .2, .3, .4))) - z <- sample(1:t0, n, replace = T, prob = rep(1 / t0, t0)) - set.seed(1683) - b.re <- rnorm(t0, 0, 2) - y <- matrix(b.re[z], nrow = n, ncol = s) + xb + err - - return(list( - y = y, x = x, b = b, gamma = gamma, z = model.matrix(~ factor(z) + 0)[, ], - b.re = b.re, Gy = G, mrfG = combnAll - )) - } -} -``` - -To simulate data with sample size $n=250$, responsible variables $s=20$ and covariates $p=300$, we can specify the corresponding parameters in the function `sim.ssur()` as follows. - -```{r} -library("BayesSUR") -library("Matrix") -n <- 250 -s <- 20 -p <- 300 -sim1 <- sim.ssur(n, s, p, seed = 1) -``` - -To simulate data from $4$ individual groups with group indicator variables following the defaul multinomial distribution $multinomial(0.1,0.2,0.3,0.4)$, we can simply add the argument `t0 = 4` in the function `sim.ssur()` as follows. - -```{r} -t0 <- 4 -sim2 <- sim.ssur(n, s, p, t0, seed = 1) # learning data -sim2.val <- sim.ssur(n, s, p, t0, seed=101) # validation data -``` - -### Run BayesSUR model with random effects - -According to the guideline of prior specification in [Zhao et al. 2023](https://doi.org/10.1093/jrsssc/qlad102), we first set the following parameters `hyperpar` and then running the BayesSUR model with random effects via `betaPrior = "reGroup"` (default `betaPrior = "independent"` with spike-and-slab priors for all coefficients). -**For illustration, we run a short MCMC** with `nIter = 300` and `burnin = 100`. -Note that here the graph used for the Markov random field prior is the true graph from the returned object of the simulation `sim2$mrfG`. - -```{r} -hyperpar <- list(mrf_d = -2, mrf_e = 1.6, a_w0 = 100, b_w0 = 500, a_w = 15, b_w = 60) -set.seed(1038) -fit2 <- BayesSUR( - data = cbind(sim2$y, sim2$z, sim2$x), - Y = 1:s, - X_0 = s + 1:t0, - X = s + t0 + 1:p, - outFilePath = "sim2_mrf_re", - hyperpar = hyperpar, - gammaInit = "0", - betaPrior = "reGroup", - nIter = 300, burnin = 100, - covariancePrior = "HIW", - standardize = F, - standardize.response = F, - gammaPrior = "MRF", - mrfG = sim2$mrfG, - output_CPO = T -) -``` - -``` -## BayesSUR -- Bayesian Seemingly Unrelated Regression Modelling -## Reading input files ... ... successfull! -## Clearing and initialising output files -## Initialising the (SUR) MCMC Chain ... ... DONE! -## Drafting the output files with the start of the chain ... DONE! -## -## Starting 2 (parallel) chain(s) for 300 iterations: -## Temperature ladder updated, new temperature ratio : 1.1 -## MCMC ends. --- Saving results and exiting -## Saved to : sim2_mrf_re1/data_SSUR_****_out.txt -## Final w : 0.148291 -## Final tau : 1.84125 w/ proposal variance: 0.408163 -## Final eta : 0.0355005 -## -- Average Omega : 0 -## Final temperature ratio : 1.1 -## -## DONE, exiting! -``` - -Check some summarized information of the results: - -```{r} -summary(fit2) -``` - -``` -## -## Call: -## BayesSUR(data = cbind(sim2$y, sim2$z, sim2$x), ...) -## -## CPOs: -## Min. 1st Qu. Median 3rd Qu. Max. -## 0.0001118321 0.0241323466 0.0349716031 0.0456556652 0.2321280902 -## -## Number of selected predictors (mPIP > 0.5): 2843 of 20x300 -## -## Top 10 predictors on average mPIP across all responses: -## X.130 X.54 X.249 X.56 X.77 X.253 X.281 X.80 -## 0.720640 0.706705 0.652730 0.650985 0.643780 0.640780 0.639045 0.636565 -## X.260 X.297 -## 0.634820 0.629595 -## -## Top 10 responses on average mPIP across all predictors: -## X.8 X.5 X.12 X.6 X.18 X.4 X.14 X.1 -## 0.4957363 0.4879933 0.4873303 0.4860670 0.4846080 0.4828333 0.4784240 0.4773090 -## X.19 X.2 -## 0.4756350 0.4742257 -## -## Expected log pointwise predictive density (elpd) estimates: -## elpd.LOO = -16437.89, elpd.WAIC = -16470.16 -## -## MCMC specification: -## iterations = 300, burn-in = 100, chains = 2 -## gamma local move sampler: bandit -## gamma initialisation: 0 -## -## Model specification: -## covariance prior: HIW -## gamma prior: MRF -## -## Hyper-parameters: -## a_w b_w nu a_tau b_tau a_eta b_eta mrf_d mrf_e a_w0 b_w0 -## 15.0 60.0 22.0 0.1 10.0 0.1 1.0 -2.0 1.6 100.0 500.0 -``` - -Compute the model performace with respect to **variable selection** - -```{r} -# compute accuracy, sensitivity, specificity of variable selection -gamma <- getEstimator(fit2) -(accuracy <- sum(data.matrix(gamma > 0.5) == sim2$gamma) / prod(dim(gamma))) -``` - -``` -## [1] 0.5358333 -``` - -```{r} -(sensitivity <- sum((data.matrix(gamma > 0.5) == 1) & (sim2$gamma == 1)) / sum(sim2$gamma == 1)) -``` - -``` -## [1] 0.5376623 -``` - -```{r} -(specificity <- sum((data.matrix(gamma > 0.5) == 0) & (sim2$gamma == 0)) / sum(sim2$gamma == 0)) -``` - -``` -## [1] 0.5355641 -``` - -Compute the model performance with respect to **response prediction** - -```{r} -# compute RMSE and RMSPE for prediction performance -beta <- getEstimator(fit2, estimator = "beta", Pmax = .5, beta.type = "conditional") -(RMSE <- sqrt(sum((sim2$y - cbind(sim2$z, sim2$x) %*% beta)^2) / prod(dim(sim2$y)))) -``` - -``` -## [1] 7.025327 -``` - -```{r} -(RMSPE <- sqrt(sum((sim2.val$y - cbind(sim2.val$z, sim2.val$x) %*% beta)^2) / prod(dim(sim2.val$y)))) -``` - -``` -## [1] 8.084381 -``` - -Compute the model performance with respect to **coefficient bias** - -```{r} -# compute bias of beta estimates -b <- sim2$b -b[sim2$gamma == 0] <- 0 -(beta.l2 <- sqrt(sum((beta[-c(1:4), ] - b)^2) / prod(dim(b)))) -``` - -``` -## [1] 0.4530592 -``` - -Compute the model performance with respect to **covariance selection** - -```{r} -g.re <- getEstimator(fit2, estimator = "Gy") -(g.accuracy <- sum((g.re > 0.5) == sim2$Gy) / prod(dim(g.re))) -``` - -``` -## [1] 0.545 -``` - -```{r} -(g.sensitivity <- sum(((g.re > 0.5) == sim2$Gy)[sim2$Gy == 1]) / sum(sim2$Gy == 1)) -``` - -``` -## [1] 0.1287129 -``` - -```{r} -(g.specificity <- sum(((g.re > 0.5) == sim2$Gy)[sim2$Gy == 0]) / sum(sim2$Gy == 0)) -``` - -``` -## [1] 0.969697 -``` - -## References - -> Leonardo Bottolo, Marco Banterle, Sylvia Richardson, Mika Ala-Korpela, Marjo-Riitta Järvelin, Alex Lewin (2021). -> A computationally efficient Bayesian seemingly unrelated regressions model for high-dimensional quantitative trait loci discovery. -> _Journal of the Royal Statistical Society: Series C (Applied Statistics)_, 70(4):886-908. DOI: [10.1111/rssc.12490](https://doi.org/10.1111/rssc.12490). - -> Zhi Zhao, Marco Banterle, Leonardo Bottolo, Sylvia Richardson, Alex Lewin, Manuela Zucknick (2021). -> BayesSUR: An R package for high-dimensional multivariate Bayesian variable and covariance selection in linear regression. -> _Journal of Statistical Software_, 100(11):1-32. DOI: [10.18637/jss.v100.i11](https://doi.org/10.18637/jss.v100.i11). - -> Zhi Zhao, Marco Banterle, Alex Lewin, Manuela Zucknick (2023). -> Multivariate Bayesian structured variable selection for pharmacogenomic studies. -> _Journal of the Royal Statistical Society: Series C (Applied Statistics)_, 73(2):420-443 qlad102. DOI: [10.1093/jrsssc/qlad102](https://doi.org/10.1093/jrsssc/qlad102). diff --git a/BayesSUR/DESCRIPTION b/DESCRIPTION similarity index 100% rename from BayesSUR/DESCRIPTION rename to DESCRIPTION diff --git a/BayesSUR/LICENSE b/LICENSE similarity index 100% rename from BayesSUR/LICENSE rename to LICENSE diff --git a/Makefile b/Makefile deleted file mode 100644 index 06e16316..00000000 --- a/Makefile +++ /dev/null @@ -1,51 +0,0 @@ -# Warning: To compile C++ manually, all Rcpp random number generation in BayesSUR/src/distr.cpp need to be adapted - -SOURCE_DIR=BayesSUR/src -VPATH=$(SOURCE_DIR) - -CC=g++ -# Uncomment clang++ and comment g++ to check compilation on both systems -#CC=clang++ -fsanitize=address,undefined -fno-sanitize=float-divide-by-zero -fno-sanitize=alignment -fno-omit-frame-pointer -g - -CFLAGS= -c -Wall -Wno-reorder -std=c++11 -I$(SOURCE_DIR)/ -DCCODE -fopenmp - -OPENLDFLAGS= -larmadillo -lpthread -lopenblas -fopenmp -NVLDFLAGS= -larmadillo -lpthread -lnvblas -fopenmp - -SOURCES_BVS=$(SOURCE_DIR)/global.cpp $(SOURCE_DIR)/utils.cpp $(SOURCE_DIR)/distr.cpp $(SOURCE_DIR)/junction_tree.cpp $(SOURCE_DIR)/HRR_Chain.cpp $(SOURCE_DIR)/SUR_Chain.cpp $(SOURCE_DIR)/drive.cpp main.cpp -#ESS_Atom.h and Parameters_type.h are interface only -OBJECTS_BVS=$(SOURCES_BVS:.cpp=.o) - -SOURCES_XML=$(SOURCE_DIR)/pugixml.cpp -OBJECTS_XML=$(SOURCES_XML:.cpp=.o) - -all:$(SOURCES_XML) $(SOURCES_BVS) BVS_NONVIDIA - -BVS_NVIDIA: OPTIM_FLAGS := -O3 -BVS_NVIDIA: $(OBJECTS_BVS) - @echo [Linking and producing executable]: - if [ -e /usr/lib/libnvblas.so -o -e /usr/lib64/libnvblas.so ] ; then $(CC) $(OBJECTS_BVS) -o BVS_Reg $(NVLDFLAGS) ; else $(CC) $(OBJECTS_BVS) -o BVS_Reg $(OPENLDFLAGS) ; fi -# I hate that this is searching only in the standard path, the correct way is with automake ./config, but that will come later... -# if you compile against nvblas run with LD_PRELOAD=/usr/lib64/libnvblas.so or similar path to libnvblas.so - -BVS_NONVIDIA: OPTIM_FLAGS := -O3 -BVS_NONVIDIA: $(OBJECTS_XML) $(OBJECTS_BVS) - @echo [Linking and producing executable]: - $(CC) $(OBJECTS_XML) $(OBJECTS_BVS) -o BVS_Reg $(OPENLDFLAGS) - -BVS_DEBUG: OPTIM_FLAGS := -O0 #maybe O1 ? -pg ? linker as well? -BVS_DEBUG: $(OBJECTS_XML) $(OBJECTS_BVS) - @echo [Linking and producing executable]: - $(CC) $(OBJECTS_XML) $(OBJECTS_BVS) -o BVS_DEBUG_Reg $(OPENLDFLAGS) -ggdb3 -g -lprofiler - -%.o: %.cpp - @echo [Compiling]: $< - $(CC) $(CFLAGS) $(OPTIM_FLAGS) -o $@ -c $< - -clean: - @echo [Cleaning: ] - rm *.o; rm $(SOURCE_DIR)/*.o; rm *_Reg; rm -rf results; - -remake: - @echo [Cleaning compilation objets only: ] - rm *.o; rm *_Reg; rm $(SOURCE_DIR)/*.o; diff --git a/BayesSUR/NAMESPACE b/NAMESPACE similarity index 100% rename from BayesSUR/NAMESPACE rename to NAMESPACE diff --git a/NEWS.md b/NEWS.md index aa20e676..6ef5b356 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,6 +1,20 @@ +# BayesSUR 2.1-9 (or BayesSUR 2.2-0 for CRAN) + +* Fixed a bug about the use of temperature parameter in `HRR_Chain.cpp`, also minor update in `SUR_Chain.cpp` + +# BayesSUR 2.1-8 + +* Added argument `beta.type` in function `plotEstimator()` to plot MPM coefficient estimates +* Fixed argument `mrfG` index issue in function `BayesSUR()` + +# BayesSUR 2.1-7 + +* Increased the threshold of predictor dimension to 100000 (previously 5000) when pre-computing XtX +* Minor changes of the parameter updates of the HRR models, especially in function `step()` + # BayesSUR 2.1-6 -* Fixed knitr VignetteBuilder +* Fixed issues of knitr VignetteBuilder and README # BayesSUR 2.1-5 diff --git a/BayesSUR/R/BayesSUR.R b/R/BayesSUR.R similarity index 100% rename from BayesSUR/R/BayesSUR.R rename to R/BayesSUR.R diff --git a/BayesSUR/R/RcppExports.R b/R/RcppExports.R similarity index 100% rename from BayesSUR/R/RcppExports.R rename to R/RcppExports.R diff --git a/BayesSUR/R/coef.BayesSUR.R b/R/coef.BayesSUR.R similarity index 100% rename from BayesSUR/R/coef.BayesSUR.R rename to R/coef.BayesSUR.R diff --git a/BayesSUR/R/datasets.R b/R/datasets.R similarity index 100% rename from BayesSUR/R/datasets.R rename to R/datasets.R diff --git a/BayesSUR/R/elpd.R b/R/elpd.R similarity index 100% rename from BayesSUR/R/elpd.R rename to R/elpd.R diff --git a/BayesSUR/R/exampleEQTL.R b/R/exampleEQTL.R similarity index 100% rename from BayesSUR/R/exampleEQTL.R rename to R/exampleEQTL.R diff --git a/BayesSUR/R/exampleGDSC.R b/R/exampleGDSC.R similarity index 100% rename from BayesSUR/R/exampleGDSC.R rename to R/exampleGDSC.R diff --git a/BayesSUR/R/fitted.BayesSUR.R b/R/fitted.BayesSUR.R similarity index 100% rename from BayesSUR/R/fitted.BayesSUR.R rename to R/fitted.BayesSUR.R diff --git a/BayesSUR/R/getEstimator.R b/R/getEstimator.R similarity index 100% rename from BayesSUR/R/getEstimator.R rename to R/getEstimator.R diff --git a/BayesSUR/R/plot.BayesSUR.R b/R/plot.BayesSUR.R similarity index 100% rename from BayesSUR/R/plot.BayesSUR.R rename to R/plot.BayesSUR.R diff --git a/BayesSUR/R/plotCPO.R b/R/plotCPO.R similarity index 100% rename from BayesSUR/R/plotCPO.R rename to R/plotCPO.R diff --git a/BayesSUR/R/plotEstimator.R b/R/plotEstimator.R similarity index 100% rename from BayesSUR/R/plotEstimator.R rename to R/plotEstimator.R diff --git a/BayesSUR/R/plotGraph.R b/R/plotGraph.R similarity index 100% rename from BayesSUR/R/plotGraph.R rename to R/plotGraph.R diff --git a/BayesSUR/R/plotMCMCdiag.R b/R/plotMCMCdiag.R similarity index 100% rename from BayesSUR/R/plotMCMCdiag.R rename to R/plotMCMCdiag.R diff --git a/BayesSUR/R/plotManhattan.R b/R/plotManhattan.R similarity index 100% rename from BayesSUR/R/plotManhattan.R rename to R/plotManhattan.R diff --git a/BayesSUR/R/plotNetwork.R b/R/plotNetwork.R similarity index 100% rename from BayesSUR/R/plotNetwork.R rename to R/plotNetwork.R diff --git a/BayesSUR/R/predict.BayesSUR.R b/R/predict.BayesSUR.R similarity index 100% rename from BayesSUR/R/predict.BayesSUR.R rename to R/predict.BayesSUR.R diff --git a/BayesSUR/R/print.BayesSUR.R b/R/print.BayesSUR.R similarity index 100% rename from BayesSUR/R/print.BayesSUR.R rename to R/print.BayesSUR.R diff --git a/BayesSUR/R/summary.BayesSUR.R b/R/summary.BayesSUR.R similarity index 100% rename from BayesSUR/R/summary.BayesSUR.R rename to R/summary.BayesSUR.R diff --git a/BayesSUR/R/targetGene.R b/R/targetGene.R similarity index 100% rename from BayesSUR/R/targetGene.R rename to R/targetGene.R diff --git a/README.md b/README.md index 63d2ea39..67db7dd0 100644 --- a/README.md +++ b/README.md @@ -3,7 +3,10 @@ -[![CRAN_Status_Badge](http://www.r-pkg.org/badges/version/BayesSUR)](https://cran.r-project.org/package=BayesSUR) +[![CRAN](http://www.r-pkg.org/badges/version/BayesSUR)](https://cran.r-project.org/package=BayesSUR) +[![r-universe](https://mbant.r-universe.dev/badges/BayesSUR)](https://mbant.r-universe.dev/BayesSUR) +[![R-CMD-check](https://github.com/zhizuio/BayesSUR/BayesSUR/workflows/R-CMD-check/badge.svg)](https://github.com/zhizuio/BayesSUR/actions) +[![License](https://img.shields.io/badge/License-MIT-yellow.svg)](https://opensource.org/licenses/MIT) @@ -38,7 +41,7 @@ Here, we show a simulation example to run the BayesSUR mdoel with random effects We design a network as the following figure (a) to construct a complex structure between $20$ response variables and $300$ predictors. It assumes that the responses are divided into six groups, and the first $120$ predictors are divided into nine groups. - + _**Figure**: True relationships between response variables and predictors. (a) Network structure between $\mathbf Y$ and $\mathbf X$. diff --git a/BayesSUR/build/vignette.rds b/build/vignette.rds similarity index 100% rename from BayesSUR/build/vignette.rds rename to build/vignette.rds diff --git a/check_and_build.R b/check_and_build.R deleted file mode 100644 index 0be68cc6..00000000 --- a/check_and_build.R +++ /dev/null @@ -1,2 +0,0 @@ -devtools::check("BayesSUR/",cran = TRUE,build_args=c('--compact-vignettes=both')) -devtools::build("BayesSUR/",args=c('--compact-vignettes=both')) diff --git a/BayesSUR/cleanup b/cleanup similarity index 100% rename from BayesSUR/cleanup rename to cleanup diff --git a/BayesSUR/configure b/configure similarity index 100% rename from BayesSUR/configure rename to configure diff --git a/BayesSUR/configure.ac b/configure.ac similarity index 100% rename from BayesSUR/configure.ac rename to configure.ac diff --git a/BayesSUR/data/datalist b/data/datalist similarity index 100% rename from BayesSUR/data/datalist rename to data/datalist diff --git a/BayesSUR/data/exampleEQTL.rda b/data/exampleEQTL.rda similarity index 100% rename from BayesSUR/data/exampleEQTL.rda rename to data/exampleEQTL.rda diff --git a/BayesSUR/data/exampleGDSC.rda b/data/exampleGDSC.rda similarity index 100% rename from BayesSUR/data/exampleGDSC.rda rename to data/exampleGDSC.rda diff --git a/BayesSUR/data/targetGene.rda b/data/targetGene.rda similarity index 100% rename from BayesSUR/data/targetGene.rda rename to data/targetGene.rda diff --git a/BayesSUR/inst/CITATION b/inst/CITATION similarity index 100% rename from BayesSUR/inst/CITATION rename to inst/CITATION diff --git a/BayesSUR/inst/doc/BayesSUR-RE.R b/inst/doc/BayesSUR-RE.R similarity index 100% rename from BayesSUR/inst/doc/BayesSUR-RE.R rename to inst/doc/BayesSUR-RE.R diff --git a/BayesSUR/inst/doc/BayesSUR-RE.Rmd b/inst/doc/BayesSUR-RE.Rmd similarity index 100% rename from BayesSUR/inst/doc/BayesSUR-RE.Rmd rename to inst/doc/BayesSUR-RE.Rmd diff --git a/BayesSUR/inst/doc/BayesSUR-RE.html b/inst/doc/BayesSUR-RE.html similarity index 100% rename from BayesSUR/inst/doc/BayesSUR-RE.html rename to inst/doc/BayesSUR-RE.html diff --git a/BayesSUR/inst/doc/BayesSUR.pdf b/inst/doc/BayesSUR.pdf similarity index 100% rename from BayesSUR/inst/doc/BayesSUR.pdf rename to inst/doc/BayesSUR.pdf diff --git a/BayesSUR/inst/doc/BayesSUR.pdf.asis b/inst/doc/BayesSUR.pdf.asis similarity index 100% rename from BayesSUR/inst/doc/BayesSUR.pdf.asis rename to inst/doc/BayesSUR.pdf.asis diff --git a/main.cpp b/main.cpp deleted file mode 100644 index 0074c99f..00000000 --- a/main.cpp +++ /dev/null @@ -1,335 +0,0 @@ -#include -#include -#include - -#ifndef CCODE -#include -#endif - -// redeclare the drive funciton -int drive( const std::string& dataFile, const std::string& mrfGFile, const std::string& blockFile, const std::string& structureGraphFile, const std::string& hyperParFile, const std::string& outFilePath, - unsigned int nIter, unsigned int burnin, unsigned int nChains, - const std::string& covariancePrior, - const std::string& gammaPrior, const std::string& gammaSampler, const std::string& gammaInit, - const std::string& betaPrior, - bool output_gamma, bool output_beta, bool output_G, bool output_sigmaRho, bool output_pi, bool output_tail, bool output_model_size, bool output_CPO, bool output_model_visit ); - -int main(int argc, char* argv[]) -{ - - unsigned int nIter = 10; // default number of iterations - unsigned int burnin = 0; - unsigned int nChains = 1; - - std::string dataFile = "data.txt"; - std::string mrfGFile = "mrfG.txt"; - std::string blockFile = "blocks.txt"; - std::string structureGraphFile = "structureGraph.txt"; - std::string hpFile = "model.xml"; - - std::string outFilePath = ""; - - std::string covariancePrior = ""; - - std::string gammaPrior = ""; -// std::string mrfGFile = ""; - std::string gammaSampler = "bandit"; - std::string gammaInit = "MLE"; - - std::string betaPrior = "independent"; - - bool out_gamma = true, out_beta = true, out_G = true, - out_sigmaRho = true, out_pi = true, out_tail = true, - out_model_size = true, out_CPO = true, out_model_visit = false; - - // ### Read and interpret command line (to put in a separate file / function?) - int na = 1; - while(na < argc) - { - if ( 0 == std::string{argv[na]}.compare(std::string{"--covariancePrior"}) ) - { - covariancePrior = std::string(argv[++na]); // use the next - - if ( covariancePrior == "sparse" || covariancePrior == "Sparse" || covariancePrior == "SPARSE" || covariancePrior == "HIW" || covariancePrior == "hiw" ) - covariancePrior = "HIW"; - else if ( covariancePrior == "dense" || covariancePrior == "Dense" || covariancePrior == "DENSE" || covariancePrior == "IW" || covariancePrior == "iw" ) - covariancePrior = "IW"; - else if ( covariancePrior == "independent" || covariancePrior == "Independent" || covariancePrior == "INDEPENDENT" || covariancePrior == "INDEP" || covariancePrior == "indep" || covariancePrior == "IG" || covariancePrior == "ig" ) - covariancePrior = "IG"; - else - { - std::cout << "Unknown covariancePrior argument: only sparse (HIW), dense(IW) or independent (IG) are available" << std::endl; - return(1); - } - - if (na+1==argc) break; // in case it's last, break - ++na; // otherwise augment counter - } - else if ( 0 == std::string{argv[na]}.compare(std::string{"--mrfGFile"}) ) - { - mrfGFile = ""+std::string(argv[++na]); // use the next - if (na+1==argc) break; // in case it's last, break - ++na; // otherwise augment counter - } - else if ( 0 == std::string{argv[na]}.compare(std::string{"--gammaPrior"}) ) - { - gammaPrior = std::string(argv[++na]); // use the next - - if ( gammaPrior == "hotspot" || gammaPrior == "HOTSPOT" || gammaPrior == "hotspots" || gammaPrior == "HOTSPOTS" || gammaPrior == "hs" || gammaPrior == "HS" ) - gammaPrior = "hotspot"; - else if ( gammaPrior == "MRF" || gammaPrior == "mrf" || gammaPrior == "markov random field" || gammaPrior == "Markov Random Field" ) - gammaPrior = "MRF"; - else if ( gammaPrior == "hierarchical" || gammaPrior == "h" || gammaPrior == "H" ) - gammaPrior = "hierarchical"; - - else - { - std::cout << "Unknown gammaPrior argument: only hotspot, MRF or hierarchical are available" << std::endl; - return(1); - } - - if (na+1==argc) break; // in case it's last, break - ++na; // otherwise augment counter - } - else if ( 0 == std::string{argv[na]}.compare(std::string{"--gammaSampler"}) ) - { - gammaSampler = std::string(argv[++na]); // use the next - - if ( gammaSampler == "MC3" || gammaSampler == "mc3" ) - gammaSampler = "MC3"; - else if ( gammaSampler == "Bandit" || gammaSampler == "bandit" || gammaSampler == "BANDIT" ) - gammaSampler = "bandit"; - else - { - std::cout << "Unknown gammaSampler method: only Bandit or MC3 are available" << std::endl; - return(1); - } - - if (na+1==argc) break; // in case it's last, break - ++na; // otherwise augment counter - } - else if ( 0 == std::string{argv[na]}.compare(std::string{"--gammaInit"}) ) - { - gammaInit = std::string(argv[++na]); // use the next - - if( gammaInit != "R" && gammaInit != "0" && gammaInit != "1" && gammaInit != "MLE") - { - std::cout << "Unknown gammaInit method: only allowed:\n\t*\tR: random init (0.5 probability)\n\t*\t0: all elements set to 0\n\t*\t1: all elements set to 1\n\t*\tMLE: computes MLE for beta and init gamma for all significant coeffs" << std::endl; - return(1); - } - - if (na+1==argc) break; // in case it's last, break - ++na; // otherwise augment counter - } - else if ( 0 == std::string{argv[na]}.compare(std::string{"--betaPrior"}) ) - { - betaPrior = std::string(argv[++na]); // use the next - - if ( betaPrior == "independent" || betaPrior == "Independent" || betaPrior == "INDEPENDENT" || betaPrior == "indep" || betaPrior == "Indep" || betaPrior == "i" || betaPrior == "I" ) - betaPrior = "independent"; - else if ( betaPrior == "gprior" || betaPrior == "gPrior" || betaPrior == "g-prior" || betaPrior == "G-Prior" || betaPrior == "GPRIOR" ) - betaPrior = "g-prior"; - else - { - std::cout << "Unknown betaPrior method: only independent is available as of yet" << std::endl; - return(1); - } - - if (na+1==argc) break; // in case it's last, break - ++na; // otherwise augment counter - } - else if ( 0 == std::string{argv[na]}.compare(std::string{"--nIter"}) ) - { - nIter = std::stoi(argv[++na]); - if (na+1==argc) break; - ++na; - } - else if ( 0 == std::string{argv[na]}.compare(std::string{"--burnin"}) ) - { - burnin = std::stoi(argv[++na]); - if (na+1==argc) break; - ++na; - } - else if ( 0 == std::string{argv[na]}.compare(std::string{"--nChains"}) ) - { - nChains = std::stoi(argv[++na]); // use the next - if (na+1==argc) break; // in case it's last, break - ++na; // otherwise augment counter - } - else if ( 0 == std::string{argv[na]}.compare(std::string{"--dataFile"}) ) - { - dataFile = ""+std::string(argv[++na]); // use the next - if (na+1==argc) break; // in case it's last, break - ++na; // otherwise augment counter - } - else if ( 0 == std::string{argv[na]}.compare(std::string{"--blockFile"}) ) - { - blockFile = ""+std::string(argv[++na]); // use the next - if (na+1==argc) break; // in case it's last, break - ++na; // otherwise augment counter - } - else if ( 0 == std::string{argv[na]}.compare(std::string{"--structureGraphFile"}) ) - { - structureGraphFile = ""+std::string(argv[++na]); // use the next - if (na+1==argc) break; // in case it's last, break - ++na; // otherwise augment counter - } - else if ( 0 == std::string{argv[na]}.compare(std::string{"--parFile"}) || 0 == std::string{argv[na]}.compare(std::string{"--hyperparFile"}) ) - { - hpFile = ""+std::string(argv[++na]); // use the next - if (na+1==argc) break; // in case it's last, break - ++na; // otherwise augment counter - } - else if ( 0 == std::string{argv[na]}.compare(std::string{"--outFilePath"}) ) - { - outFilePath = std::string(argv[++na]); // use the next - if (na+1==argc) break; // in case it's last, break - ++na; // otherwise augment counter - } - else if ( 0 == std::string{argv[na]}.compare(std::string{"--gammaOut"}) ) // From here set outputs - { - out_gamma = true; - if (na+1==argc) break; - ++na; - } - else if ( 0 == std::string{argv[na]}.compare(std::string{"--NOGammaOut"}) ) - { - out_gamma = false; - if (na+1==argc) break; - ++na; - } - else if ( 0 == std::string{argv[na]}.compare(std::string{"--betaOut"}) ) - { - out_beta = true; - if (na+1==argc) break; - ++na; - } - else if ( 0 == std::string{argv[na]}.compare(std::string{"--NOBetaOut"}) ) - { - out_beta = false; - if (na+1==argc) break; - ++na; - } - else if ( 0 == std::string{argv[na]}.compare(std::string{"--GOut"}) ) - { - out_G = true; - if (na+1==argc) break; - ++na; - } - else if ( 0 == std::string{argv[na]}.compare(std::string{"--NOGOut"}) ) - { - out_G = false; - if (na+1==argc) break; - ++na; - } - else if ( 0 == std::string{argv[na]}.compare(std::string{"--sigmaRhoOut"}) ) - { - out_sigmaRho = true; - if (na+1==argc) break; - ++na; - } - else if ( 0 == std::string{argv[na]}.compare(std::string{"--NOSigmaRhoOut"}) ) - { - out_sigmaRho = false; - if (na+1==argc) break; - ++na; - } - else if ( 0 == std::string{argv[na]}.compare(std::string{"--piOut"}) ) - { - out_pi = true; - if (na+1==argc) break; - ++na; - } - else if ( 0 == std::string{argv[na]}.compare(std::string{"--NOPiOut"}) ) - { - out_pi = false; - if (na+1==argc) break; - ++na; - } - else if ( 0 == std::string{argv[na]}.compare(std::string{"--tailOut"}) ) - { - out_tail = true; - if (na+1==argc) break; - ++na; - } - else if ( 0 == std::string{argv[na]}.compare(std::string{"--NOTailOut"}) ) - { - out_tail = false; - if (na+1==argc) break; - ++na; - } - else if ( 0 == std::string{argv[na]}.compare(std::string{"--modelSizeOut"}) ) - { - out_model_size = true; - if (na+1==argc) break; - ++na; - } - else if ( 0 == std::string{argv[na]}.compare(std::string{"--NOModelSizeOut"}) ) - { - out_model_size = false; - if (na+1==argc) break; - ++na; - } - else if ( 0 == std::string{argv[na]}.compare(std::string{"--CPOOut"}) ) - { - out_CPO = true; - if (na+1==argc) break; - ++na; - } - else if ( 0 == std::string{argv[na]}.compare(std::string{"--NOCPOOut"}) ) - { - out_CPO = false; - if (na+1==argc) break; - ++na; - } - else if ( 0 == std::string{argv[na]}.compare(std::string{"--ModelVisitOut"}) ) - { - out_model_visit = true; - if (na+1==argc) break; - ++na; - } - else if ( 0 == std::string{argv[na]}.compare(std::string{"--NOModelVisitOut"}) ) - { - out_model_visit = false; - if (na+1==argc) break; - ++na; - } - else - { - std::cout << "Unknown option: " << argv[na] << std::endl; - return(1); - } - }//end reading from command line - - if( gammaPrior == "" ) - { - if ( mrfGFile == "" ) - { - std::cout << "Using default prior for Gamma - hotspot prior" << std::endl; - gammaPrior = "hotspot"; - } - else - { - std::cout << "No value for gammaPrior was specified, but mrfG was given - choosing MRF prior" << std::endl; - gammaPrior = "MRF"; - } - - } - - int status{1}; - - try - { - status = drive(dataFile,mrfGFile,blockFile,structureGraphFile,hpFile,outFilePath, - nIter,burnin,nChains, - covariancePrior,gammaPrior,gammaSampler,gammaInit,betaPrior, - out_gamma,out_beta,out_G,out_sigmaRho,out_pi,out_tail,out_model_size,out_CPO,out_model_visit); - } - catch(const std::exception& e) - { - std::cerr << e.what() << std::endl; - } - - return status; - -} diff --git a/BayesSUR/man/BayesSUR.Rd b/man/BayesSUR.Rd similarity index 100% rename from BayesSUR/man/BayesSUR.Rd rename to man/BayesSUR.Rd diff --git a/BayesSUR/man/BayesSUR_internal.Rd b/man/BayesSUR_internal.Rd similarity index 100% rename from BayesSUR/man/BayesSUR_internal.Rd rename to man/BayesSUR_internal.Rd diff --git a/BayesSUR/man/coef.BayesSUR.Rd b/man/coef.BayesSUR.Rd similarity index 100% rename from BayesSUR/man/coef.BayesSUR.Rd rename to man/coef.BayesSUR.Rd diff --git a/BayesSUR/man/elpd.Rd b/man/elpd.Rd similarity index 100% rename from BayesSUR/man/elpd.Rd rename to man/elpd.Rd diff --git a/BayesSUR/man/exampleEQTL.Rd b/man/exampleEQTL.Rd similarity index 100% rename from BayesSUR/man/exampleEQTL.Rd rename to man/exampleEQTL.Rd diff --git a/BayesSUR/man/exampleGDSC.Rd b/man/exampleGDSC.Rd similarity index 100% rename from BayesSUR/man/exampleGDSC.Rd rename to man/exampleGDSC.Rd diff --git a/BayesSUR/man/figures/figure2.png b/man/figures/figure2.png similarity index 100% rename from BayesSUR/man/figures/figure2.png rename to man/figures/figure2.png diff --git a/BayesSUR/man/fitted.BayesSUR.Rd b/man/fitted.BayesSUR.Rd similarity index 100% rename from BayesSUR/man/fitted.BayesSUR.Rd rename to man/fitted.BayesSUR.Rd diff --git a/BayesSUR/man/getEstimator.Rd b/man/getEstimator.Rd similarity index 100% rename from BayesSUR/man/getEstimator.Rd rename to man/getEstimator.Rd diff --git a/BayesSUR/man/plot.BayesSUR.Rd b/man/plot.BayesSUR.Rd similarity index 100% rename from BayesSUR/man/plot.BayesSUR.Rd rename to man/plot.BayesSUR.Rd diff --git a/BayesSUR/man/plotCPO.Rd b/man/plotCPO.Rd similarity index 100% rename from BayesSUR/man/plotCPO.Rd rename to man/plotCPO.Rd diff --git a/BayesSUR/man/plotEstimator.Rd b/man/plotEstimator.Rd similarity index 100% rename from BayesSUR/man/plotEstimator.Rd rename to man/plotEstimator.Rd diff --git a/BayesSUR/man/plotGraph.Rd b/man/plotGraph.Rd similarity index 100% rename from BayesSUR/man/plotGraph.Rd rename to man/plotGraph.Rd diff --git a/BayesSUR/man/plotMCMCdiag.Rd b/man/plotMCMCdiag.Rd similarity index 100% rename from BayesSUR/man/plotMCMCdiag.Rd rename to man/plotMCMCdiag.Rd diff --git a/BayesSUR/man/plotManhattan.Rd b/man/plotManhattan.Rd similarity index 100% rename from BayesSUR/man/plotManhattan.Rd rename to man/plotManhattan.Rd diff --git a/BayesSUR/man/plotNetwork.Rd b/man/plotNetwork.Rd similarity index 100% rename from BayesSUR/man/plotNetwork.Rd rename to man/plotNetwork.Rd diff --git a/BayesSUR/man/predict.BayesSUR.Rd b/man/predict.BayesSUR.Rd similarity index 100% rename from BayesSUR/man/predict.BayesSUR.Rd rename to man/predict.BayesSUR.Rd diff --git a/BayesSUR/man/print.BayesSUR.Rd b/man/print.BayesSUR.Rd similarity index 100% rename from BayesSUR/man/print.BayesSUR.Rd rename to man/print.BayesSUR.Rd diff --git a/BayesSUR/man/summary.BayesSUR.Rd b/man/summary.BayesSUR.Rd similarity index 100% rename from BayesSUR/man/summary.BayesSUR.Rd rename to man/summary.BayesSUR.Rd diff --git a/BayesSUR/man/targetGene.Rd b/man/targetGene.Rd similarity index 100% rename from BayesSUR/man/targetGene.Rd rename to man/targetGene.Rd diff --git a/BayesSUR/src/BayesSUR.cpp b/src/BayesSUR.cpp similarity index 100% rename from BayesSUR/src/BayesSUR.cpp rename to src/BayesSUR.cpp diff --git a/BayesSUR/src/ESS_Atom.h b/src/ESS_Atom.h similarity index 100% rename from BayesSUR/src/ESS_Atom.h rename to src/ESS_Atom.h diff --git a/BayesSUR/src/ESS_Sampler.h b/src/ESS_Sampler.h similarity index 100% rename from BayesSUR/src/ESS_Sampler.h rename to src/ESS_Sampler.h diff --git a/BayesSUR/src/HRR_Chain.cpp b/src/HRR_Chain.cpp similarity index 100% rename from BayesSUR/src/HRR_Chain.cpp rename to src/HRR_Chain.cpp diff --git a/BayesSUR/src/HRR_Chain.h b/src/HRR_Chain.h similarity index 100% rename from BayesSUR/src/HRR_Chain.h rename to src/HRR_Chain.h diff --git a/BayesSUR/src/Makevars.in b/src/Makevars.in similarity index 100% rename from BayesSUR/src/Makevars.in rename to src/Makevars.in diff --git a/BayesSUR/src/Makevars.win b/src/Makevars.win similarity index 100% rename from BayesSUR/src/Makevars.win rename to src/Makevars.win diff --git a/BayesSUR/src/Parameter_types.h b/src/Parameter_types.h similarity index 100% rename from BayesSUR/src/Parameter_types.h rename to src/Parameter_types.h diff --git a/BayesSUR/src/RcppExports.cpp b/src/RcppExports.cpp similarity index 100% rename from BayesSUR/src/RcppExports.cpp rename to src/RcppExports.cpp diff --git a/BayesSUR/src/SUR_Chain.cpp b/src/SUR_Chain.cpp similarity index 100% rename from BayesSUR/src/SUR_Chain.cpp rename to src/SUR_Chain.cpp diff --git a/BayesSUR/src/SUR_Chain.h b/src/SUR_Chain.h similarity index 100% rename from BayesSUR/src/SUR_Chain.h rename to src/SUR_Chain.h diff --git a/BayesSUR/src/distr.cpp b/src/distr.cpp similarity index 100% rename from BayesSUR/src/distr.cpp rename to src/distr.cpp diff --git a/BayesSUR/src/distr.h b/src/distr.h similarity index 100% rename from BayesSUR/src/distr.h rename to src/distr.h diff --git a/BayesSUR/src/drive.cpp b/src/drive.cpp similarity index 100% rename from BayesSUR/src/drive.cpp rename to src/drive.cpp diff --git a/BayesSUR/src/drive.h b/src/drive.h similarity index 100% rename from BayesSUR/src/drive.h rename to src/drive.h diff --git a/BayesSUR/src/global.cpp b/src/global.cpp similarity index 100% rename from BayesSUR/src/global.cpp rename to src/global.cpp diff --git a/BayesSUR/src/global.h b/src/global.h similarity index 100% rename from BayesSUR/src/global.h rename to src/global.h diff --git a/BayesSUR/src/junction_tree.cpp b/src/junction_tree.cpp similarity index 100% rename from BayesSUR/src/junction_tree.cpp rename to src/junction_tree.cpp diff --git a/BayesSUR/src/junction_tree.h b/src/junction_tree.h similarity index 100% rename from BayesSUR/src/junction_tree.h rename to src/junction_tree.h diff --git a/BayesSUR/src/pugiconfig.hpp b/src/pugiconfig.hpp similarity index 100% rename from BayesSUR/src/pugiconfig.hpp rename to src/pugiconfig.hpp diff --git a/BayesSUR/src/pugixml.cpp b/src/pugixml.cpp similarity index 100% rename from BayesSUR/src/pugixml.cpp rename to src/pugixml.cpp diff --git a/BayesSUR/src/pugixml.hpp b/src/pugixml.hpp similarity index 100% rename from BayesSUR/src/pugixml.hpp rename to src/pugixml.hpp diff --git a/BayesSUR/src/utils.cpp b/src/utils.cpp similarity index 100% rename from BayesSUR/src/utils.cpp rename to src/utils.cpp diff --git a/BayesSUR/src/utils.h b/src/utils.h similarity index 100% rename from BayesSUR/src/utils.h rename to src/utils.h diff --git a/test.R b/test.R deleted file mode 100644 index 01496475..00000000 --- a/test.R +++ /dev/null @@ -1,52 +0,0 @@ -### Build a new version of the package -#remove.packages("BayesSUR") -#Rcpp::compileAttributes(pkgdir="BayesSUR/") -#tools::resaveRdaFiles("BayesSUR/data/example_GDSC.rda", compress="xz") -#devtools::document("BayesSUR") -#devtools::build("BayesSUR", vignettes=TRUE, args="--compact-vignettes=both") - -# $ R CMD check ./BayesSUR_1.0-2.tar.gz --as-cran -# devtools::check("BayesSUR", cran=TRUE, build_args="--compact-vignettes=both") -# $ ./BVS_Reg --covariancePrior HIW --gammaPrior MRF --dataFile data_test_GDSC.txt --blockFile blockLabels.txt -# --structureGraphFile structureGraph.txt --outFilePath results/ --nChains 2 --nIter 50 --burnin 0 --mrfGFile mrfG.txt - - -## Install the package -library(devtools) -install_github("mbant/BayesSUR/BayesSUR") -#install.packages("BayesSUR_1.0-2.tar.gz",repos = NULL,type = "source",build_vignettes = TRUE) - - -##################################################################################################### -## load the package -library(BayesSUR) -data(example_eQTL, package = "BayesSUR") -str(example_eQTL) - -# fit a SSUR model with hotspot prior -attach(example_eQTL) -fit <- BayesSUR(data = data, Y = blockList[[1]], - X = blockList[[2]], outFilePath = "results/", - nIter = 2000, nChains = 5, covariancePrior = "HIW", burnin=1000, - gammaPrior = "hotspot") - -summary(fit) - -# show the interaction of plots -plot(fit) - -# show the estimated beta, gamma and G_y -plotEstimator(fit, fig.tex=TRUE) -system(paste(getOption("pdfviewer"), "ParamEstimator.pdf")) - -# show the relationship of responses -plotResponseGraph(fit) - -# show the network representation of the associations between responses and features -plotNetwork(fit) - -# show the manhattan plot -plotManhattan(fit) - -# check the convergence of the algorithm -plotMCMCdiag(fit) diff --git a/BayesSUR/vignettes/BayesSUR-RE.Rmd b/vignettes/BayesSUR-RE.Rmd similarity index 100% rename from BayesSUR/vignettes/BayesSUR-RE.Rmd rename to vignettes/BayesSUR-RE.Rmd diff --git a/BayesSUR/vignettes/BayesSUR.pdf b/vignettes/BayesSUR.pdf similarity index 100% rename from BayesSUR/vignettes/BayesSUR.pdf rename to vignettes/BayesSUR.pdf diff --git a/BayesSUR/vignettes/BayesSUR.pdf.asis b/vignettes/BayesSUR.pdf.asis similarity index 100% rename from BayesSUR/vignettes/BayesSUR.pdf.asis rename to vignettes/BayesSUR.pdf.asis