---
title: Lineares Modell, kategoriale Prädiktoren, reine Messwiederholungsmodelle (within
  subjects designs)
author: "Peter Zezula"
subtitle: "auch als Beispiel für Multi-Level-Modelle"
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_repeated.Rmd)

Generell sind varianzanalytische Messwiederholungs- und Mixed-Designs eine Teilmenge der Multi-Level-Modelle (MLM).
Wegen der deutlich höheren Flexibilität und Mächtigkeit der MLM empfehlen wir generell deren Einsatz.

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.

Voraussetzungen: 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


Alternatives, inzwischen bevorzugtes Vorgehen

- [Multi-Level-Modelle](../mlm/mlm.html)



## Long / wide Format

Das Verständnis von long und wide als unterschiedliches Format für dieselben Daten wird vorausgesetzt.

[cf](../long_wide/long_wide.html)




# Beispiel STAI-State Normalgruppe in verschiedenen Situationen

Beschreibung des Datensatzes in den [Beispielen](/mv/example) unter STAI Exam.



## Ansatz mit ezANOVA


Ausgabe einer ANOVA Tabelle.
Mauchly's Test for Sphericity sowie die Greenhouse-Geisser und die Huyn-Feldt Korrektur werden mit ausgegeben.
Bei Messwiederholungsdesigns braucht ezANOVA das Long-Format.

```{r}
# get and prepare data
dd <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_wide.txt")
dd <- dd[dd$group == 'av',c('subj', 'neutr', 'exam1', 'exam2')]

# convert it to long format using library(reshape2)
# require(reshape2)
# dd.l <- melt(dd, id=c('subj'), variable.name="sit", value.name="anx", measured=c('neutr', 'exam1', 'exam2'))
# head(dd.l)
# rm(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 <- dd %>% tidyr::pivot_longer(names_to="sit", values_to="anx", cols=neutr:exam2)
dd.l$subj <- factor(dd.l$subj)
dd.l$sit <- factor(dd.l$sit)
head(dd.l)

#require(Rcmdr)
require(ez)

# fit model using ezANOVA
m.f <- ezANOVA(data=dd.l, dv=.(anx), wid = .(subj), within = .(sit), detailed=TRUE, type=3)
# having one variable only we can use simply the name of a variable (without quoting it)
m.f <- ezANOVA(data=dd.l, dv=anx, wid = subj, within = sit, detailed=TRUE, return_aov=T)

# parameter
m.f
summary(m.f$aov)
#ezANOVA(data=dd.l, dv=.(anx), wid = .(subj), within = .(sit), detailed=TRUE, type=3)

# ez offers ezPlot()
ezPlot(data=dd.l, dv=.(anx), wid = .(subj), within = .(sit), x=sit)

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

# create graphs
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') +
  stat_summary(fun.y=mean, geom="bar", colour='black') +
  # add error information
  #stat_summary(fun.data = mean_cl_normal, geom = "errorbar", position=position_dodge(width=0.7), width=0.2) +
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", 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)) +
  stat_summary(fun.data = mean_cl_normal, geom = "point", shape=21, size=5) +
  # add verbose axis labels
  labs(x = "Anxiety-Group", y = "STAI-Trait")

```



## ezANOVA mit detaillierteren Informationen, deskriptiven Daten und Grafiken


Über den Parameter `return_aov` bei `ezANOVA()` können wir ein aov-Objekt speichern, in dem viele Informationen mehr als in der Anova-Tabelle enthalten sind.
Beispielsweise finden sich dort die ganzen Quadratsummen, auch die verschiedenen Residual-SS.

Wir können über den Befehl `ezStats()` deskriptive Statistiken anfordern.

`ezPlot()` ist ein Wrapper für `library(ggplot2)` eignet sich, um einfach Grafiken varianzanalytischer Designs zu erstellen.

```{r}
# get and prepare data
dd <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_wide.txt")
dd <- dd[dd$group == 'av',c('subj', 'neutr', 'exam1', 'exam2')]

# convert it to long format 
#dd.l <- melt(dd, id=c('subj'), 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 <- dd %>% tidyr::pivot_longer(names_to="sit", values_to="anx", cols=neutr:exam2)
dd.l$subj <- factor(dd.l$subj)
dd.l$sit <- factor(dd.l$sit)
head(dd.l)


#require(Rcmdr)
require(ez)

# by default, ezANOVA outputs a anova table
ezANOVA(data=dd.l, dv=.(anx), wid = .(subj), within = .(sit), detailed=TRUE, type=3)

# ezANOVA can return an aov() object with flag return_aov = T
m.ez <- ezANOVA(data=dd.l, dv=.(anx), wid = .(subj), within = .(sit), detailed=TRUE, type=3, return_aov = TRUE)
m.ez
# summary shows, what parts are available
summary(m.ez)
# we can get a aov summary
summary(m.ez$aov)

# we can access special parts by number or name
# get anova table
m.ez[1]
# get aov object
m.ez['aov']
# ... or using m.ez$aov

# now we can take a look at the model matrix contrasts
head(unique(model.matrix(m.ez$aov)))

# ez also offers ezStats(). We can get the means and descriptives
ezStats(data=dd.l, dv=.(anx), wid = .(subj), within = .(sit), type=3)

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

```





