###
#
# Einfuehrung in die statistische Standardsoftware  R
#
# (C) 2005 by Gerald van den Boogaart, Greifswald
# Version 0.1
# 
# R kann von 
# http://www.cran.r-project.org
# kostenlos heruntergeladen werden und laeuft auf windows, linux und MacOS. 
#
# Teil I: R kennenlernen
# Tip: Laden Sie diese Datei mit "File/Display File"
#      Kopieren sie den Inhalt Zeile fuer Zeile in die R-Console
#      Dazu gibt es einen speziellen Knopf in der Quickbar
### Ab dem # ist es nur noch Kommentar
### Block 1: R als Taschenrechner
2+3        # Na was wohl?
1+2*3      # Punkt vor Strich
(1+2)*3    # Klammern mit runden Klammern
exp(1)     # Funktionen (z.B. die Exponentialfunktion)
? exp      # Hilfe anzeigen (am Anfang etwas kompliziert, achten Sie auf
           #  die Beispiele am Ende)
? sin      #
pi
2.5*2      # Dezimalpunkt statt Dezimalkomma
x <- 7     # x den Wert 7 zuweisen
x          # Variablen Wert verwenden
y <- 2*x   # Rechnen mit Variablen
ls()       # ls=LiSt: Zeige Liste der definierten Variablen
y          # raten sie mal anzeigen
# Aufgaben 
# 1.1 Berechnen Sie mit R den Kosinus von pi
# 1.2 Berechnen Sie mit R die Flaeche der in der folgenden Graphik
# dargestellten Wiese:
library("MASS")      # Lade neue Befehle (hier eqscplot=Equal Scale Plot)
eqscplot(c(0,110,110,0,0),c(10,10,60,40,10),type="l")
# 1.3 Auf einer kreisfoermigen Flaeche von 50cm Radius haben sie 23 Kaefer
#     gezaehlt. Wieviele Kaefer befinden sich pro Quadratmeter auf der Wiese? 

### Block 2: Rechnen mit Zahlenlisten
# In der Statistik muss man oft mit grossen Mengen on Zahlen gleichzeitig
# rechnen. 
c(1,2,3)   # Eine Zahlenliste (c = concatnate = Zusammenhaengen)
x <- c(1,2,3) # Zahlenliste in einer Variable speichern
x          # 
2*x        # mit der ganzen Liste gleichzeitig rechnen
2+x        
x*x        # mit zwei Listen rechnen
x <- 1:10  # 1 bis 10
x
x^2        
sin(x)     #
x <- seq(-3,3,by=0.1) # -3 bis 3 in 0.1 Schritten
x
y<-x^2
y
# Aufgaben:
# 2.1 Sie haben auf mehreren Wiesen Kaeferfallen aufgestellt. Manche Fallen
#     wurden von Kuehen zertreten. Es haben jeweils:
FallenZahl <- c(5,5,3,5,2,3,1,5)
# Fallen auf den verschiedenen Wiesen ausgewertet werden koennen und
# insgesamt wurden auf den Wiesen jeweils
FangSumme  <- c(33,30,18,26,8,11,4,26)
# Kaefer gefangen.
# Berechnen Sie die Anzahl der Faenge pro Falle auf jeder Wiese.

### Block 3: R als graphischer Taschenrechner
plot(x,y)           # Die Punkte in ein Koordinatensystem malen
plot(x,y,type="b")  # und verbinden (b=both)
plot(x,y,type="l")  # und die Punkte weg (l=lines), "l" ist "el" nicht "eins"
lines(x,x,col="red")# dazumalen kann man mit lines, points, arrows, segments, abline
lines(x,sin(x),col="blue")
f <- function(x) -x^2+2*x+2
plot(x,f(x))
abline(h=0)         # Horizontale bei 0
abline(v=0)         # Vertikale bei 0
abline(a=1,b=-1)    # y=a+x*b
# Aufgaben:
# 3.1 Ein fluoreszierender Marker wird an die DNA einer Hirnzelle eingebaut.
#     Mit jeder Generation halbiert sich der Anteil der markierten Nukleotide,
#     und folgt so einem exponentiellen Abfall.
#     in der ersten Generation sind circa 2% der Nukleotide markiert. 
#     Zeichnen Sie Graphik, welche den Anteil markierter Nukleotide in
#     Abhaengigkeit von der Generationenzahl darstellt.
#     Zeichnen Sie eine horzontale Linie bei einem Anteil von 0.001,
#     unterhalb dem die Markierung nicht mehr nachgewiesen werden kann.
#     Nach wievielen Generationen wird die Markierung unsichtbar? 


