library(ellipse)
attach(Sweat)
names(Sweat)

par(mfrow=c(3,3))

frame()

conf.panel(SODIUM, SWEAT, 3)
conf.panel(POTASS, SWEAT, 3)
conf.panel(SWEAT, SODIUM, 3)

frame()

conf.panel(POTASS, SODIUM, 3)
conf.panel(SWEAT, POTASS, 3)
conf.panel(SODIUM, POTASS, 3)

frame()


multi.ci <- function(v, p, alpha) {
    n   <- length(v)
    f.q <- qf(alpha, p, n-p)
    sqrt(p*(n-1)/(n-p)*f.q) * sqrt(var(v) / n)
}


bonf.ci <- function(v, p, alpha) {
    n   <- length(v)
    t.q <- qt(1-(1-alpha)/2/p, n-1)
    t.q * sqrt(var(v) / n)
}

conf.panel <- function(x, y, p) {

     n<-length(x)
     
     plot(x, y)

     lines(ellipse(cov(cbind(x, y)), t=sqrt(p*(n-1)/(n-p)*qf(0.95, p, n-p)/n),
   	         centre=c(mean(x), mean(y))))

     abline(v=mean(x)+multi.ci(x, p, 0.95))
     abline(v=mean(x)-multi.ci(x, p, 0.95))

     abline(h=mean(y)+multi.ci(y, p, 0.95))
     abline(h=mean(y)-multi.ci(y, p, 0.95))

     abline(v=mean(x)+bonf.ci(x, p, 0.95))
     abline(v=mean(x)-bonf.ci(x, p, 0.95))

     abline(h=mean(y)+bonf.ci(y, p, 0.95))
     abline(h=mean(y)-bonf.ci(y, p, 0.95))

}