# Randomization test for one way ANOVA
# Data file has column1 = group and column2 = dv
# This data file has unequal n's

data <- read.table("http://www.uvm.edu/~dhowell/StatPages/ResamplingWithR/RandomOneway/Foa.dat", header = TRUE)
attach(data)    #Using data from Foa.dat--4 groups Group and Score
names(data)
Group <- as.factor(Group)
nreps <- 5000
N <- length(Score)
n.i <- as.vector(table(Group))   # Create vector of sample sizes
k <- length(n.i)

model <- anova(lm(Score ~ Group))
obt.F <- model$"F value"[1]     # Our obtained SSbet
samp.F <- numeric(nreps)
counter <- 0
for (i in 1:nreps) {
  newScore <- sample(Score)
  newModel <- anova(lm(newScore~Group))
  samp.F[i] <- newModel$"F value"[1]
  if (samp.F[i] > obt.F) counter = counter + 1
  }
#time2 <- proc.time()
#cat(" The timing statistics are " ,(time2 - time1),"\n")
cat("The obtained value of F is ",obt.F, "\n")
pvalue <- counter/nreps

par(mfrow = c(2,1))
hist(samp.F, breaks = 50, main = "Histogram of F on Randomized Samples",
   xlab = "F value", probability = TRUE, col = "green", border = 1,
   , xlim = c(0,7), ylim = c(0,1))
 legend("topright", paste("obtained.F = ", round(obt.F, digits = 4)), col=1,  cex = 0.8)
 legend("right",paste("p-value = ",round(pvalue, digits = 4)))
 arrows( 5.5, 0.8,obt.F,0, length = .125)
 
cat("The obtained value of F of obtained means is ",obt.F, "\n")
f <- seq(0, 7,.01)
dens <- df(f,3,41)
par(new = T)
plot(f,dens, col = "red", type = "l", xlim = c(0,7), ylim = c(0,1), xlab = "", ylab = "")
#polygon(f,dens, col = "red")