## -----------------------------------------------------------------------------
## Name: 04-lm-iii.R
## Description: Script for Chapter 4 of "Predictive Modeling"
## Link: https://egarpor.github.io/PM-UC3M/
## License: https://creativecommons.org/licenses/by-nc-nd/4.0/
## Author: Eduardo García-Portugués
## Version: 6.1.0
## -----------------------------------------------------------------------------

## ----shrinkage----------------------------------------------------------------
# Load data -- baseball players statistics
data(Hitters, package = "ISLR")

# Discard NA's
Hitters <- na.omit(Hitters)

# The glmnet function works with the design matrix of predictors (without
# the ones). This can be obtained easily through model.matrix()
x <- model.matrix(Salary ~ ., data = Hitters)[, -1]
# [, -1] to remove the column of 1's associated with the intercept, since the
# intercept will be added by default in glmnet::glmnet and if we do not exclude
# it here we will end with two intercepts in the fit, one of them being NA.

# Interestingly, note that in Hitters there are two-level factors and these
# are automatically transformed into dummy variables in x -- the main advantage
# of model.matrix
head(Hitters[, 14:20])
head(x[, 14:19])

# We also need the vector of responses
y <- Hitters$Salary

## ----modelmatrix--------------------------------------------------------------
# Data with NA in the first observation and factors with two levels
data_na <- data.frame("x1" = rnorm(3), "x2" = factor(c("A", "B", "A")),
                      "x3" = factor(c("F", "F", "M")), "y" = rnorm(3))
data_na$x1[1] <- NA

# The first observation disappears!
model.matrix(y ~ ., data = data_na)

# Still removes NA's
model.matrix(y ~ ., data = data_na, na.action = "na.pass")

# Does not remove NA's
model.matrix.lm(y ~ ., data = data_na, na.action = "na.pass")

# If y ~ 0 + ., the first factor gets two dummies instead of one!
model.matrix(y ~ 0 + ., data = data_na)

## ----ridge-1------------------------------------------------------------------
# Call to the main function -- use alpha = 0 for ridge regression
library(glmnet)
ridge_mod <- glmnet(x = x, y = y, alpha = 0)
# By default, it computes the ridge solution over a set of lambdas
# automatically chosen. It also standardizes the variables by default to make
# the model fitting, since the penalization is scale-sensitive. Importantly,
# the coefficients are returned on the original scale of the predictors

# Plot of the solution path -- gives the value of the coefficients for different
# measures in xvar (penalization imposed to the model or fitness)
plot(ridge_mod, xvar = "norm", label = TRUE)
# xvar = "norm" is the default: L1 norm of the coefficients sum_j abs(beta_j)

# Versus lambda
plot(ridge_mod, label = TRUE, xvar = "lambda")

# Versus the percentage of deviance explained -- this is a generalization of the
# R^2 for generalized linear models. Since we have a linear model, this is the
# same as the R^2
plot(ridge_mod, label = TRUE, xvar = "dev")
# The maximum R^2 is slightly above 0.5

# Indeed, we can see that R^2 = 0.5461
summary(lm(Salary ~ ., data = Hitters))$r.squared

# Some persistently important predictors are 15, 14, and 19
colnames(x)[c(15, 14, 19)]
# Dummies associated to league and division are important

# What is inside glmnet's output?
names(ridge_mod)

# lambda versus R^2 -- fitness decreases when sparsity is introduced, in
# in exchange of better variable interpretation and avoidance of overfitting
plot(-log(ridge_mod$lambda), ridge_mod$dev.ratio, type = "l",
     xlab = "-log(lambda)", ylab = "R2")
ridge_mod$dev.ratio[length(ridge_mod$dev.ratio)]
# Slightly different to lm's because it compromises accuracy for speed

# The coefficients for different values of lambda are given in $a0 (intercepts)
# and $beta (slopes) or, alternatively, both in coef(ridgeMod)
length(ridge_mod$a0)
dim(ridge_mod$beta)
length(ridge_mod$lambda) # 100 lambda's were automatically chosen

# Estimated coefficients for the 50th value of lambda (includes intercept also)
coef(ridge_mod)[, 50]
ridge_mod$lambda[50]

# Zoom in path solution
plot(ridge_mod, label = TRUE, xvar = "lambda",
     xlim = -log(ridge_mod$lambda[50]) + c(-2, 2), ylim = c(-30, 10))
