## ----teengamb_seq_part---------------------------------------------------
library(faraway)
mod <- lm(gamble ~ income + status + verbal, data = teengamb)


## ----teengamb_seq_part-a-------------------------------------------------
  ## setze formel für jedes Modell zusammen
  lhs.gamble <- "gamble ~ "
  rhs.gamble <- c("1", "income", "status", "verbal")
  formulas.gamble <- sapply(seq_along(rhs.gamble), function(z) {
    paste(lhs.gamble, paste(rhs.gamble[1:z], collapse = "+"))
  })

  ## liste mit allen Modellen
  submodels.gamble <- lapply(formulas.gamble, lm, data = teengamb)
  # struktur des resultierenden Objekts:
  class(submodels.gamble)
  names(submodels.gamble)
  str(submodels.gamble, 1)

  ## SSEs der Einzelmodelle
  SSEs <- sapply(submodels.gamble, deviance)
  # alternativ/intuitiver:
  sapply(submodels.gamble, function(z) sum(residuals(z)^2))
  ## sequentielle QS als Differenz der SSEs Obermodell gegen Untermodell
  (seq.QS <- abs(diff(SSEs)))
  ## kumulative sequentielle QS
  (cum.seq.QS <- cumsum(round(seq.QS,digits=1)))


## ----teengamb_seq_part-b-------------------------------------------------
  ## SSM volles Modell
  SSM <- sum((fitted(mod) - mean(teengamb$gamble))^2)
  ## Vergleich mit Summe der sequentiellen QS
  all.equal(SSM, cumsum(seq.QS)[3])
  ## SSE volles Modell
  SSE <- sum(residuals(mod)^2)
  ## SST mit SSM und SSE aus Modell
  SST.mod <- SSM + SSE 
  ## Vergleich mit SST direkt aus Daten
  SST.data <- with(teengamb, sum((gamble - mean(gamble))^2))
  all.equal(SST.mod, SST.data)


## ----teengamb_seq_part-c-------------------------------------------------
  M2 <- mod
  M1 <- update(mod, .~. - status - verbal)
  (SSE.M1 <- deviance(M1))
  (SSE.M2 <- deviance(M2))
  (SSH <- SSE.M1 - SSE.M2)
  (F.stat <- SSH/2/SSE.M2 * (nrow(teengamb) - 3 - 1))
  ## mit linearHypothesis
  library(car)
  A <- rbind(c(0,0,1,0), c(0,0,0,1))
  linearHypothesis(M2, A, c(0,0), test = "F")


## ----buli_models---------------------------------------------------------
  load("../../daten/buli1314.Rdata")
  # library(ggplot)
  library(GGally)
  ggpairs(data = buli.1314[,-grep("Team",names(buli.1314))],
          lower=list(continuous="points",se=FALSE),
          upper=list(continuous="points",se=FALSE))
  # scatterplotMatrix(buli.1314[,-grep("Team",names(buli.1314))],smooth=FALSE)
  m1 <- lm(Points ~ Budget + AwayPointsRatio.previous + Points.previous,
           data = na.exclude(buli.1314))
  summary(m1)
  m2 <- lm(scale(Points) ~ scale(Budget) + scale(AwayPointsRatio.previous) + 
             scale(Points.previous), data = na.exclude(buli.1314))
  summary(m2)


## ----buli_graphs, fig=TRUE, fig.width=10---------------------------------
  library(sjPlot)
  ## Koeffizienten Plot: 
  sjp.lm(m1)
  ## Problem: Koeffizienten der metrischen Vars beziehen sich auf 
  # Veraenderung um eine Einheit, auch wenn eine "typische" Veraenderung 
  # z.B. eine SD sein könnte
  ## Loesung: Standardisierung der Daten
  sjp.lm(m2)

## ----dataSetup-----------------------------------------------------------
  lesen <- read.table("http://www.statistik.lmu.de/~abender/daten/lesen.txt", 
    header = TRUE)
  lm.fehler <- lm(Fehlerzahl ~ WieoftLesen, data = lesen)
  summary(lm.fehler)
  table(lesen$WieoftLesen)
  # Kategorie 5 kommt kaum vor -> mit Kat. 4 zusammenfassen
  lesen$WieoftLesen[lesen$WieoftLesen == 5] <- 4
  #als faktor umkodieren
  lesen$WieoftLesenF <- as.factor(lesen$WieoftLesen)


