---
title: "Lineares Modell, kategoriale Prädiktoren, Kontraste (geplante Kontraste)"
output:
  html_document:
    toc: true
    toc_depth: 4
    toc_float: yes
    theme: paper
header-includes:
 - \usepackage{comment}
params:
 soln: TRUE
 todo: FALSE
bibliography: ../../mv.bib
csl: ../../apa6.csl
---

### [Rmd](lm_cat_contr.Rmd)

- Kontrastkoeffizienten sind Gewichte für Mittelwertunterschiede
- Modellparameter * Ausprägung der Dummyvariablen für jeweils eine Beobachtung (ihre Faktorstufe)
- mit ihrer Hilfe können wir spezielle Hypothesen für Mittelwertunterschiede prüfen
- auch Trends, also Verläufe über die Faktor-Level bzw. Faktorstufen sind möglich
- die Parameter werden bei der Modellanpassung an die Daten geschätzt und auf Sig
- Ein einzelner Kontrast, also die Koeffizienten in einer Dummy-Variable, ist nicht an sich orthogonal oder nicht. Mehrere Kontraste können orthogonal zueinander sein. Sie sind es, wenn sie bestimmte Bedingungen erfüllen (vgl. Field et al 2012, chapter 10.4 p 414 ff). Orthogonalität betrifft also das Verhältnis von Kontrasten bzw. Dummy-Variablen zueinander.
- Eine Kontrastmatrix ist orthogonal, wenn wir die Kontrastkoeffizienten pro Zeile miteinander multiplizieren, die Zeilenergebnisse aufsummieren und das dann 0 ergibt. (Helmert-Kontrast)
(das Summenprodukt muss 0 sein)
- Eine Kontrast-Koeffizient in einer Dummy-Variable ist verknüpft mit einer bestimmten Subkategorie (Faktor-Level).  Eine 0 hier bewirkt, dass der Effekt (coefficient) auf die jeweilige Subkategorie nicht angewandt wird (Multiplikation mit 0 ergibt 0).
- Wenn ein Faktor-Level in allen Dummy-Variablen eine 0 hat, geht in die Schätzung nur noch der Intercept ein. Dies ist dann die Referenzgruppe und der Intercept ist der Mittelwert der Reaktionsvariable in dieser Gruppe.
- Kontrast-Koeffizienten werden mit den Effekten multipliziert um die Reaktionsvariable für eine spezielle Subgruppe zu berechnen.
- Eine Dummy-Variable definiert einen Effekt. -1 und +1 für zwei Subgruppen prüft den Unterschied zwischen genau diesen beiden Gruppen auf Signifikanz.  -2 für einen Level und +1 für zwei andere Level prüft den Unterschied zwischen dem Mittelwert einer Subgruppe und dem Mittelwert der Mittelwerte der beiden anderen Subgruppen auf Signifikanz.

Es gibt ein paar Standard-Kontraste, die car-Package definiert sind, darunter auch der oben verwendete `car::contr.Sum()` und `car::contr.Helmert()`.



# Vorteile und Prinzipien

Keine Alpha-Inflation - Gesamt-Risiko wird eingehalten.
Keine multiplen Vergleiche (wie bei multiplen T-Tests).
Daher auch keine Alpha Inflation, die korrigiert werden muss.
Allerdings gibt es hierzu unterschiedliche Meinungen (Vortrag Daniel Lakens, Juni 2017).

Aufteilen der erkärten Varianz.
Hypothesen vorher.
Kontraste -> dummy-Variablen -> stat. Modell
Kontraste drücken Hypothesen in Mittelwerten bzw. Mittelwertsunterschieden aus. 
Die aufgrund der Kontraste generierten Dummy-Variablen realisieren/sind das Modell.



# Standard Kontraste.

Standard-Kontraste existieren sowohl im Base-System als beispielsweise auch im `package(car)`.
im `package(car)` beginnen die Typen (Namen) der Kontraste mit einem Großbuchstaben.
`contr.sum()` vs. `contr.Sum()`

- contr.treatment bzw. car::contr.Treatment()
  default in R, Referenzkodierung
  T. steht im Dummy-Kontrast-Namen
- contr.SAS()
  wrapper für contr.treatment mit letztem Faktor-Level als Vergleichsgruppe
