################################################################################
## Einführung in R
## Modul 4: Klassische Tests
################################################################################

## copyfree, Oktober 2007
################################################################################
## Dr. Fränzi Korner-Nievergelt
## oikostat - Statistische Analysen und Beratung
## Ausserdorf 43
## CH - 6218 Ettiswil
## Tel.: +41 (0) 41 980 49 22
## fraenzi.korner@oikostat.ch
## www.oikostat.ch
################################################################################


## Inhalt
################################################################################
## 1. Daten bereitstellen
## 2. T-Test (Vergleich zweier Mittelwerte)
## 3. F-Test (Vergleich zweier Varianzen)
## 4. U-Test (Rangsummentest von Wilcoxon)
## 5. Binomialtest (Vergleich einer Proportion mit einem theoretischen Wert)
## 6. Proportionentest (Vergleich zweier Proportionen)
## 7. Chi2-Test
## 8. Fisher's Exact Test
## 9. Korrelation

## 1. Daten bereitstellen
################################################################################
## In den folgenden Codes wird der Beispielsdatensatz erzeugt. Erklärungen
## zu den Codes und zum Arbeiten mit der R-Console finden Sie in Modul 1.
nm<-12 # Anzahl Männchen, z.B. Grünfink
nw<-14 # Anzahl Weibchen
set.seed(2353)     # Zufallszahlengenerator einstellen, damit später wieder die gleichen Zufallszahlen generiert werden
fluegelm<-round(rnorm(nm, mean=77.8, sd=1.35), digits=1)
fluegelw<-round(rnorm(nw, mean=74.6, sd=1.74), digits=1) # Angabe von Mittelwert und Standardabweichung aus Cramp (1994)
dat<-data.frame(sex=c(rep("M", nm), rep("W", nw)), fluegel=c(fluegelm, fluegelw))  # Vektoren zu einem Datenfile zusammenstellen
dat$sex01<-as.numeric(dat$sex)-1
dat$Alter<-c(rep(3, 8), rep(4, 4), rep(3, 8), rep(4, 6)) # Variable Alter einfügen
dat$Alter01<-as.numeric(dat$Alter)-3
set.seed(16523)
dat$Gewicht<- -3-0.5*dat$sex01+0.2*dat$Alter01+0.3*dat$fluegel+round(rnorm(length(dat$sex), 0, 0.9), 1)

## 2. T-Test
################################################################################

# Mit dem T-Test kann ein Mittelwert mit einem theoretischen Wert oder die Mittelwerte zweier Stichproben verglichen werden.

t.test(fluegelm, fluegelw) # ist der Mittelwert in der Flügellänge zwischen M und W verschieden?
# per default verwendet R die Welch-Korrektur für unterschiedliche Varianzen in den beiden Stichproben.
# Wenn man die Welch-Korrektur nicht möchte, dann muss man sie aktiv ausschalten,
# indem man das Argument var.equal = TRUE setzt.
# Dies sollte nur gemacht werden, wenn die Varianzen in beiden Stichproben gleich sind.

t.test(fluegelm, fluegelw, var.equal=TRUE)  # klassischer T-Test


# Alternative Variante für die Spezifikation des T-Testes:
t.test(fluegel~sex, dat)  # "Formelnotation", die auch für die linearen Modelle genutzt wird.

## 3. F-Test (Vergleich von Varianzen)
################################################################################
# Sind die Varianzen in den zwei Gruppen verschieden?
var.test(fluegelm, fluegelw)
var.test(fluegel~sex, dat)  # dasselbe mit Formelnotation  (unterschiedliches Ergebnis wegen fehlendem Wert im dat)
# also dürfen wir im t-Test das Argument var.equal = TRUE setzen. Damit haben wir etwas mehr Macht.



## 4. U-Test, Wilcoxon Rangsummentest
#################################
# Vergleich zweier nicht normalverteilter Stichproben
wilcox.test(fluegel~sex, dat)

