---
title: "Lineares Modell, kategoriale Prädiktoren (linear model, categorial predictors)"
author: "Peter Zezula"
subtitle: "Varianzanalyse als Regressionsmodell"
output:
  html_document:
    toc: true
    toc_depth: 4
    toc_float: yes
    theme: paper
header-includes:
 - \usepackage{comment}
params:
 soln: TRUE
 todo: TRUE
mainfont: Helvetica
linkcolor: blue

---


### [Rmd](lm_cat.Rmd)

\[
E(y | x_1, x_2, ..., x_q) = \beta_0 + \beta_1x_1 + \beta_2x_2 + ... + \beta_qx_q
\]

bzw.

\[
y = \beta_0 + \beta_1x_1 + \beta_2x_2 + ... + \beta_qx_q + \varepsilon
\]




# Lineares Modell, kategoriale Vorhersagevariablen bzw. Prädiktoren (linear model, categorial predictors)


Dargestellt werden Varianzanalytische Designs mit Focus auf das Allgemeine Lineare Modell.

Es geht um den regressionsanalytischen Ansatz varianzanalytischer Designs.

Vorhersagevariablen sind kategorial, Reaktionsvariablen sind stetig (intervallskaliert).

In R haben kategoriale Vorhersagevariablen den Typ `factor()`


## Factor

Faktoren kodieren Gruppenzugehörigkeit.

Faktorstufen haben eine Reihenfolge, die ggf. verändert werden kann. Default ist alphabetisch nach den Namen der Faktorstufen.

Faktoren können in Zahlenwerte umgewandelt werden mit `as.numeric(factor.variable)`.

Zahlenwerte können in Faktoren umgewandelt werden mit `factor(variable)` bzw. `as.factor(variable)`

```{r}
library(tidyverse)
dd <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_wide.txt")[c('subj','gender','neutr', 'exam1', 'exam2')]
# we rename neutr to anx_neutr, it's state anxiety in neutral situation
dd <- dd %>% dplyr::rename(anx_neutr = neutr)
# add a numeric var
dd$n1.group <- c(rep(1, 90), rep(2, 90))
dd$n2.group <- c(rep(1:3, 60))

dd$fn1.group <- factor(dd$n1.group, levels=c(1,2), labels=c('one', 'two'))
# if we don't use codes, they are set to NA
dd$fn2.group.miss <- factor(dd$n2.group, levels=c(1,2), labels=c('one', 'two'))
dd$fn2.group <- factor(dd$n2.group, levels=c(1,2,3), labels=c('one', 'two', 'three'))

head(dd)
as.numeric(dd$fn2.group)
```



## Formate 'long' und 'wide'


Ein Spezifikum von R ist das sog. 'long format' beim Aufbau von Datensätzen.
Im Kern bedeutet das, dass bei R im Gegensatz zu vielen anderen bekannten Statistikpaketen Messwiederholungsdesigns (within subjects designs) keine andere Anordnung der Variablen erzwingen als Zwischengruppen-Designs (between subjects designs).

Unit long_wide 
[html](../long_wide/long_wide.html)
zum Umbau von Datensätzen.


## T-Test

Der Befehl `t.test()` ist sehr flexibel und kann sowohl mit dem Long- als auch mit dem Wide-Format umgehen.

```{r}
# comparison of different ways to calculate t.test()
dd <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_wide.txt")[c('subj','gender','neutr', 'exam1', 'exam2')]
# we rename neutr to anx_neutr, it's state anxiety in neutral situation
dd <- dd %>% dplyr::rename(anx_neutr = neutr)
# alternatives to get the same t-test
# default is t-test for independent groups (paired=FALSE)
t.test(dd$anx_neutr[dd$gender == "f"], dd$anx_neutr[dd$gender == "m"])

t.test(dd$anx_neutr ~ dd$gender)

# create wide format on the fly
dd.t <- data.frame(
  women = dd$anx_neutr[dd$gender == "f"],
  men   = dd$anx_neutr[dd$gender == "m"]
  )
head(dd.t)
# same test, now data in wide format but t.test is for independent groups
t.test(dd.t$women, dd.t$men)
# we can also get the repeated measures t-test
t.test(dd.t$women, dd.t$men, paired=TRUE)
```