abline(v = -log(ridge_mod$lambda[50]))
points(rep(-log(ridge_mod$lambda[50]), nrow(ridge_mod$beta)),
       ridge_mod$beta[, 50], pch = 16, col = 1:6)

# The squared l2-norm of the coefficients decreases as lambda increases
plot(-log(ridge_mod$lambda), sqrt(colSums(ridge_mod$beta^2)), type = "l",
     xlab = "-log(lambda)", ylab = "l2 norm")

## ----ridge-2------------------------------------------------------------------
# If we want, we can choose manually the grid of penalty parameters to explore
# The grid should be descending
ridge_mod2 <- glmnet(x = x, y = y, alpha = 0, lambda = 100:1)
plot(ridge_mod2, label = TRUE, xvar = "lambda") # Not a good choice!

# Lambda is a tuning parameter that can be chosen by cross-validation, using as
# error the MSE (other possible error can be considered for generalized models
# using the argument type.measure)

# 10-fold cross-validation. Change the seed for a different result
set.seed(12345)
kcv_ridge <- cv.glmnet(x = x, y = y, alpha = 0, nfolds = 10)

# The lambda that minimizes the CV error is
kcv_ridge$lambda.min

# Equivalent to
ind_min <- which.min(kcv_ridge$cvm)
kcv_ridge$lambda[ind_min]

# The minimum CV error
kcv_ridge$cvm[ind_min]
min(kcv_ridge$cvm)

# Potential problem! Minimum occurs at one extreme of the lambda grid in which
# CV is done. The grid was automatically selected, but can be manually inputted
range(kcv_ridge$lambda)
lambda_grid <- 10^seq(log10(kcv_ridge$lambda[1]), log10(0.1),
                      length.out = 150) # log-spaced grid
kcv_ridge2 <- cv.glmnet(x = x, y = y, nfolds = 10, alpha = 0,
                        lambda = lambda_grid)

# Much better
plot(kcv_ridge2)
kcv_ridge2$lambda.min

# But the CV curve is random, since it depends on the sample. Its variability
# can be estimated by considering the CV curves of each fold. An alternative
# approach to select lambda is to choose the largest within one standard
# deviation of the minimum error, in order to favor simplicity of the model
# around the optimal lambda value. This is known as the "one standard error rule"
kcv_ridge2$lambda.1se

# Location of both optimal lambdas in the CV loss function in dashed vertical
# lines, and lowest CV error and lowest CV error + one standard error
plot(kcv_ridge2)
ind_min2 <- which.min(kcv_ridge2$cvm)
abline(h = kcv_ridge2$cvm[ind_min2] + c(0, kcv_ridge2$cvsd[ind_min2]))
# The consideration of the one standard error rule for selecting lambda makes
# special sense when the CV function is quite flat around the minimum (hence an
# overpenalization that gives more sparsity does not affect so much the CV loss)

# Leave-one-out cross-validation. More computationally intense but completely
# objective in the choice of the fold-assignment
ncv_ridge <- cv.glmnet(x = x, y = y, alpha = 0, nfolds = nrow(Hitters),
                       lambda = lambda_grid)

# Location of both optimal lambdas in the CV loss function
plot(ncv_ridge)

# By default, cv.glmnet randomly allocates observations to folds. Passing foldid
# uses the supplied partition instead to enable running cv.glmnet several times
# on the same folds so that the CV errors are directly comparable. Must be a
# vector of length n with integer fold labels in 1:nfolds
foldid <- sample(rep(1:10, length.out = nrow(Hitters)))
kcv_ridge_foldid <- cv.glmnet(x = x, y = y, alpha = 0, foldid = foldid,
                              lambda = lambda_grid)
kcv_lasso_foldid <- cv.glmnet(x = x, y = y, alpha = 1, foldid = foldid,
                              lambda = lambda_grid)

# Minimum CV errors are comparable because they come from the same folds
min(kcv_ridge_foldid$cvm)
min(kcv_lasso_foldid$cvm)

## ----ridge-3------------------------------------------------------------------
# Inspect the best models (the glmnet fit is inside the output of cv.glmnet)
plot(kcv_ridge2$glmnet.fit, label = TRUE, xvar = "lambda")
abline(v = -log(c(kcv_ridge2$lambda.min, kcv_ridge2$lambda.1se)))

