Glmtlp Vignette
Introduction
Glmtlp is a package that makes it incredibly easy to fit a generalized linear model via penalized maximum likelihood. The regularization path is computed for the truncated lasso at a grid of values for the regularization parameter lambda. The algorithm is fast, and can exploit sparsity in the input matrix x. It fits linear, logistic and multinomial, poisson, and Cox regression models. Most supporting functions used in a popular R package glmnet can be directly used with output of glmtlp. We try our best to mimic the format of input and output in glmnet to decrease the learning curve of our package.
The authors of glmtlp are Chong Wu and Wei Pan, and is maintained by Chong Wu.
glmtlp solves the following problem
\[ \min_{\beta_0,\beta}\frac{1}{N}\sum_{i=1}^{N}w_i l(y_i,\beta_0 + \beta'x_i) + \lambda \sum_{j=1}^{p} \min(|\beta_j|/\tau,1) \]
over a grid of values of \(\lambda\) covering the entire range. Here \(l(y,\eta)\) is the negative log-likelihood contribution for observation \(i\); e.g. for the Gaussian case it is \(1/2 (y-\eta)^2\). The truncated lasso penalty is controlled by \(\tau\), and partly solves the biased issues in lasso penalty. The tuning parameter \(\lambda\) controls the overall strength of the penalty.
Installation
To install the stable version from CRAN, simply run the following from an R console (not available now):
install.packages("glmtlp")
To install the latest development builds directly from GitHub, run this instead:
if (!require("devtools"))
install.packages("devtools")
library(devtools)
devtools::install_github("ChongWu-Biostat/glmtlp")
Quick Start
The purpose of this section is to give users a general idea of the package. In terms of usage, the package is almost the same as glmtlp, however, in terms of performance, glmtlp should be better than glmnet in a wide range of scenarios. Users may have a better idea after this section.
First, we load the glmtlp package:
library(glmtlp)
We will keep most parameters as default and mainly focus on the general sense of the package. We load a set of data created beforehand for illustration. Users can either load their own data or use those saved in the package. This data is exactly same as a data set in glmnet and we use the exactly same example codes. We hope we can get some general ideas about the difference between elastic net and truncated lasso and understand why non-convex truncated lasso penalty is superior than elastic net or other convex penalties.
data("QuickStartExample")
The command loads an input matrix \(x\) and a response vector \(y\) from the package. Then we fit the model using the basic function glmTLP.
fit = glmTLP(x,y)
“fit” is an object of class glmnet that contains all the relevant information of the fitted model for further use. You can think you use a better penalty in glmnet and then get the output. Most supporting functions such as plot, print, coef and predict can be directly used for the output. For example, we can visualize the coefficients by executing the plot function:
plot(fit)
|
|
Each curve corresponds to a variable. The above figure shows the path of its coefficient against the \(l_1\)-norm of the whole coefficient vector at as \(\lambda\) varies. The axis above indicates the number of nonzero coefficients at the current \(\lambda\), which is the effective degrees of freedom (df) for the lasso. Users may also wish to annotate the curves; this can be done by setting label = TRUE in the plot command. Again, this function performs exactly the same as plot in glmnet.
A summary of the glmTLP path at each step can be printed via function print:
print(fit)
Call: glmTLP(x = x, y = y) Df %Dev Lambda [1,] 7 0.8544 1.8350000 [2,] 7 0.8623 1.6720000 [3,] 8 0.8696 1.5230000 [4,] 8 0.8761 1.3880000 [5,] 8 0.8815 1.2650000 [6,] 8 0.8859 1.1520000 [7,] 8 0.8896 1.0500000 [8,] 8 0.8927 0.9567000 [9,] 8 0.8970 0.8717000 ....
It shows the number of nonzero coefficients (Df), the percent of null deviance explained (%Dev), and the value of corresponding \(\lambda\) (Lambda). By default glmTLP calss for 100 values of lambda, where lambda is automatically determined by the internal function.
We can obtain the actual coefficients at one or more \(\lambda\)'s within the range of the sequence. Note that even though coef can return the results for \(\lambda\) out of the range, the resutls are not reliable.
coef(fit,s=0.1)
21 x 1 sparse Matrix of class "dgCMatrix" 1 (Intercept) 0.11678951 V1 1.37182915 V2 0.01835052 V3 0.75251617 V4 0.04492075 V5 -0.89167611 V6 0.60355346 V7 0.10393552 V8 0.38759151 V9 -0.01738196 V10 0.10891171 V11 0.23712445 V12 -0.05190906 V13 -0.03808770 V14 -1.14807069 V15 -0.11590021 V16 -0.03162155 V17 -0.02981730 V18 0.04185775 V19 . V20 -1.13622401
Users can also make predicitions at specific \(\lambda\)'s with new input data:
nx = matrix(rnorm(10*20),10,20) predict(fit,newx=nx,s=c(0.1,0.05))
1 2 [1,] -2.7302125 -2.7212503 [2,] -0.1261640 -0.1037789 [3,] 2.4302532 2.4263055 [4,] 4.4345692 4.4825886 [5,] -0.8703297 -0.8770905 [6,] -3.6949918 -3.6809959 [7,] -0.5279150 -0.4789763 [8,] 1.6668486 1.6533317 [9,] 0.8976858 0.9092636 [10,] -3.7293284 -3.7076979
The function glmTLP returns a sequence of models for the users to choose from. However, in many cases, users may want to the package to select one of them adaptively. We recommend use cross-validation method to select the tuning parameters.
cv.glmTLP is the main function to do cross-validation here. Again, various supporting methods such as plotting and prediction used in package glmnet can be directly applied here. We still act on the same sample data loaded before.
cvfit = cv.glmTLP(x, y,tau = 1)
cv.glmTLP returns a cv.glmnet object, which is defined in package glmnet. We encourage users the well-designed functions for potential tasks.
We can plot the object.
plot(cvfit)
|
|
The above figure shows the cross-validation curve (red dotted line), and upper and lower standard deviation curves along the \(\lambda\) sequence (error bars). Two selected \(\lambda\)'s are indicated by the vertical dotted lines.
We can view the selected \(\lambda\)'s and the corresponding coefficients. For example,
cvfit$lambda.min
[1] 0.3773556
lambda.min is the value of \(\lambda\) that gives minimum mean cross-validated error. We also save lambda.1se, which gives the most regularized model such that error within one standard error of the minimum. These two \(\lambda\) are recommended by glmnet.
Like glmTLP, we can extract coefficient via function coef.
coef(cvfit, s = "lambda.min")
21 x 1 sparse Matrix of class "dgCMatrix" 1 (Intercept) 0.146500723 V1 1.340178977 V2 . V3 0.707734699 V4 . V5 -0.847077207 V6 0.553671322 V7 0.037104906 V8 0.346252025 V9 . V10 0.007849848 V11 0.184759134 V12 . V13 . V14 -1.084353223 V15 -0.001922412 V16 . V17 . V18 . V19 . V20 -1.067892212
Like glmnet, we store the result in the sparse matrix format since the solutions along the regularization path are often sparse. If you prefer non-sparse format, use as.matrix() to convert it to matrix.
Predictions can be made based on the fitted cv.glmnet object. For example:
predict(cvfit, newx = x[1:5,], s = "lambda.min")
1
[1,] -1.3591677
[2,] 2.5760596
[3,] 0.5843585
[4,] 2.0212049
[5,] 1.5709110
newx is for the new input and s, as before, is the value(s) of \(\lambda\) at which prefictions are made.
If you are familiar with glmnet, you can find in terms of usage, glmtlp and glmnet is almost exactly the same.