-
Notifications
You must be signed in to change notification settings - Fork 1
/
precision.R
135 lines (129 loc) · 4.95 KB
/
precision.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
#' 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)
}