## Ansatz mit `lme4::lmer()`


Mit `lmer()` aus der `library(lme4)` wird ein Multi-Level-Modell angepasst.
`package(lmerTest)` baut auf `package(lme4)` auf, fügt aber Signifikanztests hinzu.

Näheres zu Multi-Level-Modellen findet sich in der Unit
[ml](../ml/ml.html).


Alternativen sind `aov()`, `car::Anova()` über spezielle Syntax zur Angabe der entsprechenden Error-Terme (für den sog. Versuchspersonen-Faktor).


```{r}
# get data
#require(Rcmdr)
require(ggplot2)
require(reshape2)
# 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")
dd <- dd[dd$group == 'av',c('subj', 'neutr', 'exam1', 'exam2')]
head(dd)

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

#psych::describe(dd)
require(Rmisc)
Rmisc::summarySE(data=dd.r, "anx", groupvars="sit")

# fit model
# random intercept model, no interaction
m.f <- lmerTest::lmer(anx ~ sit + (1|subj), data=dd.r)
# parameter
m.f
summary(m.f)
# test the model
anova(m.f)

unique(model.matrix(m.f))

# # 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), sd, na.rm=TRUE) 
#   # std. deviations
# 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)

plot(m.f)

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), width=0.2) +
	  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")

ggplot(dd.r, aes(x=sit, y=anx, fill=sit)) +
  # define colors
  scale_fill_manual(values=c("#FFFFFF", "#CCCCCC", '#AAAAAA')) +
  stat_summary(fun.y=mean, geom="bar", colour='black') +
  # add error information
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", width=0.2) +
  # put big 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-Trait")


```



## Ansatz mit `nlme::lme()`


Alternativ zum `package(lme4)` können wir auch mit dem `package(nlme)` ein Mehrebenenmodell anpassen.
Mit `lme()` aus der `library(nlme)` wird ein Multi-Level-Modell angepasst.

Näheres zu Multi-Level-Modellen findet sich in der Unit
[ml](../ml/ml.html).



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

# convert it to long format using tidyr
require(tidyr)
#dd.r <- dd %>% tidyr::gather(key=sit, value=anx, neutr, exam1, exam2)
dd.r <- dd %>% tidyr::pivot_longer(names_to="sit", values_to="anx", cols=neutr:exam2)

dd.r$subj <- factor(dd.r$subj)
dd.r$sit <- factor(dd.r$sit)
head(dd.r)

#psych::describe(dd)
require(Rmisc)
Rmisc::summarySE(data=dd.r, "anx", groupvars="sit")

# fit model
require(nlme)
# random intercept model, no interaction
m.f <- nlme::lme(fixed=anx ~ sit, random = ~1 |subj, data=dd.r)
# parameter
m.f
summary(m.f)
# test the model
anova(m.f)

#unique(model.matrix(m.f))

# # 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), sd, na.rm=TRUE) 
#   # std. deviations
# 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)

plot(m.f)

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), width=0.2) +
	  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")

ggplot(dd.r, aes(x=sit, y=anx, fill=sit)) +
  # define colors
  scale_fill_manual(values=c("#FFFFFF", "#CCCCCC", '#AAAAAA')) +
  stat_summary(fun.y=mean, geom="bar", colour='black') +
  # add error information
  stat_summary(fun.data = mean_cl_normal, geom = "errorbar", width=0.2) +
  # put big 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-Trait")