## Äquivalenter Regressionsansatz

Die Regression kann mit dem `lm()` Befehl angefordert werden.

```{r}
dd <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_classic.txt")[c('nr','gender','neutr', 'exam1', 'exam2')]
# we rename neutr to anx_neutr, it's state anxiety in neutral situation
dd <- dd %>% dplyr::rename(anx_neutr = neutr)
t.test(anx_neutr ~ gender, data=dd)
summary(lm(anx_neutr ~ gender, data=dd))

```

Das geht auch mit Dummy-Codierung.
Je nach Dummy-Codierungsart müssen die Koeffizienten anders interpretiert werden.

```{r}
df <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_classic.txt")[c('nr','gender','neutr', 'exam1', 'exam2')]
# we rename neutr to anx_neutr, it's state anxiety in neutral situation
df <- df %>% dplyr::rename(anx_neutr = neutr)
t.test(anx_neutr ~ gender, data=df)
summary(lm(anx_neutr ~ gender, data=df))
# dummy coded - reference coding - women are reference group
dummy_ref <- rep(0, length(df$gender))
dummy_ref[df$gender == 'm'] <- 1
dummy_ref
summary(lm(df$anx_neutr ~ dummy_ref))
unique(model.matrix(lm(df$anx_neutr ~ dummy_ref)))

# dummy coded, contrast coding
dummy_effect <- rep(-1, length(df$gender))
dummy_effect[df$gender == 'm'] <- 1
dummy_effect
summary(lm(df$anx_neutr ~ dummy_effect))
#head(model.matrix(lm(df$anx_neutr ~ dummy_ref)))
unique(model.matrix(lm(df$anx_neutr ~ dummy_ref)))
```

`t.test()` benutzt im Hintergrund also eine Referenzcodierung.

Je nach Codierung müssen die Parameter anders interpretiert werden.


### Abbildung dazu

```{r}
dd <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_classic.txt")[c('nr','gender','neutr', 'exam1', 'exam2')]
# we rename neutr to anx_neutr, it's state anxiety in neutral situation
dd <- dd %>% dplyr::rename(anx_neutr = neutr)
require(ggplot2)
# create a plot object, define dataframe to use, add gender to differ in colour already in base plot
#ggplot(dd, x=gender, y=av_anx_neutr, aes(nation, bmi, fill = gender)) +
ggplot(dd, x=gender, y=anx_neutr, aes(gender, anx_neutr, fill = gender)) +
  stat_summary(fun.y=mean, geom="bar", position=position_dodge()) +
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", position=position_dodge(width=0.9), width=0.2)  
```


## T-Test für abhängige Gruppen

Kleinstes varianzanalytisches Messwiederholungsdesign.
Erreicht über `paired=TRUE` in `t.test()`.

Hier im Beispiel wird die neutrale Situation mit der Examenssituation verglichen.

```{r}
# comparison of different ways to calculate t.test()
dd <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_classic.txt")[c('nr','gender','neutr', 'exam1', 'exam2')]
# we rename neutr to anx_neutr, it's state anxiety in neutral situation
dd <- dd %>% dplyr::rename(anx_neutr = neutr)
# alternatives to get the same t-test
t.test(dd$anx_neutr, dd$exam1, paired=TRUE)

# create long format on the fly
dd.t <- data.frame(
  nr     = 1:2*length(dd$nr),
  gender = factor(c(dd$gender, dd$gender), levels = c(1,2), labels = c("f", "m")),
  sit    = c(rep('anx_neutr', length(dd$nr)), rep('exam1', length(dd$nr))),
  anx    = c(dd$anx_neutr, dd$exam1)
  )
head(dd.t)
# same test, now data in long format but t.test is for dependent groups
t.test(dd.t$anx[dd.t$sit == 'anx_neutr'], dd.t$anx[dd.t$sit == 'exam1'], paired=TRUE)
t.test(anx ~ sit, data=dd.t, paired=TRUE)

# visualization
require(ggplot2)
# create a plot object, define dataframe to use, add gender to differ in colour already in base plot
#ggplot(dd, x=gender, y=av_neutr, aes(nation, bmi, fill = gender)) +
ggplot(dd.t, x=sit, y=anx, aes(sit, anx, fill = sit)) +
  stat_summary(fun.y=mean, geom="bar", position=position_dodge()) +
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", position=position_dodge(width=0.9), width=0.2)  
```



