R/enet-sim.r

In pbreheny/hdrm: High-Dimensional Regression Modeling

#' Reproduce Example 4.2
#'
#' Reproduces Example 4.2 from the book.  If you specify any options, your results may look different.
#'
#' @param N       Number of simulated realizations
#' @param rho     Vector of correlations
#' @param n       Sample size
#' @param p       Number of features
#' @param p1      Number of non-null features
#' @param b       Coefficient associated with the non-null features
#' @param alpha   Elastic net tuning parameter
#' @param corr    Correlation structure, either `cs` (compound symmetric) or `bd` (block diagonal)
#' @param seed    Seed for reproducibility
#'
#' @examples Ex4.2(N=10)
#' @export

Ex4.2 <- function(
    N = 1000,
    rho = seq(0, 0.9, 0.1),
    n = 50,
    p = 100,
    p1 = 5,
    b = 0.5,
    alpha = 0.5,
    corr = c('cs', 'bd'),
    seed=1) {
  original_seed <- .GlobalEnv$.Random.seed
  on.exit(.GlobalEnv$.Random.seed <- original_seed)
  set.seed(seed)

  corr <- match.arg(corr)
  R <- length(rho)
  alpha <- c(1, alpha)
  beta <- c(rep(b, p1), rep(0, 100-p1))

  res <- array(0, dim=c(N, R, 2), dimnames=list(1:N, rho, c("lasso", "enet")))
  pb <- txtProgressBar(0, N, style=3)
  for (i in 1:N) {
    for (j in 1:R) {
      if (corr == 'cs') {
        Data <- gen_data(n, p, rho=rho[j], beta=beta)
        pData <- gen_data(n, p, rho=rho[j], beta=beta)
      } else {
        Data <- gen_data_grp(n, J=20, K=5, rho.g=rho[j], beta=beta)
        pData <- gen_data_grp(n, J=20, K=5, rho.g=rho[j], beta=beta)
      }
      for (k in 1:2) {
        fit <- with(Data, glmnet(X, y, alpha=alpha[k]))
        P <- predict(fit, newx=pData$X)
        ind <- which.min(apply(pData$y-P, 2, crossprod))
        res[i,j,k] <- crossprod(beta - fit$beta[,ind])
      }
    }
    setTxtProgressBar(pb, i)
  }
  close(pb)
  res
}

#' Reproduce Figure 4.1
#'
#' Reproduces Figure 4.1 from the book.  If you specify any options, your results may look different.
#'
#' @param cs        Output from Ex4.2 with compound symmetric structure
#' @param bd        Output from Ex4.2 with block diagonal structure
#' @param parlist   List of arguments to pass to `par()`
#'
#' @examples
#' cs <- Ex4.2(N=5, corr='cs')
#' bd <- Ex4.2(N=5, corr='bd')
#' Fig4.1(cs, bd)
#' @export

Fig4.1 <- function(cs, bd, parlist=list(mfrow=c(1,2), mar=c(5,5,0.5,0.5), oma=c(0,0,2,0))) {
  labs <- c("Lasso", "Elastic Net")
  op <- par(parlist)
  a <- apply(cs, 2:3, mean)
  b <- apply(bd, 2:3, mean)
  matplot(as.numeric(rownames(a)), a, type='l', lty=1, lwd=3, col=pal(2), bty='n', las=1, xlab=expression(rho), ylab="MSE")
  matplot(as.numeric(rownames(b)), b, type='l', lty=1, lwd=3, col=pal(2), bty='n', las=1, xlab=expression(rho), ylab="MSE")
  toplegend(legend=labs, lwd=3, col=pal(2))
  par(op)
}