library(graper)
library(ggplot2)
Note that the implementation of logistic regression is still experimental in some parts.
Create an example data set with 4 groups, 100 train + 100 test samples and 320 features.
set.seed(123)
data <- makeExampleData(n=200, p=320, g=4,
pis=c(0.05, 0.1, 0.05, 0.1),
gammas=c(0.1, 0.1, 10, 10),
response="bernoulli")
# training data set
Xtrain <- data$X[1:100, ]
ytrain <- data$y[1:100]
# annotations of features to groups
annot <- data$annot
# test data set
Xtest <- data$X[101:200, ]
ytest <- data$y[101:200]
graper
is the main function of this package,
which allows to fit the proposed Bayesian models
with different settings on the prior (by setting spikeslab
to FALSE or TRUE)
and the variational approximation (by setting factoriseQ
to FALSE or TRUE).
By default, the model is fit with a sparsity promoting spike-and-slab prior
and a fully-factorised mean-field assumption.
The ELBO is currently not calculated in the logisitc regession framework.
fit <- graper(Xtrain, ytrain, annot, verbose=FALSE,
family="binomial", calcELB=FALSE)
## Fitting a model with 4 groups, 100 samples and 320 features.
## Fitting with random init 1
## Maximum numbers of iterations reached - no convergence or ELB not calculated
fit
## Sparse graper object for a logistic regression model with 320 predictors in 4 groups.
## Group-wise shrinkage:
## 1 2 3 4
## 170.93 4.17 297.4 329.72
## Group-wise sparsity (1 = dense, 0 = sparse):
## 1 2 3 4
## 0.35 0.1 0.29 0.25
The variational Bayes (VB) approach directly yields posterior distributions for each parameter. Note, however, that using VB these are often too concentrated and cannot be directly used for construction of confidence intervals etc. However, they can provide good point estimates.
plotPosterior(fit, "gamma", gamma0=data$gammas)
plotPosterior(fit, "pi", pi0=data$pis)
The estimated coefficients, their posterior inclusion probabilities and the intercept are contained in the result list.
# get coefficients (without the intercept)
beta <- coef(fit, include_intercept=FALSE)
# beta <- fit$EW_beta
# plot estimated versus true beta
qplot(beta, data$beta)
# get intercept
intercept <- fit$intercept
# get estimated posterior inclusion probabilities per feature
pips <- getPIPs(fit)
The function predict
can be used to make prediction on new data.
Here, we illustrate its use by predicting the response
on the test data defined above.
preds <- predict(fit, Xtest)
#SessionInfo
sessionInfo()
## R version 3.6.0 (2019-04-26)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 18.04.2 LTS
##
## Matrix products: default
## BLAS: /home/biocbuild/bbs-3.9-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.9-bioc/R/lib/libRlapack.so
##
## locale:
## [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=en_US.UTF-8 LC_COLLATE=C
## [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
## [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggplot2_3.1.1 graper_1.0.0 BiocStyle_2.12.0
##
## loaded via a namespace (and not attached):
## [1] Rcpp_1.0.1 pillar_1.3.1 compiler_3.6.0
## [4] BiocManager_1.30.4 plyr_1.8.4 tools_3.6.0
## [7] digest_0.6.18 evaluate_0.13 tibble_2.1.1
## [10] gtable_0.3.0 lattice_0.20-38 pkgconfig_2.0.2
## [13] rlang_0.3.4 Matrix_1.2-17 yaml_2.2.0
## [16] xfun_0.6 withr_2.1.2 stringr_1.4.0
## [19] dplyr_0.8.0.1 knitr_1.22 grid_3.6.0
## [22] tidyselect_0.2.5 cowplot_0.9.4 glue_1.3.1
## [25] R6_2.4.0 rmarkdown_1.12 bookdown_0.9
## [28] purrr_0.3.2 magrittr_1.5 scales_1.0.0
## [31] htmltools_0.3.6 matrixStats_0.54.0 assertthat_0.2.1
## [34] colorspace_1.4-1 labeling_0.3 stringi_1.4.3
## [37] lazyeval_0.2.2 munsell_0.5.0 crayon_1.3.4