# 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)
names(data)
data$Group <- as.factor(data$Group)
nreps <- 5000
N <- length(data$Score)
n.i <- as.vector(table(data$Group))   # Create vector of sample sizes
k <- length(n.i)

model <- anova(lm(data$Score ~ data$Group))
obt.F <- model$"F value"[1]     # Our obtained F  statistic
obt.p <- model$"Pr(>F)"[1]
cat("The obtained value of F from the standard F test is ",obt.F, "\n")
cat("This has an associated probability of ", obt.p,  "\n")
samp.F <- numeric(nreps)
counter <- 0
set.seed(1086)
for (i in 1:nreps) {
  newScore <- sample(data$Score)
  newModel <- anova(lm(newScore~data$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")
pvalue <- counter/nreps
cat("\nThe calculated value of p from randomized samples is ",pvalue, "\n")
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)

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")


# Demonstration of lmPerm library
library(lmPerm)
result <- summary(aovp(data$Score ~ as.factor(data$Group), perm = "Exact", seqs = FALSE))
cat("\n \n")
cat(" Using the lmPerm library we obtain the following results. \n")
print(result)