```



# Ausflug: Vergleich des between und within Ansatzes

Hier soll belegt werden, dass sich Effekte früher zeigen, wenn ein Within-Design möglich ist.
In den Beispielen werden dieselben Daten verwendet, es wird aber einmal so getan, 
als seien die Personen mehrfach unter verschiedenen Bedingungen gemessen worden (within), 
im anderen Fall so, als seien in den verschiedenen Bedingungen unterschiedliche Personen gemessen worden (between).


## Vergleich between, within, Ansatz mit ezANOVA


Vergleich derselben Daten, einmal behandelt als Zwischengruppenvergleich (between subjects design) und einmal als Innergruppenvergleich (within subjects design).

Jede Beobachtung bekommt einen eigenen Intercept (Versuchspersoneneneffekt).
Die Residualvarianz nimmt hierdurch ab, weil der individuelle Intercept Varianz bindet, die die Gruppenunterschiede nicht binden.
Das Verhältnis der Quadratsumme (Varianz) der Gruppe im Vergleich zur Rest-Quadratsumme (residuals) verändert sich.
Der F-Wert ist im within design wegen der Reduktion der Residualvarianz größer, was zu höherer Power führt.

Die Aufteilung der Residualvarianz ist zu sehen im Summary des Modells unter 'Random effects:'. 
Diese Anteile summieren sich zur Residualvarianz im Summary des between subjects model.
Die verschiedenen resultierenden F-Werte mit zugehörigen p-Werten sind ebenfalls zu finden.

Die Interpretation der Koeffizienten (Gewichte) und ihre Höhe ist in beiden Modellen gleich.

Achtung: Hier wird mit SS Type II gerechnet, nicht mit SS Type I wie unten (Vergleich via `aov()`. 
Bei balancierten Designs sind die Ergebnisse gleich SS Type I.

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

require(tidyr)
# dd.l <- dd %>% tidyr::gather(key=sit, value=anx, neutr, exam1, exam2)
dd.l <- dd %>% tidyr::pivot_longer(names_to="sit", values_to="anx", cols=exam1:exam2)
dd.l$subj <- factor(dd.l$subj)
dd.l$sit <- factor(dd.l$sit)
dd.l$obs <- factor(1:length(dd.l$subj))
head(dd.l)


# treat it as between subject model, subjects are 'multiplied'
ezANOVA(data=dd.l, dv=.(anx), wid = .(obs), between = .(sit), detailed=TRUE, type=2)

# repeated measures design
ezANOVA(data=dd.l, dv=.(anx), wid = .(subj), within = .(sit), detailed=TRUE, type=2)

```

Beim Vergleich der beiden Summaries sehen wir sehr deutlich die Veränderung der Quadratsummen (SSn, SSd) (Zähler- und Nenner-Quadratsumme).
Die SSn ist in beiden Modellen gleich (1682.54)
Im Between-Subjects-Modell ist die SSd, die auch die Unterschiede zwischen den Vpn enthält, viel größer (28279.65) 
als im Within-Subjects-Modell (551.45).
Daher wird auch der F-Wert viel größer bzw. der p-Wert viel kleiner. 



```{r}
# get and prepare data
dd <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/virt/v_stait_exam_wide.txt")
dd <- dd[dd$group == 'av',c('subj', 'neutr', 'exam1', 'exam2')]
# convert it to long format 
dd.l <- melt(dd, id=c('subj'), variable.name="sit", value.name="anx", measured=c('neutr', 'exam1', 'exam2'))
dd.l$obs <- 1:length(dd.l$subj)

# between subjects model
m.b <- ezANOVA(data=dd.l, dv=.(anx), wid = .(obs), between = .(sit), detailed=TRUE, type=2, return_aov = TRUE)
# within subjects model
m.r <- ezANOVA(data=dd.l, dv=.(anx), wid = .(subj), within = .(sit), detailed=TRUE, type=2, return_aov = TRUE)

# we compare the effects
m.b$ANOVA
m.r$ANOVA
# ... denominator SS ist way smaller in within design
# ... p value is way smaller
# ... nominator SS is equal

# we can get the model.matrix to check dummy contrasts
unique(model.matrix(m.b['aov']$aov))
# we can also get model.matrix for the repeated subjects design. 
# As it has subj-1 dummy columns to code effect subject, it is big
# columns are:
model.matrix(m.r['aov']$aov)[1,]
# we look at the relevant parts of the model.matrix to get an impression of subjects dummy coding
model.matrix(m.r['aov']$aov)[c(1:10, 175:180),c(1:5, 61:62)]

# the full model.matrix would be (not printed, because it is huge)
# unique(model.matrix(m.r['aov']$aov))

```

Wir können jetzt auch sehen, dass die Residualvarianz stark zurückgeht, da der 'Versuchspersonenfaktor' viel Varianz bindet.

Außerdem hilft ein Blick auf die `model.matrix()` zu verstehen, dass der Versuchspersonenfaktor tatsächlich in n-1 Dummyvariablen zerlegt wird. 
Jede Beobachtung/Vp/subj erhält also eine eigene Dummy-Variable als Faktorstufe.
Die Spalten laufen von 'Error(subj/(sit))2' ... 'Error(subj/(sit))60'.

Schön ist in der `model.matrix()` zu sehen, dass jede Versuchsperson für den 3-stufigen Faktor 'sit' je drei Mal eine 1 in ihrer Spalte hat.





## Vergleich between subjects mit within subjects Design mit `aov()`



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

#require(Rcmdr)
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.l <- melt(dd, id=c('subj'), variable.name="sit", value.name="anx", measured=c('neutr', 'exam1', 'exam2'))
#head(dd.r)
# fit between subject model
m.b <- aov(anx ~ sit, data=dd.l)
# random intercept model, no interaction
#m.f <- lmer(anx ~ sit + (subj/sit), data=dd.r)
# random intercept
m.w <- aov(anx ~ sit + Error(subj), data=dd.l)
# random intercept, random slope
m.f <- aov(anx ~ sit + Error(subj/sit), data=dd.l)
# constant in between subjects design is the same for all observations
coefficients(m.b)
coefficients(m.w)
coefficients(m.f)
# constant in within subjects design varies and depends on subject
head(coefficients(m.f)$subj)