### Block 4: Datenmatrix
data(iris)          # Beispieldaten laden
iris                # Daten anzeigen
? iris              # Informationen zum Datensatz
# Die Standarddarstellung von statistischen Daten sind sogenannte
# Datenmatrizen:
# Die Spalten enthalten verschiedene Informationen, (hier Groessen)
# Zeilen gehoeren zu je einem statistischen Individum, (hier Blueten) 
# Petal heisst Bluetenblatt
# Sepal heisst Kelchblatt 
fix(iris)           # Daten ansehen
dim(iris)           # Zeilen- und Spaltenzahl 
names(iris)         # Spaltennamen
iris$Sepal.Length   # Einzelne Daten ausgeben
iris$Sepal.Width
plot(iris$Sepal.Length,iris$Sepal.Width)

# Im Datensatz werden drei Arten verglichen:
iris$Species        # Kategorielle Daten
levels(iris$Species)# Moegliche Werte
# Ein Streudiagramm:
plot(iris$Sepal.Length,iris$Sepal.Width)
plot(iris$Sepal.Length,iris$Sepal.Width,pch=20)
# Versuchen Sie hier auch mal
# pch=10, pch="A", pch=".", cex=3, main="Super Graphik", col="yellow"
# col=2, type="l" (el!!, ohnehin kurios), xlim=c(0,10), ylim=c(0,10)
# Alle Tricks finden Sie unter ? par und ? plot

plot(iris$Sepal.Length,iris$Sepal.Width,col=c("red","green","blue")[iris$Species])
# Eine Streudiagrammmatrix (plot jede gegen jede Variable)
plot(iris)          
plot(iris,col=c("red","green","blue")[iris$Species])
#
# Oft moechte man nur einen Teildatensatz ansehen:
# Splite Datensatz nach iris$Species auf und nimm die Art "setosa":
Setosa <- split(iris,iris$Species)[["setosa"]]  # nur Setosa
Setosa
# Aufgaben
# 4.1: Wiederholen Sie Block 4 fuer den Teildatensatz Setosa
# 4.2: Speichern Sie eine der Graphiken und binden Sie in ein Word-Dokument ein.# Tip: Menue bei aktivem Graphikfenster oder ? dev.copy

# Block 5: Deskriptive Statistik
#
attach(Setosa)    # Variablen in Setosa direkt zugaenglich machen
Species
hist(Petal.Length)          # Ein Histogramm
mean(Petal.Length)          # Mittelwert
abline(v=mean(Petal.Length),col="red")
median(Petal.Length)          # Mittlerer Wert
abline(v=median(Petal.Length),col="blue")
sd(Petal.Length)            # Die Standardabweichung
# ist die Wurzel der mittleren Quadratischen Abweichung vom Mittelwert 
# Die Standardabweichung misst die Streuung
# Typischerweise liegen 60\% der Daten nicht weiter als sd vom
# Mittelwert entfernt:
abline(v=mean(Petal.Length)+c(-1,1)*sd(Petal.Length),col="green")
detach(Setosa)    # attach wieder aufheben
attach(iris)
plot(Petal.Length,Petal.Width,col=(2:4)[Species])
abline(v=mean(Setosa$Petal.Length))
abline(h=mean(Setosa$Petal.Width))

# Aufgaben
# 5.1 Geben Sie die mittleren Kelch und Bluetenblatt Groessen fuer virginica an.
# 5.2 Zeichnen Sie mit points und mit pch=20 in die Graphik  
plot(Petal.Length,Petal.Width,col=(2:4)[Species])
# die Mittelpunkte aller drei Gruppen ein. Dazu brauchen Sie drei points Befehle

# Teil II
# Block 6: Vergleich von Gruppen
#
data(iris)   # Auf ein Neues
plot(iris,col=c("red","green","blue")[iris$Species])
# Offenbar unterscheidet sich das Petal sehr deutlich
# Statistische Tests sind immer folgendermassen aufgebaut:
#
# Name: Irgend ein Name 
# Voraussetzung: z.B. Unabhaengigkeit der Beobachtungen
# Hypothese: Irgend eine Vorstellung
#  gegen
# Alternative: Das Gegenteil dieser Vorstellung
# Befehl: irgend ein R-Befehl um den Test durchzufuehren
# Der Computer berechnet jeweils einen p-Wert. Liegt der p-Wert unter 0.05,
# und sind die Voraussetzungen erfuellt, so gilt die Hypothese als
# widerlegt und die Alternative als vermutlich wahr. 0.05 ist dann eine
# Grenze fuer die Irrtumswahrscheinlichkeit. 