# The fit associated with lambda.1se (or any other lambda not included in the
# original path solution -- obtained by an interpolation) can be retrieved with
# coef() (includes intercept also)
coef(kcv_ridge2, s = kcv_ridge2$lambda.1se)

# Alternatively, one can use
predict(kcv_ridge2, type = "coefficients", s = kcv_ridge2$lambda.1se)

# Predictions for the first two observations
predict(kcv_ridge2, type = "response", s = kcv_ridge2$lambda.1se,
        newx = x[1:2, ])

# Predictions for the first observation, for all the lambdas. We can see how
# the prediction for one observation changes according to lambda
plot(-log(kcv_ridge2$lambda),
     predict(kcv_ridge2, type = "response", newx = x[1, , drop = FALSE],
             s = kcv_ridge2$lambda),
     type = "l", xlab = "-log(lambda)", ylab = " Prediction")

## ----ridge-4------------------------------------------------------------------
# Random data
p <- 5
n <- 200
beta <- seq(-1, 1, l = p)
set.seed(123124)
x <- matrix(rnorm(n * p), n, p)
y <- 1 + x %*% beta + rnorm(n)

# Mimic internal standardization of y done in glmnet, which affects the scale
# of lambda in the regularization
y <- scale(y, center = TRUE, scale = TRUE) * sqrt(n / (n - 1))

# Unrestricted fit
fit <- glmnet(x, y, alpha = 0, lambda = 0, intercept = TRUE,
              standardize = FALSE)
beta0_hat <- rbind(fit$a0, fit$beta)
beta0_hat

# Unrestricted fit matches least squares -- but recall glmnet uses an
# iterative method so it is inexact (convergence threshold thresh = 1e-7 by
# default)
X <- model.matrix(y ~ x) # A way of constructing a design matrix that is a
# data.frame and has a column of ones
solve(crossprod(X)) %*% t(X) %*% y

# Restricted fit
# glmnet considers as the regularization parameter "lambda" the value
# lambda / n (lambda being here the penalty parameter employed in the theory)
lambda <- 2
fit <- glmnet(x, y, alpha = 0, lambda = lambda / n, intercept = TRUE,
              standardize = FALSE, thresh = 1e-10)
beta_lambda_hat <- rbind(fit$a0, fit$beta)
beta_lambda_hat

# Analytical form with intercept
solve(crossprod(X) + diag(c(0, rep(lambda, p)))) %*% t(X) %*% y

## ----lasso-1------------------------------------------------------------------
# Get the Hitters data back
Hitters <- na.omit(Hitters)
x <- model.matrix(Salary ~ ., data = Hitters)[, -1]
y <- Hitters$Salary

# Call to the main function -- use alpha = 1 for lasso regression (the default)
lasso_mod <- glmnet(x = x, y = y, alpha = 1)
# Same defaults as before, same object structure

# Plot of the solution path -- now the paths are not smooth when decreasing to
# zero (they are zero exactly). This is a consequence of the l1 norm
plot(lasso_mod, xvar = "lambda", label = TRUE)

# Some persistently important predictors are 15, 14, and 19
colnames(x)[c(15, 14, 19)]
# Dummies associated to league and division are important

# Versus the R^2 -- same maximum R^2 as before
plot(lasso_mod, label = TRUE, xvar = "dev")

# Now the l1-norm of the coefficients decreases as lambda increases
plot(-log(lasso_mod$lambda), colSums(abs(lasso_mod$beta)), type = "l",
     xlab = "-log(lambda)", ylab = "l1 norm")

# 10-fold cross-validation. Change the seed for a different result
set.seed(12345)
kcv_lasso <- cv.glmnet(x = x, y = y, alpha = 1, nfolds = 10)

# The lambda that minimizes the CV error
kcv_lasso$lambda.min

# The "one standard error rule" for lambda
kcv_lasso$lambda.1se

# Location of both optimal lambdas in the CV loss function
ind_min <- which.min(kcv_lasso$cvm)
plot(kcv_lasso)
abline(h = kcv_lasso$cvm[ind_min] + c(0, kcv_lasso$cvsd[ind_min]))
# No problems now: the minimum does not occur at one extreme
# Interesting: note that the numbers on top of the figure give the number of
# coefficients *exactly* different from zero -- the number of predictors
# effectively considered in the model!
# In this case, the one standard error rule makes also sense

