---
title: "biglasso"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{biglasso}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 6, fig.height = 4)
```

# Small Data

## Linear regression

```{r}
library(biglasso)

data(colon)
X <- colon$X
y <- colon$y
dim(X)
```

```{r}
X[1:5, 1:5]
```

```{r}
# convert X to a big.matrix object
# X.bm is a pointer to the data matrix
X.bm <- as.big.matrix(X)
str(X.bm)
```

```{r}
dim(X.bm)
```

Note that this is the same as `X[1:5, 1:5]`:

```{r}
X.bm[1:5, 1:5]
```

```{r}
# fit entire solution path
fit <- biglasso(X.bm, y)
plot(fit)
```

## Cross-Validation

To run cross-validation, use `cv.biglasso()`:

```{r}
cvfit <- cv.biglasso(X.bm, y, seed = 1234, nfolds = 10)
```

After cross-validation, a few things we can do:

- plot the cross-validation plots:

```{r}
par(mfrow = c(2, 2), mar = c(3.5, 3.5, 3, 1), mgp = c(2.5, 0.5, 0))
plot(cvfit, type = "all")
```

- Summarize CV object:

```{r}
summary(cvfit)
```

- Extract non-zero coefficients at the optimal $\lambda$ value:

```{r}
coef(cvfit)[which(coef(cvfit) != 0), ]
```

## Logistic Regression

```{r}
data(Heart)
X <- Heart$X
y <- Heart$y
X.bm <- as.big.matrix(X)
fit <- biglasso(X.bm, y, family = "binomial")
plot(fit)
```

## Cox Regression

```{r}
library(survival)
X <- heart[, 4:7]
y <- Surv(heart$stop - heart$start, heart$event)
X.bm <- as.big.matrix(X)
fit <- biglasso(X.bm, y, family = "cox")
plot(fit)
```

## Multiple response Linear Regression (multitask learning)

```{r}
set.seed(10101)
n <- 300
p <- 300
m <- 5
s <- 10
b <- 1
x <- matrix(rnorm(n * p), n, p)
beta <- matrix(seq(from = -b, to = b, length.out = s * m), s, m)
y <- x[, 1:s] %*% beta + matrix(rnorm(n * m, 0, 1), n, m)
x.bm <- as.big.matrix(x)
fit <- biglasso(x.bm, y, family = "mgaussian")
plot(fit)
```

# Big Data

When the raw data file is very large, it's better to convert the raw data file
into a file-backed `big.matrix` by using a file cache. We can call function
`setupX`, which reads the raw data file and creates a backing file (.bin)
and a descriptor file (.desc) for the raw data matrix:

```{r}
# The data has 200 observations, 600 features, and 10 non-zero coefficients.
# This is not actually very big, but vignettes in R are supposed to render
# quickly, so we'll have to pretend that this can't fit in memory.
# Much larger data can be handled in the same way.
if (!file.exists("BigX.bin")) {
  X <- matrix(rnorm(1000 * 5000), 1000, 5000)
  beta <- c(-5:5)
  y <- as.numeric(X[, 1:11] %*% beta)
  write.csv(X, "BigX.csv", row.names = FALSE)
  write.csv(y, "y.csv", row.names = FALSE)
  X.bm <- setupX("BigX.csv", header = TRUE)
}
```

It's important to note that the above operation is only a one-time execution. Once
done, the data can be retrieved seamlessly by attaching its descriptor
file (`.desc`) in any new R session:

```{r, warning=F}
rm(list = c("X", "X.bm", "y")) # Pretend starting a new session
X.bm <- attach.big.matrix("BigX.desc")
y <- read.csv("y.csv")[, 1]
```

This is very appealing for big data analysis in that we don't need to "read" the
raw data again in a R session, which would be very time-consuming. The code below
again fits a lasso-penalized linear model, and runs 10-fold cross-validation:

```{r}
system.time({
  fit <- biglasso(X.bm, y)
})
```

```{r}
plot(fit)
```

For large data sets, it is often convenient to run cross validation in parallel (you can certainly use more cores on your machine, but CRAN vignettes are only allowed to use 2):

```{r}
cvfit <- cv.biglasso(X.bm, y, seed = 1234, ncores = 2, nfolds = 10)
```

```{r}
par(mfrow = c(2, 2), mar = c(3.5, 3.5, 3, 1), mgp = c(2.5, 0.5, 0))
plot(cvfit, type = "all")
```

# Links

* [biglasso on CRAN](https://cran.r-project.org/package=biglasso)
* [biglasso on GitHub](https://github.com/pbreheny/biglasso)
* [biglasso website](https://pbreheny.github.io/biglasso/index.html)
* [big.matrix manipulation](https://cran.r-project.org/package=bigmemory/index.html)