The function confp()
computes the 100(1 - alpha)%
confidence ellipsoid for a multivariate normal mean, and confidence
boxes by the univariate, Bonferroni and shadow methods. The
univariate method gives the smallest box, the shadow method gives the
largest box, just enclosing the ellipsoid.
When the data are bivariate, the ellipse and rectangles are
plotted. If plotdata = T
, the original data appear on a
scatterplot and a concentration ellipse (at the same alpha as the
confidence ellipse) is drawn.
> confp <- function (x, plotdata = F, alpha = 0.05, ...) { n <- nrow(x) p <- ncol(x) if (p >= 2) { crite <- ((p * (n - 1))/(n * (n - p))) * qf(1 - alpha, p, n - p) critm <- qt(1 - (alpha/2), n - 1)/sqrt(n) critb <- qt(1 - (alpha/(2 * p)), n - 1)/sqrt(n) smat <- var(x) mean <- apply(x, 2, mean) x.ei <- eigen(smat) hlaxe <- sqrt(x.ei$values * crite) theta <- angle(x.ei$vectors[1, 1], x.ei$vectors[2, 1]) hlbx <- sqrt(diag(smat) * crite) hlbxm <- sqrt(diag(smat)) * critm hlbxb <- sqrt(diag(smat)) * critb if (p == 2) { if (plotdata) { plot(x[, 1], x[, 2], xlab = dimnames(x)[[2]][1], ylab = dimnames(x)[[2]][2], ...) ellipse(sqrt(n) * hlaxe[1], sqrt(n) * hlaxe[2], theta, mean[1], mean[2]) } else plot(c(mean[1] - hlbx[1], mean[1] + hlbx[1]), c(mean[2] - hlbx[2], mean[2] + hlbx[2]), xlab = dimnames(x)[[2]][1], ylab = dimnames(x)[[2]][2], type = "n") points(mean[1], mean[2], pch = 3) title(main = paste(100 * (1 - alpha), "% Confidence Regions for Two Means")) rect(mean[1] - hlbx[1], mean[2] - hlbx[2], mean[1] + hlbx[1], mean[2] + hlbx[2], lty = 2) rect(mean[1] - hlbxm[1], mean[2] - hlbxm[2], mean[1] + hlbxm[1], mean[2] + hlbxm[2], lty = 4) rect(mean[1] - hlbxb[1], mean[2] - hlbxb[2], mean[1] + hlbxb[1], mean[2] + hlbxb[2], lty = 3) ellipse(hlaxe[1], hlaxe[2], theta, mean[1], mean[2]) } list(mean = mean, smat = smat, hlaxe = hlaxe, axes = x.ei$vectors, theta = theta, hlbx = hlbx, hlbxm = hlbxm, hlbxb = hlbxb) } } > confp(bones[,5:6], T, pch=19, col="forestgreen") $mean ulnadom ulna 0.70440 0.69384 $smat ulnadom ulna ulnadom 0.01156842 0.00807115 ulna 0.00807115 0.01059914 $hlaxe [1] 0.07400143 0.02926561 $axes [,1] [,2] [1,] 0.7279896 0.6855881 [2,] 0.6855881 -0.7279896 $theta [1] 0.7554113 $hlbx ulnadom ulna 0.05748732 0.05502631 $hlbxm ulnadom ulna 0.04439717 0.04249655 $hlbxb ulnadom ulna 0.05143246 0.04923066 >