---
title:  "Lineares Modell, kategoriale Prädiktoren, gemischte Modelle (mixed designs)"
author: "Peter Zezula"
subtitle: "Unit, Examples, etc"
output:
  html_document:
    toc: true
    toc_depth: 4
    toc_float: yes
    theme: paper
header-includes:
 - \usepackage{comment}
params:
 soln: True
 todo: False
mainfont: Helvetica
linkcolor: blue
editor_options: 
  chunk_output_type: console
---


### [Rmd](lm_cat_mixed.Rmd)

<!-- TODO include ANOVA Meta-Packages rstatx (pipe-friendly) and afex (Christian, lme4 inclusion, flexible Effektstärken etc.)
-->

# Gemischte varianzanalytische Designs


Gemischt werden Messwiederholungsfaktoren (einer oder mehrere) mit Zwischensubjektfaktoren (einer oder mehrere).

Voraussetzungen (wie bei Messwiederholungsdesigns): Zusätzlich zu Varianzenhomogenität über die Faktorstufen kommt noch die Homogenität der Covarianzen (Korrelationen)

Homogenität der Covarianzen wird getestet mit Mauchly's Test (bei Varianzhomogenität war es Levene-Test).
Spherizität (Gleichheit der Varianzen der Differenzen zwischen den Faktorstufen)

Bei Verletzung der Homogenität der Covarianzen: 

- Greenhouse-Geisser Korrektur
- Huyn-Feldt Korrektur


# ezANOVA aus library(ez)

## Ansatz ezANOVA: Beispiel 3x3, STAI-State von verschiedenen Gruppen in drei unterschiedlichen Situationen


Die Personen in den Gruppen gehören entweder zum unteren Drittel, zur Mitte oder zum oberen Drittel hinsichtlich ihrer mathematischen Fähigkeiten. 
Bei den verglichenen Situationen steht eine neutrale Situation zwei Prüfungssituationen gegenüber.

Im Beispiel: 

- Ein dreistufiger Between-Subjects-Faktor (Gruppe).
- Ein dreistufiger Within-Faktor (Situation).

Bei Messwiederholungs- und Mixed-Designs gibt `ezANOVA()` den
Mauchly's Test for Sphericity sowie die Greenhouse-Geisser und die Huyn-Feldt Korrektur mit aus.


```{r}
dd <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_wide.txt")
# convert it to long format 
# require(reshape2)
# dd.l <- melt(dd, id=c('subj', 'gender', 'group'), variable.name="sit", value.name="anx", measured=c('neutr', 'exam1', 'exam2'))
# head(dd.l)

# convert it to long format using tidyr
require(tidyr)
dd.l <- dd %>% tidyr::gather(key=sit, value=anx, neutr, exam1, exam2)
dd.l$subj <- factor(dd.l$subj)
dd.l$sit <- factor(dd.l$sit)
head(dd.l)


# fit model using ezANOVA, show anova table and 
require(ez)
ezANOVA(data=dd.l, dv=.(anx), wid = .(subj), between=.(group), within = .(sit), detailed=TRUE, type=2)

# descriptives - the ez way
ezStats(data=dd.l, dv=.(anx), wid = .(subj), between=.(group), within = .(sit), type=2)

# ez offers ezPlot() to visualize effects
ezPlot(data=dd.l, dv=.(anx), wid = .(subj), between=.(group), within = .(sit), x=group, split=.(sit))
ezPlot(data=dd.l, dv=.(anx), wid = .(subj), between=.(group), within = .(sit), x=sit, split=.(group))


```

'ges' im Output:  Generalized Eta-Squared measure of effect size
see: Bakeman, R. (2005). Recommended effect size statistics for repeated measures designs. 
Behavior Research Methods, 37 (3), 379-384.

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 



## Ansatz ezANOVA: Beispiel 2x3x3, STAI-State von verschiedenen Gruppen getrennt nach Geschlecht in drei unterschiedlichen Situationen


Zusätzlich zum obigen Design wird Geschlecht mit aufgenommen.

Im Beispiel: 

- Ein dreistufiger Between-Subjects-Faktor (Gruppe).
- Ein zweistufiger Between-Subjects-Faktor (Geschlecht)
- Ein dreistufiger Within-Faktor (Situation).