# Leave-one-out cross-validation
lambda_grid <- 10^seq(log10(kcv_lasso$lambda[1]), log10(0.1),
                      length.out = 150) # log-spaced grid
ncv_lasso <- cv.glmnet(x = x, y = y, alpha = 1, nfolds = nrow(Hitters),
                       lambda = lambda_grid)

# Location of both optimal lambdas in the CV loss function
plot(ncv_lasso)

## ----lasso-2------------------------------------------------------------------
# Inspect the best models
plot(kcv_lasso$glmnet.fit, label = TRUE, xvar = "lambda")
abline(v = -log(c(kcv_lasso$lambda.min, kcv_lasso$lambda.1se)))

# The model associated with lambda.min (or any other lambda not included in the
# original path solution -- obtained by an interpolation) can be retrieved with
coef(kcv_lasso, s = c(kcv_lasso$lambda.min, kcv_lasso$lambda.1se))

# Predictions for the first two observations
predict(kcv_lasso, type = "response",
        s = c(kcv_lasso$lambda.min, kcv_lasso$lambda.1se),
        newx = x[1:2, ])

## ----lasso-3------------------------------------------------------------------
# We can use lasso for model selection!
sel_preds <- coef(kcv_lasso, s = c(kcv_lasso$lambda.min,
                                   kcv_lasso$lambda.1se))[-1, ] != 0
x1 <- x[, sel_preds[, 1]]
x2 <- x[, sel_preds[, 2]]

# Least squares fit with variables selected by lasso
mod_lasso_sel1 <- lm(y ~ x1)
mod_lasso_sel2 <- lm(y ~ x2)
summary(mod_lasso_sel1)
summary(mod_lasso_sel2)

# Comparison with stepwise selection
mod_bic <- step(lm(Salary ~ ., data = Hitters), k = log(nrow(Hitters)),
                trace = 0)
summary(mod_bic)
# The lasso variable selection is similar, although the model is slightly worse
# in terms of adjusted R^2 and significance of the predictors. However, keep in
# mind that lasso is solving a constrained least squares problem, so it is
# expected to achieve better R^2 and adjusted R^2 via a selection procedure
# that employs solutions of unconstrained least squares. What is remarkable
# is the speed of lasso on selecting variables, and the fact that gives quite
# good starting points for performing further model selection

# Another interesting possibility is to run a stepwise selection starting from
# the set of predictors selected by lasso. In this search, it is important to
# use direction = "both" (default) and define the scope argument adequately.
# Work on a data frame whose columns ARE the design-matrix predictors (dummies
# included) so that lasso-selected dummy names like "LeagueN" match actual
# columns, and so that factor levels with spaces or unusual characters are
# sanitized into valid R names by data.frame()
hitters_dummies <- data.frame(Salary = y, x)
sel_names <- names(hitters_dummies)[-1][sel_preds[, 2]]
f <- reformulate(termlabels = sel_names, response = "Salary")
start <- lm(f, data = hitters_dummies) # Model with predictors selected by lasso
scope <- list(lower = ~ 1, # No predictors
              upper = terms(Salary ~ ., data = hitters_dummies)) # All preds
mod_bic_from_lasso <- step(object = start, k = log(nrow(hitters_dummies)),
                           scope = scope, trace = 0)
summary(mod_bic_from_lasso)

# Comparison in terms of BIC, slight improvement with mod_bic_from_lasso
BIC(mod_lasso_sel1, mod_lasso_sel2, mod_bic_from_lasso, mod_bic)

## ----lasso-4, fig.cap = '(ref:lasso-4-title)'---------------------------------
# Random data with predictors unrelated to the response
p <- 100
n <- 300
set.seed(123124)
x <- matrix(rnorm(n * p), n, p)
y <- 1 + rnorm(n)

# CV
lambda_grid <- exp(seq(-10, 3, l = 200))
plot(cv.glmnet(x = x, y = y, alpha = 1, nfolds = n, lambda = lambda_grid))

## ----constr-1-----------------------------------------------------------------
# Simulate data
set.seed(123456)
n <- 50
p <- 3
x1 <- rnorm(n, mean = 1)
x2 <- rnorm(n, mean = 2)
x3 <- rnorm(n, mean = 3)
eps <- rnorm(n, sd = 0.5)
y <- 1 + 2 * x1 - 3 * x2 + x3 + eps

