linreg_methods_demo.Rmd

---
title: "Exploring methods for large-scale multiple linear regression"
author: Peter Carbonetto
output: workflowr::wflow_html
---

*Add overview here.*

```{r knitr-opts, include=FALSE}
knitr::opts_chunk$set(comment = "#",collapse = TRUE,results = "hold",
                      fig.align = "center",dpi = 120)
```

Load the packages and set the seed so that we may reproduce the
results more easily.

```{r load-pkgs}
library(MASS)
library(L0Learn)
library(SSLASSO)
library(EMVS)
library(R.matlab)
library(ggplot2)
library(cowplot)
set.seed(1)
source("../R/linreg_methods_demo_functions.R")
```

Simulate a data set with correlated variables in a similar way to
Example 1 from [Zou & Hastie 2005][zh2005].

```{r sim-data}
simulate_predictors_decaying_corr <- function (n, p, s = 0.5)
  return(mvrnorm(n,rep(0,p),s^abs(outer(1:p,1:p,"-"))))
simulate_outcomes <- function (X, b, se)
  return(drop(X %*% b - 1 + rnorm(nrow(X),sd = se)))
p      <- 1000
se     <- 3
b      <- rep(0,p)
b[1:3] <- c(3,1.5,2)
Xtrain <- simulate_predictors_decaying_corr(200,p,0.5)
Xtest  <- simulate_predictors_decaying_corr(500,p,0.5)
train  <- list(X = Xtrain,y = simulate_outcomes(Xtrain,b,se))
test   <- list(X = Xtest,y = simulate_outcomes(Xtest,b,se))
btrue  <- b
```

```{r save-mat-data}
writeMat("train.mat",
         X = scale(train$X,center = TRUE,scale = FALSE),
		 y = with(train,y - mean(y)))
```

L0Learn
-------

First, let's try the [L0Learn][l0learn] from
[Hazimeh and Mazumder 2020][hm2020]. In the
[vignette][l0learn-vignette], the authors suggest using the "L0L1" or
"L0L2" penalties for better predictive performance. The package
includes an interface to automatically select the penalty parameters
$\lambda, \gamma$ that minimize the test-set error. (L0Learn has two
model fitting algorithms; the CD algorithm is faster whereas the
CDPSI can sometimes produces better fits. For expediency we'll use
the CD algorithm.)

```{r l0learn}
l0learn_cv <- with(train,L0Learn.cvfit(X,y,penalty = "L0L1",algorithm = "CD"))
i      <- which.min(sapply(l0learn_cv$cvMeans,min))
j      <- which.min(l0learn_cv$cvMeans[[i]])
gamma  <- l0learn_cv$fit$gamma[i]
lambda <- l0learn_cv$fit$lambda[[i]][j]
```

Compare the coefficient estimates against the ground truth:

```{r plot-coef-l0learn, fig.height=2.5, fig.width=2.5}
b <- as.vector(coef(l0learn_cv,gamma = gamma,lambda = lambda))
b <- b[-1]
plot_coefs(btrue,b)
```

As expected, the L0 penalty shrinks most of the coefficients to zero.

The mean squared error (MSE) summarizes the accuracy of the
predictions in the test set examples:

```{r predict-l0learn, fig.height=2.5, fig.width=2.5}
y <- as.vector(predict(l0learn_cv,newx = test$X,gamma = gamma,lambda = lambda))
plot_responses(test$y,y)
```

Let's compare this to L0Learn with the simpler L0 penalty (which is a
special case of the L0L1 penalty in which $\gamma = 0$):

```{r l0learn-l0only, fig.height=2.5, fig.width=2.5}
l0learn_cv <- with(train,L0Learn.cvfit(X,y,penalty = "L0",algorithm = "CD"))
i <- which.min(l0learn_cv$cvMeans[[1]])
lambda <- l0learn_cv$fit$lambda[[1]][i]
y <- as.vector(predict(l0learn_cv,newx = test$X,lambda = lambda))
qplot(test$y,y) +
  geom_abline(intercept = 0,slope = 1,color = "deepskyblue",
              linetype = "dashed") +
  ggtitle(sprintf("mse = %0.3f",mse(test$y,y))) +
  theme_cowplot(font_size = 12) +
  theme(plot.title = element_text(face = "plain",size = 12))
```

Indeed, in this one example at least, L0Learn with the L0L1 penalty
has better prediction performance than the L0 penalty.

Trimmed Lasso
-------------

ADD BRIEF DESCRIPTION OF TRIMMED LASSO HERE.

Use `trimmed_lasso.m` to run the method. Then load the results:

```{r trimmed-lasso}
b <- drop(readMat("trimmed_lasso_b.mat")$b)
```

Compute the MLE of the intercept:

```{r trimmed-lasso-intercept}
b0 <- with(train,mean(y - X %*% b))
```

Compare the coef estimates vs. the ground truth:

```{r plot-coef-l0learn, fig.height=2.5, fig.width=2.5}
plot_coefs(btrue,b)
```

Compare predicted Y vs. true Y:

```{r predict-trimmed-lasso, fig.height=2.5, fig.width=2.5}
y <- drop(b0 + test$X %*% b)
plot_responses(test$y,y)
```

Notes: (1) Works particularly well in this example because the true is
indeed very sparse, and we happened to choose the "correct" k (in
general this will be hard to do). (2) Cross-validation will be need to
be used to get the right k.

