R/better.subset.lm.R

In Accelerytics/mlkit: Machine Learning Kit

#' Fitting Linear Models with Better Subset Algorithm
#'
#' @description Implementation of the MM algorithm solver for the better subset
#' algorithm for a linear regression model.
#'
#' @param formula an object of class \code{\link[stats]{formula}}: a symbolic
#' description of the model to be fitted following the standard of
#' \code{\link[stats]{lm}}.
#' @param data an optional data frame, list or environment (or object coercible
#' by \code{\link[base]{as.data.frame}} to a data frame) containing the
#' variables in the model. If not found in data, the variables are taken from
#' environment (\code{formula}), typically the environment from which this
#' function is called.
#' @param m number of explanatory variables to include in the model.
#' @param intercept optional boolean indicating whether to fit an intercept. If
#' \code{TRUE}, \code{standardize} is ignored. Default is \code{FALSE}.
#' @param standardize optional boolean indicating whether to return results for
#' standardized data. If \code{intercept} is \code{TRUE}, this argument is
#' ignored. Default is \code{FALSE}.
#' \code{TRUE}, \code{standardize} is ignored. Default is \code{FALSE}.
#' @param beta.init optional initial beta parameters to use in the MM
#' algorithm. Default is \code{NULL}.
#' @param beta.tol optional absolute tolerance for rounding down parameter
#' standardized estimates. If the absolute value of a parameter estimate in the
#' standardized model is smaller than \code{beta.tol}, it is rounded down to
#' zero. Default is \code{0}, that is, no rounding.
#' @param loss.tol optional convergence tolerance on the elastic net loss in
#' the MM algorithm. Default is \code{1e-6}.
#' @param seed optional seed. Default is \code{NULL}.
#' @param verbose optional number indicating per how many iterations the
#' estimation progress is displayed. Default is \code{0}, that is, no progress
#' updates.
#'
#' @return \code{better.subset.lm} returns an object of \code{\link{class}}
#' \code{mlkit.lm.fit}. An object of class \code{mlkit.lm.fit} is a list
#' containing at least the following components:
#' \item{coefficients}{a named vector of optimal coefficients.}
#' \item{loss}{residual sum of squares for optimal coefficients.}
#' \item{r2}{coefficient of determination for optimal coefficients.}
#' \item{adj.r2}{adjusted coefficient of determination for optimal
#' coefficients.}
#'
#' @export
#'
better.subset.lm = function(formula, data, m, intercept=T, standardize=F,
  beta.init=NULL, beta.tol=0, loss.tol=1e-6, seed=NULL, verbose=0) {

  # Extract dependent variable and explanatory variables
  y = data.matrix(data[, all.vars(formula)[1]]);
  x = stats::model.matrix(formula, data); rss.tot = sum((y - mean(y)) ^ 2)

  # Store scaling parameters
  if (!standardize) {
    x.mean = colMeans(x); x.sd = apply(x, 2, stats::sd)
    y.mean = mean(y); y.sd = stats::sd(y)
  }

  # Add intercept or standardize data if necessary
  x = create.x(x, intercept, standardize)
  y = create.y(y, intercept, standardize)

  # Define constants
  Xt.X = crossprod(x); inv.lambda = 1 / eigen(Xt.X)$values[1]; set.seed(seed)
  N = nrow(x); P = ncol(x); inv.lambda.Xt.y = inv.lambda * crossprod(x, y)

  # Define helper function for displaying progress
  loss = function(beta) rss(beta, x, y)
  pline = function(i, o, n, d)
    list(c('Iteration', i), c('Loss.old', o), c('Loss.new', n), c('Delta', d))

  # Choose some inital beta_0 and compute initial loss
  b.new = initialize.beta(beta.init, x); l.new = loss(b.new)

  # Update iteration and replace old parameters by previous until convergence
  i = 0L; while (TRUE) { i = i + 1L; l.old = l.new

    # Apply MM-step and update parameters
    u = b.new - inv.lambda * (Xt.X %*% b.new) + inv.lambda.Xt.y
    abs.u = abs(u); abs.u.max = sort(abs.u, decreasing=T)[m]
    b.new = ifelse(abs.u >= abs.u.max, u, 0)
    l.new = loss(b.new); diff = l.old - l.new

    # Display progress if verbose
    if (verbose & (i %% verbose == 0))
      progress.str(pline(i, l.new, l.old, diff))

    # Break if improvement smaller than tol, that is, sufficient convergence
    if (diff / l.old < loss.tol) break
  }

  # Ensure information of last iteration is displayed
  if (verbose & (i %% verbose != 0))
    progress.str(pline(i, l.new, l.old, diff))

  # Force elements smaller than beta.tol to zero
  b.new[abs(b.new) < beta.tol] = 0

  # Descale beta if necessary
  if (intercept | standardize) beta = c(0, b.new)
  else beta = descale.beta(b.new, x.mean, x.sd, y.mean, y.sd)

  # Return mlfit object
  res = list('coefficients'=beta, 'alpha'=0, 'lambda'=0,
    'loss'=loss(b.new), 'R2'=r2(b.new, x, y), 'adj.R2'=adj.r2(b.new, x, y))
  class(res) = 'mlkit.lm.fit'
  return(res)
}