R/ridge.R

In martindurocher/floodRFA: A R-package for regional frequency of floods

#' Perform ridge regression
#'
#' Wrapper for the ridge regression function in MASS package to use with
#'  predict and residuals function. Prediction are obtained using the best
#'  regularization values.
#'
#' @param form      Formula of the model
#' @param x,xnew    Data frame of the data
#' @param lambda    List of regularization value
#' @param optim     Criteria to choose the best lambda
#'
#' @export
#'
#' @examples
#'
#'  # Example from MASS package
#' longley # not the same as the S-PLUS dataset
#' names(longley)[1] <- "y"
#'
#' fit <- FitRidge(y ~ ., longley, lambda = seq(0,0.1,0.001))
#'
#' plot(fit)
#' select(fit)
#' plot(predict(fit), residuals(fit))
#' abline( h = 0)
#'
FitRidge <- function(form, x, lambda = seq(0,0.1,0.001)){
  x <- model.frame(form,x)
  ans <- MASS::lm.ridge(form, x, lambda = lambda)
  ans$data <- x
  ans$formula <- form
  ans$optim <- which.min(ans$GCV)

  return(ans)
}

#' @export
#' @rdname FitRidge
predict.ridgelm <- function(obj, xnew = NULL, optim = 'gcv'){

  #Choose the optimal
  if(optim == 'gcv')
    oid <- obj$optim
  else if(optim == 'kHKB')
    oid <- which(abs(obj$lambda-obj$kHKB))
  else if(optim == 'kLW')
    oid <- which(abs(obj$lambda-obj$kLW))


  if(is.null(xnew))
    xnew <- model.matrix(obj$form,obj$data)
  else
    xnew <- model.matrix(obj$form,xnew)

  for(jj in 2:ncol(xnew))
    xnew[,jj] <- (xnew[,jj]-obj$xm[jj-1])/obj$scales[jj-1]

  beta <- c(obj$ym,obj$coef[,oid])
  as.numeric(xnew %*% beta)
}

#' @export
#' @rdname FitRidge
residuals.ridgelm <- function(obj, ...){
  yhat <- predict(obj,obj$data, ...)
  as.numeric(obj$data[,1] - yhat)
}