Basics

Introduction

The package implements a collection of Generalized Elastic Net (GELnet) solvers, as outlined in the following publication: https://doi.org/10.1371/journal.pcbi.1004790 This vignette covers the basic usage of the gelnet package.

Models are assembled from small, independent pieces that are combined with +, much like a ggplot2 plot is built up from geoms, scales and themes added to a base ggplot() call. A model definition always starts with gelnet(X), and is then extended with a task (what is being predicted) and one or more regularizers (how the model is penalized). The resulting definition is handed to gelnet_train() to fit the model. Nearly all of the regularization pieces have reasonable defaults, so the interface stays simple when the extra flexibility isn’t needed.

Building your first GELnet model

Let \(X\) be a samples-by-features data matrix and \(y\) be a column vector of labels. Typically, \(X\) and \(y\) are determined by the prediction task at hand, but for the purposes of this tutorial, we are going to generate them randomly:

X <- matrix( rnorm( 1000 ), 20, 50 )
y <- rnorm( 20 )

Let’s further assume that we are given a feature-feature relationship matrix \(A\). If working with genomic features, \(A\) might be the adjacency of a gene-gene interaction network. Again, we are going to generate a random matrix for the purposes of this tutorial:

A <- matrix( sample( 0:1, 50*50, repl=TRUE ), 50, 50 )
A <- A & t(A)  ## Make the matrix symmetric

We are now going to utilize the GELnet toolkit to learn a linear regression model, such that the model weights are more similar for the features that share an interaction on \(A\). As discussed in the manuscript, this can be achieved by formulating a feature-feature penalty matrix using either the graph Laplacian or \((I-D)\), where \(D\) is the graph diffusion matrix and \(I\) is the identity matrix. The gelnet package provides a function to compute the graph Laplacian from the adjacency. Here, we utilize the normalized Laplacian to keep the penalty term on the same scale as the traditional ridge regression:

library( gelnet )
L <- adj2nlapl(A)

A model definition is assembled by starting with the data matrix, adding a task (model_lin() for linear regression) and the L1- and L2-norm regularizers:

mdef <- gelnet(X) + model_lin(y) + rglz_L1(0.1) + rglz_L2(1, P = L)

where we set the L1-norm and L2-norm penalties to 0.1 and 1, respectively. The terms can be added in any order, and the definition can be built up incrementally.

The model is trained by handing the definition to gelnet_train():

model <- gelnet_train( mdef )
## Training a linear regression model
## Initial objective value: 0.433112
## Final value is 0.290151 after iteration 6

The response for new samples is computed via the dot product with the weights:

Xnew <- matrix( rnorm( 500 ), 10, 50 )
Xnew %*% model$w + model$b
##              [,1]
##  [1,] -0.07155645
##  [2,]  0.39159264
##  [3,]  0.07143397
##  [4,]  0.01555808
##  [5,]  0.50425639
##  [6,] -0.13768645
##  [7,]  0.15937134
##  [8,]  0.40538510
##  [9,]  0.02405002
## [10,]  0.53396980

Other regression problems

Linear regression is one of the three types of prediction problems supported by the package. The other two are binary logistic regression and one-class logistic regression. The latter is outlined in the following paper: https://doi.org/10.1142/9789814749411_0037

The package recognizes the problem type based on which task is added to the model definition: model_lin() for linear regression, model_blr() for binary logistic regression, and model_oclr() for one-class logistic regression. To train a binary predictor, we have to provide \(y\) as a two-level factor, where the first level is treated as the positive class.

y <- factor( y > 0, levels=c(TRUE,FALSE) )
mdef2 <- gelnet(X) + model_blr(y) + rglz_L1(0.1) + rglz_L2(1, P = L)
model2 <- gelnet_train( mdef2 )
## Training a binary logistic regression model
## Initial objective value: 0.693147
## Final value is 0.643766 after iteration 3