# die Warnung weist uns darauf hin, dass verbundene Werte (mind 2 x dieselben Werte) vorkommen
# sie muss meist nur ernst genommen werden, wenn es viele Bindungen hat.
# Die Stetigkeitskorrektur wird durchgeführt, da eine diskrete Teststatistik (U oder W) mit einer theoretischen kontinuierlichen Verteilung verglichen wird.

# Wenn die Daten jeder Stichprobe normalverteilt sind, dann hat der t-Test mehr Macht.
# Wenn die Verteilung der Daten stark von der Normalverteilung abweicht, dann hat der U-Test mehr Macht.


## 5. Binomialtest
################################################################################
# Vergleich eines Verhältnisses mit einem theoretischen Verhältnis
# Bsp. (erfunden) Ein Spatz hat Hanf- und Sonnenblumenkörner ad libitum zur Auswahl.
# Er frisst 20 Körner, wovon 5 Hanf und 15 Sonnenblumenkörner sind.
# Bevorzugt der Spatz eine der beiden Körnersorten?
# Weicht das Verhältnis 5/20 (oder 15/20) signifikant von 0.5 ab?

binom.test(5, 20)

## 6. Proportionentest
################################################################################
# Vergleich zweier Proportionen?
# Ein Fink frisst 30 Körner, wovon 20 Hanf- und 10 Sonnenblumenkörner sind.
# Unterscheiden sich Spatz und Fink in der Körnerwahl?
prop.test(c(5, 20), c(20, 30))




## 7. Chi2-Test
################################################################################
# Sind zwei kategorische Variabeln miteinander korreliert?
# Bsp. (aus Crawley 2007) Ist die Haar- und Augenfarbe miteinander korreliert?

anzahl<-matrix(c(38, 14, 11, 51), nrow=2, dimnames=list(c("blond", "dunkel"), c("blau", "braun")))
anzahl

chisq.test(anzahl)
# wenn die Yates-Stetigkeitskorrektur ausgeschalten werden möchte, dann setzt man das Argument correct = FALSE.
chisq.test(anzahl, correct=FALSE)

chitest<-chisq.test(anzahl, correct=FALSE)
str(chitest)
chitest$expected  # erwartete Anzahl Beobachtungen, falls kein Zusammenhang zwischen Augen- und Haarfarbe bestehen würde.

# Die erwarteten Zahlen sollten in allen Zellen > 5 sein. Wenn das nicht so ist, sollte der Fisher's exact Test vorgezogen werden.



## 8. Fisher's exact Test
################################################################################
# Exakter Chi-Test für kleine Stichproben.
# Bsp.(erfunden) Enthalten Steinadler-Horst auf Bäumen eher grüne Äste als Horste auf Felsen?

anzahl<-matrix(c(6, 4, 2, 8), nrow=2, dimnames=list(c("grün", "kein grün"), c("Baum", "Fels")))

fisher.test(anzahl)




## 9. Korrelation
################################################################################
plot(Gewicht~fluegel, dat)

cor.test(dat$fluegel, dat$Gewicht) # Pearsons Produkte-Moment-Korrelation (normalverteilte Daten)
cor.test(dat$fluegel, dat$Gewicht, method="spearman") # Spearman Rangkorrelation (nicht normalverteilte Daten)
# Die Warung zeigt an, dass Bindungen vorhanden sind. Wenn viele Bindungen vorhanden sind, dann sollte der p-Wert mit Vorsicht genossen werden.




## Dank
################################################################################

## Dieses Dokument entstand für einen Einführungskurs in R an der DO-G-Tagung
## vom 29. September - 3. Oktober 2007 in Giessen
## Folgenden Personen danke ich herzlich für die Unterstützung und die Werbung:
## Dr. Ommo Hüppop, Dr. Heiko Schmaljohann, Dr. Wolfgang Fiedler,
## Dr. Thomas Gottschalk, Holger Schielzeth, Dr. Roger Mundry und Dr. Boris Schröder