## ----lmLesenMittelwert---------------------------------------------------
  ## WieoftLesen Mittelwertsmodell
  lm.lesen.mw <- lm(Fehlerzahl ~ -1 + WieoftLesenF, data = lesen)
  coefs.mw <- coefficients(lm.lesen.mw)
  summary(lm.lesen.mw)


## ----groupMeans----------------------------------------------------
group.means <- tapply(lesen$Fehlerzahl, lesen$WieoftLesenF, mean)
group.means 


## ----contrastsMW---------------------------------------------------------
  ## Kontraste
  diag(1,4)
  ## Model matrix
  mm.mw <- model.matrix(lm.lesen.mw)
  # im Folgenden geben die ersten 4 Spalten die Kodierung in der Modelmatrix an, 
  # die letzte zum Vergleich die Originalvariable WieOftLesen
  head(cbind.data.frame(mm.mw, lesen$WieoftLesenF), 10)


## ----lmLesenEffekt-------------------------------------------------------
  ## WieoftLesen Effektkodierung
  lm.lesen.effect <- update(lm.lesen.mw, .~.+1, 
                            contrasts = list(WieoftLesenF = "contr.sum"))
  coefs.effect <- coefficients(lm.lesen.effect)
  summary(lm.lesen.effect)


## ----meanOfGroupMeans----------------------------------------------------
  mean.of.groupmeans<- mean(group.means)
  # Mittelwert der gruppenmittelwerte entspricht Intercept in lm.lesen.effect
  mean.of.groupmeans


## ----effectsByMeans------------------------------------------------------
  group.means - mean.of.groupmeans


## ----contrastEffektDefault-----------------------------------------------
  ## Kontraste
  contr.sum(levels(lesen$WieoftLesenF)) # "letzte" Kategorie ist "Referenz"-Kategorie
  ## model matrix
  head(cbind.data.frame(model.matrix(lm.lesen.effect), lesen$WieoftLesenF), 10)


## ----contrastsEffektCustom-----------------------------------------------
# moegliche Umkodierung der WieoftLesen Variable, Kategorie 1 als "Referenz"
lesen$WieoftLesenF2 <- factor(lesen$WieoftLesenF, levels = 4:1, labels = 4:1)
# aktualisieren des Models
lm.lesen.effect2 <- update(lm.lesen.effect, 
                           formula = .~.-WieoftLesenF + WieoftLesenF2, 
                           contrasts = list(WieoftLesenF2 = "contr.sum"))
summary(lm.lesen.effect2)


## ----lmLesenRef----------------------------------------------------------
  ## WieoftLesen Referenzkodierung
  lm.lesen.ref <- lm(Fehlerzahl ~ WieoftLesenF, data = lesen)
  coefs.ref <- coefficients(lm.lesen.ref)
  summary(lm.lesen.ref)


## ----globalMean----------------------------------------------------------
  mean(lesen$Fehlerzahl)


## ----contrastsRef--------------------------------------------------------
  ## Kontraste
  contr.treatment(levels(lesen$WieoftLesenF)) # "erste" Kategorie ist Referenz-Kategorie
  ## Model matrix
  mm.ref <- model.matrix(lm.lesen.ref)
  # im Folgenden geben die ersten 4 Spalten die Kodierung in der Modelmatrix an, 
  # die letzte zum Vergleich die Originalvariable WieOftLesen
  head(cbind.data.frame(mm.ref, lesen$WieoftLesenF), 10)


## ----coeffPlotRef, fig.width=4, fig.height=3.5, fig.keep="all"-----------
  ## Koeffizienten plots
  library(sjPlot)
  sjp.lm(lm.lesen.ref, showModelSummary=FALSE)
  
  ## oder standard plot aus termplot
  termplot(lm.lesen.ref, se=TRUE)
  
  ## oder mit effects package
  library(effects) 
  plot(allEffects(lm.lesen.ref))
  confint(lm.lesen.ref)


