library(reshape2)
library(gplots)
library(multcomp)

DF_aj = read.csv('../../Data/ASCII_Comma/Chapter_12/aj.txt', header=FALSE)
colnames(DF_aj) = c('N_Squares')
DF_aj$Species = 'AJ'

DF_c57 = read.csv('../../Data/ASCII_Comma/Chapter_12/c57.txt', header=FALSE)
colnames(DF_c57) = c('N_Squares')
DF_c57$Species = 'C57'

DF_f2 = read.csv('../../Data/ASCII_Comma/Chapter_12/f2.txt', header=FALSE)
colnames(DF_f2) = c('N_Squares')
DF_f2$Species = 'F2'

DF = rbind( DF_aj, DF_c57, DF_f2 )
DF$Species = as.factor(DF$Species)

fit = aov(N_Squares ~ Species, data=DF)
print(summary(fit))

#postscript("../../WriteUp/Graphics/Chapter12/prob_27_plotmeans.eps", onefile=FALSE, horizontal=FALSE)
plotmeans(N_Squares ~ Species, data=DF)
#dev.off()

# A Tukey comparison of means test:
#
TukeyHSD(fit)

par(las=2)
par(mar=c(5, 8, 4, 2))
plot(TukeyHSD(fit))

# Use the multcomp package:
#
par(las=2)
par(mar=c(5, 4, 6, 2))
tuk = glht(fit, linfct=mcp(Species='Tukey'))
plot(cld(tuk, level=0.05), col='lightgrey')

# A nonparametric test for one-way ANOVA:
#
print(kruskal.test(N_Squares ~ Species, data=DF))

#
# Implement the Bonferroni computed confidence intervals
#
sp = sd(DF$N_Squares) # the pooled standard deviation

group_means = aggregate(DF$N_Squares, list(DF$Species), mean)
colnames(group_means) = c('Group', 'Mean')

so = order(group_means$Mean, decreasing=TRUE)
group_means = group_means[so,]
print(group_means)

alpha = 0.05

I = 3 # the number of groups
T = table(DF$Species) # how many samples do we have of each species
print(T)
J = min(T) # a conservative number of samples per species
k = choose(I, 2)

t_crit = qt(1 - (alpha/2)/k, I*(J-1))
ci_width = (t_crit * sp) / sqrt(J/2)
print(ci_width)