## Einfaktorielle ANOVA


<!-- todo
Voraussetzungen prüfen

Modell anpassen

post hoc Kontraste

geplante Kontraste

Visualisierung

-->


## Beispiel STAI-State in Angstgruppen

Drei Gruppen miteinander vergleichen hinsichtlich der State-Angst in einer neutralen Situation (Statistik-Tutorium):

- 'normale' Studierende
- Studierende, die bereits in der Schule Statistikunterricht hatten
- Studierende mit Prüfungsangst

Kategoriale Variablen sind in R vom Variablentyp `factor`.
Umwandlung ist möglich durch den Befehl `as.factor()`.
Mehr als zweistufige Faktoren werden bei varianzanalytischen Modellen implizit in n-1 Dummy-Variablen umgewandelt.
Voreinstellung ist die Referenz-Kodierung. 
Referenz ist Faktorstufe mit dem alphabetisch ersten Namen (Klassen- bzw. Subgruppen-Label).
Mit `model.matrix()´ können wir uns die Designmatrix ansehen, die implizit generiert und verwendet wird.


```{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')]
# we rename neutr to anx_neutr, it's state anxiety in neutral situation
dd <- dd %>% dplyr::rename(anx_neutr = neutr)
head(dd)

# descriptives using summarySE
require(Rmisc)
Rmisc::summarySE(data=dd, "anx_neutr")
Rmisc::summarySE(data=dd, "anx_neutr", groupvars="group")

# grouping variables in R must be of class factor. We check that:
class(dd$group)
# aov() is a wrapper function for lm(), it's summary is more comparable to the 'usual' informations on ANOVAs
m.anx <- aov(anx_neutr ~ group, data=dd)
m.anx
# get the global F-Test
summary(m.anx)
# we can see that aov() actually is a wrapper for lm()
summary(aov(lm(anx_neutr ~ group, data=dd)))

# check means
require(effects)
effects::allEffects(m.anx)

# weights or coefficients, express the distance from reference group
m.anx$coefficients

# two coefficients for group because it has 3 levels due to dummy coding
# we can check dummy coding by inspecting the model.matrix()
unique(model.matrix(m.anx))

# visualisation via plot
plot(allEffects(m.anx), ask=FALSE)

# visualization ggplot
require(ggplot2)

ggplot(dd, x=group, y=anx_neutr, aes(group, anx_neutr, fill = group)) +
  #scale_fill_manual(values=c_colors) +
  stat_summary(fun.y=mean, geom="bar") +
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", width=0.2) +     # width refers to whiskers
  # put points on top of bars
  stat_summary(fun.data = mean_cl_normal, geom = "point", shape=21, size=5) +
  # add verbose axis labels
  labs(x = "anxiety group", y = "STAI state anxiety")

# we may also use predefined custom colors using scale_fill_manual()
c_colors <- c("#66ff66", "#ff6666", "#6666ff")
ggplot(dd, x=group, y=anx_neutr, aes(group, anx_neutr, fill = group)) +
  scale_fill_manual(values=c_colors) +
  stat_summary(fun.y=mean, geom="bar") +
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", width=0.2) +     # width refers to whiskers
  # put points on top of bars
  stat_summary(fun.data = mean_cl_normal, geom = "point", shape=21, size=5) +
  # add verbose axis labels
  labs(x = "anxiety group", y = "STAI state anxiety")

detach("package:Rmisc", unload=TRUE)
```



## Voraussetzungsüberprüfung

Der Levene-Test prüft die Homogenität der Varianzen über die Faktorstufen.
Bei Verletzung dieser Voraussetzung hilft Welch's F ratio `oneway.test()` 

```{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')]
# we rename neutr to anx_neutr, it's state anxiety in neutral situation
dd <- dd %>% dplyr::rename(anx_neutr = neutr)
# require(psych)
# psych::describe(dd$anx_neutr)
# by(dd$anx_neutr, list(dd$group), psych::describe)

# descriptives using summarySE
require(Rmisc)
Rmisc::summarySE(data=dd, "anx_neutr")
Rmisc::summarySE(data=dd, "anx_neutr", groupvars="group")

# test the differences in variability
require(car)
car::leveneTest(dd$anx_neutr, dd$group)
# if there are problems with heterogenity of variances
# compare to 'normal' test
summary(aov(anx_neutr ~ group, data=dd))
anova(aov(anx_neutr ~ group, data=dd))
# use oneway - Welch's F ratio
oneway.test(anx_neutr ~group, data=dd)

```


## Post hoc Tests

Problem multiplen Testens und damit Alpha-Inflation.
Korrektur des Alpha-Niveaus nach verschiedenen Methoden:

Bonferroni, Holm, Benjamini-Hochberg etc.
`help(p.adjust)`
Verschiedene Methoden unterscheiden sich in Bezug auf die Menge der durchzuführenden Post-Hoc-Tests, Power, wie konservativ sie sind etc.

```{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')]
# we rename neutr to anx_neutr, it's state anxiety in neutral situation
dd <- dd %>% dplyr::rename(anx_neutr = neutr)
dd$group <- factor(dd$group)
# do multiple comparisons with bonferroni correction of alpha
pairwise.t.test(dd$anx_neutr, dd$group, p.adjust.method = "bonferroni")
# Tukey and Dunnett
require(multcomp)
m.anx <- aov(anx_neutr ~ group, data=dd)
post.hocs <- multcomp::glht(m.anx, linfct = mcp(group = "Tukey"))
summary(post.hocs)
confint(post.hocs)
# Dunnett, base specifies group to compare with
# post.hocs <- glht(m.anx, linfct = mcp(group = "Dunnett"), base=1)
post.hocs <- glht(m.anx, linfct = mcp(group = "Dunnett"))
summary(post.hocs)
confint(post.hocs)
```



# Ausflug: Dummy-Kodierung und Vergleich mit dem regressionsanalytischen Ansatz

Zur Vergleichbarkeit des varianzanalytischen mit dem regressionsanalytischen Ansatzes braucht man die Dummy-Codierungen.
Es werden immer Faktorstufen - 1 Dummy-Variablen gebraucht, um **einen** Faktor zu kodieren. 


```{r}
dd <- read.delim("https://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_wide.txt")[c('subj','gender','group', 'neutr')]
# we rename neutr to anx_neutr, it's state anxiety in neutral situation
dd <- dd %>% dplyr::rename(anx_neutr = neutr)
dd.t <- dd
# do reference coding. reference group is 'normal students'. we need 2 dummy variables d.hi and d.lo
dd.t$d.hi <- 0
dd.t$d.lo <- 0
# low group contrasted in dummy var 2
dd.t$d.lo[dd.t$group == 'lo'] <- 1
# hi group contrasted in dummy var 1
dd.t$d.hi[dd.t$group == 'hi'] <- 1

# the default model
m.d1 <- aov(anx_neutr ~ group, data=dd.t)
# the reference contrast model
m.d2 <- aov(anx_neutr ~ d.hi + d.lo, data=dd.t)
# now we do it with lm()
m.d3 <- lm(anx_neutr ~ d.hi + d.lo, data=dd.t)
# get the global F-Test
summary(m.d1)
summary(m.d2)
summary(m.d3)
coefficients(m.d1)
coefficients(m.d2)
coefficients(m.d3)

```

- Die Quadratsummen von m.d2 summieren sich zu den Quadratsummen in m.d1.
- Die Residuals von m.d1 und m.d2 sind gleich.
- Der globale F-Wert von m.d1 und m.d3 sowie der p-value sind gleich.


# Ausflug: Reihenfolge von Faktorstufen

Faktoren haben in R spezielle Eigenschaften und werden in statistischen Modellen entsprechend behandelt.
Faktoren haben Stufen (`level`), also eine bestimmte Anzahl unterschiedlicher Ausprägungen, 
die sich über die Beobachtungen weg wiederholen.
Faktorstufen ordnen Beobachtungen Untergruppen zu.
`read.table()` mit seinen Unterbefehlen `read.delim()` oder `read.csv()` 
machen beim Einlesen einen "best guess", wenn sie Spalten vorfinden, die Zeichenketten (`str`) enthalten. 
Wiederholen sich Ausprägungen, wird direkt der Datentyp `factor` angenommen und gesetzt.
Tibbles, eingelesen via `read_tsv()`, `read_csv()`  aus der `library(tidyverse)`
setzen bei Spalten mit Zeichenketten immer den Datentyp `str`.
Die ggf. notwendige Umwandlung in den Typ `factor` muss explizit erfolgen.

Faktorstufen haben eine Reihenfolge.
Wenn nicht näher angegeben, ordnet R die Faktorstufen alphabetisch.
Die Reihenfolge der Faktorstufen hat Auswirkungen auf die Schätzung der Parameter in statistischen Modellen.

Die Reihenfolge kann explizit gesetzt werden über den Parameter `levels = c(...)` 
bei der Umwandlung in Faktoren. 
Mit dem Parameter `labels = c(...)` können wir die Faktorstufen auch neu benennen.
`levels` kann auch einfach einen Zahlenvektor bekommen 
und richtet sich dann nach der impliziten alphabetischen Reihenfolge.

Faktoren können auch ordered sein, was durch den Parameter `ordered = ...` gesetzt wird.
Dies gibt den Faktorstufen eine Reihenfolge, die Rangskalenniveau impliziert. 
Bei der Ausgabe eines ordered Faktors werden dann auch `<` zwischen die Faktorstufen gesetzt.

```{r warning = FALSE}
require(tidyverse)
dd <- read_tsv("https://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_wide.txt")[c('subj','gender','group', 'neutr')]
# we rename neutr to anx_neutr, it's state anxiety in neutral situation
dd <- dd %>% dplyr::rename(anx_neutr = neutr)
dd <- dd %>% arrange(desc(anx_neutr))
# we inspect group
str(dd)
head(dd$group)
# obviously, read_tsv sets gender and group to type chr
# if we want to use these in models, they have to be of type factor
# let's experiment with group
# default conversion, alphabetical order of factor levels
dd$group.fu <- factor(dd$group)
# if we want to have specific order of factor levels, we can set that using ordered
dd$group.f1 <- factor(dd$group, levels = c('lo', 'av', 'hi'))
# we can also change an existing factor's order and label it more verbos
dd$group.f2 <- factor(dd$group.f1, levels = c('av', 'lo', 'hi'), labels = c('average', 'low', 'high'))


# we take a look at the results
head(dd)
# we see <ord> as column type in our tibble
head(dd$group)
head(dd$group.f1, n=20)
head(dd$group.f2, n=10)
# the changed order is also reflected in 
as.numeric(head(dd$group.fu, n=30))
as.numeric(head(dd$group.f1, n=30))
as.numeric(head(dd$group.f2, n=30))

# if we use these factors now in linear models, they differ 
summary(lm(anx_neutr ~ group.fu, data = dd))
summary(lm(anx_neutr ~ group.f1, data = dd))
summary(lm(anx_neutr ~ group.f2, data = dd))

# coefficients reflect the different sequence, the reference group is different
lm(anx_neutr ~ group.fu, data = dd)$coefficients
lm(anx_neutr ~ group.f1, data = dd)$coefficients
lm(anx_neutr ~ group.f2, data = dd)$coefficients

# we also see the differences in the model.matrix()
head(model.matrix(lm(anx_neutr ~ group.fu, data = dd)))
head(model.matrix(lm(anx_neutr ~ group.f1, data = dd)))
head(model.matrix(lm(anx_neutr ~ group.f2, data = dd)))

```


# Ausflug: Vergleich Referenz-Kodierung vs. Effekt-Kodierung

```{r}
require(psych)
require(effects)
# 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')]
# we rename neutr to anx_neutr, it's state anxiety in neutral situation
dd <- dd %>% dplyr::rename(anx_neutr = neutr)
#psych::describe(dd)
dd.t <- dd
# do reference coding. reference group is 'normal students'. we need 2 dummy variables dr1 and dr2
dd.t$dr1 <- 0
dd.t$dr2 <- 0
# low group contrasted in dummy var 2
dd.t$dr2[dd.t$group == 'lo'] <- 1
# hi group contrasted in dummy var 1
dd.t$dr1[dd.t$group == 'hi'] <- 1
# the default model
m.anx <- aov(anx_neutr ~ group, data=dd.t)
# contrast matrix to see the dummy coding
unique(model.matrix(m.anx))
# the reference contrast model
m.anx.r <- aov(anx_neutr ~ dr1 + dr2, data=dd.t)
unique(model.matrix(m.anx.r))
# get the global F-Test
summary(m.anx)
summary(m.anx.r)
coefficients(m.anx)
coefficients(m.anx.r)
# compare coefficients with means to see their meaning
allEffects(m.anx)

# show it 
require(ggplot2)
ggplot(dd.t, x=group, y=anx_neutr, aes(group, anx_neutr, fill = group)) +
  stat_summary(fun.y=mean, geom="bar") +
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", width=0.2) +
  geom_hline(yintercept = c(coefficients(m.anx.r)[1]))


# effect coding in contrasts
dd.t$de1 <- 0
dd.t$de2 <- 0
# one group gets -1 for all dummy variables, here group 'av'
dd.t$de1[dd.t$group == 'av'] <- -1
dd.t$de2[dd.t$group == 'av'] <- -1
# low group contrasted in dummy var 2
dd.t$de2[dd.t$group == 'lo'] <- 1
# hi group contrasted in dummy var 1
dd.t$de1[dd.t$group == 'hi'] <- 1
# the common model
m.anx <- aov(anx_neutr ~ group, data=dd.t)
unique(model.matrix(m.anx))
# the reference contrast model
m.anx.e <- aov(anx_neutr ~ de1 + de2, data=dd.t)
unique(model.matrix(m.anx.e))
# get the global F-Test
summary(m.anx)
summary(m.anx.e)


# show it
ggplot(dd.t, x=group, y=anx_neutr, aes(group, anx_neutr, fill = group)) +
  stat_summary(fun.y=mean, geom="bar") +
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", width=0.2) +
  geom_hline(yintercept = c(coefficients(m.anx.e)[1]))


# the three models
coefficients(m.anx)
coefficients(m.anx.r)
coefficients(m.anx.e)
# compare coefficients of both, reference and effect dummy coding with means to see that they code means
allEffects(m.anx)
psych::describe(dd$anx_neutr)

# rm(dd.t)
```




# Ausflug: Reihenfolgeneffekt der Faktorstufen

Reference-Coding ist Default in R.
Die Referenzgruppe wird in Abhängigkeit von der Benennung der Faktorstufen gewählt (alphabetisch geordnet).
Dies hat Auswirkungen auf die Interpretation der Koeffizienten.
Die Ordnung von Faktorstufen kann aber festgelegt werden, alphabetische Sortierung ist nur Default.
Beispiele unten.

```{r}
require(psych)
require(effects)
# 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')]
# we rename neutr to anx_neutr, it's state anxiety in neutral situation
dd <- dd %>% dplyr::rename(anx_neutr = neutr)
dd$group <- factor(dd$group)
dd$n_group <- factor(c(rep('mm', 60), rep('lo',60), rep('hi',60)))
# we use the alphabetic sorting and call levels accordingly
dd$g_group <- factor(c(rep('zz', 60), rep('alo',60), rep('zhi',60))) 
m.1 <- aov(anx_neutr ~   group, data=dd)
m.2 <- aov(anx_neutr ~ n_group, data=dd)
m.3 <- aov(anx_neutr ~ g_group, data=dd)
m.1$coefficients
m.2$coefficients
m.3$coefficients
# default coding in aov() is reference coding against the first group name using alphabetical order
summary(m.1)
summary(m.2)
summary(m.3)



# show it 
require(ggplot2)
ggplot(dd, x=group, y=anx_neutr, aes(group, anx_neutr, fill = group)) +
  stat_summary(fun.y=mean, geom="bar", position=position_dodge()) +
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", position=position_dodge(width=0.9), width=0.2)  +
  geom_hline(yintercept = c(coefficients(m.1)[1]))
# n_group
ggplot(dd, x=n_group, y=anx_neutr, aes(n_group, anx_neutr, fill = n_group)) +
  stat_summary(fun.y=mean, geom="bar", position=position_dodge()) +
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", position=position_dodge(width=0.9), width=0.2)  +
  geom_hline(yintercept = c(coefficients(m.2)[1]))
# g_group
ggplot(dd, x=g_group, y=anx_neutr, aes(g_group, anx_neutr, fill = g_group)) +
  stat_summary(fun.y=mean, geom="bar", position=position_dodge()) +
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", position=position_dodge(width=0.9), width=0.2)  +
  geom_hline(yintercept = c(coefficients(m.3)[1]))

# we can redefine factor levels of existing factors to order levels due to our wishes
gr <- dd$group
# we have default alphabetical order
levels(gr)
# we can define one level to be the first by using command `relevel()`
gr <- relevel(gr, "lo")
levels(gr)
# we can also entirely define order of levels explicitly
levels(gr) <- c("hi", "lo", "av")
levels(gr)
```




# ezANOVA


Ein Paket zum komfortablen Berechnen von varianzanalytischen Modellen.
Default-Einstellungen sind ähnlich zu denen anderer Statistikprogramme, z. B. Quadratsummen Typ III.
Einfache Grafiken mit Streuungsinformationen zu den Gruppenmittelwerten.
Achtung: Im Gegensatz zu anderen Statistikprogrammen braucht ´library(ez)` für Messwiederholungsmodelle das auch ansonsten übliche Long-Format.
Siehe Unit [long_wide](../long_wide/long_wide.html)


```{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')]
# we rename neutr to anx_neutr, it's state anxiety in neutral situation
dd <- dd %>% dplyr::rename(anx_neutr = neutr)
#head(dd)
require(ez)
rr <- ezANOVA(dd, anx_neutr, subj, between=group)
# clearer with parameter names
ezANOVA(data=dd, dv=anx_neutr, wid=subj, between=group)
# we can also get a plot using
ezPlot(dd, dv=anx_neutr,  x=group, wid=subj, between = .(group))
detach("package:ez", unload=TRUE)
```

ges 	Generalized Eta-Squared measure of effect size (see in references below: Bakeman, 2005).

Effect Size Anova η2

Small | Medium | Large
--- | --- | ---
0.01 | 0.06 | 0.14 



# (geplante) Kontraste


siehe [html](./lm_cat_contr.html)



# Mehrfaktorielle ANOVA (keine Messwiederholung)


siehe lm_cat_factorial [html](./lm_cat_factorial.html)



# Covarianzanalyse

{} todo



# Reine Messwiederholungsdesigns


siehe lm_cat_repeated  [html](./lm_cat_repeated.html)



# Gemischte varianzanalytische Designs


siehe lm_cat_mixed  [html](./lm_cat_mixed.html)



# Ungleiche Zellgrößen (un-/balanced designs)


siehe lm_cat_unbalanced  [html](./lm_cat_unbalanced.html)



# Vergleich mit den Ergebnissen anderer Statistikpakete


Grundlagen finden sich in der Unit zu unbalanced designs.
[lm_cat_unbal.html]

[Link](https://myowelt.blogspot.de/2008/05/obtaining-same-anova-results-in-r-as-in.html)


# Bemerkungen Rcmdr


Rcmdr unterstützt einfaktorielle ANOVA unter 'Statistik | Mittelwertvergleiche'
'Statistik | Mittelwerte vergleichen | Einfaktorielle Varianzanalyse'

Grafische Darstellung hiervon in Rcmdr nach Auswahl des berechneten Modells:
'Modelle | wähle aktives Modell aus'
'Modelle | Grafiken | Effect plots'



# Beispiele, Übungen

Beispiele und Übungen finden sich bei den speziellen Unter-Units.

# thinking

- welche Eigenschaften unserer Beobachtungen (Versuchspersonen) werden zur Vorhersage verwendet?
- wie können die Koeffizienten (Gewichte) interpretiert werden?
  [grafische Veranschaulichung](lm_cat_factorial.html#grafische-darstellung)
- was wird vorhergesagt?
- was ist das Modell bei ANOVA-Designs?



# Links und Referenzen


Eine [Einführung](https://www.personality-project.org/r/r.guide/r.anova.html) unter Nutzung des Befehls `aov()` (base).

Viele [Anova Designs](https://www.southampton.ac.uk/~cpd/anovas/datasets/ANOVA%20in%20R.htm) incl. R-Code und Beispieldaten. 


[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](../lm_cat/lm_cat_factorial.html)
[lm_cat_unbalanced](../lm_cat/lm_cat_unbalanced.html)
[lm_cat_repeated](../lm_cat/lm_cat_repeated.html)
[lm_cat_mixed](../lm_cat/lm_cat_mixed.html)


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





# check

- continuous reaction variable, categorial explanatory variable(s)

- group membership is the predictive variable, it has to be of type `factor()`

- (Sub-)Group mean is estimated reaction variable

- parameters reflect group means or differences of group means or distances of group means to some other value f. e. grand mean

- factors and models
  - factor levels have a sequence, default is alphabetical
  - we can define factor level sequence 
    `factor(..., levels = c(...))`
  - explanatory variable of type `factor()`: in reference coding the first level is the reference group

- "classical" ANOVA designs report sums of squares (SS) and omnibus tests
  with more than 2 levels in an explanatory variable, we don't know which sub-differences are significant

- `aov()` and `car::Anova()` applied to a `lm()` result object provide ANOVA summary tables with SS and omnibus-tests for factors

- more than 2 levels of an eplanatory factor variable are split in dummy variable (levels - 1)

- using `model.matrix()`, even better combined with `unique()` shows us the dummy variables and their coding

- this is, why we have p-1 or q-1 df when we look at significance tests in categorial explanatory variables designs

- cell size differences matter
  [unbalanced designs](lm_cat_unbalanced.html)
  
- Meta-Packages oder Wrapper-Packages
  Eine komfortable Art für "klassische" varianzanalytische Messwiederholungs- und Mixed-Designs bieten  Meta-Packages wie 
  `library(ez)`, 
  die pipe-friendly `library(rstatx)` oder 
  `library(afex) mit lme4 Inclusion und flexiblen Effektstärken etc.

```{r check}
dd <- readr::read_tsv("https://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_wide.txt")
# we rename neutr to anx_neutr, it's state anxiety in neutral situation
dd <- dd %>% dplyr::rename(anx_neutr = neutr)
s.subj <- sample(dd$subj)
ggplot2::ggplot(dd, aes(x=s.subj, y=anx_neutr)) + geom_point() + geom_hline(yintercept=mean(dd$anx_neutr))
ggplot2::ggplot(dd, aes(x=s.subj, y=anx_neutr, color = group)) + geom_point() + 
  stat_summary(fun = mean, aes(x = 1, yintercept = ..y.., group = group), geom = "hline")
dd$f.group <- factor(dd$group, levels=c("lo", "av", "hi"))
mm <- lm(anx_neutr ~ f.group, data=dd)
unique(model.matrix(mm))
summary(mm)
# aov gives us the SS
aov(mm)
# package car's Anova gives us the classical SS omnibus tests and moreover Types of SS if we want
car::Anova(mm, type=3)

```



## Screencast 

- Multiple Regression - check
  [StudIP](https://studip.uni-goettingen.de/plugins.php/mediacastplugin/media/check/557da707f35cbedf27871677d1e8e982/e5b2d13c9ee2838d860b5106bbc8a5ec/126?v=lm_cat-de-2020-06-03.mp4) - 
[ownCloud](https://owncloud.gwdg.de/index.php/s/HLWrwrQOr11mHQP)


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



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

