This vignette illustrates the basic and advanced usage of the SLOPE package for R. To learn more about the theory behind SLOPE, see the SLOPE paper (Bogdan et al. 2014).

Synthetic data

For simplicity, we will use synthetic data constructed such that the true coefficient vector \(\beta\) has very few nonzero entries. The details of the construction are unimportant, but are provided below for completeness.

# Problem parameters
n = 500                         # number of observations
p = 500                         # number of variables
k = 10                          # number of variables with nonzero coefficients
amplitude = 1.2*sqrt(2*log(p))  # signal amplitude (for noise level = 1)

# Design matrix
X <- local({
  probs = runif(p, 0.1, 0.5)
  probs = t(probs) %x% matrix(1,n,2)
  X0 = matrix(rbinom(2*n*p, 1, probs), n, 2*p)
  X0 %*% (diag(p) %x% matrix(1,2,1))
})

# Coefficient and response vectors
nonzero = sample(p, k)
beta = amplitude * (1:p %in% nonzero)
y = X %*% beta + rnorm(n)

In our example, the design matrix has entries in \(\{0,1,2\}\).

1 2 3 4 5 6 7 8 9 10
0 2 1 1 0 0 1 0 1 0
2 2 1 1 0 0 0 0 1 1
1 0 1 2 1 0 0 0 1 2
1 1 1 1 1 0 0 0 0 1
2 1 2 0 0 0 0 1 2 2
1 1 2 2 2 1 0 0 1 0
0 0 1 2 0 0 1 0 0 0
2 0 0 1 1 1 0 1 1 0
1 1 0 2 0 2 1 1 1 0
2 2 2 1 0 0 1 1 0 1

Basic usage

The main entry point for the SLOPE package is the function SLOPE. To begin, we call SLOPE with all of the default options.

library(SLOPE)

result <- SLOPE(X, y)
print(result)
## 
## Call:
## SLOPE(X = X, y = y)
## 
## Selected 13 variables with coefficients:
##          91           93          122          185          198  
## 61.94609802  44.86241100  58.71531203   0.81828240  61.94609802  
##         248          331          334          340          352  
## 56.35340793  53.64336513  -0.09645369  59.81488694  59.05946534  
##         374          404          425  
## -0.20392630  52.51754483  53.61387804

The output shows the selected variables and their coefficients. The achieved false discovery proportion (FDP) is close to 0.2, the default false discovery rate (FDR).

fdp <- function(selected) sum(beta[selected] == 0) / max(1, length(selected))
fdp(result$selected)
## [1] 0.2307692

The target FDR can be adjusted using the parameter fdr:

result <- SLOPE(X, y, fdr=0.1)
fdp(result$selected)
## [1] 0.09090909

Advanced usage

Noise level

By default, the noise level is estimated from the data. If \(n\) is large enough compared to \(p\), the classical unbiased estimate of \(\sigma^2\) is used. Otherwise, the iterative SLOPE algorithm is used, as described in Section 3.2.3 of (Bogdan et al. 2014). If the noise level is known, it can be specified directly, bypassing the iterative estimation procedure. For example:

result <- SLOPE(X, y, sigma=1)
fdp(result$selected)
## [1] 0.2307692

Lambda sequence

The primary hyper-parameter in the SLOPE procedure is the decreasing sequence of \(\lambda\) values used in the sorted L1 minimization problem: \[ \min_\beta \left\{ \frac{1}{2} \Vert X\beta - y \Vert_2^2 + \sigma \cdot \sum_{i=1}^p \lambda_i |\beta|_{(i)} \right\} \] The \(\lambda\) sequence can be set using the lambda parameter. By default, we use a sequence motivated by Gaussian designs, described in Section 3.2.2 of (Bogdan et al. 2014). In most cases, this sequence should be preferred.

If you have an orthogonal design, the original Benjamini-Hochberg (BHq) sequence can also be used (see Section 1.1). We illustrate this below.

X.orth = sqrt(n) * qr.Q(qr(X))
y.orth = X.orth %*% beta + rnorm(n)
result = SLOPE(X.orth, y.orth, lambda='bhq')
fdp(result$selected)
## [1] 0.2307692

Notice that if we use the BHq sequence on a non-orthogonal design, FDR control is lost:

result = SLOPE(X, y, lambda='bhq')
fdp(result$selected)
## [1] 0.5454545

Therefore, the BHq sequence should only be used with orthogonal designs.

For completeness, we note that it is also possible to specify a custom sequence of \(\lambda\) values. (You probably don’t want to do this.)

References

Bogdan, Malgorzata, Ewout van den Berg, Chiara Sabatti, Weijie Su, and Emmanuel Candes. 2014. “SLOPE – Adaptive Variable Selection via Convex Optimization.” ArXiv:1407.3824.