## Introduction

The oem package provides estimation for various penalized linear models using the Orthogonalizing EM algorithm. Documentation for the package can be found here: oem site.

Install using the devtools package (RcppEigen must be installed first as well):

``devtools::install_github("jaredhuling/oem")``

or by cloning and building using `R CMD INSTALL`

## Citation

Xiong, S., Dai, B., Huling, J., Qian, P. Z. G. (2016) Orthogonalizing EM: A design-based least squares algorithm, Technometrics, Volume 58, Pages 285-293,
http://dx.doi.org/10.1080/00401706.2015.1054436.

Huling, J.D. and Chien, P. (2018), Fast Penalized Regression and Cross Validation for Tall Data with the OEM Package, Journal of Statistical Software, to appear, URL: https://arxiv.org/abs/1801.09661.

## Penalties

### Lasso

``````library(microbenchmark)
library(glmnet)
library(oem)
# compute the full solution path, n > p
set.seed(123)
n <- 1000000
p <- 100
m <- 25
b <- matrix(c(runif(m), rep(0, p - m)))
x <- matrix(rnorm(n * p, sd = 3), n, p)
y <- drop(x %*% b) + rnorm(n)

lambdas = oem(x, y, intercept = TRUE, standardize = FALSE, penalty = "elastic.net")\$lambda[[1]]

microbenchmark(
"glmnet[lasso]" = {res1 <- glmnet(x, y, thresh = 1e-10,
standardize = FALSE,
intercept = TRUE,
lambda = lambdas)},
"oem[lasso]"    = {res2 <- oem(x, y,
penalty = "elastic.net",
intercept = TRUE,
standardize = FALSE,
lambda = lambdas,
tol = 1e-10)},
times = 5
)``````
``````## Unit: seconds
##           expr      min       lq     mean   median       uq      max neval cld
##  glmnet[lasso] 5.325385 5.374823 5.859432 6.000302 6.292411 6.304239     5  a
##     oem[lasso] 1.539320 1.573569 1.600241 1.617136 1.631450 1.639730     5   b``````
``````# difference of results
max(abs(coef(res1) - res2\$beta[[1]]))``````
``## [1] 1.048243e-07``
``````res1 <- glmnet(x, y, thresh = 1e-12,
standardize = FALSE,
intercept = TRUE,
lambda = lambdas)

# answers are now more close once we require more precise glmnet solutions
max(abs(coef(res1) - res2\$beta[[1]]))``````
``## [1] 3.763507e-09``

### Nonconvex Penalties

``````library(sparsenet)
library(ncvreg)
# compute the full solution path, n > p
set.seed(123)
n <- 5000
p <- 200
m <- 25
b <- matrix(c(runif(m, -0.5, 0.5), rep(0, p - m)))
x <- matrix(rnorm(n * p, sd = 3), n, p)
y <- drop(x %*% b) + rnorm(n)

mcp.lam <- oem(x, y, penalty = "mcp",
gamma = 2, intercept = TRUE,
standardize = TRUE,
nlambda = 200, tol = 1e-10)\$lambda[[1]]

gamma = 4, intercept = TRUE,
standardize = TRUE,
nlambda = 200, tol = 1e-10)\$lambda[[1]]

microbenchmark(
"sparsenet[mcp]" = {res1 <- sparsenet(x, y, thresh = 1e-10,
gamma = c(2,3), #sparsenet throws an error
#if you only fit 1 value of gamma
nlambda = 200)},
"oem[mcp]"    = {res2 <- oem(x, y,
penalty = "mcp",
gamma = 2,
intercept = TRUE,
standardize = TRUE,
nlambda = 200,
tol = 1e-10)},
"ncvreg[mcp]"    = {res3 <- ncvreg(x, y,
penalty = "MCP",
gamma = 2,
lambda = mcp.lam,
eps = 1e-7)},
"oem[scad]"    = {res4 <- oem(x, y,
gamma = 4,
intercept = TRUE,
standardize = TRUE,
nlambda = 200,
tol = 1e-10)},
"ncvreg[scad]"    = {res5 <- ncvreg(x, y,
gamma = 4,
eps = 1e-7)},
times = 5
)``````
``````## Unit: milliseconds
##            expr        min         lq       mean     median         uq
##  sparsenet[mcp] 1466.54465 1492.72548 1527.32113 1517.19926 1579.70827
##        oem[mcp]   95.71381   98.09740  105.90083  105.76415  110.31668
##     ncvreg[mcp] 5196.48035 5541.69429 5669.10010 5611.31491 5865.06723
##       oem[scad]   70.74110   71.46554   80.21926   78.76494   84.25458
##    ncvreg[scad] 5289.59790 5810.69254 5801.60997 5950.84377 5964.01276
##         max neval cld
##  1580.42800     5 a
##   119.61209     5  b
##  6130.94372     5   c
##    95.87013     5  b
##  5992.90288     5   c``````
``````diffs <- array(NA, dim = c(2, 1))
colnames(diffs) <- "abs diff"
rownames(diffs) <- c("MCP:  oem and ncvreg", "SCAD: oem and ncvreg")
diffs[,1] <- c(max(abs(res2\$beta[[1]] - res3\$beta)), max(abs(res4\$beta[[1]] - res5\$beta)))
diffs``````
``````##                          abs diff
## MCP:  oem and ncvreg 1.725859e-07
## SCAD: oem and ncvreg 5.094648e-08``````

### Group Penalties