# Unser erster Test:
# Name: Kruskal Wallis Rang Summen Test
# Voraussetzungen: Unabhaengige Stichproben, stetige Verteilung 
# Hypothese: Gruppen haben im Mittel gleich grosse Werte
#  gegen
# Alternative: Verteilungen der Gruppen sind gegegeneinander versetzt
kruskal.test(Sepal.Width~Species,data=iris)

# Aehnliche Dinge testen die folgenden Tests:
# Name: Varianzanalyse
# Voraussetzungen: Unabhaengig, Normalverteilt, gleiche Streuung
# Hypothese: Gruppen haben gleichen Mittelwert
#  gegen
# Alternative: Gruppen haben unterschiedlichen Mittelwert
# Besonderheit: der p-Wert heisst Pr(>F)
anova(aov(Sepal.Width~Species,data=iris))

# Name: Zwei Stichproben t-test
# Voraussetzungen: Zwei Gruppen, unabhaengig, normalverteilt, gleiche Streuung
# Hypothese: Gruppen haben gleichen Mittelwert
#  gegen
# Alternative: Gruppen haben unterschiedlichen Mittelwert
Setosa <- split(iris,iris$Species)[["setosa"]]
Virginica <- split(iris,iris$Species)[["virginica"]]
t.test(Setosa$Sepal.Width,Virginica$Sepal.Width)

# Name: Wilcoxen Rang Summen Test
# Voraussetzungen: Zwei Gruppen, unabhaengig, stetig verteilt
# Hypothese: Gruppen haben gleichen Median
#  gegen
# Alternative: Gruppen haben unterschiedlichen Median
wilcox.test(Setosa$Sepal.Width,Virginica$Sepal.Width)

# Test fuer die Streuung:
# Name: Fligner Killean Median Test auf gleiche Varianz
# Voraussetzungen: Unabhaengige Stichproben, stetige Verteilung
# Hypothese: Gruppen haben gleiche Streuung
#  gegen
# Alternative: Gruppen haben verschiedene Streuung
fligner.test(Sepal.Width~Species,data=iris)

# Test auf Normalverteilung
# Name: Shapiro-Wilk-Test
# Voraussetzungen:  Unabhaengige Stichprobe
# Hypothese: Daten sind normalverteilt
#  gegen
# Alternative: Daten sind nicht normalverteilt
shapiro.test(iris$Petal.Length)   # Mischung aus drei Normalverteilungen
shapiro.test(Setosa$Petal.Length)
shapiro.test(Virginica$Petal.Length)

# Aufgaben:
# 6.1: Koennte die Kelchblatt-Laenge bei den einzelne Irisarten normalverteilt sein?
# 6.2: Haben die Arten virginica und versicolor unterschiedliche mittlere
#      Kelchblattlaengen

# Block 7: Regression
# Oft moechte man "je ..., desto ..." Beziehungen nachweisen
plot(iris$Petal.Length,iris$Petal.Width)
# lm = lineares Modell
MeinGeradenModell <- lm(Petal.Width~Petal.Length,data=iris)
abline(MeinGeradenModell)
# Modell: Daten = Gerade + Fehler
arrows(Petal.Length,predict(MeinGeradenModell),Petal.Length,Petal.Width)
# Test wie bei Varianzanalyse:
# Name: Test auf Regression
# Voraussetzungen: unabhaengige normalverteilte Fehler gleicher Streuung
# Hypothese: y haengt nicht von x ab gemaess y=a+Fehler
#   gegeben
# Alternative: je mehr x, desto mehr y gemaess y=a+b*x+Fehler
# Tip: p-Wert heisst hier Pr(>P)
anova(MeinGeradenModell)

# Alternativer Ansatz: Correlationstests:
# Name: Test auf Lineare (Pearsonsche) Korrelation
# Voraussetzungen: unabhaengige Zeilen, normalverteilt
# Hypothese: x,y haengen nicht linear voneinader ab
#   gegeben
# Alternative: x,y haengen linear voneinander ab
cor.test(Sepal.Length,Sepal.Width,method="pearson")
  