# Center the data and compute design matrix
x1_cen <- x1 - mean(x1)
x2_cen <- x2 - mean(x2)
x3_cen <- x3 - mean(x3)
y_cen <- y - mean(y)
X <- cbind(x1_cen, x2_cen, x3_cen)

# Linear restriction: use that
# beta_1 + beta_2 + beta_3 = 0
# beta_2 = -3
# In this case q = 2. The restriction is codified as
A <- rbind(c(1, 1, 1),
           c(0, 1, 0))
c <- c(0, -3)

# Fit model without intercept
S <- solve(crossprod(X))
beta_hat <- S %*% t(X) %*% y_cen
beta_hat

# Restricted fit enforcing A * beta = c
beta_hat_A <- beta_hat +
  S %*% t(A) %*% solve(A %*% S %*% t(A)) %*% (c - A %*% beta_hat)
beta_hat_A

# Intercept of the constrained fit
beta_hat_A_0 <- mean(y) - c(mean(x1), mean(x2), mean(x3)) %*% beta_hat_A
beta_hat_A_0

## ----multr-1------------------------------------------------------------------
# Dimensions and sample size
p <- 3
q <- 2
n <- 100

# A quick way of creating a non-diagonal (valid) covariance matrix for the
# errors
Sigma <- 3 * toeplitz(seq(1, 0.1, l = q))
set.seed(12345)
X <- mvtnorm::rmvnorm(n = n, mean = 1:p, sigma = diag(0.5, nrow = p, ncol = p))
E <- mvtnorm::rmvnorm(n = n, mean = rep(0, q), sigma = Sigma)

# Linear model
B <- matrix((-1)^(1:p) * (1:p), nrow = p, ncol = q, byrow = TRUE)
Y <- X %*% B + E

# Fitting the model (note: Y and X are matrices!)
mod <- lm(Y ~ X)
mod
# Note that the intercept is markedly different from zero -- that is because
# X is not centered

# Compare with B
B

# Summary of the model: gives q separate summaries, one for each fitted
# univariate model
summary(mod)

# Exactly equivalent to
summary(lm(Y[, 1] ~ X))
summary(lm(Y[, 2] ~ X))

## ----multr-2------------------------------------------------------------------
# When we want to add several variables of a dataset as responses through a
# formula interface, we have to use cbind() in the response. Doing
# "Petal.Width + Petal.Length ~ ..." is INCORRECT, as lm will understand
# "I(Petal.Width + Petal.Length) ~ ..." and do one single regression

# Predict Petal's measurements from Sepal's
mod_iris <- lm(cbind(Petal.Width, Petal.Length) ~
                Sepal.Length + Sepal.Width + Species, data = iris)
summary(mod_iris)

# The fitted values and residuals are now matrices
head(mod_iris$fitted.values)
head(mod_iris$residuals)

# The individual models
mod_iris1 <- lm(Petal.Width ~ Sepal.Length + Sepal.Width + Species, data = iris)
mod_iris2 <- lm(Petal.Length ~ Sepal.Length + Sepal.Width + Species, data = iris)
summary(mod_iris1)
summary(mod_iris2)

## ----multr-3------------------------------------------------------------------
# Confidence intervals for the parameters
confint(mod_iris)
# Warning! Do not confuse Petal.Width:Sepal.Length with an interaction term!
# It is meant to represent the Response:Predictor coefficient

# Prediction -- now more limited without confidence intervals implemented
predict(mod_iris, newdata = iris[1:3, ])

# MANOVA table
manova(mod_iris)

# "Same" as the "Sum Sq" and "Df" entries of
anova(mod_iris1)
anova(mod_iris2)

# anova() serves for assessing the significance of including a new predictor
# for explaining all the responses. This is based on an extension of the
# *sequential* ANOVA table briefly covered in Section 2.6. The hypothesis test
# is by default conducted with the Pillai statistic (an extension of the F-test)
anova(mod_iris)