summary(m.b)
summary(m.w)
summary(m.f)

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

# `r summary(m.b)[[1]][2,2]`
```


Im Between-Subjects-Modell `m.b` gibt es nur die SS des Effekts 'sit' 1683 und die Gesamt-Residuen 28280, in Summe 29963.

Im Within-Subjects-Modell `m.w` splitten sich die SS in die within 3105 und between 24783 und interaction 1106 und die 'sit' 968, 
in Summe die Gesamt-Residuen `3105 + 24783 + 1106 + 968 = 29962` . 

Die Residuen, gegen die die Effekte getestet werden, verkleinern sich also von 28280 (between) auf `3105 + 24783 = 27888`.






# Beispiele Übungen

## Beispiel Interesse an Schule und Arbeit in Abhängigkeit von der Altersklasse, John Quick.


two-way within subjects anova: 

[Quelle](http://www.r-bloggers.com/r-tutorial-series-two-way-repeated-measures-anova/)

[Daten](http://dl.dropbox.com/u/10246536/Web/RTutorialSeries/dataset_anova_twoWay_repeatedMeasures.csv) 
oder als
Lokale tab-delimited [Kopie](http://md.psych.bio.uni-goettingen.de/data/div/john_quick_dataset_anova_twoWay_repeatedMeasures.txt)

30 Teilnehmer, deren Interesse an Schule bzw. Arbeit in 3 verschiedenen Altersstufen gemessen wurden, mit 10, 15 und 20 Jahren.
Interesse wurde erfasst mit einer Ratingskala von 1 - 5. Je höher der Score, desto höher das Interesse.



## Interesse an Schule und Arbeit: Ansatz mit ezANOVA

```{r}
# get data
# This dataset contains a hypothetical sample of 30 participants whose interest in school and interest in work was measured at three different ages (10, 15, and 20). 
# The interest values are represented on a scale that ranges from 1 to 5 and indicate how interested each participant was in a given topic at each given age.
dd <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/div/john_quick_dataset_anova_twoWay_repeatedMeasures.txt")
# first step create age
require(reshape2)
dd.r <- melt(dd, id=c('subject'), variable.name="school_age", value.name="interest", measured=c('schoolAge10', 'schoolAge15', 'schoolAge20', 'workAge10', 'workAge15', 'workAge20'))
#require(Rcmdr)
require(ez)

require(psych)
# require(lme4) # is loaded and used by lmerTest
require(lmerTest)
# add a topic factor
##dd$topic <- c(rep('school', length(dd$schoolAge10) * 3), rep('school', length(dd$workAge10) * 3))
# convert it to long format
# add a topic factor
#dd.r$topic <- factor(c(rep('school', length(dd$schoolAge10) * 3), rep('work', length(dd$workAge10) * 3)), labels=c("school", "work"))
dd.r$topic <- factor(rep(c(1, 2), c(90, 90)), labels=c("school", "work"))
# add a age factor
dd.r$age <- factor(rep(c(1, 2, 3), c(30, 30, 30)), labels=c("a10", "a15", "a20"))

# fit model using ezANOVA
#? random intercept model, no interaction
m.f <- ezANOVA(data=dd.r, dv=.(interest), wid = .(subject), within = .(topic, age), detailed=TRUE, type=3)
# parameter
m.f

# ezANOVA can return an aov() object with flag return_aov = FALSE
m.f.aov <- ezANOVA(data=dd.r, dv=.(interest), wid = .(subject), within = .(topic, age), detailed=TRUE, type=3, return_aov = TRUE)
m.f.aov

# ez offers ezPlot() 
ezPlot(data=dd.r, dv=.(interest), wid = .(subject), within = .(topic, age), x=.(topic), split=.(age))
ezPlot(data=dd.r, dv=.(interest), wid = .(subject), within = .(topic, age), x=.(age),   split=.(topic))

# ez also offers ezStats()
ezStats(data=dd.r, dv=.(interest), wid = .(subject), within = .(topic, age), type=3)
# order matters for output
ezStats(data=dd.r, dv=.(interest), wid = .(subject), within = .(age, topic), type=3)
# ?ezStats tells us
# Fisher's Least Significant Difference is computed as sqrt(2)*qt(.975,DFd)*sqrt(MSd/N), where N is taken as the mean N per group in cases of unbalanced designs.

```




## Ansatz wie auf r-blogger Website beschrieben

[Quelle](https://dl.dropboxusercontent.com/u/10246536/Web/RTutorialSeries/example_anova_twoWay_repeatedMeasures.txt)

Der Ansatz nutzt `Anova()` der `library(car)`.

Messwiederholungsdesigns werden hier mit dem Multivariaten Ansatz gerechnet.

```{r}
#//Two-Way Repeated Measures ANOVA Example
#//Created by John M. Quick
#//http://www.johnmquick.com
#//January 22, 2011

