weighted ridge regression script

linear.R

## weighted ridge regression


linearr = function(X, y, lam = 0, weights = NULL, intercept = TRUE,
                   kernel = FALSE) {

  # checks
  n = dim(X)[1]
  p = dim(X)[2]
  if (is.null(weights)) {
    weights = rep(1, n)
  }
  if (length(weights) != n)
    stop("weights must be length ", n)
  if (length(lam) > 1)
    stop("lam must be a scalar!")
  if (lam < 0)
    stop("lam must be nonnegative!")
  if ((kernel == TRUE) & (lam == 0))
    stop("must specify lam to use kernel!")


  # initialization
  X. = as.matrix(X)
  y. = as.matrix(y)
  if (intercept == TRUE) {
    # if first column of ones, then remove it
    if (all(X.[, 1] == rep(1, n))) {
      X. = X.[, -1]
    }

    # center the data
    X_bar = (as.matrix(t(weights)) %*% X.)/sum(weights)
    y_bar = sum(weights * y.)/sum(weights)

    X. = X. - rep(1, n) %*% X_bar
    y. = y. - y_bar

  }

  root = sqrt(weights)
  X. = root * X.
  y. = root * y.


  # SVD
  svd = svd(X.)

  # adjust d vector for regularization and diagonalize
  d_adj = diag(svd$d/(svd$d^2 + lam))

  # calculate beta estimates
  betas = svd$v %*% d_adj %*% t(svd$u) %*% y.
  rownames(betas) = NULL


  # add intercept, if needed
  if (intercept == TRUE) {

    # calculate intercept
    b1 = y_bar - X_bar %*% betas
    rownames(b1) = c("intercept")
    betas = rbind(b1, betas)

  }

  returns = list(coefficients = betas)
  return(returns)
}