oem
is a package for the estimation of various penalized regression models using the oem algorithm of (Xiong et al. 2016). The focus of oem
is to provide high performance computation for big tall data. Many applications not only have a large number of variables, but a vast number of observations; oem
is designed to perform well in these settings.
The simplest way to install oem
is via the CRAN repositories as the following:
To install the development version, first install the devtools
package
install.packages("devtools", repos = "http://cran.us.rproject.org")
and then installl oem
via the install_github
function
library(devtools)
install_github("jaredhuling/oem")
First load oem
library(oem)
Simulate data
nobs < 1e4
nvars < 100
x < matrix(rnorm(nobs * nvars), ncol = nvars)
y < drop(x %*% c(0.5, 0.5, 0.5, 0.5, 1, rep(0, nvars  5))) + rnorm(nobs, sd = 4)
Fit a penalized regression model using the oem
function
fit1 < oem(x = x, y = y, penalty = "lasso")
Plot the solution path
plot(fit1)
Function Name  Functionality 


Main fitting function 

Prediction for oem objects 

Plotting for oem objects 

log Likelihood for oem objects 

Crossvalidation function 

Prediction for cv.oem objects 

Plotting for cv.oem objects 

log Likelihood for cv.oem objects 
Penalty  Option Name  Penalty Form 

Lasso 

\(\lambda \sum_{j = 1}^pd_j\beta_j\) 
Elastic Net 

\(\alpha\lambda \sum_{j = 1}^pd_j\beta_j + (1  \alpha)\lambda \sum_{j = 1}^pd_j\beta_j^2\) 
MCP 

\(\lambda \sum_{j = 1}^pd_j \int_0^{\beta_j}(1  x/(\gamma\lambda d_j))_+\mathrm{d}x\) 
SCAD 

\(\sum_{j = 1}^p p^{SCAD}_{\lambda d_j,\gamma}(\beta_j)\) 
Group Lasso 