## ----multr-4------------------------------------------------------------------
# Simulate data
n <- 500
p <- 50
q <- 10
set.seed(123456)
X <- mvtnorm::rmvnorm(n = n, mean = p:1, sigma = 5 * 0.5^toeplitz(1:p))
E <- mvtnorm::rmvnorm(n = n, mean = rep(0, q), sigma = toeplitz(q:1))
B <- 5 * (2 / (0.5 * (1:p - 15)^2 + 2) + 1 / (0.1 * (1:p - 40)^2 + 1)) %*%
  t(1 / sqrt(1:q))
Y <- X %*% B + E

# Visualize B -- dark violet is close to 0
image(1:q, 1:p, t(B), col = viridisLite::viridis(20), xlab = "q", ylab = "p")

# Lasso path fit
mfit <- glmnet(x = X, y = Y, family = "mgaussian", alpha = 1)

# A list of models for each response
str(mfit$beta, 1)

# Tuning parameter selection by 10-fold cross-validation
set.seed(12345)
kcv_lasso_m <- cv.glmnet(x = X, y = Y, family = "mgaussian", alpha = 1)
kcv_lasso_m$lambda.min
kcv_lasso_m$lambda.1se

# Location of both optimal lambdas in the CV loss function
ind_min <- which.min(kcv_lasso_m$cvm)
plot(kcv_lasso_m)
abline(h = kcv_lasso_m$cvm[ind_min] + c(0, kcv_lasso_m$cvsd[ind_min]))

# Extract the coefficients associated with some fits
coefs <- coef(kcv_lasso_m, s = c(kcv_lasso_m$lambda.min, kcv_lasso_m$lambda.1se))
str(coefs, 1)

# Predictions for the first two observations
preds <- predict(kcv_lasso_m, type = "response",
                 s = c(kcv_lasso_m$lambda.min, kcv_lasso_m$lambda.1se),
                 newx = X[1:2, ])
preds

## ----multr-5, eval = FALSE----------------------------------------------------
# manipulate::manipulate({
# 
#   # Color
#   col <- viridisLite::viridis(20)
# 
#   # Common zlim
#   zlim <- range(B) + c(-0.25, 0.25)
# 
#   # Plot true B
#   par(mfrow = c(1, 2))
#   image(1:q, 1:p, t(B), col = col, xlab = "q", ylab = "p", zlim = zlim,
#         main = "B")
# 
#   # Extract B_hat from the lasso fit, a p x q matrix
#   B_hat <- sapply(seq_along(mfit$beta), function(i) mfit$beta[i][[1]][, j])
# 
#   # Put as black rows the predictors included
#   not_zero <- abs(B_hat) > 0
#   image(1:q, 1:p, t(not_zero), breaks = c(0.5, 1),
#         col = rgb(1, 1, 1, alpha = 0.1), add = TRUE)
# 
#   # For B_hat
#   image(1:q, 1:p, t(B_hat), col = col, xlab = "q", ylab = "p", zlim = zlim,
#         main = "Bhat")
#   image(1:q, 1:p, t(not_zero), breaks = c(0.5, 1),
#         col = rgb(1, 1, 1, alpha = 0.1), add = TRUE)
# 
# }, j = manipulate::slider(min = 1, max = ncol(mfit$beta$y1), step = 1,
#                           label = "j in lambda(j)"))

## ----biglm-1------------------------------------------------------------------
# Not really "big data", but for the sake of illustration
set.seed(12345)
n <- 1e6
p <- 10
beta <- seq(-1, 1, length.out = p)^5
x1 <- matrix(rnorm(n * p), nrow = n, ncol = p)
x1[, p] <- 2 * x1[, 1] + rnorm(n, sd = 0.1) # Add some dependence to predictors
x1[, p - 1] <- 2 - x1[, 2] + rnorm(n, sd = 0.5)
y1 <- 1 + x1 %*% beta + rnorm(n)
x2 <- matrix(rnorm(100 * p), nrow = 100, ncol = p)
y2 <- 1 + x2 %*% beta + rnorm(100)
big_data1 <- data.frame("resp" = y1, "pred" = x1)
big_data2 <- data.frame("resp" = y2, "pred" = x2)

# biglm has a very similar syntax to lm -- but the formula interface does not
# work always as expected
# biglm::biglm(formula = resp ~ ., data = bigData1) # Does not work
# biglm::biglm(formula = y ~ x) # Does not work
# biglm::biglm(formula = resp ~ pred.1 + pred.2, data = bigData1) # Does work,
# but not very convenient for a large number of predictors
# Hack for automatic inclusion of all the predictors
f <- formula(paste("resp ~", paste(names(big_data1)[-1], collapse = " + ")))
biglm_mod <- biglm::biglm(formula = f, data = big_data1)