#read the dataset into an R variable using the read.csv(file) function
#dataTwoWayRepeatedMeasures <- read.csv("dataset_ANOVA_TwoWayRepeatedMeasures.csv")
dataTwoWayRepeatedMeasures <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/div/john_quick_dataset_anova_twoWay_repeatedMeasures.txt")
#display the data
#notice the atypical column arrangement for repeated measures data
head(dataTwoWayRepeatedMeasures)

#read the idata frame into an R variable
#idataTwoWayRepeatedMeasures <- read.csv("idata_ANOVA_TwoWayRepeatedMeasures.csv")
idataTwoWayRepeatedMeasures <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/div/john_quick_dataset_anova_twoWay_repeatedMeasures_idata.txt")
idataTwoWayRepeatedMeasures$Interest <- factor(idataTwoWayRepeatedMeasures$Interest)
idataTwoWayRepeatedMeasures$Age <- factor(idataTwoWayRepeatedMeasures$Age)
#display the idata frame
#notice the text values and correspondence between our idata rows and the columns in our dataset
idataTwoWayRepeatedMeasures

#create the linear model

#step 1: bind the columns
#use cbind() to bind the columns in the original dataset
interestBind <- cbind(dataTwoWayRepeatedMeasures$schoolAge10, dataTwoWayRepeatedMeasures$schoolAge15, dataTwoWayRepeatedMeasures$schoolAge20, dataTwoWayRepeatedMeasures$workAge10, dataTwoWayRepeatedMeasures$workAge15, dataTwoWayRepeatedMeasures$workAge20)

#step 2: define the model
#use lm() to generate a linear model using the bound columns from step 4
interestModel <- lm(interestBind ~ 1)

#load the car package (install first, if necessary)
library(car)
#compose the Anova(mod, idata, idesign) function
analysis <- car::Anova(interestModel, idata = idataTwoWayRepeatedMeasures, idesign = ~Interest * Age)

#use summary(object) to visualize the results of the repeated measures ANOVA
summary(analysis)

```



## Two-way Within Subjects ANOVA personality-project


Quelle: [http://personality-project.org/r/r.anova.html]

Data: [http://personality-project.org/r/datasets/R.appendix4.data] oder als lokale [Kopie](http://md.psych.bio.uni-goettingen.de/data/div/pers_proj_4.txt)

Fünf Versuchspersonen wurden gebeten, eine Liste von Wörtern zu lernen.
Es gab drei Arten von Wörtern auf der Liste.
Es gibt also einen dreifach gestuften Messwiederholungsfaktor Wortart 'valence' mit den Stufen 'neg', 'neu' und 'pos'.
Die Abrufart 'task' ist zweifach gestuft ('free', cued') und ebenfalls ein Messwiederholungsfaktor.

Dies hier dient nur der Veranschaulichung der Model-Spezifikation. 
Ein 2*3 Design mit insgesamt 5 Beobachtungen macht wenig Sinn.

```{r}
# get data
#require(Rcmdr)
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/div/pers_proj_4.txt")
# data are in long format, no conversion needed
# fit between subject model
m.b <- lm(recall ~ task * valence, data=dd)
# random intercept model, no interaction
m.u1 <- lmer(recall ~ task + valence + (1|subj), data=dd)
print(m.u1)
anova(m.u1)

# random intercept, interaction
m.u2 <- lmer(recall ~ task + valence + task:valence + (1|subj), data=dd)
print(m.u2)
anova(m.u2)

# random slope and intercept, no interaction
m.u3 <- lmer(recall ~ task + valence + (valence|subj) + (task|subj), data=dd)
print(m.u3)
anova(m.u3)


# random slope and intercept, interaction
m.u4<- lmer(recall ~ task + valence + task:valence + (valence|subj) + (task|subj), data=dd)
print(m.u4)
anova(m.u4)

# unclear what is modelled here
m.u5<- lmer(recall ~ task * valence + (1|subj) + (1|valence:subj) + (1|task:subj), data=dd)
print(m.u5)
anova(m.u5)

# model random slope and intercept for all fixed factors
# full model: random slope and intercept, interaction
m.u6<- lmer(recall ~ task * valence + (1 + valence|subj) + (1 + task|subj), data=dd)
print(m.u6)
anova(m.u6)


# compare non interaction model with full interaction model
anova(m.u1, m.u6)

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

Dasselbe mit `ezANOVA()`


```{r}
dd <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/div/pers_proj_4.txt")
require(ez)
ezANOVA(data=dd, dv=recall, wid=subj, within=.(task, valence), detailed=T, type=1 )
ezPlot(data=dd, dv=recall, x=task, wid=subj, within=.(task, valence), split=valence)
ezPlot(data=dd, dv=recall, x=valence, wid=subj, within=.(task, valence), split=task)

```



# Hintergrund: zu Formulae, auch in `nlme()`

siehe: [Mehrebenenmodelle](../multilevel/multilevel.html)



## Original-Auswertung und Einführung

Zur Anwendung kommt `aov()` und nicht `nlme()`.

Example 4. Two-Way Within-Subjects ANOVA