```{r}
dd <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_wide.txt")

# convert it to long format 
require(reshape2)
dd.l <- melt(dd, id=c('subj', 'gender', 'group'), variable.name="sit", value.name="anx", measured=c('neutr', 'exam1', 'exam2'))

# fit model using ezANOVA
require(ez)
ezANOVA(data=dd.l, dv=.(anx), wid = .(subj), between=.(gender, group), within = .(sit), detailed=TRUE, type=2)

# descriptives
ezStats(data=dd.l, dv=.(anx), wid = .(subj), between=.(gender, group), within = .(sit), type=2)
# sequence matters for output
ezStats(data=dd.l, dv=.(anx), wid = .(subj), between=.(group, gender), within = .(sit), type=2)

# graphics
ezPlot(data=dd.l, dv=.(anx), wid = .(subj), between=.(gender, group), within = .(sit), x=group, split=.(sit))
# col gives us multiple plots
ezPlot(data=dd.l, dv=.(anx), wid = .(subj), between=.(gender, group), within = .(sit), x=group, col=.(gender), split=.(sit))
ezPlot(data=dd.l, dv=.(anx), wid = .(subj), between=.(gender, group), within = .(sit), x=group, col=.(sit), split=.(gender))

```



## Ansatz ezANOVA: detaillierte Informationen über das aov-Objekt

Beispiel 2x3x3, STAI-State von verschiedenen Gruppen getrennt nach Geschlecht in drei unterschiedlichen Situationen

Im Beispiel: 

- Ein dreistufiger Between-Subjects-Faktor (Gruppe).
- Ein zweistufiger Between-Subjects-Faktor (Geschlecht)
- Ein dreistufiger Within-Faktor (Situation).

Wir generieren ein aov-Objekt über den Parameter `return_aov = TRUE`.
Dies gibt uns die Möglichkeit, Details des generierten Modells abzurufen.


```{r}
dd <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_wide.txt")
# convert it to long format 
require(reshape2)
dd.l <- melt(dd, id=c('subj', 'gender', 'group'), variable.name="sit", value.name="anx", measured=c('neutr', 'exam1', 'exam2'))

# fit model using ezANOVA
require(ez)
m.ez <- ezANOVA(data=dd.l, dv=.(anx), wid = .(subj), between=.(gender, group), within = .(sit), detailed=TRUE, type=2, return_aov = TRUE)

# summary shows list elements we can access
summary(m.ez)

# we access elements by number
m.ez[2]
# or by name
m.ez['ANOVA']

# we can get the model.matrix for our design. 
m.mx <- model.matrix(m.ez['aov']$aov)
# As it has subj-1 dummy columns to code effect subject, it is big
# columns are:
m.mx[1,]
# we can look at some significant parts to understand dummy coding
m.mx[c(2:5, 535:540),c(1:8, 185:187, 196:197)]
# ... the major part of subject coding and some of the interaction codings are cut out

# we store aov element for convenience
o.aov <- m.ez[['aov']]
# and acccess some details of effect "subj"
o.aov[['subj']]$coefficients
head(o.aov[['subj']]$residuals)
head(o.aov[['subj']]$fitted.values)
head(o.aov[['subj']]$effects)

# descriptives
ezStats(data=dd.l, dv=.(anx), wid = .(subj), between=.(gender, group), within = .(sit), type=2)

```



# library(rstatix)

`library(rstatix)` ist Pipe-friendly

```{r}
# we can do classical mixed factorial design using library(rstatix)
# we stay with omnibus tests for all effects
# library(rstatix)
library(tidyverse)
dd.w <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_wide.txt")
# convert it to long format 
# require(reshape2)
dd <- reshape2::melt(dd.w, id=c('subj', 'gender', 'group'), variable.name="sit", value.name="anx", measured=c('neutr', 'exam1', 'exam2'))

dd %>% rstatix::factorial_design(dv = anx, wid = subj, within = sit, between = c(group, gender)) -> design
design # the defined design 

r.f <- car::Anova(design$model, idata = design$idata, idesign = design$idesign, type = 3)
summary(r.f, multivariate = FALSE)

```




# library(afex)

Wrapper Package, nutzt library(lme4) bei den Multilevel-Teilen

```{r}
library(tidyverse)
dd.w <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_wide.txt")
# convert it to long format 
# require(reshape2)
dd <- reshape2::melt(dd.w, id=c('subj', 'gender', 'group'), variable.name="sit", value.name="anx", measured=c('neutr', 'exam1', 'exam2'))
# pure repeated measure design
m.1 <- afex::mixed(anx ~ sit + (1|sit) + (1|subj), data=dd)
m.1
# mixed design
m.2 <- afex::mixed(anx ~ gender + group + sit + (1|sit) + (1|subj), data=dd)
m.2
summary(m.2)$varcor
# we stay with the problem of having omnibus tests only
```




# Multilevel Ansatz mit `library(nlme)` und dem Befehl `nlme::lme()`. 

Wir können Daten mit Multi-Level Ansätzen modellieren, wenn sie hierarchisch organisiert sind.
Hier sind gruppenspezifische Parameterschätzungen möglich.
Auch Gruppen-Variablen können in die Modelle mit einbezogen werden.

