### setwd("C:/Users/Küchenhoff/Desktop/RLimo")

################################################################################
#
# Beispiel: Ascorbinsäuregehalt von Kohlköpfen
#
# Vorlesung Lineare Modelle (Helmut Küchenhoff) 20.5.2015
# Thema: Diskrete Einflussgrößen
#
# data file: kohl.txt (http://www.statistik.lmu.de/~abender/daten/kohl.txt) 
#
# contents: (1) Einlesen der Daten
#		        (2) Modellbildung
#
# last.modified: HK 19.5.2015
#
################################################################################

#### (1) Einlesen von Daten ####################################################
## aus dem m netz 
Kohl <- read.table('http://www.statistik.lmu.de/~abender/daten/kohl.txt', 
	header = TRUE)
### Kohl <- read.table('kohl.txt', 
###                     header = TRUE)

# Variablen:
# Zeit - Zeitpunkt (3 Werte als Faktorvariable) 
# gen - Genetik (2 Typen)
# gew - Gewicht
# asc - Ascorbinsäuregehalt

# Definition als Faktorvariablen 
Kohl$genf <- factor(Kohl$gen, labels=c('G39','G52'))
Kohl$Zeitf <- factor(Kohl$Zeit, labels=c('T16','T20','T21'))

Kohl
mean(Kohl$asc)
boxplot(asc~Zeitf, ylab="asc", xlab="Zeit", data=Kohl)
tapply(Kohl$asc, Kohl$Zeitf, mean) # means
tapply(Kohl$asc, Kohl$Zeitf, sd) # std. deviations

## Modell mit Referenzkodierung ##
model1 <- lm(asc ~ Zeitf, data=Kohl)
summary(model1)
#Interpreation
#Im Mittel betraegt der Ascorbinsaeuregehalt 56.4 (b0) Einheiten zum 
#Zeitpunkt T16
#Im Mittel betraegt der Ascorbinsaeuregehalt 54.15 (b0+b1) Einheiten zum 
#Zeitpunkt T20
#Im Mittel betraegt der Ascorbinsaeuregehalt 63.6 (b0+b2) Einheiten zum 
#Zeitpunkt T21
anova(model2)
#p-Wert< alpha -> H0: Glichheit der Mittelwerte  kann abgelehnt werden

## Modell mit Mittelwertcodierung ##
model2 <- lm(asc ~ Zeitf-1, data=Kohl)
summary(model2)
anova(model2)
### Vorischt: Jetzt wird test auf Alle Werte sind 0 gerechnet 
### Effektkodierung 
Kohl$D1<- (Kohl$Zeitf=="T16")-(Kohl$Zeitf=="T21")
Kohl$D2<- (Kohl$Zeitf=="T20")-(Kohl$Zeitf=="T21")
model3 <- lm(asc~D1+D2,data=Kohl)
summary(model3)
 anova(model3)
### oder direkt in R 
### Modell mit Effect-Codierung 
contrasts(Kohl$Zeitf) = contr.sum(3)
contrasts(Kohl$genf) = contr.sum(2)
contrasts(Kohl$Zeitf)
cbind(Kohl$D1,Kohl$D2)
model4<- lm(asc~Zeitf,data=Kohl)
summary(model4)
#Interpretaion wie model3

library(car)
library(effects)
plot(allEffects(model4), ask=FALSE)
#schwarze Punkte = Mittelwert in den jeweiligen Gruppen

### zeiteffekt mit referenz T20
contrasts(Kohl$Zeitf) = contr.treatment(3,base=2)
model4<- lm(asc~Zeitf,data=Kohl)
summary(model4)


### zurücksetzen 
contrasts(Kohl$Zeitf) = contrasts(Kohl$Zeitf,contrasts=F)
contrasts(Kohl$genf) = contrasts(Kohl$genf,contrasts=F)


#### (2) Modellbildung #########################################################


