remove(list=ls())
#setwd("/Volumes/MichaelsHD/Dropbox/Angewandte Statistik SoSe11/Uebungen/Blatt11/R")
setwd("/Volumes/Daten/Dropbox/Angewandte Statistik SoSe11/Uebungen/Blatt11/R")
set.seed(23514)

#Daten einlesen
mensa <- read.table("./mensa2.dat",header=T)
essen <- mensa$essen
anzahl <- mensa$anzahl
#ANOVA fuer Tukey und Scheffe
alle.aov <- aov(anzahl~essen)

#Tukey
TukeyHSD(alle.aov,"essen",ordered=TRUE)

#Scheffe
library(agricolae)
scheffe.test(alle.aov,"essen",group=TRUE)

#Teststatistiken von Hand
m <- sum(essen=="LF")
k <- 4

#Matrix fuer Tukey (wird in for-Schleife gefüllt)
tukmat <- matrix(0,ncol=4,nrow=4)
#Matrix fuer p-Werte bei Tukey (wird in for-Schleife gefüllt)
ptukmat <- tukmat
#Matrix fuer Scheffe (wird in for-Schleife gefüllt)
schefmat <- tukmat
#Vektor mit den Mittelwerten der Gruppen (also X_iPunkt)
xipunkt <- c(sum(anzahl[essen=="LF"]),sum(anzahl[essen=="GG"]),sum(anzahl[essen=="VEG"]),sum(anzahl[essen=="BIO"]))/m
xipunkt
#Sigma-Dach aus dem kompletten linearen Modell
sigmaeps <- summary(lm(anzahl~essen))$sigma
sigmaeps
#vergleichsgroesse fuer scheffe
scheff0 <- sqrt(sigmaeps*sigmaeps*qf(0.95,k-1,k*(m-1))*(k-1)/6)
scheff0
#Tukey alle
tukall <- (max(xipunkt)-min(xipunkt))*sqrt(m)/sigmaeps
tukall
#Quantil fuer Tukey
cAlphaTuk <- qtukey(0.95,k,k*(m-1))
cAlphaTuk
#P-Wert fuer alle
1-ptukey(tukall,k,k*(m-1))
#Matrizen mit einzelnen Tests
for (i in 1:3) {
	for (j in (i+1):4) {
		#Tukey
		tukmat[i,j] <- abs(xipunkt[i]-xipunkt[j])*sqrt(m)/sigmaeps
		#P-Werte fuer Tukey
		ptukmat[i,j] <- 1-ptukey(tukmat[i,j],k,k*(m-1))
		#Scheffe
		schefmat[i,j] <- (abs(xipunkt[i]-xipunkt[j]) > scheff0)
	}
}
#Reihenfolge in den Matrizen: LF, GG, VEG, BIO
tukmat
ptukmat
schefmat

#implementierte Tests und Tests aus der Vorlesung kommen zum selben Ergebnis

#Auf Blatt10 Nr.1b hatte der Vergleich LF-GG einen p-Wert von 0.05232, hier 0.198 (Tukey) und wird bei Scheffe auch nicht abgelehnt
#Der Vergleich BIO-VEG hatte einen p-Wert von 0.2585, hier 0.6465 (Tukey) und wird bei Scheffe auch nicht abgelehnt.