\(\lambda \sum_{k = 1}^Gd_k\sqrt{\sum_{j \in g_k}\beta_j^2}\) 
where \[ p_{\lambda, \gamma}^{SCAD}(\beta) = \left\{ \begin{array}{ll} \lambda\beta & \beta \leq \lambda ; \\ \frac{\beta^2  2\gamma\lambda\beta + \lambda^2}{2(\gamma  1)} & \lambda < \beta \leq \gamma\lambda ; \\ \frac{(\gamma + 1)\lambda^2}{2} & \beta > \gamma\lambda \\ \end{array} \right. \]
The following models are available currently.
Model  Option Name  Loss Form 

Linear Regression 

\(\frac{1}{2n}\sum_{i=1}^n(y_i  x_i^T\beta) ^ 2\) 
Logistic Regression 

\(\frac{1}{n}\sum_{i=1}^n\left[y_i x_i^T\beta  \log (1 + \exp\{ x_i^T\beta \} ) \right]\) 
There are plans to include support for multiple responses, binomial models (not just logistic regression), Cox’s proportional hazards model, and more if requested.
The oem algorithm is wellsuited to quickly estimate a solution path for multiple penalties simultaneously if the number of variables is not too large. The oem algorithm is only efficient for multiple penalties for linear models.
For the group lasso penalty, the groups
argument must be used. groups
should be a vector which indicates the group number for each variable.
fit2 < oem(x = x, y = y, penalty = c("lasso", "mcp", "grp.lasso"),
groups = rep(1:20, each = 5))
Plot the solution paths for all models
layout(matrix(1:3, ncol = 3))
plot(fit2, which.model = 1)
plot(fit2, which.model = 2)
plot(fit2, which.model = 3)
The following is a demonstration of oem’s efficiency for computing solutions for tuning parameter paths for multiple penalties at once.
The efficiency oem for fitting multiple penalties at once only applies to linear models. However, for linear models it is quite efficient, even for a high number of tuning parameters for many different penalties.
nobs < 1e5
nvars < 100
x2 < matrix(rnorm(nobs * nvars), ncol = nvars)
y2 < drop(x2 %*% c(0.5, 0.5, 0.5, 0.5, 1, rep(0, nvars  5))) + rnorm(nobs, sd = 4)
system.time(fit2a < oem(x = x2, y = y2, penalty = c("grp.lasso"),
groups = rep(1:20, each = 5), nlambda = 100L))
## user system elapsed
## 0.18 0.03 0.21
system.time(fit2b < oem(x = x2, y = y2, penalty = c("grp.lasso", "lasso", "mcp", "scad", "elastic.net"),
groups = rep(1:20, each = 5), nlambda = 100L))
## user system elapsed
## 0.21 0.03 0.24
system.time(fit2c < oem(x = x2, y = y2, penalty = c("grp.lasso", "lasso", "mcp", "scad", "elastic.net"),
groups = rep(1:20, each = 5), nlambda = 500L))
## user system elapsed
## 0.28 0.03 0.32
It is still more efficient to fit multiple penalties at once instead of individually for logistic regression, but the benefit is not as dramatic as for linear models.
nobs < 5e4
nvars < 100
x2 < matrix(rnorm(nobs * nvars), ncol = nvars)
y2 < rbinom(nobs, 1, prob = 1 / (1 + exp(drop(x2 %*% c(0.15, 0.15, 0.15, 0.15, 0.25, rep(0, nvars  5))))))
system.time(fit2a < oem(x = x2, y = y2, penalty = c("grp.lasso"),
family = "binomial",
groups = rep(1:20, each = 5), nlambda = 100L))
## user system elapsed
## 11.12 0.84 12.13
system.time(fit2b < oem(x = x2, y = y2, penalty = c("grp.lasso", "lasso", "mcp", "scad", "elastic.net"),
family = "binomial",
groups = rep(1:20, each = 5), nlambda = 100L))
## user system elapsed
## 47.05 4.01 52.34
Here we use the nfolds
argument to specify the number of folds for \(k\)fold cross validation
system.time(cvfit1 < cv.oem(x = x, y = y, penalty = c("lasso", "mcp", "grp.lasso"),
groups = rep(1:20, each = 5),
nfolds = 10))
## user system elapsed
## 1.14 0.05 1.21
Plot the cross validation mean squared error results for each model
layout(matrix(1:3, ncol = 3))
plot(cvfit1, which.model = 1)
plot(cvfit1, which.model = 2)
plot(cvfit1, which.model = 3)
The function xval.oem
offers accelerated cross validation for penalized linear models. In many cases is is orders of magnitude faster than cv.oem. It is only recommended for scenarios where the number of observations is larger than the number of variables. In addition to the computational gains in singlecore usage, it also benefits from parallelizaton using OpenMP (instead of using foreach, as used by cv.oem). For large enough problems, it has on a similar order of computation time as just fitting one OEM model.
nobsc < 1e5
nvarsc < 100
xc < matrix(rnorm(nobsc * nvarsc), ncol = nvarsc)
yc < drop(xc %*% c(0.5, 0.5, 0.5, 0.5, 1, rep(0, nvarsc  5))) + rnorm(nobsc, sd = 4)
system.time(cvalfit1 < cv.oem(x = xc, y = yc, penalty = "lasso",
groups = rep(1:20, each = 5),
nfolds = 10))
## user system elapsed
## 5.04 0.68 5.83
system.time(xvalfit1 < xval.oem(x = xc, y = yc, penalty = "lasso",
groups = rep(1:20, each = 5),
nfolds = 10))
## user system elapsed
## 0.83 0.03 0.94
system.time(xvalfit2 < xval.oem(x = xc, y = yc, penalty = "lasso",
groups = rep(1:20, each = 5),
nfolds = 10, ncores = 2))
## user system elapsed
## 1.04 0.05 0.64
system.time(ofit1 < oem(x = xc, y = yc, penalty = "lasso",
groups = rep(1:20, each = 5)))
## user system elapsed
## 0.20 0.01 0.21
A variety of evaluation metrics can be used for cross validation. The available metrics can be found in the table below
Model  Metric  type.measure= 

Linear Regression 
Mean squared error 

Mean absolute error 


Logistic Regression 
Deviance 

Area under the ROC curve 


Misclassification Rate 


Mean squared error of fitted mean 


Mean absolute error of fitted mean 

Consider a binary outcome setting with logistic regression.
nobs < 2e3
nvars < 20
x < matrix(runif(nobs * nvars, max = 2), ncol = nvars)
y < rbinom(nobs, 1, prob = 1 / (1 + exp(drop(x %*% c(0.25, 1, 1, 0.5, 0.5, 0.25, rep(0, nvars  6))))))
cvfit2 < cv.oem(x = x, y = y, penalty = c("lasso", "mcp", "grp.lasso"),
family = "binomial",
type.measure = "class",
groups = rep(1:10, each = 2),
nfolds = 10, standardize = FALSE)
In this case, misclassification rate is not the best indicator of performance. The classes here are imbalanced:
mean(y)
## [1] 0.069
Area under the ROC curve is an alternative classification metric to misclassification rate. It is available by setting type.measure = "auc"
.
cvfit2 < cv.oem(x = x, y = y, penalty = c("lasso", "mcp", "grp.lasso"),
family = "binomial",
type.measure = "auc",
groups = rep(1:10, each = 2),
nfolds = 10, standardize = FALSE)
With a very large dataset and computing cluster, the total size of a dataset may be very large, but if the number of variables is only moderately large (on the order of a few thousands) \(X^TX\) and \(X^TY\) may not be large and may already be available from other computations or can be computed trivially in parallel. The function oem.xtx
computes penalized linear regression models using the OEM algorithm only using \(X^TX\) and \(X^TY\). Standardization can be achieved by providing a vector of scaling factors (usually the standard deviations of the columns of x). The function is used like the following:
xtx < crossprod(xc) / nrow(xc)
xty < crossprod(xc, yc) / nrow(xc)
system.time(fit < oem(x = xc, y = yc,
penalty = c("lasso", "grp.lasso"),
standardize = FALSE, intercept = FALSE,
groups = rep(1:20, each = 5)))
## user system elapsed
## 0.18 0.02 0.20
system.time(fit.xtx < oem.xtx(xtx = xtx, xty = xty,
penalty = c("lasso", "grp.lasso"),
groups = rep(1:20, each = 5)) )
## user system elapsed
## 0.02 0.00 0.02
max(abs(fit$beta[[1]][1,]  fit.xtx$beta[[1]]))
## [1] 9.65894e15
max(abs(fit$beta[[2]][1,]  fit.xtx$beta[[2]]))
## [1] 9.769963e15
col.std < apply(xc, 2, sd)
fit.xtx.s < oem.xtx(xtx = xtx, xty = xty,
scale.factor = col.std,
penalty = c("lasso", "grp.lasso"),
groups = rep(1:20, each = 5))
The OEM package also provides functionality for ondisk computation with the big.oem
function, allowing for fitting penalized regression models on datasets too large to fit in memory. The big.oem
function uses the tools provided by the bigmemory
package, so a big.matrix object must be used for the design matrix.
set.seed(123)
nrows < 50000
ncols < 100
bkFile < "bigmat.bk"
descFile < "bigmatk.desc"
bigmat < filebacked.big.matrix(nrow=nrows, ncol=ncols, type="double",
backingfile=bkFile, backingpath=".",
descriptorfile=descFile,
dimnames=c(NULL,NULL))
# Each column value with be the column number multiplied by
# samples from a standard normal distribution.
set.seed(123)
for (i in 1:ncols) bigmat[,i] = rnorm(nrows)*i
yb < rnorm(nrows) + bigmat[,1]  bigmat[,2]
## outofmemory computation
fit < big.oem(x = bigmat, y = yb,
penalty = c("lasso", "grp.lasso"),
groups = rep(1:20, each = 5))
## fitting with inmemory computation
fit2 < oem(x = bigmat[,], y = yb,
penalty = c("lasso", "grp.lasso"),
groups = rep(1:20, each = 5))
max(abs(fit$beta[[1]]  fit2$beta[[1]]))
## [1] 1.534783e05
Computational time can be reduced a little via OpenMP parallelization of the key computational steps of the OEM algorithm. Simply use the ncores
argument to access parallelization. There is no need for the foreach package.
nobsc < 1e5
nvarsc < 500
xc < matrix(rnorm(nobsc * nvarsc), ncol = nvarsc)
yc < drop(xc %*% c(0.5, 0.5, 0.5, 0.5, 1, rep(0, nvarsc  5))) + rnorm(nobsc, sd = 4)
system.time(fit < oem(x = xc, y = yc,
penalty = c("lasso", "grp.lasso"),
standardize = FALSE, intercept = FALSE,
groups = rep(1:20, each = 25)))
## user system elapsed
## 3.30 0.07 3.40
system.time(fitp < oem(x = xc, y = yc,
penalty = c("lasso", "grp.lasso"),
standardize = FALSE, intercept = FALSE,
groups = rep(1:20, each = 25), ncores = 2))
## user system elapsed
## 3.82 0.07 2.22
If some variables should not be penalized, this can be specified through the use of the penalty.factor
argument for all penalties other than the group lasso. For the group lasso, the groupspecific weights can be modified by the group.weights
argument. penalty.factor
should be a vector of length equal to the number of columns in x
. Each element in penalty.factor
will be multiplied to the applied tuning parameter for each corresponding variable. For example, for a problem with 5 variables (ncol(x) = 5
), setting penalty.factor = c(1, 1, 1, 0, 0)
will effectively only allow penalization for the first three variables. The group.weights
argument should be a vector with length equal to the number of groups. Similarly to penalty.factor
, these weights will be multiplied to the penalty applied to each group. penalty.factor
and group.weights
can also be used to fit the adaptive lasso and adaptive group lasso, respectively.
The following example shows how to fit an adaptive lasso using oem
nobs < 1e4
nvars < 102
x < matrix(rnorm(nobs * nvars), ncol = nvars)
y < drop(x %*% c(0.5, 0.5, 0.5, 0.5, 1, 0.5, rep(0, nvars  6))) + rnorm(nobs, sd = 4)
lams < exp(seq(log(2.5), log(5e5), length.out = 100L))
ols.estimates < coef(lm.fit(y = y, x = cbind(1, x)))[1]
fit.adaptive < oem(x = x, y = y, penalty = c("lasso"),
penalty.factor = 1 / abs(ols.estimates),
lambda = lams)
group.indicators < rep(1:34, each = 3)
## norms of OLS estimates for each group
group.norms < sapply(1:34, function(idx) sqrt(sum((ols.estimates[group.indicators == idx]) ^ 2)))
fit.adaptive.grp < oem(x = x, y = y, penalty = c("grp.lasso"),
group.weights = 1 / group.norms,
groups = group.indicators,
lambda = lams)
For further information about oem
, please visit:
Xiong, Shifeng, Bin Dai, Jared Huling, and Peter Z.G. Qian. 2016. “Orthogonalizing EM: A DesignBased Least Squares Algorithm.” Technometrics, in Press.