Five subjects are asked to memorize a list of words. 
The words on this list are of three types: positive words, negative words and neutral words. 
Their recall data by word type is displayed in Appendix III. 
Note that there is a single factor (Valence ) with three levels (negative, neutral and positive). 
In addition, there is also a random factor Subject . 
Create a data file ex3 that contains this data. 
Again it is important that each observation appears on an individual line! 
Note that this is not the standard way of thinking about data. 
Example 6 will show how to transform data from the standard data table into this form.

In this example, Subject is crossed with both Task and Valence, so you must specify three different error terms: one forTask , one for Valence and one for the interaction between the two. 

Dosage  Three levels of drug were administered to 18 subjects.



In this example, Subject is crossed with both Task and Valence, so you must specify three different error terms: one forTask , one for Valence and one for the interaction between the two. Fortunately, R is smart enough to divide up the within-subjects error term properly as long as you specify it in your command. The commands are:

```{r}
# the original 
datafilename="http://personality-project.org/r/datasets/R.appendix4.data"
data.ex4=read.table(datafilename,header=T)   #read the data into a table
 data.ex4                                      #show the data
aov.ex4=aov(Recall~(Task*Valence)+Error(Subject/(Task*Valence)),data.ex4 )

summary(aov.ex4)
print(model.tables(aov.ex4,"means"),digits=3)       #report the means and the number of subjects/cell
boxplot(Recall~Task*Valence,data=data.ex4) #graphical summary of means of the 6 cells

```












# Conway

## Benefits

weniger Vpn
erhöhte Power
- durch Vpn-'Faktor'
- Error Term wird nochmals aufgeteilt in Vpn-Varianz und Restvarianz

Beispiel AB/AC gepaartes assoziatives Lernen
- Vpn lernen Liste von Begriffspaaren A-B
- Danach lernen sie eine Liste von Begriffspaaren  A-C
- Dann ist das Erinnern von B nach Hinweisreiz A gestört
- Retroaktive Interferenz

Schönes Beispiel: Auswirkung auf F-Ratio, wenn innerhalb einer Beobachtung Werte vertauscht werden und somit die unsystematische Varianz erhöht wird.

## Risks

- Reihenfolge-Effekte, Messwiederholungs-Effekte
- Lösungsansatz Counterbalancing
- missing data

Counterbalancing
- blocked design: Zufallszuweisung zu einer der möglichen Ablaufreihenfolgen
  - schon bei 3 Faktorstufen 3! = 6, bei 4 Faktorstufen 4! = 24 Kombinationsmöglichkeiten
- latin squares
  - jede Bedingung taucht mindestens einmal an jeder möglichen Postition auf
  - 3 Faktorstufen: A1, A2, A3 - A2, A3, A1 - A3, A1, A2
  - Anzahl der Kombinationen = Anzahl der Faktorstufen

Missing data
- relative Anzahl der missings: Daumenregel: weniger als 10% in einer Variablen ist akzeptabel, mehr nicht
- Muster in missings T-Test der Personen mit Missings in einer Variable gegen die ohne Missings. Sig. -> Muster
- Umgang
  - alle Beobachtungen mit Missings in relevanter Variable aus Analyse nehmen (teuer)
  - Schätzen der Missings
    - Ersetzen durch Mittelwert: Konservativ, vermindert Varianz systematisch
    - Regressionsbasierte Schätzung: besser, aber welche Variablen zur Vorhersage benutzen?
  


Bildgebende Verfahren und deren Kosten
Vpn in Reaktionszeit-Experimenten: Es gibt Langsame und Schnelle, unabhängig von der jeweiligen Aufgabe (Faktorstufe)








## aov und formulae in reinen Messwiederholungsdesigns


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

aov()

```
# wide format
subject  		storeA	storeB	storeC	storeD
lettuce			1.17	1.78	1.29	1.29
potatoes		1.77	1.98	1.99	1.99		
milk			1.49	1.69	1.79	1.59
eggs			0.65	0.99	0.69	1.09
bread			1.58	1.70	1.89	1.89
cereal			3.13	3.15	2.99	3.09
ground.beef		2.09	1.88	2.09	2.49
tomato.soup		0.62	0.65	0.65	0.69
laundry.detergent	5.89	5.99	5.99	6.99
aspirin			4.46	4.84	4.99	5.15

# long format
groceries2                                               # take a look
   price  store           subject
1   1.17 storeA           lettuce
2   1.77 storeA          potatoes
3   1.49 storeA              milk
4   0.65 storeA              eggs
5   1.58 storeA             bread
6   3.13 storeA            cereal
7   2.09 storeA       ground.beef
8   0.62 storeA       tomato.soup
9   5.89 storeA laundry.detergent
10  4.46 storeA           aspirin
11  1.78 storeB           lettuce
12  1.98 storeB          potatoes
13  1.69 storeB              milk
14  0.99 storeB              eggs
15  1.70 storeB             bread
16  3.15 storeB            cereal
17  1.88 storeB       ground.beef
18  0.65 storeB       tomato.soup
19  5.99 storeB laundry.detergent
20  4.84 storeB           aspirin
21  1.29 storeC           lettuce
22  1.99 storeC          potatoes
23  1.79 storeC              milk
24  0.69 storeC              eggs
25  1.89 storeC             bread
26  2.99 storeC            cereal
27  2.09 storeC       ground.beef
28  0.65 storeC       tomato.soup
29  5.99 storeC laundry.detergent
30  4.99 storeC           aspirin
31  1.29 storeD           lettuce
32  1.99 storeD          potatoes
33  1.59 storeD              milk
34  1.09 storeD              eggs
35  1.89 storeD             bread
36  3.09 storeD            cereal
37  2.49 storeD       ground.beef
38  0.69 storeD       tomato.soup
39  6.99 storeD laundry.detergent
40  5.15 storeD           aspirin

aov.out = aov(price ~ store + Error(subject/store), data=groceries2)

```

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))




