| Title: | Adaptive penalization in high-dimensional regression and classification with external covariates using variational Bayes |
|---|---|
| Description: | This package enables regression and classification on high-dimensional data with different relative strengths of penalization for different feature groups, such as different assays or omic types. The optimal relative strengths are chosen adaptively. Optimisation is performed using a variational Bayes approach. |
| Authors: | Britta Velten [aut, cre], Wolfgang Huber [aut] |
| Maintainer: | Britta Velten <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 1.21.0 |
| Built: | 2024-09-16 05:18:13 UTC |
| Source: | https://github.com/bioc/graper |
Function to obtain estimated coefficients from a fitted graper model.
## S3 method for class 'graper' coef(object, include_intercept = TRUE, ...)## S3 method for class 'graper' coef(object, include_intercept = TRUE, ...)
object |
fitted graper model as obtained from |
include_intercept |
whether to include the estimated intercept value in the output |
... |
other arguments |
1-Column matrix of estimated coefficients.
# create data dat <- makeExampleData() # fit the graper model to the data fit <- graper(dat$X, dat$y, dat$annot) # extract the model coefficients coef(fit)# create data dat <- makeExampleData() # fit the graper model to the data fit <- graper(dat$X, dat$y, dat$annot) # extract the model coefficients coef(fit)
Function to obtain estimated posterior inclusion probabilities per feature from a fitted graper model.
getPIPs(object)getPIPs(object)
object |
fitted graper model as obtained from |
1-Column matrix of estimated posterior inclusion probabilities.
# create data dat <- makeExampleData() # fit the graper model to the data fit <- graper(dat$X, dat$y, dat$annot) # extract the posterior inclusion probabilities from the fitted model getPIPs(fit)# create data dat <- makeExampleData() # fit the graper model to the data fit <- graper(dat$X, dat$y, dat$annot) # extract the posterior inclusion probabilities from the fitted model getPIPs(fit)
Fit a regression model with graper given a matrix
of predictors (X), a response vector (y) and
a vector of group memberships for each predictor
in X (annot).
For each group a different strength of penalization
is determined adaptively.
graper(X, y, annot, factoriseQ = TRUE, spikeslab = TRUE, intercept = TRUE, family = "gaussian", standardize = TRUE, n_rep = 1, max_iter = 3000, th = 0.01, d_tau = 0.001, r_tau = 0.001, d_gamma = 0.001, r_gamma = 0.001, r_pi = 1, d_pi = 1, calcELB = TRUE, verbose = TRUE, freqELB = 1, nogamma = FALSE, init_psi = 1)graper(X, y, annot, factoriseQ = TRUE, spikeslab = TRUE, intercept = TRUE, family = "gaussian", standardize = TRUE, n_rep = 1, max_iter = 3000, th = 0.01, d_tau = 0.001, r_tau = 0.001, d_gamma = 0.001, r_gamma = 0.001, r_pi = 1, d_pi = 1, calcELB = TRUE, verbose = TRUE, freqELB = 1, nogamma = FALSE, init_psi = 1)
X |
design matrix of size n (samples) x p (features) |
y |
response vector of size n |
annot |
factor of length p indicating group membership of each feature (column) in X |
factoriseQ |
if set to TRUE, the variational distribution is assumed to fully factorize across features (faster, default). If FALSE, a multivariate variational distribution is used. |
spikeslab |
if set to TRUE, a spike and slab prior on the coefficients (default). |
intercept |
whether to include an intercept into the model |
family |
Likelihood model for the response, either "gaussian" for linear regression or "binomial" for logistic regression |
standardize |
whether to standardize the predictors to unit variance |
n_rep |
number of repetitions with different random initializations to be fit |
max_iter |
maximum number of iterations |
th |
convergence threshold for the evidence lower bound (ELB) |
d_tau |
hyper-parameters for prior of tau (noise precision) |
r_tau |
hyper-parameters for prior of tau (noise precision) |
d_gamma |
hyper-parameters for prior of gamma (coefficients' prior precision) |
r_gamma |
hyper-parameters for prior of gamma (coefficients' prior precision) |
r_pi |
hyper-parameters for Beta prior of the mixture probabilities in the spike and slab prior |
d_pi |
hyper-parameters for Beta prior of the mixture probabilities in the spike and slab prior |
calcELB |
whether to calculate the evidence lower bound (ELB) |
verbose |
whether to print out intermediate messages during fitting |
freqELB |
frequency at which the evidence lower bound (ELB) is to be calculated, i.e. each freqELB-th iteration |
nogamma |
if TRUE, the normal prior will have same variance for all groups (only relevant for spikeslab = TRUE) |
init_psi |
initial value for the spike variables |
The function trains the graper model given
a matrix of predictors (X), a response vector (y)
and a vector of group memberships for each predictor in X
(annot). For each feature group as specified in
annot a penalty factor and sparsity level is learnt.
By default it uses a Spike-and-Slab prior on the coefficients
and uses a fully factorized variational distribution
in the inference. This provides a fast way to train the model.
Using spikeslab=FALSE a ridge regression like model can
be fitted using a normal instead of the spike and slab prior.
Setting factoriseQ = FALSE gives a more exact inference
scheme based on a multivariate variational distribution,
but can be much slower.
As the optimization is non-convex is can be helpful to
use multiple random initializations by setting n_rep
to a value larger 1. The returned model is then chosen
as the optimal fit with respect to the evidence lower bound (ELB).
Depending on the response vector a linear regression model
(family = "gaussian") or a logistic regression model
(family = "binomial") is fitted.
Note, that the implementation of logistic regression is still
experimental.
A graper object containing
estimated model coefficients in liner/logistic regression
estimated posterior-inclusion probabilities for each feature
estimated intercept term
annotation vector of features to the groups as
specified when calling graper
estimated penalty factor per group
estimated sparsity level per group (from 1 (dense) to 0 (sparse))
estimated noise precision
parameters of the variational distributions of beta, gamma, tau and pi
final value of the evidence lower bound
values of the evidence lower bound for all iterations
other options used when calling graper
# create data dat <- makeExampleData() # fit a sparse model with spike and slab prior fit <- graper(dat$X, dat$y, dat$annot) fit # print fitted object beta <- coef(fit, include_intercept=FALSE) # model coeffients pips <- getPIPs(fit) # posterior inclusion probabilities pf <- fit$EW_gamma # penalty factors per group sparsities <- fit$EW_pi # sparsity levels per group # fit a dense model without spike and slab prior fit <- graper(dat$X, dat$y, dat$annot, spikeslab=FALSE) # fit a dense model using a multivariate variational distribution fit <- graper(dat$X, dat$y, dat$annot, factoriseQ=TRUE, spikeslab=FALSE)# create data dat <- makeExampleData() # fit a sparse model with spike and slab prior fit <- graper(dat$X, dat$y, dat$annot) fit # print fitted object beta <- coef(fit, include_intercept=FALSE) # model coeffients pips <- getPIPs(fit) # posterior inclusion probabilities pf <- fit$EW_gamma # penalty factors per group sparsities <- fit$EW_pi # sparsity levels per group # fit a dense model without spike and slab prior fit <- graper(dat$X, dat$y, dat$annot, spikeslab=FALSE) # fit a dense model using a multivariate variational distribution fit <- graper(dat$X, dat$y, dat$annot, factoriseQ=TRUE, spikeslab=FALSE)
Simulate data from the graper model with groups of equal size and pre-specified parameters gamma, pi and tau.
makeExampleData(n = 100, p = 200, g = 4, gammas = c(0.1, 1, 10, 100), pis = c(0.5, 0.5, 0.5, 0.5), tau = 1, rho = 0, response = "gaussian", intercept = 0)makeExampleData(n = 100, p = 200, g = 4, gammas = c(0.1, 1, 10, 100), pis = c(0.5, 0.5, 0.5, 0.5), tau = 1, rho = 0, response = "gaussian", intercept = 0)
n |
number of samples |
p |
number of features |
g |
number of groups |
gammas |
vector of length g, specifying the slab precision of the prior on beta per group |
pis |
vector of length g, specifying the probability of s to be 1 (slab) |
tau |
noise precision |
rho |
correlation of design matrix (Toeplitz structure) |
response |
"gaussian" for continuous response from a linear regression model, "bernoulli" for a binary response from a logistic regression model. |
intercept |
model intercept (default: 0) |
list containing the design matrix X,
the response y, the feature annotation to
groups annot as well as the different parameters
in the Bayesian model and the correlation strength rho
dat <- makeExampleData()dat <- makeExampleData()
Simulate data from the graper model with groups of unequal size and pre-specified parameters gamma, pi and tau.
makeExampleDataWithUnequalGroups(n = 100, pg = c(100, 100, 10, 10), gammas = c(0.1, 10, 0.1, 10), pis = c(0.5, 0.5, 0.5, 0.5), tau = 1, rho = 0, response = "gaussian", intercept = 0)makeExampleDataWithUnequalGroups(n = 100, pg = c(100, 100, 10, 10), gammas = c(0.1, 10, 0.1, 10), pis = c(0.5, 0.5, 0.5, 0.5), tau = 1, rho = 0, response = "gaussian", intercept = 0)
n |
number of samples |
pg |
vector of length g (desired number of groups) with number of features per group |
gammas |
vector of length g, specifying the slab precision of the prior on beta per group |
pis |
vector of length g, specifying the probability of s to be 1 (slab) |
tau |
noise precision (only relevant for gaussian response) |
rho |
correlation of design matrix (Toeplitz structure) |
response |
"gaussian" for continuous response from a linear regression model, "bernoulli" for a binary response from a logistic regression model. |
intercept |
model intercept (default: 0) |
list containin the design matrix X,
the response y,
the feature annotation to groups annot as well as
the different parameters in the Bayesian model
and the correlation strength rho
dat <- makeExampleDataWithUnequalGroups()dat <- makeExampleDataWithUnequalGroups()
Function to plot the evidence lower bound (ELBO) over iterations to monitor the convergence of the algorithm.
plotELBO(fit)plotELBO(fit)
fit |
fit as produced by |
a ggplot object
dat <- makeExampleData() fit <- graper(dat$X, dat$y, dat$annot) plotELBO(fit)dat <- makeExampleData() fit <- graper(dat$X, dat$y, dat$annot) plotELBO(fit)
Function to plot the group-wise penalty factors (gamma) and sparsity levels.
plotGroupPenalties(fit)plotGroupPenalties(fit)
fit |
fit as produced by |
a ggplot object
dat <- makeExampleData() fit <- graper(dat$X, dat$y, dat$annot) plotGroupPenalties(fit)dat <- makeExampleData() fit <- graper(dat$X, dat$y, dat$annot) plotGroupPenalties(fit)
Function to plot the posterior of the model parameters obtained by graper from the variational inference framework.
plotPosterior(fit, param2plot, beta0 = NULL, gamma0 = NULL, tau0 = NULL, pi0 = NULL, s0 = NULL, jmax = 2, range = NULL)plotPosterior(fit, param2plot, beta0 = NULL, gamma0 = NULL, tau0 = NULL, pi0 = NULL, s0 = NULL, jmax = 2, range = NULL)
fit |
fit as produced by |
param2plot |
which parameter of the graper model to plot (gamma, beta, tau or s) |
beta0 |
true beta (if known) |
gamma0 |
true gamma (if known) |
tau0 |
true tau (if known) |
pi0 |
true pi (if known) |
s0 |
true s (if known) |
jmax |
maximal number of components per group to plot (for beta and s) |
range |
plotting range (x-axis) |
a ggplot object
# create data dat <- makeExampleData() # fit the graper model fit <- graper(dat$X, dat$y, dat$annot) # plot posterior distribution of the gamma parameter plotPosterior(fit, param2plot="gamma")# create data dat <- makeExampleData() # fit the graper model fit <- graper(dat$X, dat$y, dat$annot) # plot posterior distribution of the gamma parameter plotPosterior(fit, param2plot="gamma")
Function to predict the response on a new data set using a fitted graper model.
## S3 method for class 'graper' predict(object, newX, type = c("inRange", "response", "link"), ...)## S3 method for class 'graper' predict(object, newX, type = c("inRange", "response", "link"), ...)
object |
fitted graper model as obtained from |
newX |
Predictor matrix of size n_test
(number of new test samples) x p (number of predictors)
(same feature structure as used in |
type |
type of prediction returned, either:
|
... |
other arguments |
A vector with predictions.
# create data dat <- makeExampleData() # split data into train and test sets of equal size ntrain <- dat$n / 2 # fit the model to the train data fit <- graper(dat$X[seq_len(ntrain), ], dat$y[seq_len(ntrain)], dat$annot) # make predictions on the test data ypred <- predict(fit, dat$X[seq_len(ntrain) + dat$n / 2, ]) # create data for logistic regression dat <- makeExampleData(response="bernoulli") # split data into train and test sets of equal size ntrain <- dat$n / 2 # fit the graper model for a logistic model fit <- graper(dat$X[seq_len(ntrain), ], dat$y[seq_len(ntrain)], dat$annot, family="binomial") # make predictions on the test data ypred <- predict(fit, dat$X[seq_len(ntrain) + dat$n / 2, ], type = "inRange")# create data dat <- makeExampleData() # split data into train and test sets of equal size ntrain <- dat$n / 2 # fit the model to the train data fit <- graper(dat$X[seq_len(ntrain), ], dat$y[seq_len(ntrain)], dat$annot) # make predictions on the test data ypred <- predict(fit, dat$X[seq_len(ntrain) + dat$n / 2, ]) # create data for logistic regression dat <- makeExampleData(response="bernoulli") # split data into train and test sets of equal size ntrain <- dat$n / 2 # fit the graper model for a logistic model fit <- graper(dat$X[seq_len(ntrain), ], dat$y[seq_len(ntrain)], dat$annot, family="binomial") # make predictions on the test data ypred <- predict(fit, dat$X[seq_len(ntrain) + dat$n / 2, ], type = "inRange")
Function to print a fitted graper model.
## S3 method for class 'graper' print(x, ...)## S3 method for class 'graper' print(x, ...)
x |
fitted graper model as obtained from |
... |
additional print arguments |
Print output.
# create data dat <- makeExampleData() # fit the graper model fit <- graper(dat$X, dat$y, dat$annot) # print a summary of the fitted model print(fit)# create data dat <- makeExampleData() # fit the graper model fit <- graper(dat$X, dat$y, dat$annot) # print a summary of the fitted model print(fit)