#' 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)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.