# Ausflug: Kontrast Koeffizienten in Messwiederholungsdesigns

Aus dem Buch Rutherford stammt das folgende Beispiel.
8 Personen
Faktor A (Lerntechnik, zweistufig) als Messwiederholungsfaktor (within subjects factor)
Faktor B (Lernzeit, dreistufig) als Messwiederholungsfaktor (within subjects factor)

Bemerkung: Nur eine Reihenfolge sinnvoll (memorize -> mnemonics). Dadurch unkontrollierte Reihenfolgeeffekte.

Datensatz long mit Effekt-Kodierung.
Die Dummy-Kodierung drückt Effekte aus als Abweichungen vom Gesamtmittelwert (grand mean) und nicht als Abweichungen von einer Referenzgruppe.
Die Koeffizienten müssen also anders interpretiert werden, die Anova-Tabelle bleibt gleich.

[data](http://md.psych.bio.uni-goettingen.de/data/div/rutherford_2x3_rep_long_effectcoded.txt)

Berechnung der Effekte unter Nutzung der Dummy-Variablen, die im Datensatz hinterlegt sind.

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

unique(dd[,3:length(names(dd))])
#mm <- aov(anx ~ sit + Error(subj/sit), data=dd.l)
mm <- aov(y ~ Lernbed + Error(vp/Lernbed), data=dd)
mm
summary(mm)

# todo: to be completed

```







# Beispiele, Übungen / examples and exercises


## Bushtucker

Beispiel Bushtucker von Andy Field.

Vergleich der Anwendung verschiedener Pakete bei einem einfaktoriellen Messwiederholungsmodell:
- `aov()` mit Errorterm
- `nlme::lme()`
- `lme4::lmer()`

[example_bushtucker.html](example_bushtucker.html)


## Ratten

Das Ratten-Beispiel von Everitt (2010)
[link](../example/rats/rats.html)


## Moore

Beispiel Moore
[link](../example/moore/moore.html)


todo: posthoc test
http://stats.stackexchange.com/questions/49108/post-hoc-test-after-2-factor-repeated-measures-anova-in-r



## Beispiel Hays


[Quelle](http://www.psych.upenn.edu/~baron/rpsych/rpsych.html#htoc60)


```{r}
require(psych)
require(effects)
# across subjects then conditions
data1<-c(
49,47,46,47,48,47,41,46,43,47,46,45,
48,46,47,45,49,44,44,45,42,45,45,40,
49,46,47,45,49,45,41,43,44,46,45,40,
45,43,44,45,48,46,40,45,40,45,47,40)
# create data frame and arm factors
matrix(data1, 
  ncol= 4, 
  dimnames = list(paste("subj", 1:12), c("Shape1.Color1", "Shape2.Color1", "Shape1.Color2", "Shape2.Color2")))
Hays.df <- data.frame(rt = data1, 
subj = factor(rep(paste("subj", 1:12, sep=""), 4)),
shape = factor(rep(rep(c("shape1", "shape2"), c(12, 12)), 2)),
color = factor(rep(c("color1", "color2"), c(24, 24))))

head(Hays.df)

m.r <- aov(rt ~ shape * color + Error(subj/(shape * color)), data=Hays.df)
m.r
summary(m.r)


# ez way
require(ez)
ezANOVA(data=Hays.df, dv=rt, wid=subj, within=.(shape, color), detailed=T, type=1 )
ezPlot(data=Hays.df, dv=rt, x=shape, wid=subj, within=.(shape, color), split=color)
#ezANOVA(data=Hays.df, dv=rt, wid=subj, within=.(shape, color), detailed=T, type=3 )
#ezANOVA(data=dd.l, dv=.(anx), wid = .(subj), within = .(sit), detailed=TRUE, type=1)
rm(Hays.df)
```



## Beispiel mit den Daten von Rutherford

[Beispiel](ex_rutherford.html)  mit den Daten von Rutherford (2002).








# Links und Referenzen

[Site](http://www.r-statistics.com/tag/ez/)
mit Links auf relevante Sites über Messwiederholungsdesigns mit R

Bates Vortrag mit sleepstudy data
[http://lme4.r-forge.r-project.org/slides/2011-03-16-Amsterdam/2Longitudinal.pdf]
{}
[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"
[everitt]: http://www.crcpress.com/product/isbn/9781439807699 "Brian S. Everitt (2010). Multivariable Modeling and Multivariate Analysis for the Behavioral Sciences"




# Check

## 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.

## zwei MLM-Packages, lme4 und lmer

## lmertest: Signifikanztests in lme4

## Wie kommen wir zu den Formeln?

Nehmen wir an, wir interessieren uns für die Erklärkraft der intervallskalierten Variable $year$ für das ebenfalls intervallskalierte Outcome-Maß  $gdp$ in einem hierarchischen Datensatz, gruppiert nach Variable $country$,

Dabei haben wir ein Modell mit random Intercept und random Slope im Kopf:

$$ \widehat{gdp}_{ij} = fixed_{00} + random_{0i} + (fixed_{1i} + random_{1i}) \cdot year_{ij} + \varepsilon_{ij} $$

oder

$$ \widehat{gdp}_{ij} = fixed\_intercept + fixed\_slope\_year \cdot year_{ij} + random\_intercept\_country_{i} + random\_slope\_country_{i} \cdot year_{ij} + \varepsilon_{ij} $$
i country
j year



## Example GDP

It ist basically the same example that you can find in our sheet_repeated_measure.

```{r}
require(tidyverse)
gdp_wide <- read_delim("https://md.psych.bio.uni-goettingen.de/mv/data/div/gdp.dat", "\t")
# gdp <- gather(gdp_wide, country, gdp, austria:usa)
gdp <- gdp_wide %>% tidyr::pivot_longer(names_to="country", values_to="gdp", cols=austria:usa)
gdp <- gdp %>% mutate(country_fac = factor(country)) %>% 
  mutate(year_fac = factor(year)) %>%
  mutate(year_nr = year - 1949)
require(nlme)
m1a <- lm(gdp ~ year, data = gdp)
# we better use interval scaled year_nr, so we get a slope and not a lot of dummy variable effects
m1 <- gls(gdp ~ year_nr, data = gdp, method = 'ML')
m2 <- lme(gdp ~ year_nr, data = gdp,
          random = ~1 | country_fac, 
          method = "ML")
m3 <- lme(gdp ~ year_nr, data  = gdp, 
          random = ~1 + year_nr | country_fac,
          method ="ML")
head(gdp[,c("year", "year_nr")])
summary(m1)
summary(m2)
summary(m3)

# we can see that models go better AIC and BIC go down
cat(AIC(m1), AIC(m2), AIC(m3))
cat(BIC(m1), BIC(m2), BIC(m3))

# also residuals are becoming smaller
cat(sd(m1$residuals), sd(m2$residuals), sd(m3$residuals))

# do random slopes matter?
anova(m2, m3) # ... yes

head(unique(model.matrix(m3, data=gdp)))

# we take a closer look at what the model offers
coef(m3)
fixef(m3)
ranef(m3)

gdp$pred_m1 <- predict(m1)
gdp$pred_m2 <- predict(m2)
gdp$pred_m3 <- predict(m3)

ncountries <- length(levels(gdp$country_fac))
nyears <- 34

# we calculate predicted value of first year of Austria manually using coefficients
fixef(m3)[1]  + ranef(m3)[[1]][1]  + (fixef(m3)[2] + ranef(m3)[[2]][1]) * gdp$year_nr[1]

gdp_plot_1 <- ggplot(data = gdp, aes(year, gdp,
                                     group = country_fac,
                                     #shape= country_fac,
                                     color = country_fac
                                     )) +
  # we scale x back to years
  scale_x_continuous(breaks = c(1, 11, 21, 31), label = c("1950", "1960", "1970", "1980")) +
  geom_line(legend=FALSE) +
  geom_point(legend=FALSE) +
  labs(x = "Year", y = "GDP per capita", color = "Country")

( gdp_plot_m1 <- gdp_plot_1 + geom_line(aes(y = pred_m1), size=2) )

( gdp_plot_m2 <- gdp_plot_1 + geom_line(aes(y = pred_m2)) )

( gdp_plot_m3 <- gdp_plot_1 + geom_line(aes(y = pred_m3)) )

( gdp_plot_m <- gdp_plot_1 + geom_line(aes(y = pred_m3))  + geom_line(aes(y = pred_m1), group=1, size=2, color="black", linetype="dashed") )


# details on variances in random effects
VarCorr(m3)

# importance of level variable, here country
#require(ICC)	# also returns confidence intervall
#ICCest(x=country_fac, y=gdp, data=gdp)

```




### Screencast 

A Screencast can be found in [lm_cat_mixed](./lm_cat_mixed.html). It explains not only pure repeated measure ANOVA models, but also.

- 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)

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

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

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