- contr.helmert() bzw car::contr.Helmert()
  Außer der letzten wird jede Kategorie (Faktorstufe) mit dem mittleren Effekt der folgenden Faktorstufen verglichen
- contr.sum() bzw car::contr.Sum()
  uses ‘sum to zero contrasts’. Gives orthogonal contrasts where you compare every level to the overall mean.
- contr.poly()
  gibt es nur im Base-System, nicht in `library(car)`

Design Matrix: Die Koeffizienten kann man sich für einen bestimmten Faktor ausgeben lassen. Intern wird eine entsprechende Menge an Dummy-Variablen gebildet. 

Kontraste können als Eigenschaft(property) an eine Variable vom Typ `factor()` angehängt werden.
Beispiel: im Data-Frame die vierstufige Variable V1 (Typ factor) mit einem referenzkodierten Kontrast versehen, 
dabei die Stufe 3 als Referenzkategorie wählen.
`contrasts(dataframe$V1) <- contr.Treatment(4, base=3)` 
Danach `aov()` aufrufen.

**Achtung** in obigem Beispiel gibt es Bezüge auf die Reihenfolge der Faktorstufen. 
Diese Reihenfolge richtet sich nach der **alphabetischen** Reihenfolge der Faktorstufen-Bezeichnungen (levels).
Die Reihenfolge kann man sich über den Befehl `levels()` ausgeben lassen. 
Für obigen Faktor V1 würde der Befehl lauten: `levels(dataframe$V1)`.
Transparenter wird das Ganze, wenn man beides miteinander verbindet:
`contrasts(dataframe$V1) <- contr.Treatment(levels(dataframe$V1), base=3)`.
Hier wird der Bezug auf den dritten Level des factor V1 unmittelbar klar.
 
Kontraste werden in Faktorstufen -1 Dummy-Variablen umgesetzt, die man sich ausgeben lassen kann `model.matrix()`.
Hierüber hat man eine Kontrollmöglichkeit über das Vorgehen der Berechnung ("Welche Faktorstufe ist Referenzgruppe?").

Eigendefinierte Spezialkontraste können über den Weg der Speicherung als Eigenschaft des Faktors ebenfalls spezifiziert werden. 
`contrasts(dataframe$V1) <- cbind(c(-1,0,-1,+2), c(-1,2,-2,0), c(-1,1,-1,1))`.
Kontrolliert werden kann der mit einem Faktor verknüpfte Kontrast über die Ausgabe der Variablen bzw. `contrasts(FaktorVariablenName)`

Je nach Art des Kontrastes müssen die Koeffizienten/Gewichte (coefficients) unterschiedlich interpretiert werden.
Die varianzanalytischen aggregierten Tests sind identisch.



# Beispiel STAI-State in einer neutralen Situation in verschiedenen Gruppen (Geschlecht, Angstgruppe)


Es werden hier verschiedene Arten von Kontrasten demonstriert, die nicht unbedingt alle inhaltlich Sinn machen.

-Default Einstellung ist contr.Treatment mit der ersten Faktorstufe als Referenzstufe. 
-Hier ist das die Gruppe mit dem Label 'av' (average student), die eher zufällig auch inhaltlich die adäquate Vergleichsgruppe ist.