#### Modell mit zwei Haupteffeklten 
Kohl_M0 <- lm(asc  ~ genf  + Zeitf , data=Kohl)
summary(Kohl_M0)
#Interpretaion
#Im Mittel betraegt die Ascorbinsaeuregehalt 69.750 in den 
#Referenzkategorien d.h. zum Zeitpunkt T21 und Gen52
#Der Ascorbinsaeuregehalt betraegt in Gruppe Gen39 im Mittel 12.900 
#Einheiten weniger als in der Gruppe Gen52, bei konstanthalten des 
#Zeitpunktes (ceteris paribus).
#Der Ascorbinsaeuregehalt betraegt im Mittel zum Zeitpunkt T16 6.9 
#Eineiten weniger als zum Zeitpunkt T21, bei konstanthalten des Gens
#Der Ascorbinsaeuregehalt betraegt im Mittel zum Zeitpunkt T20 9.15 
#Eineiten weniger als zum Zeitpunkt T21. c.p.


plot(allEffects(Kohl_M0), ask=FALSE)
anova(Kohl_M01)
#sowohl genf als auch Zeitf haben signifikanten Einfluss, siehe p-Wert
#####Orthogonales Design fuehrt zu gleichen Schaetzern wie das einfache Modell ###
Kohl_M01 <- lm(asc ~ Zeitf+genf  + genf *Zeitf , data=Kohl)
summary(Kohl_M01)
#Im Mittel betraegt der asc 71.8 (b0) Einheiten zum Zeitpunkt T21 und Gen52
#Im Vergleich zu Gen52 betraegt die asc in Gen39 17 Einheiten weniger,
#zum Zeitpunkt T21.

#In G52:
#Im Vergleich zu T21 betraegt die asc in T16 (T20) im Mittel 9.3=b1 
#(12.9=b2)Einheiten weniger (in G52).

#In G39:
#Im Vergleich zu T21 betraegt die asc in T16 (T20) im Mittel 4.5(=b1+b4) 
#(5.4=(b2+b5)) Einheiten weniger (in G39, addition Interaktionsterm).

Anova(Kohl_M01)
#Interkation nicht signifikant.
#Effekt von Gen haengt nicht signifiaknt von der Zeitvariable ab
#bzw Effekt von Zeitvariable haengt nicht von der Auspraegung von Gen ab.

plot(allEffects(Kohl_M01), ask=FALSE)

###Effektkodierung liefert interprtierbare betas  
contrasts(Kohl$Zeitf) = contr.sum(3)
contrasts(Kohl$genf) = contr.sum(2)
Kohl_M02 <- lm(asc ~ Zeitf+genf  + genf *Zeitf , data=Kohl)
summary(Kohl_M02)

### (a) Modell mit Haupteffekten Genetik und Gewicht ###########################

with(Kohl,plot(gew,asc))

Kohl_M1 <- lm(asc ~ genf + gew, data=Kohl)
summary(Kohl_M1)
Kohl$gewb<-Kohl$gew-mean(Kohl$gew)

# Test mit Partiellen Quadratsummen  (Sum Sq)
Anova(Kohl_M1, type="II")

# Darstellung der Effekte 
##trellis.device(theme="col.whitebg")
plot(allEffects(Kohl_M1), ask=FALSE)

# Zusammenhang von gew und asc als 2 Parallle Geraden 
plot(Kohl$gew,Kohl_M1$fitted.values)


### (b) Modell mit Haupteffekten  Gentik, Zeit und Gewicht #####################

Kohl_M2 <- lm(asc ~ genf+Zeitf+gew, data=Kohl)
summary(Kohl_M2)

# Testen der Effekte mit partiellen Quadartsumme 
Anova(Kohl_M2, type="II")
# alle Einfluesse signifikant

# Darstellung der Effekte 
plot(allEffects(Kohl_M2), ask=FALSE)

#  Zusammenhang von gew und asc als 6 Parallle Geraden 

plot(Kohl$gew,Kohl_M2$fitted.values,col=Kohl$Zeit)
# -> Verschiedene Farben für die unterschiedlichen Zeiten 

### (c) Modell mit Genetik, Gewicht und Interaktion ############################
Kohl_M3 <- lm(asc ~ genf+gew+gew*genf, data=Kohl)
summary(Kohl_M3)
#Interpretation der Interaktionen 
#(unterschiedlische Steigung in den verschiedenen Gen-Gruppen)
#Steigt gew um eine Einheit in der Gruppe gen52 so steigt asc im Mittel
#um 1279.48 (b2)
#Steigt gew um eine Einheit in der Gruppe gen39 so steigt asc im Mittel
#um 1324.31 (b2+b3)

