Confidence Regions for a Multivariate Normal Mean

2003-11-26

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 
 
> 

Statistics 4M03/6M03