---
title: "Beispiel Bushtucker"
author: "Peter Zezula"
subtitle: "Beispiel für einfaktorielles Messwiederholungsmodell, Ansatz mit verschiedenen Paketen"
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](example_bushtucker.Rmd)


Aus Andy Fields, Exploring statistics using R (2012) chapter 13.4.7.2 

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


```{r bushtucker}
library(tidyverse)
dd.u <- read_tsv("http://md.psych.bio.uni-goettingen.de/mv/data/af/af-bushtucker-universal.txt")[1:8,1:4]
dd.u$subj <- 1:nrow(dd.u)
dd.l <- dd.u %>% gather(key = "animal", value="seconds", stick:witchetty)
dd.l$f.animal <- factor(dd.l$animal)
dd.l$f.subj   <- factor(dd.l$subj)
dd.l <- dd.l %>% group_by(f.subj) %>% mutate(subj.ms = mean(seconds, na.rm=T))
dd <- dd.l[c("f.subj", "f.animal", "seconds", "subj.ms")]

# independent group ANOVA
# lm(seconds ~ f.animal,                          data=dd) # would also show dummy effects in summary
ma.0 <- aov(   seconds ~ f.animal,                          data=dd)
# repeated measures ANOVA
ma.r <- aov(   seconds ~ f.animal + Error(f.subj/f.animal), data=dd)
summary(ma.0)
summary(ma.r)
# ... residuals go down and sig goes up

# using lme()
library(nlme)
mn.0 <- lme(seconds ~ 1,        random = ~1|f.subj, data = dd, method = "ML")   # random intercept of subjects already included  
mn.a <- lme(seconds ~ f.animal, random = ~1|f.subj, data = dd, method = "ML")   # 
anova(mn.a)
anova(mn.0, mn.a)

# using lme4() or better lmerTest
library(lmerTest)
m4.0 <- lmer(seconds ~            (1|f.subj), data=dd, REML=F)   # random intercept of subjects already included  
m4.1 <- lmer(seconds ~ f.animal + (1|f.subj), data=dd, REML=F)   #
anova(m4.1)
anova(m4.0, m4.1)

# lme() model proposed in Field (2012) is not adequate
mf.n <- lme(seconds ~ f.animal, random = ~1|f.subj/f.animal, data = dd, method = "ML")
summary(mf.n)
# this can be clearly seen by comparing
anova(mn.a, mf.n)
# ... the two models are nested by allowing more random effects, in fact too much

# confusion might be a result of doing it previously using aov() like shown above
mf.a <- aov(seconds ~ f.animal + Error(f.subj/f.animal), data = dd)
summary(mf.a)

```

Model `mf.n` ist so, wie von Andy Field in Kapitel 13.4.7.2 auf Buchseite 574 vorgeschlagen. Es macht aber keinen Sinn, da mit `subj/f.animal` die gesamte Varianz aufgeklärt ist. Die Residuals sinken auf 0, der Model-Vergleich wird nicht signifikant, weil es viel zu viele Parameter gibt. Die Informationskriterien AIC und BIC sinken nicht, sondern steigen.


# Extension "drink"

A second second level variable 

```{r}
# dd data object above has to be available
# we generate some data and add a drink variable with 4 levels
psych::describe(dd$subj.ms)
set.seed(17)
# cor(dd$subj.ms, dd$subj.ms + rnorm(nrow(dd)))

dd$x <- dd$subj.ms + rnorm(nrow(dd))
x <- sort(dd$x)
gar <- split(x, ceiling(seq_along(x)/(length(x)/4)))
# dd$drink <- split(x, ceiling(seq_along(x)/8))
dd$drink <- NA
dd$drink[dd$x %in% gar[[1]]] <- "a" 
dd$drink[dd$x %in% gar[[2]]] <- "b" 
dd$drink[dd$x %in% gar[[3]]] <- "c" 
dd$drink[dd$x %in% gar[[4]]] <- "d" 
dd$f.drink <- factor(dd$drink)

mn.ad <- nlme::lme(seconds ~ f.animal + f.drink, random = ~1|f.subj, data = dd, method = "ML") 
anova(mn.a, mn.ad)
summary(mn.ad)

# we do it using lmerTest
m4.ad <- lmer(seconds ~ f.animal + f.drink + (1|f.subj), data=dd, REML=F)   #
anova(m4.1, m4.ad)
summary(m4.ad)
```



# Literature

Pinheiro & Bates (2000)

Comparison of repeated measures designs using different packages [source](http://www.dwoll.de/rexrepos/posts/anovaMixed.html)

https://ethz.ch/content/dam/ethz/special-interest/math/statistics/sfs/Education/Advanced%20Studies%20in%20Applied%20Statistics/course-material/repeatedMeasures/skript_RM.pdf

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