In addition to the group lasso, the oem package offers computation for the group MCP, group SCAD, and group sparse lasso penalties. All aforementioned penalties can also be augmented with a ridge penalty.

``````library(gglasso)
library(grpreg)
library(grplasso)
# compute the full solution path, n > p
set.seed(123)
n <- 10000
p <- 200
m <- 25
b <- matrix(c(runif(m, -0.5, 0.5), rep(0, p - m)))
x <- matrix(rnorm(n * p, sd = 3), n, p)
y <- drop(x %*% b) + rnorm(n)
groups <- rep(1:floor(p/10), each = 10)

grp.lam <- oem(x, y, penalty = "grp.lasso",
groups = groups,
nlambda = 100, tol = 1e-10)\$lambda[[1]]

microbenchmark(
"gglasso[grp.lasso]" = {res1 <- gglasso(x, y, group = groups,
lambda = grp.lam,
intercept = FALSE,
eps = 1e-8)},
"oem[grp.lasso]"    = {res2 <- oem(x, y,
groups = groups,
intercept = FALSE,
standardize = FALSE,
penalty = "grp.lasso",
lambda = grp.lam,
tol = 1e-10)},
"grplasso[grp.lasso]"    = {res3 <- grplasso(x=x, y=y,
index = groups,
standardize = FALSE,
center = FALSE, model = LinReg(),
lambda = grp.lam * n * 2,
control = grpl.control(trace = 0, tol = 1e-10))},
"grpreg[grp.lasso]"    = {res4 <- grpreg(x, y, group = groups,
eps = 1e-10, lambda = grp.lam)},
times = 5
)``````
``````## Unit: milliseconds
##                 expr        min         lq       mean     median         uq
##   gglasso[grp.lasso]  679.59049  724.16350  858.99280  801.79179  865.83580
##       oem[grp.lasso]   59.84769   62.23879   64.11779   63.36026   64.30146
##  grplasso[grp.lasso] 3714.92601 3753.18663 4322.32431 4537.50185 4786.80867
##    grpreg[grp.lasso] 1216.21794 1248.84647 1270.46132 1269.71047 1287.75969
##         max neval cld
##  1223.58241     5 a
##    70.84075     5  b
##  4819.19839     5   c
##  1329.77201     5 a``````
``````diffs <- array(NA, dim = c(2, 1))
colnames(diffs) <- "abs diff"
rownames(diffs) <- c("oem and gglasso", "oem and grplasso")
diffs[,1] <- c(  max(abs(res2\$beta[[1]][-1,] - res1\$beta)), max(abs(res2\$beta[[1]][-1,] - res3\$coefficients))  )
diffs``````
``````##                      abs diff
## oem and gglasso  1.303906e-06
## oem and grplasso 1.645871e-08``````

#### Bigger Group Lasso Example

``````set.seed(123)
n <- 500000
p <- 200
m <- 25
b <- matrix(c(runif(m, -0.5, 0.5), rep(0, p - m)))
x <- matrix(rnorm(n * p, sd = 3), n, p)
y <- drop(x %*% b) + rnorm(n)
groups <- rep(1:floor(p/10), each = 10)

# fit all group penalties at once
"sparse.group.lasso")
system.time(res <- oem(x, y,
penalty = grp.penalties,
groups  = groups,
alpha   = 0.25, # mixing param for l2 penalties
tau     = 0.5)) # mixing param for sparse grp lasso ``````
``````##    user  system elapsed
##   2.043   0.222   2.267``````

### Fitting Multiple Penalties

The oem algorithm is quite efficient at fitting multiple penalties simultaneously when p is not too big.

``````set.seed(123)
n <- 100000
p <- 100
m <- 15
b <- matrix(c(runif(m, -0.25, 0.25), rep(0, p - m)))
x <- matrix(rnorm(n * p, sd = 3), n, p)
y <- drop(x %*% b) + rnorm(n)

microbenchmark(
"oem[lasso]"    = {res1 <- oem(x, y,
penalty = "elastic.net",
intercept = TRUE,
standardize = TRUE,
tol = 1e-10)},
"oem[lasso/mcp/scad/ols]"    = {res2 <- oem(x, y,
penalty = c("elastic.net", "mcp",
"grp.mcp", "sparse.grp.lasso",
"grp.mcp.net", "mcp.net"),
gamma = 4,
tau = 0.5,
alpha = 0.25,
groups = rep(1:10, each = 10),
intercept = TRUE,
standardize = TRUE,
tol = 1e-10)},
times = 5
)``````
``````## Unit: milliseconds
##                     expr      min       lq     mean   median       uq      max
##               oem[lasso] 125.6408 126.7870 130.3534 127.3374 133.4962 138.5055
##  oem[lasso/mcp/scad/ols] 148.3162 152.1743 153.0176 152.4529 154.4300 157.7144
##  neval cld
##      5  a
##      5   b``````
``````#png("../mcp_path.png", width = 3000, height = 3000, res = 400);par(mar=c(5.1,5.1,4.1,2.1));plot(res2, which.model = 2, main = "mcp",lwd = 3,cex.axis=2.0, cex.lab=2.0, cex.main=2.0, cex.sub=2.0);dev.off()
#

layout(matrix(1:4, ncol=2, byrow = TRUE))
plot(res2, which.model = 1, lwd = 2,
xvar = "lambda")
plot(res2, which.model = 2, lwd = 2,
xvar = "lambda")
plot(res2, which.model = 4, lwd = 2,
xvar = "lambda")
plot(res2, which.model = 7, lwd = 2,
xvar = "lambda")``````