Nothing
#' Measures of Precision and Shrinkage for Ridge Regression
#'
#' @name precision
#' @aliases precision precision.ridge precision.lm
#'
#' @description
#' Calculates measures of precision based on the size of the estimated
#' covariance matrices of the parameters and shrinkage of the parameters in a
#' ridge regression model. %% ~~ A concise (1-5 lines) description of what the
#' function does. ~~
#'
#' Three measures of (inverse) precision based on the \dQuote{size} of the
#' covariance matrix of the parameters are calculated. Let \eqn{V_k} be the
#' covariance matrix for a given ridge constant, and let \eqn{\lambda_i , i= 1,
#' \dots p} be its eigenvalues
#' \enumerate{
#' \item \eqn{\log | V_k | = \log \prod \lambda} or \eqn{|V_k|^{1/p} =(\prod \lambda)^{1/p}}
#' measures the linearized
#' volume of the covariance ellipsoid and corresponds conceptually to Wilks'
#' Lambda criterion
#' \item \eqn{ trace( V_k ) = \sum \lambda} corresponds conceptually to Pillai's trace criterion
#' \item \eqn{ \lambda_1 = max (\lambda)} corresponds to Roy's largest root criterion.
#' }
#'
#' @param object An object of class \code{ridge} or \code{lm}
#' @param det.fun Function to be applied to the determinants of the covariance
#' matrices, one of \code{c("log","root")}.
#' @param normalize If \code{TRUE} the length of the coefficient vector is
#' normalized to a maximum of 1.0.
#' @param \dots Other arguments (currently unused)
#'
#' @return A data.frame with the following columns
#' \item{lambda}{The ridge constant}
#' \item{df}{The equivalent effective degrees of freedom}
#' \item{det}{The \code{det.fun} function of the determinant of the covariance matrix}
#' \item{trace}{The trace of the covariance matrix}
#' \item{max.eig}{Maximum eigen value of the covariance matrix}
#' \item{norm.beta}{The root mean square of the estimated coefficients, possibly normalized}
#'
#' @note Models fit by \code{lm} and \code{ridge} use a different scaling for
#' the predictors, so the results of \code{precision} for an \code{lm} model
#' will not correspond to those for \code{ridge} with ridge constant = 0.
#'
#' @author Michael Friendly
#' @export
#' @seealso \code{\link{ridge}},
#' @keywords regression models
#' @examples
#'
#' longley.y <- longley[, "Employed"]
#' longley.X <- data.matrix(longley[, c(2:6,1)])
#'
#' lambda <- c(0, 0.005, 0.01, 0.02, 0.04, 0.08)
#' lridge <- ridge(longley.y, longley.X, lambda=lambda)
#' clr <- c("black", rainbow(length(lambda)-1, start=.6, end=.1))
#' coef(lridge)
#'
#' (pdat <- precision(lridge))
#' # plot log |Var(b)| vs. length(beta)
#' with(pdat, {
#' plot(norm.beta, det, type="b",
#' cex.lab=1.25, pch=16, cex=1.5, col=clr, lwd=2,
#' xlab='shrinkage: ||b|| / max(||b||)',
#' ylab='variance: log |Var(b)|')
#' text(norm.beta, det, lambda, cex=1.25, pos=c(rep(2,length(lambda)-1),4))
#' text(min(norm.beta), max(det), "Variance vs. Shrinkage", cex=1.5, pos=4)
#' })
#'
#' # plot trace[Var(b)] vs. length(beta)
#' with(pdat, {
#' plot(norm.beta, trace, type="b",
#' cex.lab=1.25, pch=16, cex=1.5, col=clr, lwd=2,
#' xlab='shrinkage: ||b|| / max(||b||)',
#' ylab='variance: trace [Var(b)]')
#' text(norm.beta, trace, lambda, cex=1.25, pos=c(2, rep(4,length(lambda)-1)))
#' # text(min(norm.beta), max(det), "Variance vs. Shrinkage", cex=1.5, pos=4)
#' })
#'
#'
precision <- function(object, det.fun, normalize, ...) {
UseMethod("precision")
}
# DONE: allow choice of log.det or det()^{1/p}
#' @exportS3Method
precision.ridge <- function(object,
det.fun=c("log","root"),
normalize=TRUE, ...) {
tr <- function(x) sum(diag(x))
maxeig <- function(x) max(eigen(x)$values)
V <- object$cov
p <- ncol(coef(object))
det.fun <- match.arg(det.fun)
ldet <- unlist(lapply(V, det))
ldet <- if(det.fun == "log") log(ldet) else ldet^(1/p)
trace <- unlist(lapply(V, tr))
meig <- unlist(lapply(V, maxeig))
norm <- sqrt(rowMeans(coef(object)^2))
if (normalize) norm <- norm / max(norm)
data.frame(lambda=object$lambda,
df=object$df,
det=ldet,
trace=trace,
max.eig=meig,
norm.beta=norm)
}
#' @exportS3Method
precision.lm <- function(object, det.fun=c("log","root"), normalize=TRUE, ...) {
V <- vcov(object)
beta <- coefficients(object)
if (names(beta[1]) == "(Intercept)") {
V <- V[-1, -1]
beta <- beta[-1]
}
# else warning("No intercept: precision may not be sensible.")
p <- length(beta)
det.fun <- match.arg(det.fun)
ldet <- det(V)
ldet <- if(det.fun == "log") log(ldet) else ldet^(1/p)
trace <- sum(diag(V))
meig <- max(eigen(V)$values)
norm <- sqrt(mean(beta^2))
if (normalize) norm <- norm / max(norm)
res <- list(df=length(beta),
det=ldet,
trace=trace,
max.eig=meig,
norm.beta=norm)
unlist(res)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.