tabletten <- read.table("tabletten.dat", header=T)
m <- 10; k <- 7
tab.aov <- aov(Data~Labs, data=tabletten)
cat("\nDie R-interne Varianzanalyse ergibt:\n")

print(summary(tab.aov))

cat("\nDie Hypothese, dass keine systematischen Unterschiede bestehen, kann\n")
cat("also zu praktisch jedem beliebigen Niveau verworfen werden.\n\n")

# Konfidenzintervalle für die Differenzen der Erwartungswerte - wegen der
# Symmetrie (das Konfidenzintervall für (j,l) ist gleich dem an 0
# gespiegelten von (l,j)) genügt es, die unteren Schranken zu betrachten.

# zuerst hole ich mir die empirischen Mittelwerte...
means <- model.tables(tab.aov, type="means")$tables$Labs

#und berechne die Differenzen.
meandiff <- matrix(rep(means,7)-rep(means, rep(7,7)), nrow=7)


dotest <- function(alpha){
	# Für die Konfidenzintervalle benötige ich nun noch c und s
	#alpha <- 0.05
	c <- qtukey(1-alpha, k, m*k-k)
	s <- sqrt(summary(tab.aov)[[1]]$"Mean Sq"[2])
	
	# Nun kann ich eine Matrix mit den unteren Intervallgrenzen erzeugen...
	d <- c*s/sqrt(m)
	mlower <- meandiff-d
	
	# und ausgeben, welche Paare "aus der Reihe fallen"
	test <- ifelse(mlower>0, 1,0)
	
	cat("Die folgenden Paare unterscheiden sich signifikant (Niveau:",alpha)
	cat("):\n")
	for (i in 1:7) {for (j in 1:7) {if (test[i,j]==1) cat("(",i,j,")")}}
	cat("\n")
}
dotest(0.05)
cat("\nOder wem das lieber ist:\n")  
dotest(0.01)

#rm(tabletten, tab.aov)