## ----Zusammenfassung-----------------------------------------------------
## Mittelwertmodell
coefs.mw
## Referenzkodierung
coefs.ref
# erwartungswerte
c(coefs.ref[1], coefs.ref[1] + coefs.ref[-1])
## effekt kodierung
coefs.effect
# erwartungswerte
c(coefs.effect[1] + coefs.effect[-1], sum(coefs.effect[1], -coefs.effect[-1]))


## ----linHyp--------------------------------------------------------------
  summary(lm.lesen.ref)
  A <- rbind(c(0,2,-1,0),c(0,3,0,-1))
  c <- c(0,0)
  linearHypothesis(lm.lesen.ref,A,rhs=c,test="F")


## ----Aufgabe3------------------------------------------------------------
## einlesen der Daten
weight <- read.table("http://www.statistik.lmu.de/~abender/unterlagen/limo14/Datensatz4.txt", 
  header = TRUE)

# Berechnung von Mittelwert und Standardabweichung fuer jede Kombination
tapply(weight$weightgain, list(weight$source, weight$amount), mean)
tapply(weight$weightgain, list(weight$source, weight$amount), sd)

# graphische Darstellung: Abweichung der Mittelwerte vom Gesamtmittelwert
# (Betrachtung der einzelnen Faktoren, nicht ihrer Kombination)
plot.design(weight)
mean(weight$weightgain)
tapply(weight$weightgain, weight$source, mean)
tapply(weight$weightgain, weight$amount, mean)


# zweifaktorielle Anova (mit Interaktion)
wg.aov <- aov(weightgain ~ source * amount, data=weight)

# Designmatrix
head(model.matrix(wg.aov))
cbind(model.matrix(wg.aov), weight[, 1:2])

summary(wg.aov)

# graphische Darstellung: Interaktion
interaction.plot(weight$amount, weight$source, weight$weightgain, legend=FALSE, 
  las = 1)
legend(1.5, 95, legend = levels(weight$source), title = "weight$source", 
  lty = c(2,1), bty="n")

# Parameterschaetzer
coef(wg.aov)
coef(lm.sa <- lm(weightgain ~ source * amount, data = weight))
## alternative visualisierung
plot(allEffects(lm.sa))



### Modell mit Effektkodierung
wg.aov.effect <- aov(weightgain ~ source + amount + source:amount, data=weight, 
    contrasts=list(source=contr.sum, amount=contr.sum))

# Designmatrix
head(model.matrix(wg.aov.effect))

summary(wg.aov.effect)
# (keine Veraenderung)

# Parameterschaetzer
coef(wg.aov.effect)


## ----contrastsSplit------------------------------------------------------
# definiere Kontraste-Matrix
contr.split <- lower.tri(matrix(1, nrow = 4, ncol = 4), diag = TRUE) * 1
contr.split
# man kann die selbst definierte Kontrastematrix der lm Funktion wie ueblich 
# ueber das contrasts Argument uebergeben 
#(Achtung: erste Spalte der Kontrastematrix weglassen)
lm.lesen.split <- update(lm.lesen.ref, 
        contrasts = list(WieoftLesenF = contr.split[, -1]))
coefs.split <- coefficients(lm.lesen.split)
summary(lm.lesen.split)


## ----splitGmeans---------------------------------------------------------
group.means# Gruppenmittelwerte
cumsum(coefs.split)# kumulativ aufsummierte Koeffizienten


## ----modelMatrixSplit----------------------------------------------------
# siehe Uebung fuer das Schema zur Kodierung der zu erzeugenden Dummy-Variablen
# hier wird in einem Schritt die Modelmatrix erzeugt:
mm.split <- (matrix(1:4, nrow = nrow(lesen), ncol = 4, byrow = TRUE) <= lesen$WieoftLesen) * 1
# und die Spaltennamen angepasst
colnames(mm.split) <- paste0("WieoftLesen", 1:4)
lesen.neu <- cbind.data.frame(Fehlerzahl = lesen$Fehlerzahl, mm.split)
# passe Modell mit allen Variablen im Datensatz an, lasse Intercept weg
lm.lesen.split2 <- lm(Fehlerzahl ~ . - 1, data = lesen.neu)
summary(lm.lesen.split2)

# pruefe ob koeffizientenschaetzer identisch
all(coefficients(lm.lesen.split2) == coefficients(lm.lesen.split))