| Title: | Variable Selection for Heterogeneous Populations |
|---|---|
| Description: | Provides variable selection and estimation routines for models with main effects stratified on multiple binary factors. The 'vennLasso' package is an implementation of the method introduced in Huling, et al. (2017) <doi:10.1111/biom.12769>. |
| Authors: | Jared Huling [aut, cre]
|
| Maintainer: | Jared Huling <[email protected]> |
| License: | GPL (>= 2) |
| Version: | 0.1.6 |
| Built: | 2026-11-14 10:30:36 UTC |
| Source: | https://github.com/cran/vennLasso |
Cross Validation for the vennLasso
cv.vennLasso( x, y, groups, lambda = NULL, compute.se = FALSE, conf.int = NULL, type.measure = c("mse", "deviance", "class", "auc", "mae", "brier"), nfolds = 10, foldid, grouped = TRUE, keep = FALSE, parallel = FALSE, ... )cv.vennLasso( x, y, groups, lambda = NULL, compute.se = FALSE, conf.int = NULL, type.measure = c("mse", "deviance", "class", "auc", "mae", "brier"), nfolds = 10, foldid, grouped = TRUE, keep = FALSE, parallel = FALSE, ... )
x |
input matrix or SparseMatrix of dimension nobs x nvars. Each row is an observation, each column corresponds to a covariate |
y |
numeric response vector of length nobs |
groups |
A list of length equal to the number of groups containing vectors of integers indicating the variable IDs for each group. For example, groups=list(c(1,2), c(2,3), c(3,4,5)) specifies that Group 1 contains variables 1 and 2, Group 2 contains variables 2 and 3, and Group 3 contains variables 3, 4, and 5. Can also be a matrix of 0s and 1s with the number of columns equal to the number of groups and the number of rows equal to the number of variables. A value of 1 in row i and column j indicates that variable i is in group j and 0 indicates that variable i is not in group j. |
lambda |
A user-specified sequence of lambda values. Left unspecified, the a sequence of lambda values is automatically computed, ranging uniformly on the log scale over the relevant range of lambda values. |
compute.se |
logical flag. If |
conf.int |
value between 0 and 1 indicating the level of the confidence intervals to be computed. For example
if |
type.measure |
One of |
nfolds |
number of folds for cross-validation. default is 10. 3 is smallest value allowed. |
foldid |
an optional vector of values between 1 and nfold specifying which fold each observation belongs to. |
grouped |
Like in glmnet, this is an experimental argument, with default |
keep |
If |
parallel |
If TRUE, use parallel foreach to fit each fold. Must register parallel before hand, such as doMC. |
... |
parameters to be passed to vennLasso |
An object with S3 class "cv.vennLasso"
library(Matrix) set.seed(123) n.obs <- 150 n.vars <- 25 true.beta.mat <- array(NA, dim = c(3, n.vars)) true.beta.mat[1,] <- c(-0.5, -1, 0, 0, 2, rep(0, n.vars - 5)) true.beta.mat[2,] <- c(0.5, 0.5, -0.5, -0.5, 1, -1, rep(0, n.vars - 6)) true.beta.mat[3,] <- c(0, 0, 1, 1, -1, rep(0, n.vars - 5)) rownames(true.beta.mat) <- c("1,0", "1,1", "0,1") true.beta <- as.vector(t(true.beta.mat)) x.sub1 <- matrix(rnorm(n.obs * n.vars), n.obs, n.vars) x.sub2 <- matrix(rnorm(n.obs * n.vars), n.obs, n.vars) x.sub3 <- matrix(rnorm(n.obs * n.vars), n.obs, n.vars) x <- as.matrix(rbind(x.sub1, x.sub2, x.sub3)) conditions <- as.matrix(cbind(c(rep(1, 2 * n.obs), rep(0, n.obs)), c(rep(0, n.obs), rep(1, 2 * n.obs)))) y <- rnorm(n.obs * 3, sd = 3) + drop(as.matrix(bdiag(x.sub1, x.sub2, x.sub3)) %*% true.beta) fit <- cv.vennLasso(x = x, y = y, groups = conditions, nfolds = 3) fitted.coef <- predict(fit$vennLasso.fit, type = "coefficients", s = fit$lambda.min) (true.coef <- true.beta.mat[match(dimnames(fit$vennLasso.fit$beta)[[1]], rownames(true.beta.mat)),]) round(fitted.coef, 2) ## effects smaller for logistic regression ## Not run: true.beta.mat <- true.beta.mat / 2 true.beta <- true.beta / 2 # logistic regression example#' y <- rbinom(n.obs * 3, 1, prob = 1 / (1 + exp(-drop(as.matrix(bdiag(x.sub1, x.sub2, x.sub3)) %*% true.beta)))) bfit <- cv.vennLasso(x = x, y = y, groups = conditions, family = "binomial", nfolds = 3) fitted.coef <- predict(bfit$vennLasso.fit, type = "coefficients", s = bfit$lambda.min) (true.coef <- true.beta.mat[match(dimnames(bfit$vennLasso.fit$beta)[[1]], rownames(true.beta.mat)),]) round(fitted.coef, 2) ## End(Not run)library(Matrix) set.seed(123) n.obs <- 150 n.vars <- 25 true.beta.mat <- array(NA, dim = c(3, n.vars)) true.beta.mat[1,] <- c(-0.5, -1, 0, 0, 2, rep(0, n.vars - 5)) true.beta.mat[2,] <- c(0.5, 0.5, -0.5, -0.5, 1, -1, rep(0, n.vars - 6)) true.beta.mat[3,] <- c(0, 0, 1, 1, -1, rep(0, n.vars - 5)) rownames(true.beta.mat) <- c("1,0", "1,1", "0,1") true.beta <- as.vector(t(true.beta.mat)) x.sub1 <- matrix(rnorm(n.obs * n.vars), n.obs, n.vars) x.sub2 <- matrix(rnorm(n.obs * n.vars), n.obs, n.vars) x.sub3 <- matrix(rnorm(n.obs * n.vars), n.obs, n.vars) x <- as.matrix(rbind(x.sub1, x.sub2, x.sub3)) conditions <- as.matrix(cbind(c(rep(1, 2 * n.obs), rep(0, n.obs)), c(rep(0, n.obs), rep(1, 2 * n.obs)))) y <- rnorm(n.obs * 3, sd = 3) + drop(as.matrix(bdiag(x.sub1, x.sub2, x.sub3)) %*% true.beta) fit <- cv.vennLasso(x = x, y = y, groups = conditions, nfolds = 3) fitted.coef <- predict(fit$vennLasso.fit, type = "coefficients", s = fit$lambda.min) (true.coef <- true.beta.mat[match(dimnames(fit$vennLasso.fit$beta)[[1]], rownames(true.beta.mat)),]) round(fitted.coef, 2) ## effects smaller for logistic regression ## Not run: true.beta.mat <- true.beta.mat / 2 true.beta <- true.beta / 2 # logistic regression example#' y <- rbinom(n.obs * 3, 1, prob = 1 / (1 + exp(-drop(as.matrix(bdiag(x.sub1, x.sub2, x.sub3)) %*% true.beta)))) bfit <- cv.vennLasso(x = x, y = y, groups = conditions, family = "binomial", nfolds = 3) fitted.coef <- predict(bfit$vennLasso.fit, type = "coefficients", s = bfit$lambda.min) (true.coef <- true.beta.mat[match(dimnames(bfit$vennLasso.fit$beta)[[1]], rownames(true.beta.mat)),]) round(fitted.coef, 2) ## End(Not run)
function to estimate the hierarchical sparsity parameter for a desired level of sparsity for simulated hierarchical coefficients
estimate.hier.sparsity.param( ncats, nvars, avg.hier.zeros = 0.3, nsims = 150, effect.size.max = 0.5, misspecification.prop = 0 )estimate.hier.sparsity.param( ncats, nvars, avg.hier.zeros = 0.3, nsims = 150, effect.size.max = 0.5, misspecification.prop = 0 )
ncats |
number of categories to stratify on |
nvars |
number of variables |
avg.hier.zeros |
desired percent of zero variables among the variables with hierarchical zero patterns. |
nsims |
number of simulations to estimate the average sparsity. A larger number will be more accurate but take much longer. |
effect.size.max |
maximum magnitude of the true effect sizes |
misspecification.prop |
proportion of variables with hierarchical missingness misspecified |
set.seed(123) # estimate hier.sparsity.param for 0.15 total proportion of nonzero variables # among vars with hierarchical zero patterns ## Not run: hsp <- estimate.hier.sparsity.param(ncats = 3, nvars = 25, avg.hier.zeros = 0.15, nsims = 100) ## End(Not run) # the above results in the following value hsp <- 0.6341772 # check that this does indeed achieve the desired level of sparsity mean(replicate(100, mean(genHierSparseBeta(ncats = 3, nvars = 25, hier.sparsity.param = hsp) != 0) )) sparseBeta <- genHierSparseBeta(ncats = 3, nvars = 25, hier.sparsity.param = hsp) ## Not run: hsp2 <- estimate.hier.sparsity.param(ncats = 2, nvars = 100, avg.hier.zeros = 0.30, nsims = 50) # 0.5778425 hsp3 <- estimate.hier.sparsity.param(ncats = 3, nvars = 100, avg.hier.zeros = 0.30, nsims = 50) # 0.4336312 hsp4 <- estimate.hier.sparsity.param(ncats = 4, nvars = 100, avg.hier.zeros = 0.30, nsims = 50) # 0.2670061 hsp5 <- estimate.hier.sparsity.param(ncats = 5, nvars = 100, avg.hier.zeros = 0.30, nsims = 50) # 0.146682 ## End(Not run) # 0.07551241 for hsp6set.seed(123) # estimate hier.sparsity.param for 0.15 total proportion of nonzero variables # among vars with hierarchical zero patterns ## Not run: hsp <- estimate.hier.sparsity.param(ncats = 3, nvars = 25, avg.hier.zeros = 0.15, nsims = 100) ## End(Not run) # the above results in the following value hsp <- 0.6341772 # check that this does indeed achieve the desired level of sparsity mean(replicate(100, mean(genHierSparseBeta(ncats = 3, nvars = 25, hier.sparsity.param = hsp) != 0) )) sparseBeta <- genHierSparseBeta(ncats = 3, nvars = 25, hier.sparsity.param = hsp) ## Not run: hsp2 <- estimate.hier.sparsity.param(ncats = 2, nvars = 100, avg.hier.zeros = 0.30, nsims = 50) # 0.5778425 hsp3 <- estimate.hier.sparsity.param(ncats = 3, nvars = 100, avg.hier.zeros = 0.30, nsims = 50) # 0.4336312 hsp4 <- estimate.hier.sparsity.param(ncats = 4, nvars = 100, avg.hier.zeros = 0.30, nsims = 50) # 0.2670061 hsp5 <- estimate.hier.sparsity.param(ncats = 5, nvars = 100, avg.hier.zeros = 0.30, nsims = 50) # 0.146682 ## End(Not run) # 0.07551241 for hsp6
function to generate coefficient matrix with hierarchical sparsity
genHierSparseBeta( ncats, nvars, hier.sparsity.param = 0.5, avg.hier.zeros = NULL, effect.size.max = 0.5, misspecification.prop = 0 )genHierSparseBeta( ncats, nvars, hier.sparsity.param = 0.5, avg.hier.zeros = NULL, effect.size.max = 0.5, misspecification.prop = 0 )
ncats |
number of categories to stratify on |
nvars |
number of variables |
hier.sparsity.param |
parameter between 0 and 1 which determines how much hierarchical sparsity there is. To achieve a desired total level of sparsity among the variables with hierarchical sparsity, this parameter can be estimated using the function 'estimate.hier.sparsity.param' |
avg.hier.zeros |
desired percent of zero variables among the variables with hierarchical zero patterns. If this is specified, it will override the given hier.sparsity.param value and estimate it. This takes a while |
effect.size.max |
maximum magnitude of the true effect sizes |
misspecification.prop |
proportion of variables with hierarchical missingness misspecified |
set.seed(123) # estimate hier.sparsity.param for 0.15 total proportion of nonzero variables # among vars with hierarchical zero patterns # NOT RUN: Takes a long time # hsp <- estimate.hier.sparsity.param(ncats = 3, nvars = 25, avg.hier.zeros = 0.15, nsims = 100) # the above results in the following value hsp <- 0.6341772 # check that this does indeed achieve the desired level of sparsity mean(replicate(100, mean(genHierSparseBeta(ncats = 3, nvars = 25, hier.sparsity.param = hsp) != 0) )) sparseBeta <- genHierSparseBeta(ncats = 3, nvars = 25, hier.sparsity.param = hsp)set.seed(123) # estimate hier.sparsity.param for 0.15 total proportion of nonzero variables # among vars with hierarchical zero patterns # NOT RUN: Takes a long time # hsp <- estimate.hier.sparsity.param(ncats = 3, nvars = 25, avg.hier.zeros = 0.15, nsims = 100) # the above results in the following value hsp <- 0.6341772 # check that this does indeed achieve the desired level of sparsity mean(replicate(100, mean(genHierSparseBeta(ncats = 3, nvars = 25, hier.sparsity.param = hsp) != 0) )) sparseBeta <- genHierSparseBeta(ncats = 3, nvars = 25, hier.sparsity.param = hsp)
function to generate data with hierarchical sparsity
genHierSparseData( ncats, nvars, nobs, nobs.test = 100, hier.sparsity.param = 0.5, avg.hier.zeros = NULL, prop.zero.vars = 0.5, effect.size.max = 0.5, misspecification.prop = 0, family = c("gaussian", "binomial", "coxph"), sd = 1, snr = NULL, beta = NULL, tau = 10, covar = 0 )genHierSparseData( ncats, nvars, nobs, nobs.test = 100, hier.sparsity.param = 0.5, avg.hier.zeros = NULL, prop.zero.vars = 0.5, effect.size.max = 0.5, misspecification.prop = 0, family = c("gaussian", "binomial", "coxph"), sd = 1, snr = NULL, beta = NULL, tau = 10, covar = 0 )
ncats |
number of categories to stratify on |
nvars |
number of variables |
nobs |
number of observations per strata to simulate |
nobs.test |
number of independent test observations per strata to simulate |
hier.sparsity.param |
parameter between 0 and 1 which determines how much hierarchical sparsity there is. To achieve a desired total level of sparsity among the variables with hierarchical sparsity, this parameter can be estimated using the function 'estimate.hier.sparsity.param' |
avg.hier.zeros |
desired percent of zero variables among the variables with hierarchical zero patterns. If this is specified, it will override the given hier.sparsity.param value and estimate it. This takes a while |
prop.zero.vars |
proportion of all variables that will be zero across all strata |
effect.size.max |
maximum magnitude of the true effect sizes |
misspecification.prop |
proportion of variables with hierarchical missingness misspecified |
family |
family for the response variable |
sd |
standard devation for gaussian simulations |
snr |
signal-to-noise ratio (only used for |
beta |
a matrix of true beta values. If given, then no beta will be created and data will be simulated from the given beta |
tau |
rate parameter for |
covar |
scalar, pairwise covariance term for covariates |
set.seed(123) dat.sim <- genHierSparseData(ncats = 3, nvars = 100, nobs = 200) # estimate hier.sparsity.param for 0.15 total proportion of nonzero variables # among vars with hierarchical zero patterns ## Not run: hsp <- estimate.hier.sparsity.param(ncats = 3, nvars = 50, avg.hier.zeros = 0.15, nsims = 100) ## End(Not run) # the above results in the following value hsp <- 0.6270698 # check that this does indeed achieve the desired level of sparsity mean(replicate(50, mean(genHierSparseBeta(ncats = 3, nvars = 50, hier.sparsity.param = hsp) != 0) )) dat.sim2 <- genHierSparseData(ncats = 3, nvars = 100, nobs = 200, hier.sparsity.param = hsp) sparseBeta <- genHierSparseBeta(ncats = 3, nvars = 100, hier.sparsity.param = hsp) ## generate data with already generated beta dat.sim3 <- genHierSparseData(ncats = 3, nvars = 100, nobs = 200, beta = sparseBeta) ## complete example: ## 50% sparsity: hsp <- 0.2626451 dat.sim <- genHierSparseData(ncats = 3, nvars = 25, nobs = 150, nobs.test = 1000, hier.sparsity.param = hsp, prop.zero.vars = 0.5, effect.size.max = 0.25, family = "gaussian") x <- dat.sim$x x.test <- dat.sim$x.test y <- dat.sim$y y.test <- dat.sim$y.test grp <- dat.sim$group.ind grp.test <- dat.sim$group.ind.test fit.adapt <- cv.vennLasso(x, y, grp, adaptive.lasso = TRUE, nlambda = 25, family = "gaussian", abs.tol = 1e-5, rel.tol = 1e-5, maxit = 1000, irls.maxit = 15L, gamma = 0.2, standardize = FALSE, intercept = TRUE, nfolds = 3, model.matrix = TRUE) preds.a <- predict(fit.adapt$vennLasso.fit, x.test, grp.test, s = fit.adapt$lambda.min, type = 'response')set.seed(123) dat.sim <- genHierSparseData(ncats = 3, nvars = 100, nobs = 200) # estimate hier.sparsity.param for 0.15 total proportion of nonzero variables # among vars with hierarchical zero patterns ## Not run: hsp <- estimate.hier.sparsity.param(ncats = 3, nvars = 50, avg.hier.zeros = 0.15, nsims = 100) ## End(Not run) # the above results in the following value hsp <- 0.6270698 # check that this does indeed achieve the desired level of sparsity mean(replicate(50, mean(genHierSparseBeta(ncats = 3, nvars = 50, hier.sparsity.param = hsp) != 0) )) dat.sim2 <- genHierSparseData(ncats = 3, nvars = 100, nobs = 200, hier.sparsity.param = hsp) sparseBeta <- genHierSparseBeta(ncats = 3, nvars = 100, hier.sparsity.param = hsp) ## generate data with already generated beta dat.sim3 <- genHierSparseData(ncats = 3, nvars = 100, nobs = 200, beta = sparseBeta) ## complete example: ## 50% sparsity: hsp <- 0.2626451 dat.sim <- genHierSparseData(ncats = 3, nvars = 25, nobs = 150, nobs.test = 1000, hier.sparsity.param = hsp, prop.zero.vars = 0.5, effect.size.max = 0.25, family = "gaussian") x <- dat.sim$x x.test <- dat.sim$x.test y <- dat.sim$y y.test <- dat.sim$y.test grp <- dat.sim$group.ind grp.test <- dat.sim$group.ind.test fit.adapt <- cv.vennLasso(x, y, grp, adaptive.lasso = TRUE, nlambda = 25, family = "gaussian", abs.tol = 1e-5, rel.tol = 1e-5, maxit = 1000, irls.maxit = 15L, gamma = 0.2, standardize = FALSE, intercept = TRUE, nfolds = 3, model.matrix = TRUE) preds.a <- predict(fit.adapt$vennLasso.fit, x.test, grp.test, s = fit.adapt$lambda.min, type = 'response')
log likelihood function for fitted vennLasso objects
## S3 method for class 'vennLasso' logLik(object, ...)## S3 method for class 'vennLasso' logLik(object, ...)
object |
fitted "vennLasso" model object. |
... |
not used |
library(Matrix) set.seed(123) n.obs <- 200 n.vars <- 50 true.beta.mat <- array(NA, dim = c(3, n.vars)) true.beta.mat[1,] <- c(-0.5, -1, 0, 0, 2, rep(0, n.vars - 5)) true.beta.mat[2,] <- c(0.5, 0.5, -0.5, -0.5, 1, -1, rep(0, n.vars - 6)) true.beta.mat[3,] <- c(0, 0, 1, 1, -1, rep(0, n.vars - 5)) rownames(true.beta.mat) <- c("1,0", "1,1", "0,1") true.beta <- as.vector(t(true.beta.mat)) x.sub1 <- matrix(rnorm(n.obs * n.vars), n.obs, n.vars) x.sub2 <- matrix(rnorm(n.obs * n.vars), n.obs, n.vars) x.sub3 <- matrix(rnorm(n.obs * n.vars), n.obs, n.vars) x <- as.matrix(rbind(x.sub1, x.sub2, x.sub3)) conditions <- as.matrix(cbind(c(rep(1, 2 * n.obs), rep(0, n.obs)), c(rep(0, n.obs), rep(1, 2 * n.obs)))) y <- rnorm(n.obs * 3, sd = 3) + drop(as.matrix(bdiag(x.sub1, x.sub2, x.sub3)) %*% true.beta) fit <- vennLasso(x = x, y = y, groups = conditions) logLik(fit)library(Matrix) set.seed(123) n.obs <- 200 n.vars <- 50 true.beta.mat <- array(NA, dim = c(3, n.vars)) true.beta.mat[1,] <- c(-0.5, -1, 0, 0, 2, rep(0, n.vars - 5)) true.beta.mat[2,] <- c(0.5, 0.5, -0.5, -0.5, 1, -1, rep(0, n.vars - 6)) true.beta.mat[3,] <- c(0, 0, 1, 1, -1, rep(0, n.vars - 5)) rownames(true.beta.mat) <- c("1,0", "1,1", "0,1") true.beta <- as.vector(t(true.beta.mat)) x.sub1 <- matrix(rnorm(n.obs * n.vars), n.obs, n.vars) x.sub2 <- matrix(rnorm(n.obs * n.vars), n.obs, n.vars) x.sub3 <- matrix(rnorm(n.obs * n.vars), n.obs, n.vars) x <- as.matrix(rbind(x.sub1, x.sub2, x.sub3)) conditions <- as.matrix(cbind(c(rep(1, 2 * n.obs), rep(0, n.obs)), c(rep(0, n.obs), rep(1, 2 * n.obs)))) y <- rnorm(n.obs * 3, sd = 3) + drop(as.matrix(bdiag(x.sub1, x.sub2, x.sub3)) %*% true.beta) fit <- vennLasso(x = x, y = y, groups = conditions) logLik(fit)
Overlapping Group Lasso (OGLasso)
oglasso( x, y, delta = NULL, group, fused = NULL, family = c("gaussian", "binomial", "coxph"), nlambda = 100L, lambda = NULL, lambda.min.ratio = NULL, lambda.fused = 0, alpha = NULL, group.weights = NULL, adaptive.lasso = FALSE, adaptive.fused = FALSE, penalty.factor = NULL, penalty.factor.fused = NULL, gamma = 1, standardize = TRUE, intercept = TRUE, compute.se = FALSE, rho = NULL, dynamic.rho = TRUE, maxit = 500L, abs.tol = 1e-05, rel.tol = 1e-05, irls.tol = 1e-05, irls.maxit = 100L )oglasso( x, y, delta = NULL, group, fused = NULL, family = c("gaussian", "binomial", "coxph"), nlambda = 100L, lambda = NULL, lambda.min.ratio = NULL, lambda.fused = 0, alpha = NULL, group.weights = NULL, adaptive.lasso = FALSE, adaptive.fused = FALSE, penalty.factor = NULL, penalty.factor.fused = NULL, gamma = 1, standardize = TRUE, intercept = TRUE, compute.se = FALSE, rho = NULL, dynamic.rho = TRUE, maxit = 500L, abs.tol = 1e-05, rel.tol = 1e-05, irls.tol = 1e-05, irls.maxit = 100L )
x |
input matrix of dimension nobs by nvars. Each row is an observation, each column corresponds to a covariate |
y |
numeric response vector of length nobs |
delta |
vector of length equal to the number of observations with values in 1 and 0, where a 1 indicates the observed time is a death and a 0 indicates the observed time is a censoring event |
group |
A list of length equal to the number of groups containing vectors of integers
indicating the variable IDs for each group. For example, |
fused |
matrix specifying generalized lasso penalty formulation. Each column corresponds to each variable and each row corresponds to a new penalty term, ie if row 1 has the first entry of 1 and the second entry of -1, then the penalty term lambda.fused * |beta_1 - beta_2| will be added. Not available now |
family |
"gaussian" for least squares problems, "binomial" for binary response |
nlambda |
The number of lambda values. Default is 100. |
lambda |
A user-specified sequence of lambda values. Left unspecified, the a sequence of lambda values is automatically computed, ranging uniformly on the log scale over the relevant range of lambda values. |
lambda.min.ratio |
Smallest value for lambda, as a fraction of lambda.max, the (data derived) entry
value (i.e. the smallest value for which all parameter estimates are zero). The default
depends on the sample size nobs relative to the number of variables nvars. If
|
lambda.fused |
tuning parameter for fused (generalized) lasso penalty |
alpha |
currently not used. Will be used later for fused lasso |
group.weights |
A vector of values representing multiplicative factors by which each group's penalty is to be multiplied. Often, this is a function (such as the square root) of the number of predictors in each group. The default is to use the square root of group size for the group selection methods. |
adaptive.lasso |
Flag indicating whether or not to use adaptive lasso weights. If set to |
adaptive.fused |
Flag indicating whether or not to use adaptive fused lasso weights. |
penalty.factor |
vector of weights to be multiplied to the tuning parameter for the group lasso penalty. A vector of length equal to the number of groups |
penalty.factor.fused |
vector of weights to be multiplied to the tuning parameter for the fused lasso penalty. A vector of length equal to the number of variables. mostly for internal usage |
gamma |
power to raise the MLE estimated weights by for the adaptive lasso. defaults to 1 |
standardize |
Logical flag for |
intercept |
Should intercept(s) be fitted ( |
compute.se |
Should standard errors be computed? If |
rho |
ADMM parameter. must be a strictly positive value. By default, an appropriate value is automatically chosen |
dynamic.rho |
|
maxit |
integer. Maximum number of ADMM iterations. Default is 500. |
abs.tol |
absolute convergence tolerance for ADMM iterations for the relative dual and primal residuals.
Default is |
rel.tol |
relative convergence tolerance for ADMM iterations for the relative dual and primal residuals.
Default is |
irls.tol |
convergence tolerance for IRLS iterations. Only used if |
irls.maxit |
integer. Maximum number of IRLS iterations. Only used if |
An object with S3 class "oglasso"
library(vennLasso) set.seed(123) n.obs <- 1e3 n.vars <- 50 true.beta <- c(rep(0,2), 1, -1, rep(0, 8), 0.5, -0.5, 1, rep(0, 35)) x <- matrix(rnorm(n.obs * n.vars), n.obs, n.vars) y <- rnorm(n.obs, sd = 3) + drop(x %*% true.beta) groups <- c(list(c(1,2), c(2,3), c(3,4,5), 5:10, 6:12, 7:15), lapply(16:50, function(x) x)) ## Not run: fit <- oglasso(x = x, y = y, group = groups) ## End(Not run)library(vennLasso) set.seed(123) n.obs <- 1e3 n.vars <- 50 true.beta <- c(rep(0,2), 1, -1, rep(0, 8), 0.5, -0.5, 1, rep(0, 35)) x <- matrix(rnorm(n.obs * n.vars), n.obs, n.vars) y <- rnorm(n.obs, sd = 3) + drop(x %*% true.beta) groups <- c(list(c(1,2), c(2,3), c(3,4,5), 5:10, 6:12, 7:15), lapply(16:50, function(x) x)) ## Not run: fit <- oglasso(x = x, y = y, group = groups) ## End(Not run)
Plot method for cv.vennLasso fitted objects
Plotting method for vennLasso fitted objects
## S3 method for class 'cv.vennLasso' plot(x, sign.lambda = 1, ...) ## S3 method for class 'vennLasso' plot( x, which.subpop = 1, xvar = c("norm", "lambda", "loglambda", "dev"), xlab = iname, ylab = "Coefficients", ... )## S3 method for class 'cv.vennLasso' plot(x, sign.lambda = 1, ...) ## S3 method for class 'vennLasso' plot( x, which.subpop = 1, xvar = c("norm", "lambda", "loglambda", "dev"), xlab = iname, ylab = "Coefficients", ... )
x |
fitted |
sign.lambda |
Either plot against log(lambda) (default) or its negative if |
... |
other graphical parameters for the plot |
which.subpop |
which row in the coefficient matrix should be plotting? Each row corresponds to a particular combination of the specified stratifying variables |
xvar |
What is on the X-axis. "norm" plots against the L1-norm of the coefficients, "lambda" against the log-lambda sequence, and "dev" against the percent deviance explained. |
xlab |
character value supplied for x-axis label |
ylab |
character value supplied for y-axis label |
set.seed(123) dat.sim <- genHierSparseData(ncats = 3, nvars = 25, nobs = 100, hier.sparsity.param = 0.5, prop.zero.vars = 0.5, effect.size.max = 0.25, family = "gaussian") x <- dat.sim$x x.test <- dat.sim$x.test y <- dat.sim$y y.test <- dat.sim$y.test grp <- dat.sim$group.ind grp.test <- dat.sim$group.ind.test fit.adapt <- cv.vennLasso(x, y, grp, adaptive.lasso = TRUE, nlambda = 25, nfolds = 4) plot(fit.adapt) library(Matrix) set.seed(123) n.obs <- 200 n.vars <- 50 true.beta.mat <- array(NA, dim = c(3, n.vars)) true.beta.mat[1,] <- c(-0.5, -1, 0, 0, 2, rep(0, n.vars - 5)) true.beta.mat[2,] <- c(0.5, 0.5, -0.5, -0.5, 1, -1, rep(0, n.vars - 6)) true.beta.mat[3,] <- c(0, 0, 1, 1, -1, rep(0, n.vars - 5)) rownames(true.beta.mat) <- c("1,0", "1,1", "0,1") true.beta <- as.vector(t(true.beta.mat)) x.sub1 <- matrix(rnorm(n.obs * n.vars), n.obs, n.vars) x.sub2 <- matrix(rnorm(n.obs * n.vars), n.obs, n.vars) x.sub3 <- matrix(rnorm(n.obs * n.vars), n.obs, n.vars) x <- as.matrix(rbind(x.sub1, x.sub2, x.sub3)) conditions <- as.matrix(cbind(c(rep(1, 2 * n.obs), rep(0, n.obs)), c(rep(0, n.obs), rep(1, 2 * n.obs)))) y <- rnorm(n.obs * 3, sd = 3) + drop(as.matrix(bdiag(x.sub1, x.sub2, x.sub3)) %*% true.beta) fit <- vennLasso(x = x, y = y, groups = conditions) layout(matrix(1:3, ncol = 3)) plot(fit, which.subpop = 1) plot(fit, which.subpop = 2) plot(fit, which.subpop = 3)set.seed(123) dat.sim <- genHierSparseData(ncats = 3, nvars = 25, nobs = 100, hier.sparsity.param = 0.5, prop.zero.vars = 0.5, effect.size.max = 0.25, family = "gaussian") x <- dat.sim$x x.test <- dat.sim$x.test y <- dat.sim$y y.test <- dat.sim$y.test grp <- dat.sim$group.ind grp.test <- dat.sim$group.ind.test fit.adapt <- cv.vennLasso(x, y, grp, adaptive.lasso = TRUE, nlambda = 25, nfolds = 4) plot(fit.adapt) library(Matrix) set.seed(123) n.obs <- 200 n.vars <- 50 true.beta.mat <- array(NA, dim = c(3, n.vars)) true.beta.mat[1,] <- c(-0.5, -1, 0, 0, 2, rep(0, n.vars - 5)) true.beta.mat[2,] <- c(0.5, 0.5, -0.5, -0.5, 1, -1, rep(0, n.vars - 6)) true.beta.mat[3,] <- c(0, 0, 1, 1, -1, rep(0, n.vars - 5)) rownames(true.beta.mat) <- c("1,0", "1,1", "0,1") true.beta <- as.vector(t(true.beta.mat)) x.sub1 <- matrix(rnorm(n.obs * n.vars), n.obs, n.vars) x.sub2 <- matrix(rnorm(n.obs * n.vars), n.obs, n.vars) x.sub3 <- matrix(rnorm(n.obs * n.vars), n.obs, n.vars) x <- as.matrix(rbind(x.sub1, x.sub2, x.sub3)) conditions <- as.matrix(cbind(c(rep(1, 2 * n.obs), rep(0, n.obs)), c(rep(0, n.obs), rep(1, 2 * n.obs)))) y <- rnorm(n.obs * 3, sd = 3) + drop(as.matrix(bdiag(x.sub1, x.sub2, x.sub3)) %*% true.beta) fit <- vennLasso(x = x, y = y, groups = conditions) layout(matrix(1:3, ncol = 3)) plot(fit, which.subpop = 1) plot(fit, which.subpop = 2) plot(fit, which.subpop = 3)
plotting function to investigate estimated coefficients
plotCoefs(object, s = NULL, ...)plotCoefs(object, s = NULL, ...)
object |
fitted vennLasso object |
s |
lambda value for the predictions. Only one can be specified at a time |
... |
other graphical parameters for the plot |
set.seed(123) dat.sim <- genHierSparseData(ncats = 3, nvars = 25, nobs = 200) fit <- vennLasso(x = dat.sim$x, y = dat.sim$y, groups = dat.sim$group.ind) plotCoefs(fit, s = fit$lambda[22])set.seed(123) dat.sim <- genHierSparseData(ncats = 3, nvars = 25, nobs = 200) fit <- vennLasso(x = dat.sim$x, y = dat.sim$y, groups = dat.sim$group.ind) plotCoefs(fit, s = fit$lambda[22])
plotting function to investigate hierarchical structure of selection
plotSelections(object, s = NULL, type = c("d3.tree"), ...)plotSelections(object, s = NULL, type = c("d3.tree"), ...)
object |
fitted vennLasso object |
s |
lambda value for the predictions. Only one can be specified at a time |
type |
type of plot to make. Currently only "d3.tree" and "igraph.tree" available |
... |
other graphical parameters for the plot |
set.seed(123) dat.sim <- genHierSparseData(ncats = 3, nvars = 25, nobs = 200) fit <- vennLasso(x = dat.sim$x, y = dat.sim$y, groups = dat.sim$group.ind) plotSelections(fit, s = fit$lambda[32])set.seed(123) dat.sim <- genHierSparseData(ncats = 3, nvars = 25, nobs = 200) fit <- vennLasso(x = dat.sim$x, y = dat.sim$y, groups = dat.sim$group.ind) plotSelections(fit, s = fit$lambda[32])
plotting function for venn diagrams of overlapping conditions
plotVenn( conditions, condition.names = NULL, lty = "blank", fill.colors = c("royalblue1", "goldenrod1", "mediumvioletred", "turquoise3", "firebrick1"), ... )plotVenn( conditions, condition.names = NULL, lty = "blank", fill.colors = c("royalblue1", "goldenrod1", "mediumvioletred", "turquoise3", "firebrick1"), ... )
conditions |
condition matrix such as the one given to vennLasso() function. It can have up to 5 conditions |
condition.names |
names of the conditions (equal to the number of columns of conditions) |
lty |
standard 'lty' graphical parameter for line type around circles. default is no lines |
fill.colors |
vector of colors for plotting. Set fill.colors = NULL for no colors |
... |
other graphical parameters for the plot |
library(Matrix) set.seed(123) n.obs <- 200 n.vars <- 50 true.beta.mat <- array(NA, dim = c(3, n.vars)) true.beta.mat[1,] <- c(-0.5, -1, 0, 0, 2, rep(0, n.vars - 5)) true.beta.mat[2,] <- c(0.5, 0.5, -0.5, -0.5, 1, -1, rep(0, n.vars - 6)) true.beta.mat[3,] <- c(0, 0, 1, 1, -1, rep(0, n.vars - 5)) rownames(true.beta.mat) <- c("1,0", "1,1", "0,1") true.beta <- as.vector(t(true.beta.mat)) x.sub1 <- matrix(rnorm(n.obs * n.vars), n.obs, n.vars) x.sub2 <- matrix(rnorm(n.obs * n.vars), n.obs, n.vars) x.sub3 <- matrix(rnorm(n.obs * n.vars), n.obs, n.vars) x <- as.matrix(rbind(x.sub1, x.sub2, x.sub3)) conditions <- as.matrix(cbind(c(rep(1, 2 * n.obs), rep(0, n.obs)), c(rep(0, n.obs), rep(1, 2 * n.obs)))) y <- rnorm(n.obs * 3, sd = 3) + drop(as.matrix(bdiag(x.sub1, x.sub2, x.sub3)) %*% true.beta) fit <- vennLasso(x = x, y = y, groups = conditions) vennobj <- plotVenn(conditions)library(Matrix) set.seed(123) n.obs <- 200 n.vars <- 50 true.beta.mat <- array(NA, dim = c(3, n.vars)) true.beta.mat[1,] <- c(-0.5, -1, 0, 0, 2, rep(0, n.vars - 5)) true.beta.mat[2,] <- c(0.5, 0.5, -0.5, -0.5, 1, -1, rep(0, n.vars - 6)) true.beta.mat[3,] <- c(0, 0, 1, 1, -1, rep(0, n.vars - 5)) rownames(true.beta.mat) <- c("1,0", "1,1", "0,1") true.beta <- as.vector(t(true.beta.mat)) x.sub1 <- matrix(rnorm(n.obs * n.vars), n.obs, n.vars) x.sub2 <- matrix(rnorm(n.obs * n.vars), n.obs, n.vars) x.sub3 <- matrix(rnorm(n.obs * n.vars), n.obs, n.vars) x <- as.matrix(rbind(x.sub1, x.sub2, x.sub3)) conditions <- as.matrix(cbind(c(rep(1, 2 * n.obs), rep(0, n.obs)), c(rep(0, n.obs), rep(1, 2 * n.obs)))) y <- rnorm(n.obs * 3, sd = 3) + drop(as.matrix(bdiag(x.sub1, x.sub2, x.sub3)) %*% true.beta) fit <- vennLasso(x = x, y = y, groups = conditions) vennobj <- plotVenn(conditions)
Prediction for Cross Validation Hierarchical Lasso Object
## S3 method for class 'cv.vennLasso' predict(object, newx, group.mat, s = c("lambda.min"), use.refit = FALSE, ...)## S3 method for class 'cv.vennLasso' predict(object, newx, group.mat, s = c("lambda.min"), use.refit = FALSE, ...)
object |
fitted cv.vennLasso object |
newx |
new matrix for predictions |
group.mat |
A matrix of the group memberships for now. Ignore the rest: A list of length equal to the number of groups containing vectors of integers indicating the variable IDs for each group. For example, groups=list(c(1,2), c(2,3), c(3,4,5)) specifies that Group 1 contains variables 1 and 2, Group 2 contains variables 2 and 3, and Group 3 contains variables 3, 4, and 5. Can also be a matrix of 0s and 1s with the number of columns equal to the number of groups and the number of rows equal to the number of variables. A value of 1 in row i and column j indicates that variable i is in group j and 0 indicates that variable i is not in group j. |
s |
lambda value for the predictions. defaults to all values computed in the vennLasso object |
use.refit |
Should the refitted beta estimates be used for prediction? Defaults to FALSE. If TRUE then the beta estimates from the model refit on just the selected covariates are used |
... |
parameters to be passed to predict.vennLasso |
predictions or coefficients
Prediction for Hierarchical Shared Lasso
## S3 method for class 'vennLasso' predict( object, newx, group.mat, s = NULL, use.refit = FALSE, type = c("link", "response", "coefficients", "nonzero", "class", "nvars", "median", "survival"), ... )## S3 method for class 'vennLasso' predict( object, newx, group.mat, s = NULL, use.refit = FALSE, type = c("link", "response", "coefficients", "nonzero", "class", "nvars", "median", "survival"), ... )
object |
fitted vennLasso object |
newx |
new matrix for predictions |
group.mat |
A matrix of the group memberships for now. Ignore the rest: A list of length equal to the number of groups containing vectors of integers indicating the variable IDs for each group. For example, groups=list(c(1,2), c(2,3), c(3,4,5)) specifies that Group 1 contains variables 1 and 2, Group 2 contains variables 2 and 3, and Group 3 contains variables 3, 4, and 5. Can also be a matrix of 0s and 1s with the number of columns equal to the number of groups and the number of rows equal to the number of variables. A value of 1 in row i and column j indicates that variable i is in group j and 0 indicates that variable i is not in group j. |
s |
lambda value for the predictions. defaults to all values computed in the vennLasso object |
use.refit |
Should the refitted beta estimates be used for prediction? Defaults to FALSE. If TRUE then the beta estimates from the model refit on just the selected covariates are used |
type |
type of predictions to be made. |
... |
parameters to be passed to vennLasso |
predictions or coefficients
Fitting vennLasso models
vennLasso( x, y, groups, family = c("gaussian", "binomial"), nlambda = 100L, lambda = NULL, lambda.min.ratio = NULL, lambda.fused = NULL, penalty.factor = NULL, group.weights = NULL, adaptive.lasso = FALSE, adaptive.fused = FALSE, gamma = 1, standardize = FALSE, intercept = TRUE, one.intercept = FALSE, compute.se = FALSE, conf.int = NULL, rho = NULL, dynamic.rho = TRUE, maxit = 500L, abs.tol = 1e-05, rel.tol = 1e-05, irls.tol = 1e-05, irls.maxit = 100L, model.matrix = FALSE, ... )vennLasso( x, y, groups, family = c("gaussian", "binomial"), nlambda = 100L, lambda = NULL, lambda.min.ratio = NULL, lambda.fused = NULL, penalty.factor = NULL, group.weights = NULL, adaptive.lasso = FALSE, adaptive.fused = FALSE, gamma = 1, standardize = FALSE, intercept = TRUE, one.intercept = FALSE, compute.se = FALSE, conf.int = NULL, rho = NULL, dynamic.rho = TRUE, maxit = 500L, abs.tol = 1e-05, rel.tol = 1e-05, irls.tol = 1e-05, irls.maxit = 100L, model.matrix = FALSE, ... )
x |
input matrix of dimension nobs by nvars. Each row is an observation, each column corresponds to a covariate |
y |
numeric response vector of length nobs |
groups |
A list of length equal to the number of groups containing vectors of integers
indicating the variable IDs for each group. For example, |
family |
|
nlambda |
The number of lambda values. Default is 100. |
lambda |
A user-specified sequence of lambda values. Left unspecified, the a sequence of lambda values is automatically computed, ranging uniformly on the log scale over the relevant range of lambda values. |
lambda.min.ratio |
Smallest value for lambda, as a fraction of |
lambda.fused |
tuning parameter for fused lasso penalty |
penalty.factor |
vector of weights to be multiplied to the tuning parameter for the group lasso penalty. A vector of length equal to the number of groups |
group.weights |
A vector of values representing multiplicative factors by which each group's penalty is to be multiplied. Often, this is a function (such as the square root) of the number of predictors in each group. The default is to use the square root of group size for the group selection methods. |
adaptive.lasso |
Flag indicating whether or not to use adaptive lasso weights. If set to |
adaptive.fused |
Flag indicating whether or not to use adaptive fused lasso weights. |
gamma |
power to raise the MLE estimated weights by for the adaptive lasso. defaults to 1 |
standardize |
Should the data be standardized? Defaults to |
intercept |
Should an intercept be fit? Defaults to |
one.intercept |
Should a single intercept be fit for all subpopulations instead of one
for each subpopulation? Defaults to |
compute.se |
Should standard errors be computed? If |
conf.int |
level for confidence intervals. Defaults to |
rho |
ADMM parameter. must be a strictly positive value. By default, an appropriate value is automatically chosen |
dynamic.rho |
|
maxit |
integer. Maximum number of ADMM iterations. Default is 500. |
abs.tol |
absolute convergence tolerance for ADMM iterations for the relative dual and primal residuals.
Default is |
rel.tol |
relative convergence tolerance for ADMM iterations for the relative dual and primal residuals.
Default is |
irls.tol |
convergence tolerance for IRLS iterations. Only used if |
irls.maxit |
integer. Maximum number of IRLS iterations. Only used if |
model.matrix |
logical flag. Should the design matrix used be returned? |
... |
not used |
An object with S3 class "vennLasso"
library(Matrix) # first simulate heterogeneous data using # genHierSparseData set.seed(123) dat.sim <- genHierSparseData(ncats = 2, nvars = 25, nobs = 200, hier.sparsity.param = 0.5, prop.zero.vars = 0.5, family = "gaussian") x <- dat.sim$x conditions <- dat.sim$group.ind y <- dat.sim$y true.beta.mat <- dat.sim$beta.mat fit <- vennLasso(x = x, y = y, groups = conditions) (true.coef <- true.beta.mat[match(dimnames(fit$beta)[[1]], rownames(true.beta.mat)),]) round(fit$beta[,,21], 2) ## fit adaptive version and compute confidence intervals afit <- vennLasso(x = x, y = y, groups = conditions, conf.int = 0.95, adaptive.lasso = TRUE) (true.coef <- true.beta.mat[match(dimnames(fit$beta)[[1]], rownames(true.beta.mat)),])[,1:10] round(afit$beta[,1:10,28], 2) round(afit$lower.ci[,1:10,28], 2) round(afit$upper.ci[,1:10,28], 2) aic.idx <- which.min(afit$aic) bic.idx <- which.min(afit$bic) # actual coverage # actual coverage mean(true.coef[afit$beta[,-1,aic.idx] != 0] >= afit$lower.ci[,-1,aic.idx][afit$beta[,-1,aic.idx] != 0] & true.coef[afit$beta[,-1,aic.idx] != 0] <= afit$upper.ci[,-1,aic.idx][afit$beta[,-1,aic.idx] != 0]) (covered <- true.coef >= afit$lower.ci[,-1,aic.idx] & true.coef <= afit$upper.ci[,-1,aic.idx]) mean(covered) # logistic regression example ## Not run: set.seed(123) dat.sim <- genHierSparseData(ncats = 2, nvars = 25, nobs = 200, hier.sparsity.param = 0.5, prop.zero.vars = 0.5, family = "binomial", effect.size.max = 0.5) # don't make any # coefficients too big x <- dat.sim$x conditions <- dat.sim$group.ind y <- dat.sim$y true.beta.b <- dat.sim$beta.mat bfit <- vennLasso(x = x, y = y, groups = conditions, family = "binomial") (true.coef.b <- -true.beta.b[match(dimnames(fit$beta)[[1]], rownames(true.beta.b)),]) round(bfit$beta[,,20], 2) ## End(Not run)library(Matrix) # first simulate heterogeneous data using # genHierSparseData set.seed(123) dat.sim <- genHierSparseData(ncats = 2, nvars = 25, nobs = 200, hier.sparsity.param = 0.5, prop.zero.vars = 0.5, family = "gaussian") x <- dat.sim$x conditions <- dat.sim$group.ind y <- dat.sim$y true.beta.mat <- dat.sim$beta.mat fit <- vennLasso(x = x, y = y, groups = conditions) (true.coef <- true.beta.mat[match(dimnames(fit$beta)[[1]], rownames(true.beta.mat)),]) round(fit$beta[,,21], 2) ## fit adaptive version and compute confidence intervals afit <- vennLasso(x = x, y = y, groups = conditions, conf.int = 0.95, adaptive.lasso = TRUE) (true.coef <- true.beta.mat[match(dimnames(fit$beta)[[1]], rownames(true.beta.mat)),])[,1:10] round(afit$beta[,1:10,28], 2) round(afit$lower.ci[,1:10,28], 2) round(afit$upper.ci[,1:10,28], 2) aic.idx <- which.min(afit$aic) bic.idx <- which.min(afit$bic) # actual coverage # actual coverage mean(true.coef[afit$beta[,-1,aic.idx] != 0] >= afit$lower.ci[,-1,aic.idx][afit$beta[,-1,aic.idx] != 0] & true.coef[afit$beta[,-1,aic.idx] != 0] <= afit$upper.ci[,-1,aic.idx][afit$beta[,-1,aic.idx] != 0]) (covered <- true.coef >= afit$lower.ci[,-1,aic.idx] & true.coef <= afit$upper.ci[,-1,aic.idx]) mean(covered) # logistic regression example ## Not run: set.seed(123) dat.sim <- genHierSparseData(ncats = 2, nvars = 25, nobs = 200, hier.sparsity.param = 0.5, prop.zero.vars = 0.5, family = "binomial", effect.size.max = 0.5) # don't make any # coefficients too big x <- dat.sim$x conditions <- dat.sim$group.ind y <- dat.sim$y true.beta.b <- dat.sim$beta.mat bfit <- vennLasso(x = x, y = y, groups = conditions, family = "binomial") (true.coef.b <- -true.beta.b[match(dimnames(fit$beta)[[1]], rownames(true.beta.b)),]) round(bfit$beta[,,20], 2) ## End(Not run)