# Name: Test auf Spearman Rang Korrelation
# Voraussetzungen: unabhaengigkeit der Zeilen, stetige Verteilung
# Hypothese: x,y sind unabhaengig
#   gegeben
# Alternative: je mehr x, desto mehr/oder weniger y
cor.test(Sepal.Length,Sepal.Width,method="spearman")

# Aufgabe:
# 7.1: Untersuchen Sie den Zusammenhang von Kelchblatt und Bluetenblattlaenge
#      fuer die Art setosa.


# Block 8: Eigene Daten einlesen
# Eigene Daten einzulesen ist in R machmal etwas schwierig
# Wir ueben das an einem besonders haarigen Beispiel:
# 1) Laden Sie die Datei http://... aus dem Internet herunter
# 2) Oeffnen Sie die Datei mit Notepad und loeschen Sie
#    alles ausser den Daten und den Variablenamen
# 3) Laden Sie mit den folgenden Befehlen die Daten ein.
#    Modifizieren Sie den Dateinamen und Pfadnamen entsprechend.
#    Verwenden Sie / statt \
#    sep="\t" heisst Tab as Trennzeichen verwenden
#    header=TRUE heisst Variablennamen mit einlesen
daten <- read.table("c:/meinverzeichnis/DatenDatei.html",header=TRUE,sep="\t")
# Tritt hier ein fehler auf?
# -> Ist der Dateiname richtig?
# -> Sind die Daten gut als Tabelle strukturier
# -> Ist ueberfluessiger Text in der Datei oder ueberfluessige Zeilen
daten
# Test: Sieht das ordentlich aus?
# Dann pruefen:
names(daten)
# Stimmen die Namen? Falls nicht:
# -> Formatfehler in Datei suchen
sapply(daten,class) 
# Stimmen die angegebene Datentypen? Falls nicht:
# -> Punkte statt , verwenden
# -> Fehlende Werte auf NA setzen
# -> Formatfehler in Datei suchen

# Aufgabe
# 8.1: Stellen Sie die Daten graphisch dar
# 8.2: Beantworten Sie die in der Beschreibung der Daten aufgeworfene
#      statistische Frage 


# Block 9: Kategorielle Daten
table(iris$Species)
barplot(table(iris$Species))

# Name: Chi-Quadrat-Test auf diskrete Gleichverteilung
# Voraussetzungen: Genuegend unabhaengige Daten
# Hypothese: Gleiche Wahrscheinlichkeit fuer Gruppenzugehoerigkeit
#   gegen
# Alternative: Ungleiche Wahrscheinlichkeit fuer Gruppenzugehoerigkeit
chisq.test(table(iris$Species))

data(Titanic)
? Titanic
mosaicplot(Titanic, main = "Survival on the Titanic")
## Higher survival rates in children?
apply(Titanic, c(3, 4), sum)
## Higher survival rates in females?
apply(Titanic, c(2, 4), sum)
## Use loglm() in package 'MASS' for further analysis ...

# Name: Fishers exact Test
# Voraussetzungen: Kontingenztafel
# Hypothese: Einteilung in Spalten und Zeilenmerkmal sind unabhaengig
#   gegen
# Alternative: Einteilungen in Spalten und Zeilenmerkmal sind abhaengig
fisher.test(apply(Titanic, c(2, 4), sum))

# Name: Chi-Quadrat-Test fuer Kontingenztafeln
# Voraussetzungen: Kontingenztafel
# Hypothese: Einteilung in Spalten und Zeilenmerkmal sind unabhaengig
#   gegen
# Alternative: Einteilungen in Spalten und Zeilenmerkmal sind abhaengig
chisq.test(apply(Titanic, c(2, 4), sum))

# Aufgabe:
# 9.1 Galt "Frauen und Kinder zuerst" auf der Titanic?
# 9.2 Halten Sie das mit den 50/50/50 fuer Zufall?

# Block 10: Ausblick
# R kann mehr: Hier ist ein kurzer Ausblick
boxplot(Sepal.Length~Species,data=iris,notch=T) 
help.search("regression") # Hilfeseite zu Thema suchen
help.start()   # Hilfe, weitere Doku und Links ins Internet
# Falls es nicht geht can man den Browser auch angeben:
# help.start(browser="C:/Programme/Netscape") 
demo()         # Verfuegbare Demonstrationen
demo(graphics) # Demonstration zur Graphiken
demo(persp)