From s200b@stat.Berkeley.EDU Tue Feb 22 13:48:27 2000
Date: Tue, 22 Feb 2000 13:47:36 -0800 (PST)
From: Ben Bolstad <s200b@stat.Berkeley.EDU>
To: bolstad@stat.Berkeley.EDU

outside <- 0
outsideScheffe <- 0
for (i in 1:10000){
norms <- matrix(rnorm(30),10,3)
smean <- apply(norms,2,mean)
svars <- apply(norms,2,var)

s <- sqrt((9*svars[1] + 9*svars[2] + 9*svars[3])/(30-3))

#first use t
# mu1 - mu2

ll1 <- (smean[1] - smean[2]) - qt(0.975,27)*s*sqrt(1/10 + 1/10)
ul1 <- (smean[1] - smean[2]) + qt(0.975,27)*s*sqrt(1/10 + 1/10)

# mu1 - mu3
ll2 <- (smean[1] - smean[3]) - qt(0.975,27)*s*sqrt(1/10 + 1/10)
ul2 <- (smean[1] - smean[3]) + qt(0.975,27)*s*sqrt(1/10 + 1/10)

# mu2 - mu3

ll3 <- (smean[2] - smean[3]) - qt(0.975,27)*s*sqrt(1/10 + 1/10)
ul3 <- (smean[2] - smean[3]) + qt(0.975,27)*s*sqrt(1/10 + 1/10)

#Now Scheffe multiple comparison

ll1a <- (smean[1] - smean[2]) - s*sqrt((3-1)*qf(0.95,2,27)*(1/10+1/10))
ul1a <- (smean[1] - smean[2]) + s*sqrt((3-1)*qf(0.95,2,27)*(1/10+1/10))


ll2a <- (smean[1] - smean[3]) - s*sqrt((3-1)*qf(0.95,2,27)*(1/10+1/10))
ul2a <- (smean[1] - smean[3]) + s*sqrt((3-1)*qf(0.95,2,27)*(1/10+1/10))

ll3a <- (smean[2] - smean[3]) - s*sqrt((3-1)*qf(0.95,2,27)*(1/10+1/10))
ul3a <- (smean[2] - smean[3]) + s*sqrt((3-1)*qf(0.95,2,27)*(1/10+1/10))

outside <- outside + sum(sum(0 < c(ll1,ll2,ll3)) + sum(0 > c(ul1,ul2,ul3)) > 0)
outsideScheffe <- outsideScheffe + sum(sum(0 < c(ll1a,ll2a,ll3a)) + sum(0 > c(ul1a,ul2a,ul3a)) >0)

}
print(outside/10000)
print(outsideScheffe/10000)

#look at some example intervals
print(matrix(c(ll1,ll2,ll3,ul1,ul2,ul3), 3,2))
print(matrix(c(ll1a,ll2a,ll3a,ul1a,ul2a,ul3a), 3,2))

# here are the results of running this simulation

> source("code.S")
[1] 0.1159
[1] 0.0394

> print(matrix(c(ll1,ll2,ll3,ul1,ul2,ul3), 3,2))
           [,1]     [,2]
[1,] -0.9656419 1.141835
[2,] -0.9316393 1.175838
[3,] -1.0197360 1.087741
> print(matrix(c(ll1a,ll2a,ll3a,ul1a,ul2a,ul3a), 3,2))
          [,1]     [,2]
[1,] -1.242040 1.418234
[2,] -1.208038 1.452237
[3,] -1.296135 1.364140