An introduction to bestridge

2.1 Regression: linear model

We apply the methods to the data set reported in Scheetz et al. (2006), which contains gene expression levels of 18975 probes obtained from 120 rats. We are interested in finding probes that are related to that of gene TRIM32 using linear regression. This gene had been known to cause Bardet-Biedl syndrome, which is a genetic disease of multiple organ systems including the retina. For simplicity, the number of probes is reduced to 500 by picking out maximum 500 marginal correalation to Trim32 gene.

library(bestridge)
data("trim32", package = "bestridge")

Supposed we wants to fit a best subset ridge regression model with λ taking values in 100 grid points from 100 to 0.001 and the model size between max {p, n/log (n)} (the decimal portion is rounded) and 1, which is the default model size range where p denotes the number of total variables and n the sample size, via Generalized Information Criterion. We call the bsrr function as follows:

y <- trim32[, 1]
x <- as.matrix(trim32[, -1])
lm.bsrr <- bsrr(x, y)

The fitted coefficients can be extracted by running the coef function.

coef(lm.bsrr, sparse = TRUE)

To make a prediction on the new data, a predict function can be used as follows:

predict.bsrr <- predict(lm.bsrr, newx = x)

If the option newx is omitted, the fitted linear predictors are used.

2.2 Classification: logistic regression

We now turn to the classification task. We use the data set duke (West et al. 2001). This data set contains microarray experiment for 86 breast cancer patients. We are interested in classify the patients into estrogen receptor-positive (Status = 0) and estrogen receptor-negative (Status = 1) based on the expression level of the considered genes.

data("duke")
y <- duke$y
x <- as.matrix(duke[, -1])

Supposed we wants to do this classification based on logistic regression with a BSRR model. This can be realized through calling the bsrr function with the parameter family = "binomial as follows:

logi.bsrr <- bsrr(x, y, family = "binomial", method = "sequential")
#> BSS with sequential path

Calling the plot routine on an bsrr object fitted by type = "bsrr" will provide a heatmap of GIC of the fitted path.

plot(logi.bsrr)

3.1 Censored response: Cox proportional hazard model

The bestridge package also supports studying the relationship between predictors variables and survival time based on the Cox proportional hazard model. We now demonstrate the usage on a real data set, Alizadeh, Eisen, Davis, Ma, Lossos, Rosenwald, Boldrick, Sabet, Tran, Yu et al. (2000)(Alizadeh et al. 2000): gene-expression data in lymphoma patients. We illustrate the usage on a subset of the full data set. There were 50 patients with measurements on 1000 probes.

data(patient.data)

x <- patient.data$x
y <- patient.data$time
status <- patient.data$status

To implement the best subset ridge regression on the Cox proportional hazard model, call the bestridge function with family specified to be "cox".

cox.bsrr <- bsrr(x, cbind(y, status), family = "cox")

We use the summary function to draw a summary of the fitted bsrr object. 10 probes among the total 1000 are selected according to the result.

summary(cox.bsrr)
#> -------------------------------------------------------------------------------------------
#>     Penalized Primal-dual active algorithm with tuning parameter determined by 
#>     powell method using golden section method for line search 
#> 
#>     Best model with k = 10 lambda = 0.002287302 includes predictors: 
#> 
#>        x4        x5      x114      x195      x394      x533      x604      x656 
#> 14.285682 -8.229564 -6.160176  2.227762 11.106649 -4.282352  5.125808 -4.680172 
#>      x879      x913 
#> -8.494358  3.594189 
#> 
#>     log-likelihood:   -37.16859 
#>     deviance:          74.33718 
#>     GIC:               168.5587
#> -------------------------------------------------------------------------------------------

3.2 Criterion for tuning parameters selection

+ Information criterion: AIC, BIC, GIC, EBIC

+ Cross-validation

So far, we have been stick to the default Generalized Information Criterion for tuning parameter selections. In this package, we provide a bunch of criteria including the Akaike information criterion (Akaike 1974) and Bayesian information criterion (Schwarz et al. 1978), Generalized Information Criterion (Konishi and Kitagawa 1996), and extended BIC (Chen and Chen 2008, 2012), as well as cross-validation. By default, bsrr selects the tuning parameters according to the Generalized Information Criterion.

To choose one of the information criteria to select the optimal values of model size s and shrinkage λ, the input value tune in the bsrr function needs to be specified:

lm.bsrr.ebic <- bsrr(x, y, tune = "ebic")

To use cross-validation for parameter selections, set the input value of tune = "cv". By default, 5-fold cross-validation will be conducted.

lm.bsrr.cv <- bsrr(x, y, tune = "cv", nfolds = 5)

3.3 Paths for tuning parameters selection

+ Sequential method

+ Powell's method

We shall give a more explicit description of the parameters tuning paths. As mentioned in the introduction, we either apply a method = "sequential" path for examining each combination of the two tuning parameters s and λ sequentially and chose the optimal one according to certain information criteria or cv, or a less computational burdensome Powell method (method = "psequential" or method = "pgsection"). The initial starting point x₀ of the Powell method is set to be the (s_min, λ_min), and the initial search vectors v₁, v₂ are the coordinate unit vectors. The method minimizes the value of chosen information criterion or cross-validation error by a bi-directional search along each search vector on the 2-dimensional tuning parameter space composed of s and λ, in turn. Callers can choose to find the optimal combination along each search vector sequentially or conducting a Golden-section search by setting method = "psequential" or method = "pgsection". If warm.start = TRUE, we use the estimate from the last iteration in the PDAS algorithm as a warm start.