SSLASSO
-------

Next let's try both the "adaptive" and "separable" variants of the
[Spike-and-Slab LASSO][sslasso]:

```{r sslasso-1}
sslasso_adapt <-
  with(train,SSLASSO(X,y,penalty = "adaptive",variance = "unknown"))
sslasso_sep <-
  with(train,SSLASSO(X,y,penalty = "separable",variance = "unknown"))
```

Let's take a look at the SSLASSO with separable penalty first. As far
as I can tell, the SSLASSO does not provide an automated way to select
the spike and slab penalties, so I'll choose one by hand:

```{r plot-coef-sslasso-sep, fig.height=2.5, fig.width=2.5}
b <- sslasso_sep$beta[,16]
qplot(btrue,b) +
  geom_abline(intercept = 0,slope = 1,color = "deepskyblue",
              linetype = "dashed") +
  theme_cowplot(font_size = 12)	 
```

Compare these estimates to the estimates obtained by the adaptive
SSLASSO:

```{r plot-coef-sslasso-adapt, fig.height=2.5, fig.width=2.5}
b <- sslasso_adapt$beta[,16]
qplot(btrue,b) +
  geom_abline(intercept = 0,slope = 1,color = "deepskyblue",
              linetype = "dashed") +
  theme_cowplot(font_size = 12)	 
```

Another annoying thing about the SSLASSO package is that it does not
provide a "predict" method. So we will compute the predictions by
hand after extracting the :

```{r predict-sslasso, fig.height=2.5, fig.width=5}
b0 <- sslasso_sep$intercept[16]
b  <- sslasso_sep$beta[,16]
y  <- drop(b0 + test$X %*% b)
p1 <- qplot(test$y,y) +
  geom_abline(intercept = 0,slope = 1,color = "deepskyblue",
              linetype = "dashed") +
  ggtitle(sprintf("separable (mse = %0.3f)",mse(test$y,y))) +
  theme_cowplot(font_size = 12)
b0 <- sslasso_adapt$intercept[16]
b  <- sslasso_adapt$beta[,16]
y  <- drop(b0 + test$X %*% b)
p2 <- qplot(test$y,y) +
  geom_abline(intercept = 0,slope = 1,color = "deepskyblue",
              linetype = "dashed") +
  ggtitle(sprintf("adaptive (mse = %0.3f)",mse(test$y,y))) +
  theme_cowplot(font_size = 12)
plot_grid(p1,p2)
```

EMVS
----

The [EMVS package][emvs] seems to be better documented (unfortunately,
it was removed from CRAN). It also has two variants with different
priors, the "independent" prior (which is recommended by the authors)
and the "conjugate" prior. Let's compare their performance in this
simulated example.

```{r emvs, results="hide"}
emvs_conj = with(train,
  EMVS(y,X,v0 = seq(0.1,2,length.out = 20),v1 = 10,independent = FALSE))
emvs_ind = with(train,
  EMVS(y,X,v0 = exp(seq(-18,-1,length.out = 20)),v1 = 1,independent = TRUE))
```

*Add text here.*

```{r emvs-conj, fig.height=2.5, fig.width=2.5}
i  <- which.max(emvs_conj$log_g_function)
b  <- with(emvs_conj,betas[i,] * prob_inclusion[i,])
qplot(btrue,b) +
  geom_abline(intercept = 0,slope = 1,color = "deepskyblue",
              linetype = "dashed") +
  theme_cowplot(font_size = 12)
```

*Add text here.*

```{r emvs-ind, fig.height=2.5, fig.width=2.5}
b <- with(emvs_ind,betas[14,] * prob_inclusion[14,])
qplot(btrue,b) +
  geom_abline(intercept = 0,slope = 1,color = "deepskyblue",
              linetype = "dashed") +
  theme_cowplot(font_size = 12)
```

*Add text here.*

```{r predict-emvs, fig.height=2.5, fig.width=5}
b0 <- emvs_conj$intersects[i]
b  <- emvs_conj$betas[i,]
y  <- drop(b0 + test$X %*% b)
p1 <- qplot(test$y,y) +
  geom_abline(intercept = 0,slope = 1,color = "deepskyblue",
              linetype = "dashed") +
  ggtitle(sprintf("conjugate (mse = %0.3f)",mse(test$y,y))) +
  theme_cowplot(font_size = 12)
b0 <- emvs_ind$intersects[14]
b  <- emvs_ind$beta[14,]
y  <- drop(b0 + test$X %*% b)
p2 <- qplot(test$y,y) +
  geom_abline(intercept = 0,slope = 1,color = "deepskyblue",
              linetype = "dashed") +
  ggtitle(sprintf("independent (mse = %0.3f)",mse(test$y,y))) +
  theme_cowplot(font_size = 12)
plot_grid(p1,p2)
```

[emvs]: https://cran.r-project.org/package=EMVS
[sslasso]: https://cran.r-project.org/package=SSLASSO
[zh2005]: https://doi.org/10.1111/j.1467-9868.2005.00503.x
[hm2020]: https://doi.org/10.1287/opre.2019.1919
[l0learn]: https://cran.r-project.org/package=L0Learn
[l0learn-vignette]: https://cran.r-project.org/web/packages/L0Learn/vignettes/L0Learn-vignette.html