```{r}
# get data
# dd <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_wide.txt")[c('subj','gender','group', 'neutr')]
dd <- readr::read_tsv("https://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_wide.txt")[c('subj','gender','group', 'neutr')]
dd$gender <- factor(dd$gender)
dd$group <- factor(dd$group)
# although existing in the base system in car package more transparent
require(car)
# levels of factor group
levels(dd$group)
# show design matrices for some standard contrasts

# standard contrasts
contr.Treatment(levels(dd$group))

# refernce group is 'lo' group, default is first level (group 1)
contr.Treatment(levels(dd$group), base=3)

# compare each group with last group, like SAS does/did
contr.SAS(levels(dd$group))

# compare each group with mean effect of all subsequent categories, except last group
contr.Helmert(levels(dd$group))

# contr.Sum effect coding
contr.Sum(levels(dd$group))

# polynomial contrasts
contr.poly(levels(dd$group))


# do an analysis with a special contrast
unique(model.matrix(aov(neutr ~ group, data=dd, contrasts = list(group = contr.Treatment))))
coefficients(aov(neutr ~ group, data=dd, contrasts = list(group = contr.Treatment)))

unique(model.matrix(aov(neutr ~ group, data=dd, contrasts = list(group = contr.Sum))))
coefficients(aov(neutr ~ group, data=dd, contrasts = list(group = contr.Sum)))

unique(model.matrix(aov(neutr ~ group, data=dd, contrasts = list(group = contr.poly))))
coefficients(aov(neutr ~ group, data=dd, contrasts = list(group = contr.poly)))

unique(model.matrix(aov(neutr ~ group, data=dd)))
coefficients(aov(neutr ~ group, data=dd))

require(Rmisc)
Rmisc::summarySE(data=dd, "neutr")
Rmisc::summarySE(data=dd, "neutr", groupvars="group")

# you may save a specific contrast for a factor
contrasts(dd$group) <- contr.Treatment(3, base=2)
# and control the current settings
dd$group['contrasts']


# you may define your own contrasts also
# compare all groups to the lo goup
levels(dd$group)
c1 <- c(1, 1, -2)
c2 <- c(-1, 1, 0)
contrasts(dd$group) <- cbind(c1, c2)
# check the design matrix
dd$group
dd$group['contrasts']
dd$group[2][1]
model.matrix(aov(neutr ~ group, data=dd))
coefficients(aov(neutr ~ group, data=dd))
summary(aov(neutr ~ group, data=dd))

```




# Aus Kontrasten resultierende Dummy-Variablen speichern/ausgeben

In R können die für ein Modell erzeugten Kontrastkoeffizienten bzw. Dummy-Variablen mit dem Befehl `model.matrix()` generiert werden.
Die Design-Matrix (contr.()) zeigt das Prinzip hinter der Dummy-Codierung.
Dadurch hat man Zugriff hierauf und kann sie beispielsweise zu einem Data-Frame hinzufügen oder speichern.
Analyse der Dummy-Matrix bzw. die Kontrolle der contrast-settings eines factors sind eine sichere Möglichkeit, das berechnete Modell nachzuvollziehen.

```{r}
# get data
dd <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_wide.txt")[c('subj','gender','group', 'neutr')]

m.anx <- aov(neutr ~ group, data=dd)

dd <- cbind(dd, model.matrix(neutr ~ group, data=dd))
head(dd)

# generate effect coded contrast matrix
mm.e <- model.matrix(aov(neutr ~ group, data=dd, contrasts = list(group = contr.Sum)))
```

R has a "sub-language" to translate formulas into design matrix, and in the spirit of the language you can take advantage of it. 



# Voreingestellte Kontraste (default contrasts)

Kontraste können auch als Default gesetzt sein und verändert werden

```{r}
# get data
#dd <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_wide.txt")[c('subj','gender','group', 'neutr')]
dd <- readr::read_tsv("https://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_wide.txt")[c('subj','gender','group', 'neutr')]
dd$gender <- factor(dd$gender)
dd$group <- factor(dd$group)

default.options <- options('contrasts')
default.options
coefficients(aov(neutr ~ group, data=dd))
# with changed default options the same call to aov() produces different results
options(contrasts=c('contr.sum','contr.poly'))
coefficients(aov(neutr ~ group, data=dd))

# reset to old default
options(contrasts=default.options[['contrasts']])
```





# 2x2 oder 1x4 Design mit Kontrasten

Im folgenden Beispiel wird ein eigentlich 2*2 varianzanalytisches Design auch über einen Faktor mit 4 Faktorstufen gerechnet, wobei mit Kontrasten gearbeitet wird. 
Bei diesem Beispiel werden auch spezifische, nicht orthogonale Kontraste angegeben.