Wir können "klassische" varianzanalytische Messwiederholungsdesigns und Mixed-Designs als Multi-Level-Designs betrachten.
Sie sind dann ein Spezialfall einer Multi-Level-Analyse, in der nur random Intercepts betrachtet werden.


Multi-Level-Ansatz wie in Unit [ml](../ml/ml.html)

Benutzt wird `lme()` aus dem `package(nlme)`.


Ein dreistufiger Between-Subjects-Faktor (Gruppe).
Ein zweistufiger Between-Subjects-Faktor (Geschlecht).
Ein dreistufiger Within-Faktor (Situation).



```{r}
dd <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_wide.txt")
# convert it to long format 
# require(reshape2)
# dd.l <- melt(dd, id=c('subj', 'gender', 'group'), variable.name="sit", value.name="anx", measured=c('neutr', 'exam1', 'exam2'))

# convert it to long format using tidyr
require(tidyr)
dd.l <- dd %>% tidyr::gather(key=sit, value=anx, neutr, exam1, exam2)
dd.l$subj <- factor(dd.l$subj)
dd.l$sit <- factor(dd.l$sit)
head(dd.l)


# fit model using nlme::lme
require(nlme)

# m.1 <- lme(fixed = anx ~ 1 + gender * group, random = ~1 + sit | subj, data=dd.l, method='ML')

# null model to compare with
m.0 <- lme(fixed= anx ~ 1, random = ~1 | subj, data=dd.l, method='ML')
# random intercept model, this is equivqlent to ezANOVA example above
m.1 <- lme(fixed = anx ~ 1 + gender * group * sit, random = ~1 | subj, data=dd.l, method='ML')
# we may also adapt a random intercept random slope model:
m.2 <- lme(fixed = anx ~ 1 + gender * group * sit, random = ~1 + sit | subj, data=dd.l, method='ML')

# compare models
anova(m.0, m.1)
# ... the ANOVA equvalent model differs from the zero model
anova(m.1, m.2)
# ... integrating random slope in m.2 is even better and significant different from m.1

# see details
summary(m.1)
summary(m.2)

# we might want to compare this with ezANOVA results (commented out)
# m.ez <- ezANOVA(data=dd.l, dv=.(anx), wid = .(subj), between=.(gender, group), within = .(sit), detailed=TRUE, type=2, return_aov = TRUE)
# summary(m.ez)
# m.ez

# we look at group descriptives
require(Rmisc)
summarySE(data=dd.l, "anx", c("group", "gender", "sit"))
detach("package:Rmisc", unload=TRUE)

# we can get at fixed effects
fixef(m.1)

# and at random effects, only the first few here because we have many subjects
head(ranef(m.1))

# we can predict individual values of level 1
head(predict(m.1, level=1))
dd.l$pred.vals <- predict(m.1, level=1) # Vorhersage auf Level 1

library(ICC)
ICCest(x = subj, y = anx, data = dd.l)


```




# Multilevel Ansatz mit `lmer4::lmer()` bzw. `library(lmerTest)`


Wie vorgestellt in 
[Field (2012, p. 573)][field] **Repeated-Measures Designs(GLM 4)**

Benutzt wird `lmer()` aus dem `package(lmer4)`.
Wegen der Diskussion um die korrekten P-Werte bei lmer4 Designs haben die Autoren von `library(lmer4)` beschlossen, _keine_ p-Werte mehr auszugeben.
Unter bestimmten Umständen (classical design, balanced, nested, etc.) ist die Berechnung von p-Werten weniger strittig.
Wir können sie erhalten, wenn wir `library(lmerTest)` benutzen, die ihrerseits auf `library(lmer4)` aufbaut.

- Ein dreistufiger Between-Subjects-Faktor (Gruppe).
- Ein zweistufiger Between-Subjects-Faktor (Geschlecht).
- Ein dreiftufiger Within-Subjects-Faktor (Situation).



```{r}
# # get data
# #require(reshape2)
# #require(psych)
# # require(lme4) # is loaded and used by lmerTest
# # require(lmerTest)
# dd <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_wide.txt")
# # we rename situation neutr to become reference group. This might be closer to our research design.
# # releveling has no impact in lme summaries.
# dd$aneutr <- dd$neutr
# dd$neutr <- NULL
# 
# # convert it to long format 
# # dd.l <- melt(dd, id=c('subj', 'gender', 'group'), variable.name="sit", value.name="anx", measured=c('neutr', 'exam1', 'exam2'))
# 
# # convert it to long format using tidyr
# require(tidyr)
# dd.l <- dd %>% tidyr::gather(key=sit, value=anx, aneutr, exam1, exam2)
# dd.l$subj <- factor(dd.l$subj)
# dd.l$sit <- factor(dd.l$sit)
# head(dd.l)


dd <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_wide.txt")
# convert it to long format using tidyr
require(tidyr)
dd.l <- dd %>% tidyr::gather(key=sit, value=anx, neutr, exam1, exam2)
dd.l$subj <- factor(dd.l$subj)
dd.l$sit <- factor(dd.l$sit)
head(dd.l)

#require(lme4)
require(lmerTest)
# m.f <- lme4::lmer(anx ~ group + gender + (1 | subj) + (1 | sit), data=dd.l)
# m.f
# summary(m.f)

m.t <- lmerTest::lmer(anx ~ group + gender + (1 | subj) + (1 | sit), data=dd.l)
m.t
summary(m.t)

```



