stat.glmnet_coefdiff: Importance statistics based on a GLM with cross-validation

Description

Fits a generalized linear model via penalized maximum likelihood and cross-validation. Then, compute the difference statistic $$W_j = |Z_j| - |\tilde{Z}_j|$$ where $Z_j$ and $\tilde{Z}_j$ are the coefficient estimates for the jth variable and its knockoff, respectively. The value of the regularization parameter $\lambda$ is selected by cross-validation and computed with glmnet.

Usage

stat.glmnet_coefdiff(X, X_k, y, family = "gaussian", cores = 2, ...)

Value

A vector of statistics $W$ of length p.

Arguments

X: n-by-p matrix of original variables.
X_k: n-by-p matrix of knockoff variables.
y: vector of length n, containing the response variables. Quantitative for family="gaussian", or family="poisson" (non-negative counts). For family="binomial" should be either a factor with two levels, or a two-column matrix of counts or proportions (the second column is treated as the target class; for a factor, the last level in alphabetical order is the target class). For family="multinomial", can be a nc>=2 level factor, or a matrix with nc columns of counts or proportions. For either "binomial" or "multinomial", if y is presented as a vector, it will be coerced into a factor. For family="cox", y should be a two-column matrix with columns named 'time' and 'status'. The latter is a binary variable, with '1' indicating death, and '0' indicating right censored. The function Surv() in package survival produces such a matrix. For family="mgaussian", y is a matrix of quantitative responses.
family: response type (see above).
cores: Number of cores used to compute the statistics by running cv.glmnet. Unless otherwise specified, the number of cores is set equal to two (if available).
...: additional arguments specific to glmnet (see Details).

Details

This function uses the glmnet package to fit a generalized linear model via penalized maximum likelihood.

The statistics $W_j$ are constructed by taking the difference between the coefficient of the j-th variable and its knockoff.

By default, the value of the regularization parameter is chosen by 10-fold cross-validation.

The default response family is 'gaussian', for a linear regression model. Different response families (e.g. 'binomial') can be specified by passing an optional parameter 'family'.

The optional nlambda parameter can be used to control the granularity of the grid of $\lambda$'s. The default value of nlambda is 500, where p is the number of columns of X.

If the family is 'binomial' and a lambda sequence is not provided by the user, this function generates it on a log-linear scale before calling 'glmnet'.

For a complete list of the available additional arguments, see cv.glmnet and glmnet.

Examples

Run this code

set.seed(2022)
p=200; n=100; k=15
mu = rep(0,p); Sigma = diag(p)
X = matrix(rnorm(n*p),n)
nonzero = sample(p, k)
beta = 3.5 * (1:p %in% nonzero)
y = X %*% beta + rnorm(n)
knockoffs = function(X) create.gaussian(X, mu, Sigma)

# Basic usage with default arguments
result = knockoff.filter(X, y, knockoffs=knockoffs, 
                           statistic=stat.glmnet_coefdiff)
print(result$selected)

# Advanced usage with custom arguments
foo = stat.glmnet_coefdiff
k_stat = function(X, X_k, y) foo(X, X_k, y, nlambda=200)
result = knockoff.filter(X, y, knockoffs=knockoffs, statistic=k_stat)
print(result$selected)