By default, the tuning parameters are determined by the Powell method through a golden-section line search and we exploit warm starts. The minimum λ value is set to be 0.01 and the maximum 100. The maximum model size is the smaller one of the total number of variables p and n/log (n). Advanced users of this toolkit can change this default behavior and supply their own tuning path.

To conduct parameter selections sequentially, users should specify the method to be “sequential” and provide a list of λ to the lambda.list as well as an increasing list of model size to the s.list.

my.lambda.list <- exp(seq(log(10), log(0.01), length.out = 10))
my.s.list <- 1:10

lm.bsrr.seq <- bsrr(x, y, method = "sequential", s.list = my.s.list,
                   lambda.list = my.lambda.list)

To conduct parameter selections using the Powell method on a user-supplied tuning parameter space, users should assign values to s.min, s.max, lambda.min, and lambda.max. 100 values of λ will be generated decreasing from lambda.max to lambda.min on the log scale. Here we do the line search sequentially.

my.s.min <- 1
my.s.max <- 10
my.lambda.min <- 0.01
my.lambda.max <- 10

lm.bsrr.powell <- bsrr(x, y, method = "pgsection", 
                      s.min = my.s.min, s.max = my.s.max, 
                      lambda.min = my.lambda.min, lambda.max = my.lambda.max)

3.4 Sure independence screening

We also provide feature screening option to deal with ultrahigh dimensional data for computational and statistical efficiency. Users can apply the sure independence screening method (Saldana and Feng 2018) based on maximum marginal likelihood estimators, to pre-exclude some irrelevant variables before fitting a model by pressing a integer to screening.num. The SIS will filter a set of variables with size equals to screening.num. Then the active set updates are restricted on this set of variables.

lm.bsrr.screening <- bsrr(x, y, screening.num = round(nrow(x)/log(nrow(x))))

4.1 Settings

In this section, we will conduct a monte carol study on our proposed best subset ridge regression method and make a comparison with other commonly used variable selection methods in three aspects. The first aspect is the predictive performance on a held-out testing data of the same size of training data. For linear regression, this is defined as $\frac{\Vert X \beta^\dagger - X \hat{\beta}\Vert _2}{\Vert X \beta^\dagger\Vert_2}$, where β^† is the true coefficient and β̂ is an estimator. For logistic regression, we calculate the classification accuracy by the area under the ROC curve (AUC). The second aspect is the coefficient estimation performance defined as ‖β^† − β̂‖, where β^† denote the underlying true coefficients and β^† is out estimation. The third aspect is the selection performance in terms of true positive rate (TPR), false positive rate (FPR), the Matthews correlation coefficient (MCC) score, and the model size.

We generate a multivariate Gaussian data matrix X_n × p ∼ MVN(0, Σ) and consider the variables to have exponential correlation. n denotes the sample size and p the number of total variables, and we set n = 200, p = 2000. We consider the following instance of Σ := ((σ_ij)), setting each entry σ_ij = 0.5^{|i − j|}. We use a sparse coefficient vector β^† with 20 equi-spaced nonzero entries, each set to 1. Then the response vector y is generated with gaussian noise added. For linear regression, y = Xβ^† + ϵ. For logistic regression, y = Bernoulli(Pr(Y = 1)), where Pr(Y = 1) = exp (x^Tβ^† + ϵ)/(1 + exp (x^Tβ^† + ϵ)). ϵ ∼ N(0, σ²) is independent of X. We define the signal-to-noise ratio (SNR) by SNR = $\frac{Var(X\beta^\dagger)}{Var(\epsilon)} = \frac{\beta^{\dagger T} \Sigma \beta^\dagger}{\sigma^2}$.

we compare the performance of the following methods:

1. (BSRR methods): Tuning parameters selected through a sequential path and the Powell method. Denoted as Grid-BSRR and Powell-BSRR.

2. (BSS method): The best subset selection method. We consider the primal-dual active subset selection method and use our own implementation in bestridge package.

3. (L₁ method): Lasso estimator. We use the implementation of glmnet.

4. (L₁L₂ method): Elastic net estimator. This uses a combination of the L₁ and L₂ regularization. We use the implementation of glmnet.

For BSRR estimators, the 2D grid of tuning parameters has λ taking 100 values between λ_max and λ_min on a log scale. The λ_max and λ_min will be specified. For the BSS method, the model size parameter takes value in {1, …, 30}. The parameter combined the L₁ penalty and the L₂ penalty in the Elastic net is fixed at 0.5. And for L₁ and L₁L₂ methods, we use the default settings. The regularization parameters are chosen by 5-fold cross-validation.

4.2 Linear regression

Results for linear regression are reported in Figure 1. When SNR is low, BSS has the largest estimation error and selects the smallest model size compared with other methods, while the extra L₂ penalty in BSRR has effectively lowered the estimation error of BSS. For high SNR, the performance of BSS and BSRR is similar to each other. In terms of prediction error, all methods are good in the high SNR scenario, but it costs the Lasso and the Elastic Net a very dense supports while still cannot catch up with BSS and BSRR. Notice that for SNR 100, the selected model size of the Elastic Net is almost 7 times the true model size. This trend is no better for the Lasso. Of all measures across the whole SNR range, BSRR generally exhibits excellent performance—accurate in variables selection and prediction.

4.3 Logistic regression

Results for linear regression are reported in Figure 2. The result shows that as the SNR rises, the extra shrinkage in BSRR methods helps the best subset ridge regression make impressive improvement in terms of prediction accuracy, estimation ability, and variable selection, where it outperforms the state-of-the-art variable selection methods Lasso and Elastic Net.