# Testen der Effekte mit partiellen Quadartsumme 
Anova(Kohl_M3, type="II")
# -> Beachte: Interaktion wird bei Betrachtung der Haupteffekte nicht in 
# das Modell einbezogen

Anova(Kohl_M1, type="II")
# -> QS der Haupteffekte mit Modell M1 identisch
# Nur Partielle QS nicht sinnvoll bei Vorhandensein von Interaktionen 
Anova(Kohl_M3, type="III")

# Darstellung der Effekte 
plot(allEffects(Kohl_M3), ask=FALSE)

# Zusammenhang von gew und asc als 2 Geraden für die beiden Gruppen der 
# Genetik mit unterschiedlicher Steigung  
plot(Kohl$gew,Kohl_M3$fitted.values,col=Kohl$gen)
# -> Verschiedene Farben für die unterschiedliche Genetik 


### (d) Jetzt zusätzlich Zeit ins Modell #######################################

Kohl_M4 <- lm(asc ~ genf + Zeitf + gew + genf*gew, data=Kohl)
summary(Kohl_M4)

# Testen der Effekte mit partiellen Quadartsumme 
Anova(Kohl_M4, type="II")
# -> Interaktion wird nicht bei Betrachtung der Haupteffekte nicht in 
# das Modell einbezogen (Identitsch mit M1)
# -> Beachte Interaktion gew*gen nicht signifikant

# Darstellung der Effekte 
plot(allEffects(Kohl_M4), ask=FALSE)

# Zusammenhang von gew und asc als 6 Gerden 
plot(Kohl$gew,Kohl_M4$fitted.values)
# -> Verschiedene Farben für die unterschiedliche Genetik 
# -> paralell-Verschiebung nach Zeit (einfacher Effekt)   
# -> Unterschiedliche Steigungen nach Genetik (Interaktion gen*gew) kaum 
# erkennbar (Erinnerung Interaktion nicht signifikant)

### (e) Modell mit Haupteffekten und Interaktion gew und Zeit ##################

Kohl_M5 <- lm(asc ~ genf+Zeitf+gew+Zeitf*gew, data=Kohl)
summary(Kohl_M5)

# Testen der Effekte mit partiellen Quadartsumme 
Anova(Kohl_M5, type="II")
# -> Interaktion Zeit*gew nicht signifikant! 

# Darstellung der Effekte 
plot(allEffects(Kohl_M5), ask=FALSE)

#  Zusammenhang von gew und asc als 6 Gerden 
plot(Kohl$gew,Kohl_M5$fitted.values,col=Kohl$Zeit)
# -> Verschiedene Farben für die unterschiedliche Zeit  
# -> Unterschiedliche Steigungen nach Zeit (Interaktion Zeit*gew) 
# -> paralell-Verschiebung nach Zeit (einfacher Effekt)   
# -> Unterschied in der Steigung besonders in Zeitgruppe 16 (Schwarz zu den 
# anderen beiden Gruppen   

#### (f) Modell mit allen Interaktionen ########################################

Kohl_M6 <- lm(asc ~ Zeitf*gew*genf, data=Kohl)
# -> Beachte hier verkürzte Schreibwiese: Wenn die 3er Interaktion im Modell ist,
# werden automatisch alle Haupteffekte und die in der 3er Interaktion enthaltenen 
# Effekte in das Modell einbezogen 
# Kohl_M6 <- lm(asc ~ gen+Zeit+gew+Zeit*gew+ Zeit*gen+gew*gen + Zeit*gew*gen , data=Kohl)  
# liefert identisches Ergebnis 

summary(Kohl_M6)

# Testen der Effekte mit partiellen Quadartsumme 
Anova(Kohl_M6, type="II")
# -> Alle Interaktion nicht signifikant! 
# -> Einfaches Modell wohl ausreichend 

# Darstellung der Effekte
plot(allEffects(Kohl_M6), ask=FALSE)

#  Zusammenhang von gew und asc als 6 Geraden mit unterschiedlichen Steigungen  
plot(Kohl$gew,Kohl_M6$fitted.values,col=Kohl$Zeit)
# -> Verschiedene Farben für die unterschiedliche Genetik