```{r}
dd <- read.delim(header=F, text="
1	a1	b1	a1b1	10
2	a1	b1	a1b1	11
3	a1	b1	a1b1	9
4	a1	b1	a1b1	12
5	a1	b1	a1b1	10
6	a1	b2	a1b2	15
7	a1	b2	a1b2	17
8	a1	b2	a1b2	13
9	a1	b2	a1b2	14
10	a1	b2	a1b2	16
11	a2	b1	a2b1	9
12	a2	b1	a2b1	8
13	a2	b1	a2b1	10
14	a2	b1	a2b1	11
15	a2	b1	a2b1	12
16	a2	b2	a2b2	5
17	a2	b2	a2b2	7
18	a2	b2	a2b2	4
19	a2	b2	a2b2	6
20	a2	b2	a2b2	4")
colnames(dd) <- c('subj', 'A', 'B', 'AB','y')
head(dd)

# we set type factor for A, B, and AB
dd$A  <- factor(dd$A)
dd$B  <- factor(dd$B)
dd$AB <- factor(dd$AB)

# graph it
require(ggplot2)
# create a plot object, define dataframe to use, add gender to differ in colour already in base plot
ggplot(dd, x=B, y=y, aes(B, y, fill = A)) +
  stat_summary(fun=mean, geom="bar", position=position_dodge()) +
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", position=position_dodge(width=0.9), width=0.2)  

# the 'normal' 2*2 way, main effects only
summary(aov(y ~ A + B, data=dd))
# main effect B doesn't seem to exist

# the 'normal' way 2*2 Design full model including interaction
summary(aov(y ~ A * B, data=dd))
# we find a significant interaction


# alternatively put as a single factor design
# create a plot object, define dataframe to use, add gender to differ in colour already in base plot
ggplot(dd, x=AB, y=y, aes(AB, y)) +
  stat_summary(fun=mean, geom="bar", position=position_dodge()) +
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", position=position_dodge(width=0.9), width=0.2)  
require(ez)
# descriptives - the ez way
ezStats(data=dd, dv=.(y), wid=(subj), between=.(AB))
# grand mean and more descriptives of y
require(psych)
psych::describe(dd$y)

require(car)
# standard contrasts may be called directly in aov()
# reference coding, reference group is 1 (a1b1), coefficients are differences to reference group
coefficients(aov(y ~ AB, data=dd, contrasts = list(AB = contr.Treatment)))
# effect coding, coefficients are differences referred to grand mean
coefficients(aov(y ~ AB, data=dd, contrasts = list(AB = contr.Sum)))
# polynomial contrasts
coefficients(aov(y ~ AB, data=dd, contrasts = list(AB = contr.poly)))

# we may look at the contrast coefficients
contr.treatment(levels(dd$AB), base=3)
# contrast names in package(car) are a bit more verbose
contr.Treatment(levels(dd$AB), base=3)


# another way to specify contrasts is to make them properties of factor variables
# specify a reference contrast, take group number 3 (a2b1) as reference group, coefficients are differences to reference group
contrasts(dd$AB) <- contr.treatment(4, base=3)
# or more verbose
contrasts(dd$AB) <- contr.Treatment(4, base=3)

coefficients(aov(y ~ AB, data=dd))
summary.lm(aov(y ~ AB, data=dd))
             
# special contrasts work in a similar way, defining a contrasts property of a factor
# non orthogonal contrast
# both b1 (a1b1, a2b1) groups vs. group 4 (a2b2)
# both b1 (a1b1, a2b1) groups vs. group 3 (a1b2)
# both b1 (a1b1, a2b1) groups vs. both b2 groups (a1b2, a2b2)
contrasts(dd$AB) <- cbind(c(-1,0,-1,+2), c(-1,2,-2,0), c(-1,1,-1,1))
# check the design matrix
contrasts(dd$AB)
# check the dummy variables
cbind(dd, model.matrix(aov(y ~ AB, data=dd)) )
# contrasts specified express an effect contrast and coefficients express differences from grand mean
coefficients(aov(y ~ AB, data=dd))
summary.lm(aov(y ~ AB, data=dd))

# check reference contrasts and compare first group (a1b1) with the other ones, non orthogonal contrasts
contrasts(dd$AB) <- cbind(c(0,1,0,0), c(0,0,1,0), c(0,0,0,1))
# check the dummy variables
cbind(dd, model.matrix(aov(y ~ AB, data=dd)) )
# contrasts specified express an effect contrast and coefficients express differences from reference group (here a1b1)
coefficients(aov(y ~ AB, data=dd))
summary.lm(aov(y ~ AB, data=dd))


```



# Planned Contrasts

We extend the [sheet_anova]() a bit to understand.

subgroups defined by a factor variable ar ordered.
contrast coefficients apply to the subgroups in this order.
one contrast is one dummy variable.
the dummy variables can be put into the regression formula to estimate the reaction variable for a specific subgroup (subgroup combination in factorial designs).
a contrast vector applied to a factor variable using the command `contrasts()` is similar to a column of `model.matrix()` but ordered due to the factor levels.

