### R code from vignette source 'Blatt1_Lsg.Rnw'
### Encoding: ISO8859-1


###################################################
### code chunk number 2: dataSetup
###################################################
lesen <- read.table("lesen.txt",header = TRUE)
str(lesen)
summary(lesen)

table(lesen$WieoftLesen)# Kategorie 5 kommt kaum vor -> mit Kat. 4 zusammenfassen
# Kat1 = oft, Kat5 = fast nie
lesen$WieoftLesen[lesen$WieoftLesen == 5] <- 4
#als faktor umkodieren
lesen$WieoftLesenF <- as.factor(lesen$WieoftLesen)


###################################################
### code chunk number 3: lmLesenDummy
###################################################
# WieoftLesen dummy
# R kodiert (ungeordnete) faktoren per default dummy (mit erstem Level als Referenz)
lm.lesen.dummy <- lm(Fehlerzahl ~ WieoftLesenF, data = lesen)
coefs.dummy <- coefficients(lm.lesen.dummy)
summary(lm.lesen.dummy)
# Intercept: E-Wert in gruppe 1 (schueler die oft lesen)


###################################################
### code chunk number 4: contrastsDummy
###################################################
# contraste koennen auch manuell erzeugt werden
# dummy kodierung per contr.treatment
contr.treatment(levels(lesen$WieoftLesenF))
# extrahieren der modematrix
mm.dummy <- model.matrix(lm.lesen.dummy)
dim(mm.dummy)
head(cbind.data.frame(mm.dummy[, grep("^WieoftLesenF.", colnames(mm.dummy))], 
                lesen$WieoftLesenF), 10)


## Mittelwerte
global.mean <- mean(lesen$Fehlerzahl)
group.means <- tapply(lesen$Fehlerzahl, lesen$WieoftLesenF, mean)
group.means
# Koeffizienten geben Differenz zum Erwartungswert der Referenzkategorie an:
coefs.dummy[1] + coefs.dummy[-1]
group.means - group.means[1]
# mittlerer Erwartungswert
mean.groupmeans <- mean(group.means)
mean.groupmeans  # beta0 bei Effektkodierungs sein


###################################################
### code chunk number 5: contrastEffektDefault
###################################################
# WieoftLesen effektkodierung (per contr.sum)
contr.sum(levels(lesen$WieoftLesenF))# "falsche" Referenzkategorie
lm.lesen.effect <- update(lm.lesen.dummy, 
        contrasts = list(WieoftLesenF = "contr.sum"))
coefs.effect <- coefficients(lm.lesen.effect)
summary(lm.lesen.effect)
# Intercept: Mittelwert der Gruppenmittelwerte 
# (wenn keine weiteren kovariablen im Modell)



###################################################
### code chunk number 7: effectsByMeans
###################################################
group.means - mean.groupmeans
coefs.effect


###################################################
### code chunk number 8: contrastsEffektCustom
###################################################
# moegliche umkodierungen
lesen$WieoftLesenF2 <- factor(lesen$WieoftLesenF, levels = 4:1, labels = 4:1)

# kategorie 1 als "Referenz"
lm.lesen.effect2 <- update(lm.lesen.dummy, 
        formula = .~.-WieoftLesenF + WieoftLesenF2, 
        contrasts = list(WieoftLesenF2 = "contr.sum"))
summary(lm.lesen.effect2)
# Intercept und betas bleiben gleich (aber Vorsicht bei Interpretation)
# Belabelung bezieht sich nicht auf die belabelung des urspruenglichen Faktors



# Split-Kodierung:
# Kontraste in R nicht vorimplementiert -> selbst definieren: 
###################################################
### code chunk number 9: contrastsSplit
###################################################
# definiere Kontraste-Matrix
contr.split <- lower.tri(matrix(1, nrow = 4, ncol = 4), diag = TRUE) * 1
contr.split
lm.lesen.split <- update(lm.lesen.dummy, 
        contrasts = list(WieoftLesenF = contr.split[, -1]))
coefs.split <- coefficients(lm.lesen.split)
summary(lm.lesen.split)


###################################################
### code chunk number 10: splitGmeans
###################################################
group.means
cumsum(coefs.split)
# Intercept: erwartungswert in gruppe 1 (wie dummy kodierung)
# Betas: veraenderung bei uebergang von gruppe j-1 zu gruppe j
# (d.h. oft -> haeufig, haeufig -> regelmaessig, regelmaessig -> selten/fast nie)


###################################################
### code chunk number 11: modelMatrixSplit
###################################################
# model matrix (by hand)
mm.split <- (matrix(1:4, nrow = nrow(lesen), ncol = 4, byrow = TRUE) <= 
            lesen$WieoftLesen) * 1
colnames(mm.split) <- paste0("WieoftLesen", 1:4)
head(cbind(mm.split, WieoftLesen = lesen$WieoftLesenF), 10)

lm.lesen.split2 <- lm(Fehlerzahl ~ . - 1,# herausnehmen des Intercept
        data = cbind.data.frame(Fehlerzahl = lesen$Fehlerzahl, mm.split))
summary(lm.lesen.split2)


## modell mit weiteren kovariablen: 
lm.lesen <- lm(Fehlerzahl ~ sex + Jahrgang + Lesezeitmin + WieoftLesenF, 
        data = lesen)
summary(lm.lesen)
# beta0 entspricht nicht mehr dem Mittelwert von WieoftLesenF1
# beta0 = E(y|WieoftLesenF = 1, sex = 0, Jahrgang = 0, Lesezeitmin = 0)


###################################################
## code chunk number 12: modelAll
###################################################

### mit Dummykodierung ###
lm.lesen.dummy.all <- lm(Fehlerzahl ~ sex+Jahrgang+Lesezeitmin+WieoftLesenF, data = lesen)
coefs.dummy.all <- coefficients(lm.lesen.dummy.all)
summary(lm.lesen.dummy.all)

# Bsp: Erwartete Fehlerzahl für eine weibliche Person aus Gruppe 4 (selten/fast nie lesen), 
# der 4. Klasse, mit Lesezeitmin 30: 
coefs.dummy.all[1] + 30*coefs.dummy.all[4] + coefs.dummy.all[7]


### mit Effektkodierung ###
lm.lesen.effect.all <- lm(Fehlerzahl ~ sex+Jahrgang+Lesezeitmin+WieoftLesenF, data = lesen,
                          contrasts = list(WieoftLesenF = "contr.sum"))
coefs.effect.all <- coefficients(lm.lesen.effect.all)
summary(lm.lesen.effect.all)

# Bsp: Erwartete Fehlerzahl für eine weibliche Person aus Gruppe 4 (selten/fast nie lesen), 
# der 4. Klasse, mit Lesezeitmin 30: 
coefs.effect.all[1] + 30*coefs.effect.all[4] - sum(coefs.effect.all[5:7]) 
# Achtung: Referenz ist hier Kategorie 4! 


### mit Splitkodierung ###
lm.lesen.split.all <- lm(Fehlerzahl ~ sex+Jahrgang+Lesezeitmin+WieoftLesenF, data = lesen,
                         contrasts = list(WieoftLesenF = contr.split[, -1]))
coefs.split.all <- coefficients(lm.lesen.split.all)
summary(lm.lesen.split.all)

# Bsp: Erwartete Fehlerzahl für eine weibliche Person aus Gruppe 4 (selten/fast nie lesen), 
# der 4. Klasse, mit Lesezeitmin 30: 
coefs.split.all[1] + 30*coefs.split.all[4] + sum(coefs.split.all[5:7])