R/renet.R

In CorReg: Linear Regression Based on Linear Structure Between Variables

# ' correction for enet not working when X has only one variable
renet <- function(x, y, lambda = 0, max.steps, normalize = TRUE, intercept = TRUE,
                  trace = FALSE, eps = .Machine$double.eps)
{
   call <- match.call()
   nm <- dim(x)
   n <- nm[1]
   m <- nm[2]
   im <- seq(m)
   one <- rep(1, n)
   vn <- dimnames(x)[[2]]
   meanx <- drop(one %*% x) / n
   if (intercept == FALSE) {
      meanx <- rep(0, m)
   }
   x <- scale(x, meanx, FALSE)
   normx <- sqrt(drop(one %*% (x^2)))
   if (normalize == FALSE) {
      normx <- rep(1, m)
   }
   if (any(normx < eps * sqrt(n)))
      stop("Some of the columns of x have zero variance")
   names(normx) <- NULL
   x <- scale(x, FALSE, normx)
   mu <- mean(y)
   if (intercept == FALSE) {
      mu <- 0
   }
   y <- drop(y - mu)
   d1 <- sqrt(lambda)
   d2 <- 1 / sqrt(1 + lambda)
   Cvec <- drop(t(y) %*% x) * d2
   ssy <- sum(y^2)
   residuals <- c(y, rep(0, m))
   if (lambda > 0) {
      maxvars <- m
   }
   if (lambda == 0) {
      maxvars <- min(m, n - 1)
   }
   if (missing(max.steps)) {
      max.steps <- 50 * maxvars
   }
   L1norm <- 0
   penalty <- max(abs(Cvec))
   beta <- rep(0, m)
   betactive <- list(NULL)
   first.in <- integer(m)
   active <- NULL
   Actset <- list(NULL)
   df <- 0
   if (lambda != 0) {
      Cp <- ssy
   }
   ignores <- NULL
   actions <- as.list(seq(max.steps))
   drops <- FALSE
   Sign <- NULL
   R <- NULL
   k <- 0
   while ((k < max.steps) & (length(active) < maxvars)) {
      action <- NULL
      k <- k + 1
      inactive <- if (k == 1) {
         im
      } else {im[-c(active, ignores)]}
      C <- Cvec[inactive]
      Cmax <- max(abs(C))
      if (!any(drops)) {
         new <- abs(C) == Cmax
         C <- C[!new]
         new <- inactive[new]
         for (inew in new) {
            R <- updateRR(x[, inew], R, x[, active], lambda)
            if (attr(R, "rank") == length(active)) {
               nR <- seq(length(active))
               R <- R[nR, nR, drop = FALSE]
               attr(R, "rank") <- length(active)
               ignores <- c(ignores, inew)
               action <- c(action, -inew)
               if (trace)
                  message("LARS-EN Step", k, ":\t Variable", inew,
                          "\tcollinear; dropped for good\n")
            } else {
               if (first.in[inew] == 0)
                  first.in[inew] <- k
               active <- c(active, inew)
               Sign <- c(Sign, sign(Cvec[inew]))
               action <- c(action, inew)
               if (trace)
                  message("LARS-EN Step", k, ":\t Variable", inew,
                          "\tadded\n")
            }
         }
      } else action <- -dropid
      Gi1 <- backsolve(R, backsolvet(R, Sign))
      A <- 1 / sqrt(sum(Gi1 * Sign))
      w <- A * Gi1
      u1 <- drop(x[, active, drop = FALSE] %*% w * d2)
      u2 <- rep(0, m)
      u2[active] <- d1 * d2 * w
      u <- c(u1, u2)
      if (lambda > 0) {
         maxvars <- m - length(ignores)
      }
      if (lambda == 0) {
         maxvars <- min(m - length(ignores), n - 1)
      }
      if (length(active) >= maxvars) {
         gamhat <- Cmax / A
      } else {
         a <- (drop(u1 %*% x[, -c(active, ignores)]) + d1 *
                  u2[-c(active, ignores)]) * d2
         gam <- c((Cmax - C) / (A - a), (Cmax + C) / (A + a))
         gamhat <- min(gam[gam > eps], Cmax / A)
         Cdrop <- c(C - gamhat * a, -C + gamhat * a) - (Cmax -
                                                           gamhat * A)
      }
      dropid <- NULL
      b1 <- beta[active]
      z1 <- -b1 / w
      zmin <- min(z1[z1 > eps], gamhat)
      if (zmin < gamhat) {
         gamhat <- zmin
         drops <- z1 == zmin
      } else drops <- FALSE
      beta[active] <- beta[active] + gamhat * w
      betactive[[k]] <- beta[active]
      Actset[[k]] <- active
      residuals <- residuals - (gamhat * u)
      Cvec <- (drop(t(residuals[1:n]) %*% x) + d1 * residuals[-(1:n)]) *
         d2
      L1norm <- c(L1norm, sum(abs(beta[active])) / d2)
      penalty <- c(penalty, penalty[k] - abs(gamhat * A))
      if (any(drops)) {
         dropid <- seq(drops)[drops]
         for (id in rev(dropid)) {
            if (trace)
               message("LARS-EN Step", k, ":\t Variable", active[id],
                       "\tdropped\n")
            R <- downdateR(R, id)
         }
         dropid <- active[drops]
         beta[dropid] <- 0
         active <- active[!drops]
         Sign <- Sign[!drops]
      }
      if (!is.null(vn))
         names(action) <- vn[abs(action)]
      actions[[k]] <- action
   }
   allset <- Actset[[1]]
   if (k >= 2) {
      for (i in 2:k) {
         allset <- union(allset, Actset[[i]])
      }
   }
   
   
   allset <- sort(allset)
   max.p <- length(allset)
   beta.pure <- matrix(0, k + 1, max.p)
   for (i in 2:(k + 1)) {
      for (j in 1:length(Actset[[i - 1]])) {
         l <- c(1:max.p)[allset == Actset[[i - 1]][j]]
         beta.pure[i, l] <- betactive[[i - 1]][j]
      }
   }
   beta.pure <- beta.pure / d2
   dimnames(beta.pure) <- list(paste(0:k), vn[allset])
   k <- dim(beta.pure)[1]
   df <- 1:k
   for (i in 1:k) {
      a <- drop(beta.pure[i, ])
      df[i] <- 1 + length(a[a != 0])
   }
   residuals <- y - x[, allset, drop = FALSE] %*% t(beta.pure)
   beta.pure <- scale(beta.pure, FALSE, normx[allset])
   RSS <- apply(residuals^2, 2, sum)
   R2 <- 1 - RSS / RSS[1]
   Cp <- ((n - m - 1) * RSS) / rev(RSS)[1] - n + 2 * df
   object <- list(call = call, actions = actions[seq(k)], allset = allset,
                  beta.pure = beta.pure, vn = vn, mu = mu, normx = normx[allset],
                  meanx = meanx[allset], lambda = lambda, L1norm = L1norm,
                  penalty = penalty * 2 / d2, df = df, Cp = Cp, sigma2 = rev(RSS)[1] / (n -
                                                                                           m - 1))
   class(object) <- "enet"
   object
}