Nothing
#' Ellipses Plot - Visual Correlation Matrix Comparison
#'
#' @description When there are 3 or more variables in the data, this function produces
#' a matrix with ellipses drawn in the upper triangle. The ellipse in cell
#' \eqn{i,j} of the plot is drawn to be a contour of a standard bivariate
#' normal with correlation \eqn{\rho_{ij}}. One ellipse is drawn in each
#' cell for each model in the \code{covfm} object. When there are 2
#' variables in the data, this function produces a scatter plot of the data
#' with an overlaid 95\% confidence ellipse for each model in the
#' \code{covfm} object.
#'
#' @param x a \code{"covfm"} object.
#' @param \dots additional arguments are ignored.
#'
#' @return x is invisibly returned.
#'
#' @importFrom graphics legend lines plot polygon strheight strwidth text
#'
#' @export
ellipsesPlot.covfm <- function(x, ...)
{
n.models <- length(x)
mod.names <- names(x)
p <- dim(vcov(x[[1]]))[1]
## if p == 2 plot data with overlaid ellipse ##
if(p == 2) {
old.par <- par(pty = "s")
on.exit(par(old.par))
ellipse <- function(loc, A)
{
detA <- A[1, 1] * A[2, 2] - A[1, 2]^2
dist <- sqrt(qchisq(0.95, 2))
ylimit <- sqrt(A[2, 2]) * dist
y <- seq(-ylimit, ylimit, 0.01 * ylimit)
sqrt.discr <- detA/A[2,2]^2 * (A[2,2] * dist^2 - y^2)
sqrt.discr[c(1, length(sqrt.discr))] <- 0.0
sqrt.discr <- sqrt(sqrt.discr)
b <- loc[1] + A[1, 2] / A[2, 2] * y
x1 <- b - sqrt.discr
x2 <- b + sqrt.discr
y <- loc[2] + y
rbind(cbind(x1, y), cbind(rev(x2), rev(y)))
}
z <- list()
x.min <- Inf
x.max <- -Inf
y.min <- Inf
y.max <- -Inf
for(i in 1:n.models) {
z[[i]] <- ellipse(center(x[[i]]), vcov(x[[i]]))
x.min <- min(x.min, z[[i]][,1])
x.max <- max(x.max, z[[i]][,1])
y.min <- min(y.min, z[[i]][,2])
y.max <- max(y.max, z[[i]][,2])
}
X <- try(eval(x[[1]]$call$data, envir = sys.parent(2)), silent = TRUE)
if(!is.null(X) && !inherits(X, "try-error")) {
X <- as.matrix(X)
x.min <- min(x.min, X[,1])
x.max <- max(x.max, X[,1])
y.min <- min(y.min, X[,2])
y.max <- max(y.max, X[,2])
} else {
X <- matrix(NA, 1, 2)
}
center <- c(mean(c(x.min, x.max)), mean(c(y.min, y.max)))
s.range <- max(abs(c(center[1] - x.min, x.max - center[1],
center[2] - y.min, y.max - center[2])))
if(n.models == 1)
header <- "95% Tolerance Ellipse"
else
header <- "95% Tolerance Ellipses"
plot(X,
xlim = c(center[1] - s.range, center[1] + s.range),
ylim = c(center[2] - s.range, center[2] + s.range),
xlab = names(center(x[[1]]))[1],
ylab = names(center(x[[1]]))[2],
main = header,
col = "lightgray",
pch = 16,
asp = 1)
for(i in 1:length(z))
polygon(z[[i]], density = 0, lty = i, col = i, lwd = 2)
pos <- ifelse(vcov(x[[1]])[1,2] > 0, "topleft", "topright")
legend(x = pos, legend = mod.names, col = 1:n.models, lty = 1:n.models,
lwd = 2, bty = "n")
}
## if p > 2 plot matrix of ellipses ##
else {
old.par <- par(mar = rep(0.1, 4), pty = "s")
on.exit(par(old.par))
labels <- dimnames(vcov(x[[1]]))[[1]]
tl.margin <- max(strwidth(paste("WW", labels, sep = ""), units = "inches"))
tl.margin <- 1.05*tl.margin / par()$fin[1]
if(tl.margin < 0.05) {
tl.margin <- 0.05
cex.labels <- 1.0
} else if(tl.margin > 0.2) {
cex.labels <- 0.2 / tl.margin
tl.margin <- 0.2
} else {
cex.labels <- 1.0
}
br.margin <- 0.025*(n.models - 1)
plot(NA, NA,
xlim = c(0.5 - p*tl.margin, 0.5 + (1 + br.margin)*p),
ylim = c(0.5 - p*br.margin, (1 + tl.margin)*p + 0.5),
type = "n",
axes = FALSE,
xlab = "",
ylab = "")
ht.corr <- (1.5 * n.models - 0.5) * strheight("8", units = "user")
wt.corr <- strwidth("-0.00", units = "user")
cex.corr <- min(c(0.75 / max(c(ht.corr, wt.corr)), 1.25))
vert <- (1:n.models) * 1.6
vert <- vert - mean(vert) + 0.5
for(k in 1:n.models) {
s <- sqrt(diag(vcov(x[[k]])))
X <- vcov(x[[k]]) / (s %o% s)
xc <- col(X)
yc <- row(X)[p:1,]
ut <- row(X) < col(X)
lt <- row(X) > col(X)
pts <- c(seq(0.0, 2*pi, length.out = 181), NA)
xs <- sapply(X[ut], function(u, v) cos(v + acos(u)/2), v = pts)
xs <- 0.475 * xs + rep(xc[ut], each = 182)
ys <- sapply(X[ut], function(u, v) cos(v - acos(u)/2), v = pts)
ys <- 0.475 * ys + rep(yc[ut], each = 182)
polygon(x = as.vector(xs),
y = as.vector(ys),
density = 0,
lwd = n.models - k + 1,
col = k,
lty = k)
corr <- X[lt]
corr[corr > 0.99 & corr < 1.00] <- 0.99
corr[corr < -0.99 & corr > -1.00] <- -0.99
corr <- format(round(corr, digits = 2))
text(xc[lt] + 0.5*strwidth("0.00", cex = cex.corr),
yc[lt],
labels = corr,
adj = c(1.0, vert[k]),
col = k,
cex = cex.corr)
}
lines(c(1, p), c(p, 1), col = "gray")
text(1:p, p+0.5,
paste(" ", labels, sep = ""),
adj = 0,
srt = 90,
cex = cex.labels)
text(0.5, p:1,
paste(labels, " ", sep = ""),
adj = 1,
cex = cex.labels)
legend(x = "bottomright",
legend = mod.names,
lwd = n.models:1,
col = 1:n.models,
lty = 1:n.models,
bty = "n",
horiz = TRUE)
}
invisible(x)
}
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.