# STATISTIK I
# LOESUNGSVORSCHLAG BLATT 9
################################

##############
# AUFGABE 2 ##
##############

# a)

library(MASS)
schlaf <- sleep

# nicht parametrischer Test
wilcox.test(schlaf[1:10,1],schlaf[11:20,1],paired=TRUE)

# per Hand: Man bildet die Raenge der Betraege der Differenzen
# und addiert diejenigen die positiven sind. Hier ist die Teststatistik somit 0,
# da Medikament 2 immer groessr oder gleich dem Medikament 1 ist.
# Diese Teststatistik ist unter H_0 approximativ normalverteilt mit 
# Erwartungswert 27.5 und Varianz 96.25. Somit ergibt sich der p-Wert approximativ:
2*pnorm(0,27.5,sqrt(96.25))
# 0.005062032
# Somit ist statistisch belegt, dass die Medikamente eine unterschiedliche Wirkung haben.


# parametrischer Test. Wichtig ist hierbei die Annahme einer Normalverteilung der Differenzen
# da die Stichprobe sehr klein ist (n=10)!
t.test(schlaf[1:10,1],schlaf[11:20,1],paired=TRUE)

# per Hand:
# bilde die Differenzen
dif <- schlaf[1:10,1]-schlaf[11:20,1]

# Annahme: Differenzen sind normalverteilt. Da Varianz unbekannt muss diese geschaetzt werden.
sa_dif <- sd(dif)

# die Teststatistik ist hier T= mean(dif)/(sd(diff)/sqrt(10)) und diese ist unter H_0
# t-verteilt mit 9 Freiheitsgraden
sqrt(10)*mean(dif)/sd(dif)
# Wert der Teststatistik: -4.062128

# Somit ergibt sich ein p-Wert von:
2*pt(sqrt(10)*mean(dif)/sd(dif),df=9)
# 0.00283289

# b)
agento <- read.table("http://rosuda.org/lehre/SS12/Datensaetze/AgentOrange.txt",header=T,sep="",quote="")

# erstmal einen Überblick über die Daten verschaffen, z.B. so
table(agento$Dioxin,agento$Veteran)

# die beiden Ausreißer sind beide in der Gruppe der Vietnam-Veteranen, also in der deutlich
# groesseren Gruppe und haben somit vermutlich nicht einen so erheblichen Einfluss, jedoch muesste man
# das genauer ueberpruefen (einmal mit den beiden "Ausreissern" einmal ohne)

# Varianzen gleich?
var(agento$Dioxin[agento$Veteran=="Vietnam"])
# Vietnam ungefähr 7
var(agento$Dioxin[agento$Veteran=="Other"])
# Other ungefähr 5

# Man koennte an dieser Stelle einen Test durchfuehren, ob die Varianzen signifikant verschieden sind, z.B. so
var.test(agento$Dioxin[agento$Veteran=="Vietnam"],agento$Dioxin[agento$Veteran=="Other"])
# nicht signifikant verschieden. Solche Tests wurden allerdings in der Vorlesung nicht besprochen und sind
# somit natuerlich nicht klausurrelevant. Aber es gibt sie.

t.test(agento$Dioxin[agento$Veteran=="Vietnam"],agento$Dioxin[agento$Veteran=="Other"],paired=FALSE)
# p-Wert: 0.7713
# nicht signifikant, keine Unterschiede



##############
# AUFGABE 3 ##
##############

# Chi-Quadrat-Unabhaengigkeits per Hand

# beobachtete Werte:
ds <- matrix(c(2350,121,2337,197),nrow=2)
rownames(ds) <- c("Schulreif","Nicht schulreif")
colnames(ds) <- c("Maedchen","Jungen")


# erwartete Werte bei Unabhaengigkeit
ds <- matrix(c(2314,157,2373,161),nrow=2)
rownames(ds) <- c("Schulreif","Nicht schulreif")
colnames(ds) <- c("Maedchen","Jungen")


# Wert der Teststatistik:
((2314-2350)^2/2314)+((2373-2337)^2/2373)+((157-121)^2/157)+((161-197)^2/161)

# Unter der Annahme von Unabhaengigkeit ist diese Teststatistik chi-quadrat verteilt
# mit einem Freiheitsgrad (Anzahl Zeilen - 1)*(Anzahl Spalten - 1) 

# p-Wert
1-pchisq(17.41068,df=1)

# also hoch signifikant, das heisst die Variablen sind nicht unabhaengig.

# Funktion in R:
chisq.test(ds,correct=FALSE)

# Visualisierung mit einem Mosaicplot:
mosaicplot(ds)

# sieht nicht unabhaengig aus! Dies bestaetigt der Test.

##############
# AUFGABE 4 ##
##############

# beobachtete Werte:
sono <- matrix(c(0.0702,0.0999,0.5814,0.0728,0.0182,0.087,0.0663,0.0041),ncol=2)

sum_col <- colSums(sono)
sum_row <- rowSums(sono)

# erwartete Werte:
erw <- outer(sum_row,sum_col)

# Die Funktion "outer()" berechnet das outer product. 

n <- 4323

# Teststatistik:
n*sum((sono-erw)^2/erw)
# 609.5442

# p-Wert
1-pchisq(609.5442,3)
# = 0

# Dies bedeutet, dass H0 sicher verworfen werden kann, das heisst die
# Variablen sind nicht unabhaengig

# b)
n <- 100

# Teststatistik:
n*sum((sono-erw)^2/erw)
# X^2 = 14.1

# p-Wert
1-pchisq(14.1,3)

# Ja auch dieser p-Wert ist immer noch hoch signifikant

# c)
# Ab welchem n wird dieser Test bei geg. Proportionen signifikant?
# n*sum((sono-erw)^2/erw) >= qchisq(0.95,4)
qchisq(0.95,3)/sum((sono-erw)^2/erw)
# das heisst fuer n >= 56

# => Hat man also 2 Tabellen mit unterschiedlicher Beobachtungszahl n,
# so kann man diese Verfahrensweise einsetzen, um die Staerke der Abweichung
# von der Unabhaenigkeit zu vergleichen.