# lm's call
lm_mod <- lm(formula = resp ~ ., data = big_data1)

# The reduction in size of the resulting object is more than notable
print(object.size(biglm_mod), units = "KB")
print(object.size(lm_mod), units = "MB")

# Summaries
s1 <- summary(biglm_mod)
s2 <- summary(lm_mod)
s1
s2

# Further information
s1$mat # Coefficients and their inferences
s1$rsq # R^2
s1$nullrss # SST (as in Section 2.6)

# Extract coefficients
coef(biglm_mod)

# Prediction works as usual
predict(biglm_mod, newdata = big_data2[1:5, ])

# Must contain a column for the response
# predict(biglmMod, newdata = bigData2[1:5, -1]) # Error

# Update the model with training data
update(biglm_mod, moredata = big_data2)

# AIC and BIC
AIC(biglm_mod, k = 2)
AIC(biglm_mod, k = log(n))

# Features not immediately available for biglm objects: stepwise selection by
# step, residuals, variance of the error, model diagnostics, and vifs

# Workaround for obtaining hat(sigma)^2 = SSE / (n - p - 1), SSE = SST * (1 - R^2)
(s1$nullrss * (1 - s1$rsq)) / s1$obj$df.resid
s2$sigma^2

## ----biglm-2, fig.cap = '(ref:biglm-2-title)', fig.margin = FALSE-------------
# Model selection adapted to big data models
reg <- leaps::regsubsets(biglm_mod, nvmax = p, method = "exhaustive")
plot(reg) # Plot best model (top row) to worst model (bottom row)

# Summarize (otherwise regsubsets's output is hard to decipher)
subs <- summary(reg)
subs

# Lots of useful information
str(subs, 1)

# Get the model with lowest BIC
subs$which
subs$bic
subs$which[which.min(subs$bic), ]

# Show the display in Figure 4.6
subs$which[order(subs$bic), ]

# It also works with ordinary linear models and it is much faster and
# informative than step
reg <- leaps::regsubsets(resp ~ ., data = big_data1, nvmax = p,
                         method = "backward")
subs <- summary(reg)
subs$bic
subs$which[which.min(subs$bic), ]

# Compare it with step
step(lm(resp ~ ., data = big_data1), trace = 0,
     direction = "backward", k = log(n))

## ----biglm-3------------------------------------------------------------------
# Size of the response
print(object.size(rnorm(1e6)) * 1e2, units = "GB")

# Size of the predictors
print(object.size(rnorm(1e6)) * 1e2 * 10, units = "GB")

## ----biglm-4------------------------------------------------------------------
# Linear regression with n = 10^8 and p = 10
n <- 10^8
p <- 10
beta <- seq(-1, 1, length.out = p)^5

# Number of chunks for splitting the dataset
n_chunks <- 1e3
n_small <- n / n_chunks

# Simulates reading the first chunk of data
set.seed(12345)
x <- matrix(rnorm(n_small * p), nrow = n_small, ncol = p)
x[, p] <- 2 * x[, 1] + rnorm(n_small, sd = 0.1)
x[, p - 1] <- 2 - x[, 2] + rnorm(n_small, sd = 0.5)
y <- 1 + x %*% beta + rnorm(n_small)

# First fit
big_mod <- biglm::biglm(y ~ x, data = data.frame(y, x))

# Update fit
# pb <- txtProgressBar(style = 3)
for (i in 2:n_chunks) {

  # Simulates reading the i-th chunk of data
  set.seed(12345 + i)
  x <- matrix(rnorm(n_small * p), nrow = n_small, ncol = p)
  x[, p] <- 2 * x[, 1] + rnorm(n_small, sd = 0.1)
  x[, p - 1] <- 2 - x[, 2] + rnorm(n_small, sd = 0.5)
  y <- 1 + x %*% beta + rnorm(n_small)

  # Update the fit
  big_mod <- update(big_mod, moredata = data.frame(y, x))

  # Progress
  # setTxtProgressBar(pb = pb, value = i / nChunks)

}

# Final model
summary(big_mod)
print(object.size(big_mod), units = "KB")