If we were to score the training data using this model, we can observe that the positive samples are receiving higher scores than the negative ones

data.frame( scores= X %*% model2$w + model2$b, labels= y )
##         scores labels
## 1   0.40905286   TRUE
## 2   0.01897029  FALSE
## 3  -0.05150661  FALSE
## 4   0.78077524   TRUE
## 5  -0.19938677  FALSE
## 6   0.14124765   TRUE
## 7   0.61981896   TRUE
## 8   0.30265754   TRUE
## 9   0.21585262   TRUE
## 10  0.71104158   TRUE
## 11 -0.25601148  FALSE
## 12  0.05138601  FALSE
## 13  0.62110636   TRUE
## 14 -0.29814396  FALSE
## 15 -0.14075878  FALSE
## 16 -0.27964162  FALSE
## 17  0.62046173   TRUE
## 18 -0.08534166  FALSE
## 19  0.64124516   TRUE
## 20  0.32721525   TRUE

However, if there is class imbalance, the scores will tend to be skewed towards the class with more samples. This can be addressed by using the balanced argument of model_blr():

mdef2bal <- gelnet(X) + model_blr(y, balanced=TRUE) + rglz_L1(0.1) + rglz_L2(1, P = L)
model2bal <- gelnet_train( mdef2bal )
## Training a binary logistic regression model
## Initial objective value: 0.693147
## Final value is 0.647321 after iteration 3
data.frame( scores= X %*% model2bal$w + model2bal$b, labels= y )
##         scores labels
## 1   0.23932585   TRUE
## 2  -0.16018540  FALSE
## 3  -0.22850485  FALSE
## 4   0.61663355   TRUE
## 5  -0.37963015  FALSE
## 6  -0.03271190   TRUE
## 7   0.45084906   TRUE
## 8   0.12681418   TRUE
## 9   0.04003994   TRUE
## 10  0.54365911   TRUE
## 11 -0.43814507  FALSE
## 12 -0.12795722  FALSE
## 13  0.45168603   TRUE
## 14 -0.47969599  FALSE
## 15 -0.32031123  FALSE
## 16 -0.46129722  FALSE
## 17  0.45222823   TRUE
## 18 -0.26564079  FALSE
## 19  0.47517577   TRUE
## 20  0.15656425   TRUE

Traditionally, the loss function for logistic regression is averaged over \(n\), the number of samples. This causes every sample to make the same contribution to the loss, which is what causes the skew towards the larger class. By using the balanced flag, the problem is reformulated slightly such that the loss is averaged over the positive and negative samples separately, and then the mean of both averages is used as the overall loss.

Finally, we can build a one-class logistic regression model using just the positive samples. To train a one-class model we add model_oclr(), which requires no labels:

j <- which( y == TRUE )
mdef1 <- gelnet(X[j,]) + model_oclr() + rglz_L1(0.1) + rglz_L2(1, P = L)
model1 <- gelnet_train( mdef1 )
## Training a one-class model
## Initial objective value: 0.693147

The model can now be used as a detector that recognizes the positive samples

data.frame( scores= X %*% model1$w, labels= y )
##         scores labels
## 1   0.28372681   TRUE
## 2  -0.06447266  FALSE
## 3   0.05392755  FALSE
## 4   0.58725516   TRUE
## 5   0.15966498  FALSE
## 6   0.34038067   TRUE
## 7   0.32707188   TRUE
## 8  -0.06816523   TRUE
## 9   0.43508772   TRUE
## 10  0.44229930   TRUE
## 11 -0.31340585  FALSE
## 12 -0.10261685  FALSE
## 13  0.32982234   TRUE
## 14 -0.21052052  FALSE
## 15  0.02635640  FALSE
## 16 -0.29992964  FALSE
## 17  0.39294222   TRUE
## 18 -0.28778595  FALSE
## 19  0.44788940   TRUE
## 20  0.38330742   TRUE