# Grafik mit ggplot2

Hier ein expliziter Weg.

`library(ez)` erzeugt mit `ezPlot()` ebenfalls ggplot-Grafiken.


```{r}
require(reshape2)
dd <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_wide.txt")
# convert it to long format 
dd.l <- melt(dd, id=c('subj', 'gender', 'group'), variable.name="sit", value.name="anx", measured=c('neutr', 'exam1', 'exam2'))


require(ggplot2)
ggplot(dd.l, aes(x=sit, y=anx)) +
  # define colors
  scale_fill_manual(values=c("#FFFFFF", "#CCCCCC", '#AAAAAA')) +
  # create white background etc
  #theme_bw() +
  geom_point() +
  # add verbose axis labels
  labs(x = "Anxiety-Group", y = "STAI-Trait")


ggplot(dd.l, aes(x=sit, y=anx, fill=sit)) +
  # define colors
  scale_fill_manual(values=c("#FFFFFF", "#CCCCCC", '#AAAAAA')) +
  # create white background etc
  #theme_bw() +
  stat_summary(fun.y=mean, geom="bar", position=position_dodge(width=0.7), colour='black') +
  # add error information
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", position=position_dodge(width=0.7), width=0.2) +
  # put big points on top of bars
  stat_summary(fun.data = mean_cl_normal, geom = "point", shape=21, size=5, position=position_dodge(width=0.7)) +
  # add verbose axis labels
  labs(x = "Anxiety-Group", y = "STAI-Trait")


```



# lmer formulae und Erklärungen

