## pdf("E:\\blue_drive\\Rmisc\\graphs\\meniscus.pdf") meniscus <- read.table("http://www.stat.ufl.edu/~winner/data/meniscus.dat",header=F, col.names=c("method","loadfail","displace","stiff")) attach(meniscus) set.seed(13579) r <- 3 n_T <- 18 (ybarall <- mean(loadfail)) (ybar1 <- mean(loadfail[1:6])); (var1 <- var(loadfail[1:6])) (ybar2 <- mean(loadfail[7:12])); (var2 <- var(loadfail[7:12])) (ybar3 <- mean(loadfail[13:18])); (var3 <- var(loadfail[13:18])) (sstr <- 6*((ybar1-ybarall)^2 + (ybar2-ybarall)^2 + (ybar3-ybarall)^2)) (sse <- (6-1)*(var1 + var2 + var3)) (f_obs <- (sstr/(r-1))/(sse/(n_T-r))) num_samp <- 10^5-1 f_rand <- rep(0,num_samp) for (i1 in 1:num_samp) { randgrp <- sample(n_T,size=n_T,replace=FALSE) ybar1 <- mean(loadfail[randgrp[1:6]]); var1 <- var(loadfail[randgrp[1:6]]) ybar2 <- mean(loadfail[randgrp[7:12]]); var2 <- var(loadfail[randgrp[7:12]]) ybar3 <- mean(loadfail[randgrp[13:18]]); var3 <- var(loadfail[randgrp[13:18]]) sstr <- 6*((ybar1-ybarall)^2 + (ybar2-ybarall)^2 + (ybar3-ybarall)^2) sse <- (6-1)*(var1 + var2 + var3) f_rand[i1] <- (sstr/(r-1))/(sse/(n_T-r)) } summary(f_rand) (p_rand <- (sum(f_rand >= f_obs)+1)/(num_samp+1)) hist(f_rand,xlab="F*",main="Randomization Distribution for F-test") abline(v=f_obs)