Here we use a factor with 4 levels.
$$\hat{y}_i = b_0 + b_1 \cdot D_1 + b_2 \cdot D_2 + b_3 \cdot D_3 + b_4 \cdot D_4$$

```{r}
library(tidyverse)
dd <- read_delim("https://pzezula.pages.gwdg.de/data/Superhero.dat", delim = "\t")
# we define a new variable hero.f with factorized categorial levels and label them also
dd$hero.f <- factor(dd$hero,
  levels = c(1:4), labels =
    c("Spiderman", "Superman", "Hulk", "Ninja Turtle")
)
# all default
m.0 <- lm(injury ~ hero.f, dd)
summary.lm(m.0)
unique(model.matrix(m.0))

mm <- m.0
g.l <- "Superman"
# mm["(Intercept)"] + 


# the contrast of task 9 of anova sheet, only one contrast defined, the other two are set automatically and orthogonal to the first one
contrasts(dd$hero.f) <- c(1, -2, 1, 0)
m.1 <- lm(injury ~ hero.f, dd)
unique(model.matrix(m.1))

# a contrast
contrasts(dd$hero.f) <- cbind(c(3,-1,-1,-1), c(0, 2,-1,-1), c(0, 0, 1, -1) )
m.2 <- lm(injury ~ hero.f, dd)
unique(model.matrix(m.2))

# we double the coefficients
contrasts(dd$hero.f) <- cbind(c(6,-2,-2,-2), c(0, 4,-2,-2), c(0, 0, 2, -2) )
m.3 <- lm(injury ~ hero.f, dd)
unique(model.matrix(m.3))

# simply doubling the contrast-coefficients divides the model coefficients by 2
m.2$coefficients
m.3$coefficients

# the significance tests stay the same
summary(m.2)
summary(m.3)

# a standard contrast of package car
contrasts(dd$hero.f) <- car::contr.Sum(levels(dd$hero.f))
m.4 <- lm(injury ~ hero.f, dd)
unique(model.matrix(m.4))

# orthogonal contrast Helmert
contrasts(dd$hero.f) <- car::contr.Helmert(levels(dd$hero.f))
m.5 <- lm(injury ~ hero.f, dd)
unique(model.matrix(m.5))

require(Rmisc)
Rmisc::summarySE(data=dd, "injury", groupvars="hero.f")

```



# Beispiele Übungen / Exercises

{}todo




# Links und Referenzen


Package library(memisc): Convenience Methods for Setting Contrasts

Erklärung und Verstehbeispiel zu contr.sum, contr.helmert, contr.treatment
[link](https://www.clayford.net/statistics/tag/sum-contrasts/)

[Field (2012)][field] **Kapitel comparing means**

<!-- this does not appear in html-page. It's used only to arm links and references. -->
[field]: https://www.uk.sagepub.com/dsur/main.htm "Field et al (2012). Discovering Statistics using R"

[lm_cat_factorial]:  https://md.psych.bio.uni-goettingen.de/mv/unit/lm_cat/lm_cat_factorial.html
[lm_cat_unbalanced]: https://md.psych.bio.uni-goettingen.de/mv/unit/lm_cat/lm_cat_unbalanced.html
[lm_cat_repeated]:   https://md.psych.bio.uni-goettingen.de/mv/unit/lm_cat/lm_cat_repeated.html
[lm_cat_mixed]:      https://md.psych.bio.uni-goettingen.de/mv/unit/lm_cat/lm_cat_mixed.html

Sammlung R und Anova, Einstiegspunkte [https://www.r-statistics.com/tag/repeated-measures-anova/]


# check

- we can apply contrasts to factors

- various standard contrasts are predefined

- if contrasts are orthogonal, they account for independent parts of variance

- we can also define and use our own contrasts


```{r check}
```

`r if(!params$todo) {"\\begin{comment}"}`
<!--
todo: 
unbalanced einfaktoriell: Problem sind ggf. ungleiche Varianzen in den ungleich großen Gruppen. Weighted means or grand mean.
Ansonsten sind ungleiche Zellbesetzungen bei mehrfaktoriellen Designs ein Problem.
todo: include [this](https://www.clayford.net/statistics/tag/sum-contrasts/)
-->

`r if(!params$todo) {"\\end{comment}"}`


Version: `r format(Sys.time(), '%d %B, %Y %H:%M')`