[Quelle](http://stats.stackexchange.com/questions/13166/rs-lmer-cheat-sheet)



    What's the difference between (~1 +....) and (1 | ...) and (0 | ...) etc.?

Say you have variable V1 predicted by categorical variable V2, which is treated as a random effect, and continuous variable V3, which is treated as a linear fixed effect. Using lmer syntax, simplest model (M1) is:

V1 ~ (1|V2) + V3

This model will estimate:

P1: A global intercept

P2: Random effect intercepts for V2 (i.e. for each level of V2, that level's intercept's deviation from the global intercept)

P3: A single global estimate for the effect (slope) of V3

The next most complex model (M2) is:

V1 ~ (1|V2) + V3 + (0+V3|V2)

This model estimates all the parameters from M1, but will additionally estimate:

P4: The effect of V3 within each level of V2 (more specifically, the degree to which the V3 effect within a given level deviates from the global effect of V3), while enforcing a zero correlation between the intercept deviations and V3 effect deviations across levels of V2.

This latter restriction is relaxed in a final most complex model (M3):

V1 ~ (1+V3|V2) + V3

In which all parameters from M2 are estimated while allowing correlation between the intercept deviations and V3 effect deviations within levels of V2. Thus, in M3, an additional parameter is estimated:

P5: The correlation between intercept deviations and V3 deviations across levels of V2

Usually model pairs like M2 and M3 are computed then compared to evaluate the evidence for correlations between fixed effects (including the global intercept).

Now consider adding another fixed effect predictor, V4. The model:

V1 ~ (1+V3*V4|V2) + V3*V4

would estimate:

P1: A global intercept

P2: A single global estimate for the effect of V3

P3: A single global estimate for the effect of V4

P4: A single global estimate for the interaction between V3 and V4

P5: Deviations of the intercept from P1 in each level of V2

P6: Deviations of the V3 effect from P2 in each level of V2

P7: Deviations of the V4 effect from P3 in each level of V2

P8: Deviations of the V3-by-V4 interaction from P4 in each level of V2

P9 Correlation between P5 and P6 across levels of V2

P10 Correlation between P5 and P7 across levels of V2

P11 Correlation between P5 and P8 across levels of V2

P12 Correlation between P6 and P7 across levels of V2

P13 Correlation between P6 and P8 across levels of V2

P14 Correlation between P7 and P8 across levels of V2

Phew, That's a lot of parameters! And I didn't even bother to list the variance parameters estimated by the model. 
What's more, if you have a categorical variable with more than 2 levels that you want to model as a fixed effect, 
instead of a single effect for that variable you will always be estimating k-1 effects (where k is the number of levels), 
thereby exploding the number of parameters to be estimated by the model even further.



# Ein paar unfertige Versuche


## pitch and attractiveness


[http://seriousstats.wordpress.com/tag/repeated-measures-anova/]



```{r}
# get data
dd <- read.csv('http://www2.ntupsychology.net/seriousstats/pitch.csv')
head(dd)
library(lme4)
pitch.me <- lmer(pitch ~ base + attract + (1|Face) + (1|Participant), data=dd)
pitch.me


```




## Example: Politeness and pitch

http://www.bodowinter.com/tutorial/bw_LME_tutorial.pdf

The difference in politeness level is represented in the column called “attitude”.
In that column, “pol” stands for polite and “inf” for informal. 
Sex is represented as “F” and “M” in the column “gender”. 
The dependent measure is “frequency”,
which is the voice pitch measured in Hertz (Hz).
To remind you, higher values mean higher pitch.

The interesting random effects for us are in the column “subject” and “scenario”,
the latter being the name of the item column (remember the different scenarios
like “asking for a favor”?).

```{r eval=FALSE}
# get data
# dd = read.csv("http://www.bodowinter.com/tutorial/politeness_data.csv")
dd <- readr::read_csv("http://md.psych.bio.uni-goettingen.de/mv/data/div/politeness_data.csv")
dd$subject <- factor(dd$subject)
dd$gender <- factor(dd$gender)
dd$attitude <- factor(dd$attitude)

require(psych)
psych::describe(dd)
# missing data?
which(is.na(dd)==T)
# a first boxplot
boxplot(frequency ~ attitude*gender,
col=c("white","lightgray"),dd)


# pitch.dat <- read.csv('http://www2.ntupsychology.net/seriousstats/pitch.csv')
# head(pitch.dat)
# library(lme4)
# pitch.me <- lmer(pitch ~ base + attract + (1|Face) + (1|Participant), data=pitch.dat)
# pitch.me


```


`lmer()` muss random effect haben, sonst läuft es nicht.
Modell mit Vp-Faktor und Faktor scenario:

```{r eval=FALSE}
# dd = read.csv("http://www.bodowinter.com/tutorial/politeness_data.csv")
dd <- readr::read_csv("http://md.psych.bio.uni-goettingen.de/mv/data/div/politeness_data.csv")
dd$subject <- factor(dd$subject)
dd$gender <- factor(dd$gender)
dd$attitude <- factor(dd$attitude)

# fit model
politeness.model = lmer(frequency ~ attitude + (1|subject) + (1|scenario), data=dd)
# look at the model
politeness.model
```

Note that the lmer() function (just like the lm() function in tutorial 1) took
whatever comes first in the alphabet to be the reference level. “inf” comes before
“pol”, so the slope represents the change from “inf” to “pol”. 


Aufnahme von Geschlecht ins Modell.
Geschlecht ist ein 'fixed effekt'. 
Es gibt zwei Stufen (Ausprägungen) und wir erwarten im Mittel eine höhere Stimme bei Frauen.

```{r eval=FALSE}
# dd = read.csv("http://www.bodowinter.com/tutorial/politeness_data.csv")
dd <- readr::read_csv("http://md.psych.bio.uni-goettingen.de/mv/data/div/politeness_data.csv")
dd$subject <- factor(dd$subject)
dd$gender <- factor(dd$gender)
dd$attitude <- factor(dd$attitude)

# fit model
politeness.model = lmer(frequency ~ attitude + gender + (1|subject) + (1|scenario), data=dd)
# look at the model
politeness.model
```

Geschlecht bindet einen Teil der Variation, die vorher im subject-Faktor enthalten war.

Signifikanztests via Modellvergleich.
Wegen der Diskussion um p-Werte in Mixed-Models hier mit Maximum-Likelihood Schätzung (REML-Flag).

```{r eval=FALSE}
# dd = read.csv("http://www.bodowinter.com/tutorial/politeness_data.csv")
dd <- readr::read_csv("http://md.psych.bio.uni-goettingen.de/mv/data/div/politeness_data.csv")
dd$subject <- factor(dd$subject)
dd$gender <- factor(dd$gender)
dd$attitude <- factor(dd$attitude)

# fit model
politeness.null  = lmer(frequency ~            gender + (1|subject) + (1|scenario), data=dd, REML=FALSE)
politeness.model = lmer(frequency ~ attitude + gender + (1|subject) + (1|scenario), data=dd, REML=FALSE)
# perform likelihood ratio test of the difference
anova(politeness.null,politeness.model)
# random intercept model, so coefficients include subject and situation parameters
coef(politeness.model)
```


random slope model (interaction)

```{r eval=FALSE}
# dd = read.csv("http://www.bodowinter.com/tutorial/politeness_data.csv")
# dd = read.csv("http://www.bodowinter.com/uploads/1/2/9/3/129362560/politeness_data.csv")
dd <- readr::read_csv("http://md.psych.bio.uni-goettingen.de/mv/data/div/politeness_data.csv")
dd$subject <- factor(dd$subject)
dd$gender <- factor(dd$gender)
dd$attitude <- factor(dd$attitude)

# fit model
politeness.model = lmer(frequency ~ attitude + gender + (1+attitude|subject) + (1+attitude|scenario), data=dd, REML=FALSE)
coef(politeness.model)
# test the effect of attitude, interaction included
politeness.null = lmer(frequency ~             gender + (1+attitude|subject) + (1+attitude|scenario), data=dd, REML=FALSE)
anova(politeness.null,politeness.model)
```


Note that the null model needs to have the same random effects structure. 
So, if your full model is a random slope model, your null model also needs to be a random slope model.

influencial cases
dfbeta() funktioniert nicht bei mixed models
library(influence.ME) ansehen

```{r eval=FALSE}
require(lmerTest)
# politeness = read.csv("http://www.bodowinter.com/tutorial/politeness_data.csv")
dd <- readr::read_csv("http://md.psych.bio.uni-goettingen.de/mv/data/div/politeness_data.csv")
dd$subject <- factor(dd$subject)
dd$gender <- factor(dd$gender)
dd$attitude <- factor(dd$attitude)

# self written 
#{} check somenumber
somenumber <- 5
all.res=numeric(nrow(politeness))
for(i in 1:nrow(politeness)){
  #myfullmodel=lmer(frequency~attitude+(1+attitude|subject), data=politeness,POP[-i,])
  myfullmodel=lmer(frequency~attitude+(1+attitude|subject), data=dd,REML=FALSE)
  all.res[i]=fixef(myfullmodel)[somenumber]
  }

```




## aov und formulae in reinen Messwiederholungsdesigns


[Quelle](http://ww2.coastal.edu/kingw/statistics/R-tutorials/repeated.html)

aov()

Two factor design with repeated measures on both factors:

DV ~ IV1 * IV2 + Error(subject/(IV1*IV2))

Three factor mixed design with repeated measures on IV2 and IV3:

DV ~ IV1 * IV2 * IV3 + Error(subject/(IV2*IV3))





# Übungen / Exercises

cf thinking in [lm_cat](lm_cat.html#thinking)


## Priming Beispiel

Beispiel Prime
[link](../example/prime/prime.html)

<!-- todo: share example -->

<!-- todo: Wollschläger example -->


# Links und Referenzen

http://dwoll.de/rexrepos/posts/anovaMixed.html

Wollschläger
http://www.uni-kiel.de/psychologie/rexrepos/posts/anovaRBFpq.html

Bodo Winter: Mixed Models
Tutorial and data seem to have moved. Some links:
[http://www.bodowinter.com/tutorial/bw_LME_tutorial.pdf]
[http://www.bodowinter.com/uploads/1/2/9/3/129362560/bw_lme_tutorial1.pdf]
Data
`dd = read.csv("http://www.bodowinter.com/tutorial/politeness_data.csv")`
`dd = read.csv("http://www.bodowinter.com/uploads/1/2/9/3/129362560/politeness_data.csv")`


Hao Zhang
[http://www.stat.purdue.edu/~zhanghao/MAS/handout/anova%20with%20mixed%20effects.pdf]

Motor Control Lab
http://gribblelab.wordpress.com/2009/03/09/repeated-measures-anova-using-r/

Missing values
[Mixed-Models-for-Repeated-Measures](http://www.uvm.edu/~dhowell/StatPages/More_Stuff/Mixed-Models-Repeated/Mixed-Models-for-Repeated-Measures1.html)

Ein Video zum Unterschied fixed effects models vs. random effects models
[here](https://www.youtube.com/watch?v=ySTb5Nrhc8g)

Zur 'Correlation of Fixed Effects' in lmer
[Erklärungsansatz](http://stats.stackexchange.com/questions/57240/how-do-i-interpret-the-correlations-of-fixed-effects-in-my-glmer-output)
[Formel](https://stat.ethz.ch/pipermail/r-sig-mixed-models/2009q1/001926.html)

Zur Diskussion p-Values und Vorgehen
(statexchange Beitrag)[http://stats.stackexchange.com/questions/63464/is-this-an-acceptable-way-to-analyse-mixed-effect-models-with-lme4-in-r]
{}
[Field (2012, p. 245)][field] **Kapitel 7: Regression**

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



# gar


```



ggplot(dd.r, aes(x=sit, y=anx)) +
  # define colors
  scale_fill_manual(values=c("#FFFFFF", "#CCCCCC", '#AAAAAA')) +
  # create white background etc
  #theme_bw() +
  geom_point() +
  # add verbose axis labels
  labs(x = "Anxiety-Group", y = "STAI-Trait")


ggplot(dd.r, aes(x=sit, y=anx, fill=sit)) +
  # define colors
  scale_fill_manual(values=c("#FFFFFF", "#CCCCCC", '#AAAAAA')) +
  # create white background etc
  #theme_bw() +
  stat_summary(fun.y=mean, geom="bar", position=position_dodge(width=0.7), colour='black') +
  # add error information
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", position=position_dodge(width=0.7), width=0.2) +
  # put big points on top of bars
  stat_summary(fun.data = mean_cl_normal, geom = "point", shape=21, size=5, position=position_dodge(width=0.7)) +
  # add verbose axis labels
  labs(x = "Anxiety-Group", y = "STAI-Trait")
 ```
 
 
## Multi-Level-Ansatz: Beispiel STAI-State verschiedene Gruppen in verschiedenen Situationen

Ein dreistufiger Between-Subjects-Faktor (Gruppe).
Ein dreiftufiger Within-Faktor (Situation).

```{r}
dd <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_wide.txt")
require("reshape2")
require(psych)
# require(lme4) # is loaded and used by lmerTest
require(lmerTest)
head(dd)
psych::describe(dd)

# convert it to long format 
dd.r <- reshape2::melt(dd, id=c('subj', 'gender', 'group'), variable.name='sit', value.name='anx', measured=c('neutr', 'exam1', 'exam2'))

# random intercept model (global intercept for subjects, not specific for fixed factor level)
m.1 <- lmer(anx ~ group + sit + (1|subj), data=dd.r)
m.1
summary(m.1)
anova(m.1)

# descriptives - the Rcmdr way
tapply(dd.r$anx, list(situation=dd.r$sit), mean, na.rm=TRUE) # means
tapply(dd.r$anx, list(situation=dd.r$sit), function(x) sum(!is.na(x))) # counts
# the psych package way
by(dd.r$anx, list(dd.r$sit), psych::describe)

# visualization
require(ggplot2)
ggplot(dd, x=group, y=neutr, aes(group, 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)  

ggplot(dd.r, aes(x=sit, y=anx)) +
  # define colors
  scale_fill_manual(values=c("#FFFFFF", "#CCCCCC", '#AAAAAA')) +
  # create white background etc
  #theme_bw() +
  geom_point() +
  # add verbose axis labels
  labs(x = "Anxiety-Group", y = "STAI-Trait")


```

# Check

Check, whether you have understood the following statements

- In between factors we have different observations in every sub-category.
- In within factors we have the same observations over the sub-cateogories of the factor. This implies dependencies. 
  Several measures of the same observation are usually more similar than several measures of different observations.
- In Anova designs the terms between-subjects-design and within-subjects-design can be found. Subjects are often our observations, but not always.
  Therefore we prefer the term within or between factor.
- Mixed models mix between factors and within factors. 
- When we have continuous explanatory variables in our design, this is often called ANCOVA (Analysis of Covariance)
- We recommend using multi level techniques (packages) to analyse mixed designs usually using the `package(nlme)` 
  but you might also consider using `package(lme4)`, which uses a different syntax
- When we apply MLM to repeated measures designs we consider the observations to be a level above level 1.
  Measuring them repeatedly implies a hierarchy in our data.
  We allow random effects for our observations, f. e. random intercepts due to individual characteristics of them.
- Random effects and their coefficients apply only to the adequate higher level subgroup whereas fixed effects are fix for all observations.
  So random coefficients are specific additional effects (corrections) over fixed effects in members of sub-levels.
- It's always a good idea to take a look at the coefficients of the adapted MLM to understand, what is happening.
  This also helps for visualizations of the effects.
- All the aspects of categorial explanatory variables apply also (dummy-coding, contrasts, unbalanced designs, type of SS, etc)
- Prediction includes random effects also. 
  If we do not have the same observations again, this makes no sense, as the characteristics of a specific observation are used. 
  But we can limit the prediction to levels used. `level=0` uses fixed effects only.

- Example
  
```{r}
# convert it to long format using tidyr
library(tidyverse)
dd.w <- read_tsv("http://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_wide.txt")

dd <- dd.w %>% tidyr::gather(key=sit, value=anx, neutr, exam1, exam2)
dd$subj.f <- factor(dd$subj)
dd$group.f <- factor(dd$group, levels=c("av", "lo", "hi"))      # tendency to show anxiety, actual anxiety in variable dd$anx
dd$sit.f <- factor(dd$sit, levels=c("neutr", "exam1", "exam2")) # situation 
head(dd)

library(nlme)
m.00  <- gls(anx ~ 1,                   data=dd, method="ML")
m.0   <- lme(anx ~ 1, random = ~1|subj, data=dd, method="ML")
m.g   <- lme(anx ~ group.f, random = ~1|subj, data=dd, method="ML")
m.s   <- lme(anx ~ sit.f, random = ~1|subj, data=dd, method="ML")
m.g.s <- lme(anx ~ group.f + sit.f, random = ~1|subj, data=dd, method="ML")
m.gs  <- lme(anx ~ group.f * sit.f, random = ~1|subj, data=dd, method="ML")

summary(m.00)
ggplot(data=dd) + geom_point(aes(x=sit.f, y=anx, color=group.f), size=2) + 
    geom_point(aes(x=sit.f, y=predict(m.00)), color="grey") +
    geom_point(aes(x=sit.f, y=m.00$coefficients[[1]]), size = 4, color = "grey") +
    geom_line(aes(x=sit.f, y=m.00$coefficients[[1]]), color = "grey", linetype = "dashed") +
    facet_grid(.~ group.f)  


summary(m.0)
fixef(m.0)
head(ranef(m.0))
ggplot(data=dd) + geom_point(aes(x=sit.f, y=anx, color=group.f), size=2) + 
    geom_point(aes(x=sit.f, y=predict(m.0)), color="grey") +
    geom_point(aes(x=sit.f, y=m.0$coefficients[[1]]), size = 4, color = "grey") +
    geom_line(aes(x=sit.f, y=m.0$coefficients[[1]]), color = "grey", linetype = "dashed") +
    facet_grid(.~ group.f)  

fixef(m.g)
head(ranef(m.g))
summary(m.g)

fixef(m.s)
head(ranef(m.s))
(pl <- ggplot(data=dd) + geom_point(aes(x=sit.f, y=anx, color=group.f), size=2) + 
    geom_point(aes(x=sit.f, y=predict(m.s)), color="grey") +
    stat_summary(fun = mean, aes(x=sit.f, y=predict(m.s)), geom="point", size = 4, color = "grey") +
    stat_summary(fun = mean, aes(x=sit.f, y=predict(m.s), group=subj), geom="line",color = "grey", linetype = "dashed") +
    facet_grid(.~ group.f)  
  )
summary(m.s)

fixef(m.gs)
head(ranef(m.gs))
ggplot(data=dd) + geom_point(aes(x=sit.f, y=anx, group=group.f, color=group.f), size=2) + 
    geom_point(aes(x=sit.f, y=predict(m.gs))) +
    stat_summary(fun = mean, aes(x=sit.f, y=predict(m.gs)), geom="point", size = 4, color = "grey") +
    # stat_summary(fun = mean, aes(x=sit.f, y=predict(m.gs)), geom="line", color = "grey", linetype = "dashed") +
     stat_summary(fun = mean, aes(x=sit.f, y=predict(m.gs), group=subj), geom="line",color = "grey", linetype = "dashed") +
    facet_grid(.~ group.f)
summary(m.gs)

# our prediction is a combination of fixed and random effects
head(model.matrix(m.gs, data=dd))
head(unique(model.matrix(m.gs, data=dd)))
# as.numeric(rownames(head(unique(model.matrix(m.gs, data=dd)))))

head(m.gs$fitted)
m.gs$fitted[as.numeric(rownames(head(unique(model.matrix(m.gs, data=dd))))),]

head(predict(m.gs))
dd.n <- tibble(group.f = factor(c("av","av","av","lo","lo","lo","hi","hi","hi")), 
               sit.f   = factor(c("neutr","exam1", "exam2","neutr","exam1", "exam2","neutr","exam1", "exam2")),
               subj    = c(180,179,177,181,182,183,184,185,186)          
               )
# subj is 1:180, predict needs all variables in newdata that occur in the model behind, random variables also
predict(m.gs, dd.n)
# prediction only on base of fixed effects or on base of a specific max. level can be done using the level parameter
predict(m.gs, dd.n, level=0)  # here only the fixed effects are used to do prediction

```


### Screencast 


- Mixed Anova Models - Check
  [StudIP](https://studip.uni-goettingen.de/plugins.php/mediacastplugin/media/check/557da707f35cbedf27871677d1e8e982/e9e10aeeb74d3726e6818694beb7db82/126?v=lm_cat_mixed_mlm-de-2020-06-.mp4) - 
[ownCloud](https://owncloud.gwdg.de/index.php/s/seCMpuLzkHIs8rD)


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