---
title: "Multi-Level-Modelle"
subtitle: "auch als Alternative zu varianzanalytischen Messwiederholungsmodellen"
author: "Peter Zezula"
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

---

### [Rmd](mlm.Rmd)

`r if(!params$todo) {"\\begin{comment}"}`
todo:
guter profunder vergleichender Überblick mit lme4 und differenzierter grafischer Aufbereitung, hauptsächlich orthodont dataset
[linear mixed modeling](http://sia.webpopix.org/lme.html)

[this](https://cran.r-project.org/doc/contrib/Bliese_Multilevel.pdf)
and [this](http://www.stats.ox.ac.uk/~snijders/MLB_new_S.pdf)
and [this](http://jaredknowles.com/journal/2015/8/12/announcing-mertools)


ein Beispiel mit Ländern und länderspezifischen Strukturen, kulturellen Unterschieden, politischen Systemen, Korruption, finanzieller Ausstattung, ...

Messwiederholungsdesigns als Multi-Level-Modell mit nlme
http://www.jason-french.com/tutorials/repeatedmeasures.html

Multi-Level-Modelle berücksichtigen Abhängigkeiten in Daten, die durch hierarchische Datenstruktur entstehen.

Klassische Ansätze Aggregation und Disaggregation
todo: was ist das
Können mit bestimmten Problemen nicht adäquat umgehen.

Likelihood ratio test: Die Likelihoods werden logarithmiert (und mit -2 multipliziert (-2LL)).
Bei zwei Modellen, die hierarchisch sind und für die man je eine Likelihood hat, ist die Differenz der beiden LogLiklihood Chi^2 verteilt.
Dies nutzen wir beim Likelihood-Ratio-Test zum statistischen Vergleich zweier Modelle.

nested designs

Level 1 enthält die höchste Auflösungsstufe von Vorhersagevariablen.
Die Reaktionsvariable ist immer eine Level 1 Variable.

Bei individuellen Random-Coefficients können nur die fixed effects zur Vorhersage von Reaktionsvariablen bei unbekannten/neuen Beobachtungen benutzt werden.
Allerdings geben die random effects einen Eindruck der Größenordnung, die durch eventuelle andere Prädiktoren aufgeklärt werden könnte.


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





# Theorie: Multi-Level-Modelle als Ersatz für Repeated-Measures-ANOVA-Designs


## Einführende Präsentation von Thomas Schultze-Gerlach und Christan Treffenstädt

Die Präsentation als [rmd](http://www.e-scientifics.de/content/r-space/mlm/gemi_mlm_1.rmd)

Chapter 18 of Serious Stats introduces multilevel models by considering them as an extension of repeated measures ANOVA models that can cope with missing outcomes, time-varying covariates and can relax the sphericity assumption of conventional repeated measures ANOVA. They can also deal with other – less well known – problems such as having stimuli that are random factor (e.g., see this post on my Psychological Statistics blog). Last, but not least, multilevel generalised linear models allow you to have discrete and bounded outcomes (e.g., dichotomous, ordinal or count data) rather than be constrained by as assuming a continuous response with normal errors.
[Quelle](https://seriousstats.wordpress.com/2013/04/18/using-multilevel-models-to-get-accurate-inferences-for-repeated-measures-anova-designs/)



## Ein paar Punkte

unconditional model: Modell, in dem nur die Ebenen definiert werden. Null-Modell mit Ebenen. Auf jeder Ebene eine Konstante als Vorhersagewert.

Ebenenspezifische Modelle, also für jede Ausprägung von Ebene eine einzelne Analyse. Generalisierbarkeit fraglich.





## Ein paar praktische Tips


Darstellungen hier mit `nlme::lme()`.
Eine gängige Alternative für Multi-Level-Modelle in R ist `lme4::lmer()` oder besser `lmerTest::lmer()`, ein Wrapper für `lme4::lmer()` der auch Signifikanztests mit ausgibt.

Bei Convergenzfehlern, die nicht auf falschen Modellen beruhen: 
[post](http://stats.stackexchange.com/questions/40647/lme-error-iteration-limit-reached)

Correlation in Ausgabe bei Random Effects
Da wir für jede Beobachtung auf Ebenen über Ebene 1 Parameterschätzungen für jede einzelne Beobachtung erhalten, 
können wir diese Parameter miteinander korrelieren. 
Dies ist Teil der Ausgabe.


## Hintergrund: Formulae in `lme4::lmer()`

[Quelle](http://stats.stackexchange.com/questions/18428/formula-symbols-for-mixed-model-using-lme4)

The general trick is, as mentioned in another answer, is that the formula follows the form dependent ~ independent | grouping. 
The grouping is generally a random factor, you can include fixed factors without any grouping 
and you can have additional random factors without any fixed factor (an intercept-only model). 
A + between factors indicates no interaction, a * indicates interaction.

For random factors, you have three basic variants:

-    Intercepts only by random factor: (1 | random.factor)
-    Slopes only by random factor: (0 + fixed.factor | random.factor)
-    Intercepts and slopes by random factor: (1 + fixed.factor | random.factor)

Note that variant 3 has the slope and the intercept calculated in the same grouping, i.e. at the same time. 
If we want the slope and the intercept calculated independently, i.e. without any assumed correlation between the two, we need a fourth variant:

    Intercept and slope, separately, by random factor: (1 | random.factor) + (0 + fixed.factor | random.factor)

There's also a nice summary in response to a related question that you should look at.

If you're up to digging into the math a bit, Barr et al. (2013) summarize the lmer syntax quite nicely in their Table 1, adapted here to meet the constraints of tableless markdown. That paper dealt with psycholinguistic data, so the two random effects are Subject and Item.

Models and equivalent lme4 formula syntax:

Nr  | Formula                                          | Syntax
--- | -----------------------------------------------  | -------------------------------------------
1.  | $Y_{si}=β_0+β_1X_i+e_{si}$                       | n/a (Not a mixed-effects model)
2.  | $Y_{si}=β_0+S_{0s}+β_1X_i+e_{si}$                | `Y ∼ X+(1∣Subject)`
3.  | $Y_{si}=β_0+S_{0s}+(β_1+S_{1s})Xi+e_{si}$        | `Y ∼ X+(1 + X∣Subject)`
4.  | $Y_{si}=β_0+S_{0s}+I_{0i}+(β_1+S_{1s})Xi+e_{si}$ | `Y ∼ X+(1 + X∣Subject)+(1∣Item)`
5.  | $Y_{si}=β_0+S_{0s}+I_{0i}+β_1X_i+e_{si}$         | `Y ∼ X+(1∣Subject)+(1∣Item)`
6.  | As (4), but S0s, S1s independent                 | `Y ∼ X+(1∣Subject)+(0 + X∣ Subject)+(1∣Item)`
7.  | $Y_{si}=β_0+I_{0i}+(β_1+S_{1s})X_i+e_{si}$       | `Y ∼ X+(0 + X∣Subject)+(1∣Item)`


References:

Barr, Dale J, R. Levy, C. Scheepers und H. J. Tily (2013). Random effects structure for confirmatory hypothesis testing: Keep it maximal. Journal of Memory and Language, 68:255– 278.




# Hintergrund: MR als Multi-Level Modell mit Maximum Likelihood Ratio Test


Einfache Regression Körpergewicht und Körpergröße, ML-Parameterschätzung mit lme aus `library(nlme)`.
Modelltest eines Nullmodells gegen ein Modell mit Slope.

```{r}
ddf <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/virt/v_bmi_preprocessed.txt", fileEncoding = "UTF-8")

require(nlme)
m.0 <- lme(fixed= weight ~ 1, random = ~1 | subj, data=ddf, method='ML')
summary(m.0)
m.1 <- lme(fixed= weight ~ 1 + height, random = ~1 | subj, data=ddf, method='ML')
summary(m.1)

anova(m.0, m.1)

# compare it with least square linear model
lm.0 <- lm(ddf$weight ~ 1)
summary(lm.0)
mean(ddf$weight)
lm.1 <- lm(ddf$weight ~ ddf$height)
summary(lm.1)
anova(lm.0, lm.1)

# get a plot for it
plot(ddf$height, ddf$weight)
abline(lm.1, col = "red")

```



# Intraclass Correlation Coefficient (ICC)

proportion of the variation that is being contributed by different levels of random effects.
[source](http://tolstoy.newcastle.edu.au/R/help/06/05/26408.html)

Wenn in höherem Level nur noch wenig Residualvarianz übrig ist, ist der ICC hoch. Die Kontext-Einflüsse (der Level-Ausprägung) sind also stark bzw. haben großen Einfluss auf die Ausprägung der AV.

Je größer der Varianzanteil zwischen den Gruppen (Stufen auf einem Level) ist, desto weniger Varianz besteht anteilig innerhalb der Gruppen. Das bedeutet, dass sich die Werte innerhalb der Gruppen ähnlicher sind, also stärker miteinander korellieren.

Partitionierung der Residualvarianz:

- die Residualvarianz ist die Varianz, die durch das (Gesamt-)Modell nicht aufgeklärt wird (variance of the residuals)
- Varianz des Intercept für einen Faktor ist die Varianz um den Grand Mean, also die Residualvarianz des Grand Mean
- der Varianzanteil des Faktors an der Gesamt-Residualvarianz ist der ICC dieses Faktors
- Residualvarianz(Faktor) / Residualvarianz(Faktor) + Residualvarianz(Modell)
- je kleiner der Anteil der Residualvarianz(Modell) wird, desto näher an 1 ist die ICC und desto mehr lohnt sich die Aufnahme des Faktors als Level
- je größer der Anteil der Residualvarianz(Modell), desto näher an 0 ist die ICC, desto weniger lohnt sich die Aufnahme des Faktors als Level
- je weiter die Faktorstufen auseinanderliegen (große Residualvarianz des Grand Mean) und 
  je weniger Restvarianz innerhalb der Faktorstufe übrig ist (Residualvarianz des Modells), desto größer ist der ICC


todo: Grafik hierzu

ICC-Paper of Andy Field (2012) on companon site.
[paper](https://studysites.uk.sagepub.com/dsur/study/DSUR%20Smart%20Alex-Labcoat%20Leni-Self%20Test%20Answers/DSUR%20Chapter%2019%20Web%20Material.pdf)

ICC and VarCorr

- ICC tells us, whether MLM is necessary? Variance bound by level (group) membership should be considerable (recommendation: >10%) to use MLM instead of GLM
- We use `VarCorr()` to get the necessary variances to calculate ICC
- in `nlme::lme()` we have to co `as.numeric()` on the `VarCorr()` object, because content type is string.
- `VarCorr()` is different in `nlme::lme()` and `lme4::lmer()`

Example in "check" below.




# `lme()` und missing data

[source](http://r.789695.n4.nabble.com/Missing-data-and-LME-models-and-diagnostic-plots-td896351.html)

Mixed models based on likelihood methods can often handle missing
observations within subjects, but they not do well with missing
individual elements in the design matrices (think unit nonresponse vs
item nonresponse in the survey world).  Continuing with the example I
recently sent to you


```{r}

set.seed(0112358)
library(nlme)

# we create a dataframe withount missings
school <- factor(rep(1:20, each=4))
year <- rep(2006:2009, 20)
year.c <- year - mean(year)
tmt <- sample(0:1, 20, replace = TRUE)[school]
math <- rnorm(80, 2 + tmt + .001*year + .0001*tmt*year, 1.5) + rnorm(20)[school]
tmt <- factor(tmt)
df.f <- data.frame(math, school, tmt, year, year.c)
rm(math, school, year, year.c, tmt)

# Now create missing observations,
# -- note just a single obs for the first school:
df.m <- df.f[-c(2:6, sample(75, 5)), ]
table(df.m$school)
#
# school:  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20
# numobs:  1  2  4  3  4  4  4  4  4  3  3  4  4  4  3  4  4  4  4  4

m.f <- lme(math ~ year.c*tmt, data = df.f, random=~1|school)
m.m <- lme(math ~ year.c*tmt, data = df.m, random=~1|school)

intervals(m.f)
intervals(m.m)

```

Note that even with the missing obs, in this instance we got pretty
good estimates of the variance components (parameters: within-group sd
= 1.5; between-school sd = 1).  However, I had to change the seed from
the one I used in the last email because, as is not uncommon with
unbalanced data and relatively small sample sizes, there were
convergence problems.  If you were using classical MoM estimators
fitting a mixed model with unbalanced data would lead to very
precarious inferences (not that the likelihood inference isn't without
it share of pitfalls, but that's another email...)



# Mehrebenenmodelle mit mehr als zwei Ebenen


formula                                 | description
--------------------------------------- | ------------------------------------------------------------
$Y_{ijk} = \gamma_{000} + v_{00k} + u_{0jk} + \varepsilon_{ijk}$                                     | intercept
$+ \gamma_{100} \cdot X_{1_{ijk}} + \gamma_{010} \cdot X_{2_{0jk}} + \gamma_{001} \cdot X_{3_{00k}}$ | Fixed Effects
$+ v_{10k} \cdot X_{1_{ijk}} + u_{1jk} \cdot X_{1_{ijk}}$                                            | random slopes des Level-1-Prädiktors
$+ v_{20k} \cdot X_{2_{0jk}}$                                                                        | random slope des Level-2-Prädiktors
$+ \gamma_{110} \cdot X_{1_{ijk}} \cdot X_{2_{0jk}}$ <br> $+ \gamma_{101} \cdot X_{1_{ijk}} \cdot X_{3_{00k}}$ <br> $+ \gamma_{011} \cdot X_{2_{0jk}} \cdot X_{3_{00k}}$                                                                                         | cross-level-2-fach-Interaktionen
$+ \gamma_{111} \cdot X_{1_{ijk}} \cdot X_{2_{0jk}} \cdot X_{3_{00k}}$                               | cross-level-3-fach-Interaktion




# Hintergrund: Correlation in Summaries von `nlme::gls()` bzw. `nlme::lme()`

Im Summary von `nlme::gls()` bzw. `nlme::lme()` taucht immer auch eine Korrelationstabelle mit auf, sowohl bei den fixed, als auch bei den random Effects.

Der Koeffizient zeigt die Korrelation zwischen Intercept und Slope, wenn der Slope variieren würde, also z. B. bei Wiederholung von Erhebungen. 
Für die einzelne untersuchte Stichprobe gibt es keine berechenbare Korrelation zwischen Intercept und Slope, weil sie für die Stichprobe konstant sind.

Eine Erklärung und [Visualisierung](https://stats.stackexchange.com/questions/171125/correlation-between-ols-estimators-for-intercept-and-slope)
des Intercept in der Correlation Tabelle von `nlme::gls()` bzw. `nlme::lme()`

![Visualisation](bPp50.jpg)

Nach Zentrieren der Daten verschwindet die Korrelation, sie ist also abhängig von der Skalierung der Erklärvariablen.
Meist ist diese Korrelation für die Interpretation irrelevant.


<!--

An example:

```{r}
library(tidyverse)
library(nlme)
# data come from library(mlmRev)
dd <- mlmRev::Exam
dd$normexam.c <- dd$normexam - mean(dd$normexam)
dd$standLRT.c <- dd$standLRT - mean(dd$standLRT)

# effect of centerin in fixed effect
mg <- gls(normexam ~ standLRT, data=dd, metho="ML")
summary(mg)
mg.c <- gls(normexam.c ~ standLRT.c, data=dd, metho="ML")
summary(mg.c)

# here follow some tests {}todo
# categorial predictor school and continuous predictor reading test: standLRT
mg.2 <- gls(normexam ~ school + standLRT, data=dd, method="ML", correlation=NULL)
# summary(mg.2)
# using MLM 
m.2 <- lme(fixed = normexam ~ standLRT,
       random = ~ standLRT | school,
       data=dd,
       method="ML")

names(summary(mg.2))
summary(mg.2)$corBeta[1:6,1:6]

names(summary(m.2))
summary(m.2)$corFixed

# effect of centering

m.c <- lme(fixed = normexam.c ~ standLRT.c,
       random = ~ standLRT.c | school,
       data=dd,
       method="ML")

summary(m.c)
summary(m.c)$corFixed


dd %>% dplyr::filter(school == "14" | school == "5") %>% 
  ggplot(aes(x=standLRT, y=normexam, color = school)) + geom_point() + geom_smooth(method=lm, se=FALSE)

```




The GLS covariance of estimated fixed effects (including the 
intercept) is (X' \hat{\Sigma}^-1 X)^-1, where X is the model matrix 
and \Sigma is the response/error covariance matrix.  This will 
generally have non-zero off-diagonals, implying the fixed effects 
estimates are correlated.  The values you're inquiring about are the 
scaled estimates of those off-diagonals. 

[from](https://stackoverflow.com/questions/45515595/summary-extract-correlation-coefficient)

```{r eval = FALSE}
# The full variance covariance matrix of the coefficient estimates is:
fm <- lm(c1 ~ c2)
vcov(fm)
# and in particular sqrt(diag(vcov(fm))) equals coef(summary(fm))[, 2]
# The corresponding correlation matrix is:
cov2cor(vcov(fm))
# The correlation between the coefficient estimates is:
cov2cor(vcov(fm))[1, 2]

# example
dd <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/div/sat.txt")
mm <- lm(univ_GPA ~ math_SAT + verb_SAT, data=dd)
head(model.matrix(mm))
# we can get the crossprod of the model.matrix(mm) by
t(model.matrix(mm)) %*% model.matrix(mm)
crossprod(model.matrix(mm)) -> mmm
# we call this SSCP (Sums of Squares and Cross Products) Matrix
1/(nrow(model.matrix(mm))-1) * t(model.matrix(mm)) %*% model.matrix(mm)
# this is called the variance covariance matrix

# to get the inverse of this matrix
solve(mmm)
sss <- solve(mmm, diag(ncol(mmm)))  # the resulting matrix is as close as possible to uncorrelated explanatory variables ??
sss
solve(t(model.matrix(mm)) %*% model.matrix(mm))  # covariance matrix
cov2cor(solve(t(model.matrix(mm)) %*% model.matrix(mm)))  # correlation matrix









summary(mm)$sigma   # Residual standard error
summary(mm)$sigma^2 # Residual variance
mm$df.residual      # Residual df

# variance-covariance matrix of model estimates
sss                     * summary(mm)$sigma^2
solve(mmm)              * summary(mm)$sigma^2
solve(model.matrix(mm)) * summary(mm)$sigma^2
# is the same as
vcov(mm)

# use variances from variance-covariance matrix 
# to to calculate the standard error of the estimates (col. "Std.Error" in summary)
std.err <- sqrt(diag(vcov(mm)))



t(model.matrix(mm)) %*% diag(dim(model.matrix(mm))[2])
crossprod((model.matrix(mm)), diag(dim(model.matrix(mm))[2]))
# mm <- lm(univ_GPA ~ math_SAT, data=dd)

# unscaled covariance matrix of the extimated effects (?)
summary(mm)$cov.unscaled

cov2cor(summary(mm)$cov.unscaled)

summary(mm, correlation=TRUE)

head(model.matrix(mm))
coef(summary(mm)) # part of summary output

sqrt(diag(vcov(mm)))
coef(summary(mm))[, 2]
cov2cor(vcov(mm))
cov2cor(vcov(mm))[1, 2]

# wrong: sd(mm$residuals)  # sd of residuals is not exactly standard error of residuals, df differ 

# standard error of estimates is Residual standard error and can be calculated by
sqrt(sum(mm$residuals^2)/(nrow(dd)-2))
sqrt(sum(mm$residuals^2)/mm$df.residual)

# Constant
sqrt(sum((mm$fitted.values - mm$coefficients[1])^2)/mm$df.residual)



sqrt(diag(vcov(mm)) * sd(mm$residuals)^2)

matrix(c(1,2,4,3,0,-1), nrow=3)
A <- matrix(c(1,4,0,  2,3,-1), nrow = 3)
B <- matrix(c(4,2,1,0,0,1), nrow = 2)
A
B
A %*% B
```

### Standard error of the estimate

sqrt(sum(hat(y) - y)^2 / df-residuals)
df-residuals = 2 when 1 explanatory and 1 reaction variable

[source](https://stats.stackexchange.com/questions/324260/manually-calculate-the-parameters-std-error-of-lm-output-in-r/324651)
Std. Error of the Intercept:
sqrt (( 1/n ]) + (x_mean)^2 / Σ(xi - x_mean)^2)



[source](https://www.youtube.com/watch?v=r-txC-dpI-E)

Sol https://www.datascienceblog.net/post/machine-learning/linear_models/
https://statistikguru.de/lexikon/sscp-matrix.html

-->

# Parameter Corr in `nlme::lme()` {#corr}

Inspired by [StackExchange](https://stats.stackexchange.com/questions/213719/r-default-correlation-in-nlmelme)
[https://stats.stackexchange.com/questions/367957/autoregression-in-nlme-undestanding-how-to-specify-corar1]

Daten sind simuliert. 
Wir können uns vorstellen, `es`value` seien die Stimmungs- bzw. Befindlichkeitswerte als Rating zwischen -10 und 10 im Laufe einer Woche.

```{r}
set.seed(100)
n.subj <- 20
n.times <- 6
# subj.c <-rep(NA, n.subj)
subj.c <- rnorm(n.subj, sd=2)
subj.s <- rnorm(n.subj, sd=1)
mat = matrix(ncol = n.times, nrow = n.subj)
# converting the matrix to data frame
dd=data.frame(mat)
# line <- rep(subj.c[1], n.times)
# line <- line + rnorm(n.subj, sd=0.2)
for (subj in 1:n.subj) {
  resid <- rnorm(n.times, 0, sd=2)
  for (tt in 1:n.times) {
    dd[subj,tt] <- subj.c[subj] +  (tt - n.times/2) * subj.s[subj] + resid[tt]
  }
}
for (cc in 1:ncol(dd)) {colnames(dd)[cc] <- paste('t', as.character(cc), sep="_")}
# for (cc in 1:ncol(dd)) {colnames(dd)[cc] <- as.numeric(cc)}
dd <- cbind(1:nrow(dd), dd)
colnames(dd)[1] <- "subj"
dd.w <- dd
dd <- dd.w %>% tidyr::pivot_longer(cols = 2:ncol(dd), names_to = "time.s", values_to = "value")
dd <- dd %>% dplyr::rowwise() %>% dplyr::mutate(time = as.numeric(strsplit(time.s, split="_")[[1]][2]))
# dd <- dd %>% dplyr::mutate(time.n = as.numeric(time))
dd %>% ggplot(aes(x=time, y=value, group=subj, color=factor(subj))) +
  geom_point() + 
  geom_line()


#value = rnorm(100)
#time = rep(1:5, each = 20)
#id = rep(letters[1:20],5)
#dd = data.frame(value, time, id)

fit.0 <- dd %>% nlme::lme(value ~ time, random=~ 1 | subj,
              correlation = NULL,
              data=.)
fit.1 <- dd %>% nlme::lme(value ~ time, random=~ 1 | subj,
              correlation=corAR1(),
              data=.)

# random intercept and slope model without correlation between time steps
fit.2a <- dd %>% nlme::lme(value ~ time, random=~ 1 + time | subj,
              correlation=NULL,
              data=.)
fit.2b <- dd %>% nlme::lme(value ~ time, random=~ 1 + time | subj,
              correlation=corAR1(form=~1), # time is equally spaced (one daily value measured)
              # correlation=corAR1(form=~1|subj), # time is equally spaced and subj are uncorrelated
              data=.)
# just to show syntax and idea, corCAR1
# random intercept and slope model with time as continuous variable for correlation structure per subj
fit.2c <- dd %>% nlme::lme(value ~ time, random=~ 1 + time | subj,
              correlation=corCAR1(form=~time),  # this would also account for varying time lags
              data=.)


anova(fit.2a, fit.2b)

anova(fit.0, fit.1, fit.2a)

summary(fit.0)
summary(fit.1)
summary(fit.2b)



```



# Beispiele


## Beispiel OxBoys

[source](https://rpubs.com/haj3/multimodel)


This dataset contains three features from 26 students in Oxford: student id, (measured) height, and (measurement) year.


```{r}
df <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/div/oxboys.txt", fileEncoding = "UTF-8")
head(df)

# Plot the relationship between height and year, and draw a linearly regressed line ignoring the id variable.
plot(df$occasion, df$height, col="blue", main="Oxford Students: Height vs. Year", cex.lab = 1.5, cex.main = 1.5, cex.axis = 1.5, xlim = c(0,10))
abline(lm(height~occasion, data=df), col="red", lwd=5)
library(nlme)

# Plot the relationship between height and year, but this time, fit a different linear regression for each individual. In ggplot, you can do it by specifying the group variable: ggplot(oxboys, aes(x=year, y=height, group=id)) + stat_smooth(method=“lm”) Other plotting packages also should have some such facility. Comment on differences between (a) and (b).
library(ggplot2)

ggplot(data=df, aes(x=occasion, y=height, group=subj)) + 
  stat_smooth(method="lm", color="red") + 
  ggtitle("Oxford Boys: Height vs. Occasion") +
  expand_limits(x=c(min(df$occasion),max(df$occasion))) +
  scale_x_continuous(breaks=seq(min(df$occasion), max(df$occasion), 1)) + 
  theme(axis.title = element_text(face="bold", colour="blue", size=20), 
    plot.title = element_text(family = "Trebuchet MS", color="blue", face="bold", size=22))

```

The first graph plots the relationship between year and height for all students. The abline represents the line of best fit for a global MLR model. In contrast, the ggplot graph shows the relationship between year and height with respect to each student (group=id). This is a graphic representation of a local model (i.e., 26 different linear regressions).



## Example Pitch


[Quelle](https://seriousstats.wordpress.com/2013/04/18/using-multilevel-models-to-get-accurate-inferences-for-repeated-measures-anova-designs/)
or [see](http://seriousstats.wordpress.com/tag/repeated-measures-anova/)

Pitch hier Stimmlage

**Description**: This experiment (from a collaboration with Tim Wells and Andrew Dunn) 
involves looking at the pitch of male voices making attractiveness ratings with respect to female faces. 
The effect of interest (for this example) is whether average pitch goes up or done for higher ratings (and if so, by how much). 
A conventional ANOVA is problematic because this is a design with two fully crossed random factors – each participant (n = 30) sees each face (n = 32) 
and any conclusions ought to generalise both to other participants and (crucially) to other faces. 
Furthermore, there is a time-varying covariate – the baseline pitch of the numerical rating when no face is presented.


### Anwendung von `lme4::lmer()`

```{r}
# we could get the original data via
# pitch.dat <- read.csv('http://www2.ntupsychology.net/seriousstats/pitch.csv')
# we use a local copy for convenience
# write.table(pitch.dat, "pitch.txt", quote=F, sep="\t", row.names=F)
df <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/div/pitch.txt")

head(df)

require(psych)
describe(df)
require(Rmisc)
summarySE(data=df, "pitch")
summarySE(data=df, "pitch", c("attract"))
summarySE(data=df, "pitch", c("Context"))
summarySE(data=df, "pitch", c("attract", "Context"))

require(ggplot2)
ggplot(data=df, aes(x=attract, y=pitch)) + geom_point() + geom_smooth(method = lm, se=F)
# regression line indicates a slightly growing pitch with growing attractivity

# require(ez)
# ezANOVA(data=)

library(lme4)

pitch.me <- lmer(pitch ~ base + attract + (1|Face) + (1|Participant), data=df)

pitch.me

coef(pitch.me)
fixef(pitch.me)
ranef(pitch.me)

head(predict(pitch.me))


```



## Beispiel Auswirkung von Schlafdeprivation auf Reaktionszeiten

Der Datensatz ist Teil der `library(lme4)`.
Beschreibung des Datensatzes findet sich unter `help(sleepstudy)`

```{r}
require(lme4)
df <- sleepstudy
head(sleepstudy)

# some descriptives
require(psych)
psych::describe(sleepstudy)
require(Rmisc)
summarySE(data=df, "Reaction")
summarySE(data=df, "Reaction", "Days")
summarySE(data=df, "Reaction", "Subject")

require(ggplot2)
ggplot(data=df, aes(x=Days, y=Reaction), color=Subject) + geom_point(aes(color=Subject)) + geom_smooth(method = lm, se=F)

# we convert Days to type factor
df$f.Days <- factor(df$Days)

# the classical anova approach
require(ez)
m.ez <- ezANOVA(data=df, dv=.(Reaction), wid = .(Subject), within = .(f.Days), detailed=TRUE, type=2, return_aov = T)
m.ez
# descriptives
ezStats(data=df, dv=.(Reaction), wid = .(Subject), within = .(f.Days))
# and a plot
ezPlot(data=df,  dv=.(Reaction), wid = .(Subject), within = .(f.Days), x=.(f.Days))
       
# we store aov element for convenience
o.aov <- m.ez[['aov']][['Subject:f.Days']]
# we calculate predicted values and store them for further use
# reference group is mean at day zero 256.6518
mean.day.0 <- 256.6518
df$pred.ez <- NA
adds.one <- c(0, o.aov[['coefficients']])
adds <- rep(adds.one, length(unique(df$Subject)))
# predicted values are group means
# we get to group means by adding the coefficients to the mean of the reference group
df$pred.ez <-adds + mean.day.0 

# correlation of predicted values and reaction variable 
cor(df$Reaction, df$pred.ez)

# standard deviation of residuals
sd(df$Reaction - df$pred.ez)

require(lattice)
xyplot(Reaction ~ Days | Subject, sleepstudy, type = c("g","p","r"),
       index = function(x,y) coef(lm(y ~ x))[1],
       xlab = "Days of sleep deprivation",
       ylab = "Average reaction time (ms)", aspect = "xy")
(fm1 <- lmer(Reaction ~ Days + (Days|Subject), sleepstudy))
(fm2 <- lmer(Reaction ~ Days + (1|Subject) + (0+Days|Subject), sleepstudy))

```




```{r}
require(lme4)
df <- sleepstudy
# detach package lme4, we only need it's sleepstudy data 
# detach("package:lme4", unload=TRUE)

require(nlme)
m.0 <- lme(fixed= Reaction ~ 1, random = ~1 | Subject, data=df, method='ML')
summary(m.0)
m.1 <- lme(fixed= Reaction ~ 1 + Days, random = ~1 | Subject, data=df, method='ML')
summary(m.1)

m.2 <- lme(fixed= Reaction ~ 1 + Days, random = ~1 + Days | Subject, data=df, method='ML')
summary(m.2)

anova(m.0, m.1, m.2)

# an lme object "understands" predict 
df$pred.m.1 <- predict(m.1)
cor(df$Reaction, df$pred.m.1)
df$pred.m.2 <- predict(m.2)
cor(df$Reaction, df$pred.m.2)

# variance of residuals
sd(m.1$residuals)
# ... is smaller than in classical analysis

# we can get at the coefficients via
fixef(m.1)
ranef(m.1)

# we extend the data with effects of model m.1
df$inter <- NA
df$slope <- NA
require(dplyr)
df <- df %>% dplyr::group_by(Subject) %>% dplyr::mutate(inter=fixef(m.1)[1] + ranef(m.1)[Subject,1])
df <- df %>% dplyr::mutate(slope=fixef(m.1)[2])

# we visualize the random intercept fixed slope model
require(ggplot2)
pplot <- ggplot(data=df, aes(x=Days, y=Reaction), color=Subject) + 
  geom_point(aes(color=Subject)) +
  geom_abline(aes(intercept=inter, slope=slope, color = Subject)) +
  geom_smooth(method = lm, se=F)
pplot

# we extend the data with effects of model m.2
df$inter <- NA
df$slope <- NA
require(dplyr)
df <- df %>% dplyr::group_by(Subject) %>% dplyr::mutate(inter=fixef(m.2)[1] + ranef(m.2)[Subject,1])
df <- df %>% dplyr::group_by(Subject) %>% dplyr::mutate(slope=fixef(m.2)[2] + ranef(m.2)[Subject,2])


# we visualize the random intercept random slope model
pplot <- ggplot(data=df, aes(x=Days, y=Reaction), color=Subject) + 
  geom_point(aes(color=Subject)) +
  geom_abline(aes(intercept=inter, slope=slope, color = Subject)) +
  geom_smooth(method = lm, se=F, size=3) 
pplot








# pplot <- ggplot(data=df, aes(x=Days, y=Reaction), color=Subject) + geom_point(aes(color=Subject)) 
# pplot + geom_smooth(method = lm, se=F)
# 
# # data for ggplot have to be in one dataframe
# require(dplyr)
# df.t <- dplyr::data_frame(subj = 1:length(ranef(m.2)[,1]), r.int = ranef(m.2)[,1], r.slope = ranef(m.2)[,2], f.int <- fixef(m.2)[1], f.slope <- fixef(m.2)[2])
# names(df.t) <- c("subj", "r.int", "r.slope", "f.int", "f.slope")
# df.t <- df.t %>% dplyr::mutate(intercept = f.int + r.int, slope = f.slope + r.slope)
# df.t <- df.t %>% dplyr::mutate(x0 = 0, y0 = intercept, x9 = 9, y9 = intercept + 9 * slope, Days = 0, Reaction=300)
# head(df.t)
# 
# # df.t$subj <- 1:length(ranef(m.2)[1])
# # df.t <- cbind(df.t, ranef(m.2))
# # df.t$fix.intercept <- fixef(m.2)[1]
# # df.t$fix.slope     <- fixef(m.2)[2]
# # names(df.t) <- c("subj", "r.int", "r.slope", "f.int", "f.slope")
# # df.t <- df.t %>% dplyr::mutate(intercept = f.int + r.int, slope = f.slope + f.int)
# # head(df.t)
# 
# # we extend the data with effects
# df$inter <- NA
# df$slope <- NA
# for (ss in rownames(ranef(m.2))){
#   df$inter[df$Subject == ss] <- fixef(m.2)[1] + ranef(m.2)[ss,1]
#   df$slope[df$Subject == ss] <- fixef(m.2)[2] + ranef(m.2)[ss,2]
# }
# 
# #pplot <- pplot + geom_line(data=df.t, aes(intercept, slope), color=subj)
# #pplot <- ggplot(data=df, aes(x=Days, y=Reaction), color=Subject, show.legend=F) + geom_point(aes(color=Subject), show.legend=F) 
# pplot <- ggplot(data=df, aes(x=Days, y=Reaction), color=Subject) + geom_point(aes(color=Subject)) 
# pplot <- pplot + geom_abline(aes(intercept=inter, slope=slope, color = Subject)) 
# #pplot <- pplot + geom_line() 
# #pplot <- pplot + geom_segment(data=df.t, mapping=aes(x=x0, y=y0, xend=x9, yend=y9))
# #pplot <- pplot + geom_segment(data=df.t, mapping=aes(x=x0, y=y0, xend=x9, yend=y9, color=factor(subj)))
# pplot <- pplot + geom_smooth(method = lm, se=F)
# pplot
# 
# #pplot <- pplot + geom_line(data=df.t, aes(mapping = NULL, x=x0, y=intercept), color=subj)
# #pplot <- pplot + geom_line(data=df.t, aes(intercept = intercept, slope = slope, color=factor(subj)))
# pplot
# 
# 
# 
# 
# pplot + geom_abline(intercept = fixef(m.2)[1], slope = fixef(m.2)[2])              
# for (nn in 1:length(ranef(m.2)[,1])){
#   #cat(ranef(m.2)[nn])
#   c.int   <- ranef(m.2)[nn,1] + fixef(m.2)[1]
#   #cat(c.int)
#   c.slope <- ranef(m.2)[nn,2] + fixef(m.2)[2]
#   #cat(c.slope)
#   pplot <- pplot +  geom_abline(intercept = c.int, slope = c.slope)      
# }
# pplot <- pplot + geom_abline(intercept = fixef(m.2)[1], slope = fixef(m.2)[2], size=2)              
# pplot 



```



## Beispiel Schule, Abschlussnote und Schuleingangstest-Niveau


[Quelle](https://www.thetransfertest.com/schools-intake-information)

Im britischen Schulsystem gibt es beim Übergang von Primary School zu Secondary School einen standardisierten Test, 
in Nordirland beispielsweise der "Common Entrance Assessment (CEA)" durchgeführt von der AQE (Association for Quality Education). 
Auf Basis der Ergebnisse dieses Tests können Eltern ihre Kinder in den Highschools anmelden. 
Die Highschools haben unterschiedliche Minimal-Punktwerte, also Cut-Off-Werte, um Schüler zu akzeptieren.
Die Schuleingangstests haben einen Mittelwert von 100.

Ein Beispiel: Belfast High School. Hier wurden 2016 nur überdurchschnittliche Schüler aufgenommen.
```
AQE Highest 2016: 124
AQE Lowest 2016: 101
AQE Applications 2016: 220
AQE 1st Preference 2016: 179
AQE Admitted 2016: 136
```

Die Abschlussnote von Schülern soll vorhergesagt werden `normexam`.
Die Schüler sind an verschiedenen Schulen.
Vorhersagevariable auf Level Schüler ist `standLRT` (Literarisches Textverständnis).
Vorhersagevariable auf Schul-Level ist `schavg` (Durchschnittliches Aufnahmeniveau der jeweiligen Schule).

Die Schul-Level Variablen erkennen wir daran, dass ihre Ausprägung für die jeweilige Schule gleich ist, also innerhalb einer Schule für alle Schülerinnen und Schüler gleich bleiben (sich wiederholen).

Der verwendete Datensatz beschreibt 65 Schulen in London von insgesamt 4059 Schülern.
Infos hierzu gibt es mit `?mlmRev::Exam`


A data frame with 4059 observations on the following 9 variables.

var       | description
--------- | -----------------------------------------------------------------------------------------------------------------------
school    | School ID - a factor.
normexam  | Normalized exam score.
schgend   | School gender - a factor. Levels are mixed, boys, and girls.
schavg    | School average of intake score.
vr        | Student level Verbal Reasoning (VR) score band at intake - a factor. Levels are bottom 25%, mid 50%, and top 25%.
intake    | Band of student's intake score - a factor. Levels are bottom 25%, mid 50% and top 25%./
standLRT  | Standardised LR test score (London Reading Test - a reading test taken when they were 11 years old)
sex       | Sex of the student - levels are F and M.
type      | School type - levels are Mxd and Sngl.
student   | Student id (within school) - a factor




```{r}
require(nlme)
require(mlmRev)
# mlmRev hat data Exam
dd <- mlmRev::Exam
names(dd) 
head(dd)
require(Rmisc)
summarySE(dd, "normexam", c("school"))

require(psych)
psych::describe (dd)

# we specify the null model
m.0 <- lme(fixed = normexam ~ 1,   # 1, a constant, predicts and results in mean
         data = dd,
         random = ~ 1 | school,
         method="ML")      # levels are defined, school defines level two
m.0
# ... we find mean in fixed (Intercept)
# ... we find sd of residuals for (Intercept) level 2 and level 1 Residual


# we add a level 1 predictor
# random intercept, fixed predictor in individual level
m.1 <- lme(fixed = normexam ~ standLRT, random = ~ 1 | school, data = dd, method="ML") 
m.1
# ... sd of Intercept goes down
# ... sd of Residuals goes down
fixef(m.1)
head(ranef(m.1))


# we add a level 2 predictor
# random intercept, random slope
m.2 <- lme(fixed = normexam ~ standLRT, random = ~ standLRT | school, data = dd, method="ML")
m.2
fixef(m.2)
head(ranef(m.2))


# we can add a level 2 predictor on level 1 also
# random intercept, individual and group level predictor. There is a random slope also.
m.3 <- lme(fixed = normexam ~ standLRT + schavg, random = ~ standLRT | school, data = dd, method="ML")
m.3
fixef(m.3)
head(ranef(m.3))

# we can also add cross-level interaction
# random intercept, cross-level interaction
m.4 <- lme(fixed = normexam ~ standLRT * schavg, random = ~ standLRT | school, data=dd, method="ML") 
fixef(m.4)
head(ranef(m.4))


# check model (Null-Model is automatically generated, deviance of null-model can be found in model$null.deviance)
deviance(m.0)
deviance(m.1)
anova(m.0)
anova(m.1, test="F")
# which is default
anova(m.1)


```




## Beispiel: Schultyp und Zukunftsperspektive


Ecologic fallacy: Systematiken auf höheren Ebenen werden als Effekt auf level 1 fehlinterpretiert.

`y` ist Reaktionsvariable
`x` ist Vorhersagevariable auf Level 1
`treatment` ist Gruppierungsvariable für Level 2
`g.eff` ist Variable auf Level 2
ev. `subj` kann benutzt werden als zweite Level 2 Variable (Betrachtung als Messwiederholungsdesign)


Szenario:

- x ist pot. kogn. Leistungsfähigkeit
- y ist Zukunftsaussichten, Zuversicht in die Zukunft
- Gruppen sind Schulformen: Haupt- Realschule, Gymnasium
- Effekt: Das Wissen um die Aussichtslosigkeit in Hauptschulen


Repeated-Measures-ANOVA würde Mittelwertsunterschiede zwischen den Schultypen prüfen. 
Schultyp würde als kategoriale Variable aufgefasst.


(Es handelt sich um simulierte Daten)

```{r}
# we read data like generated above
df <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/div/schooltype_future.txt")
df$schooltype <- factor(df$schooltype)
# correlation is
cor(df$x, df$y)

# complete scatterplot seems to visualize this correlation, we might notice some asymmetry in distribution
require(ggplot2)
# we visualize the random intercept random slope model
pplot <- ggplot(data=df, aes(x=x, y=y)) + 
  #geom_point(aes(color=Subject)) +
  geom_point() + 
  #geom_abline(aes(intercept=inter, slope=slope, color = Subject)) +
  geom_smooth(method = lm, se=F, size=1) +
	xlab('cognitive abilities') + 
	ylab('future perspectives')
pplot

# but if we look at the groups ...
pplot <- ggplot(data=df, aes(x=x, y=y, color=schooltype, shape=schooltype)) + 
  #geom_point(aes(color=Subject)) +
  geom_point() +
  #geom_smooth(method="lm", fill=NA)
  #geom_abline(aes(intercept=inter, slope=slope, color = Subject)) +
  #geom_smooth(method = lm, se=F, size=2, ) 
	xlab('cognitive abilities') + 
	ylab('future perspectives')
pplot

# group specific correlations show huge differences
cor(df[df$schooltype == 'a', 'x'], df[df$schooltype=='a', 'y'])
cor(df[df$schooltype == 'b', 'x'], df[df$schooltype=='b', 'y'])
cor(df[df$schooltype == 'c', 'x'], df[df$schooltype=='c', 'y'])

# and the descriptives
require(Rmisc)
summarySE(data=df, "x", c("schooltype"))
summarySE(data=df, "y", c("schooltype"))
summarySE(data=df, "y")

# require(psych)
# psych::describeBy(df, group="schooltype")
# psych::describe(df)

# we fit a repeated measures anova
require(ez)
m.ez <- ezANOVA(data=data.frame(df), dv=y, wid=subj, within=schooltype, detailed=T, return_aov = T)
m.ez

require(nlme)
# for demonstration reasons we fit a level 1 only zero model
df.t <- cbind(df, 1)
names(df.t) <- c("subj", "schooltype", "g.eff", "x", "y", 'const')
m.00 <- lme(fixed = y ~ 1, random = ~ 1 | const, data=df.t, method='ML')
# fitted value is the same for all subjects
unique(predict(m.00))
# correlation cannot be computed due to the missing variance in the predicted values

# zero multi level model
m.0 <- lme(fixed = y ~ 1, random = ~1 | schooltype, data=df, method='ML')
summary(m.0)
# get colors
g <- ggplot_build(pplot)
c.colors <- unique(g$data[[1]]["colour"])
# fitted value is the same for all subjects within a level2 group
unique(predict(m.0))

# ... random intercept of schooltype 'explains' variance and reduces residual variance
pplot + 
  #geom_hline(yintercept=m.0$coefficients$random[['schooltype']][1], color=c1) 
  geom_hline(yintercept=m.0$coefficients$random[['schooltype']][1], color=c.colors[1,1]) +
  geom_hline(yintercept=m.0$coefficients$random[['schooltype']][2], color=c.colors[2,1]) +
  geom_hline(yintercept=m.0$coefficients$random[['schooltype']][3], color=c.colors[3,1])

# level 1 predictor, random intercept
m.1 <- lme(fixed = y ~ 1 + x, random = ~1 | schooltype, data=df, method='ML')
summary(m.1)

pplot + 
  geom_abline(intercept = fixef(m.1)[1] + ranef(m.1)[1,1], slope = fixef(m.1)[2], color=c.colors[1,1]) + 
  geom_abline(intercept = fixef(m.1)[1] + ranef(m.1)[2,1], slope = fixef(m.1)[2], color=c.colors[2,1]) + 
  geom_abline(intercept = fixef(m.1)[1] + ranef(m.1)[3,1], slope = fixef(m.1)[2], color=c.colors[3,1]) + 
	geom_abline(intercept = fixef(m.1)[1], slope = fixef(m.1)[2], size = 2)


# we don't want to run into an optimization error ("iteration limit reached without convergence")
ctrl <- lmeControl(opt='optim');
# flow.lme <- lme(rate ~ nozzle, error= nozzle|operator, control=ctrl, data=Flow);
# level 2 variable, random intercept and random slope
m.2 <- lme(fixed = y ~ 1 + x, random = ~ 1 + x | schooltype, data=df, control=ctrl, method='ML')
summary(m.2)

# we visualize the group specific random intercept and random slope
pplot + 
  geom_abline(intercept = fixef(m.2)[1] + ranef(m.2)[1,1], slope = fixef(m.2)[2] + ranef(m.2)[1,2], color=c.colors[1,1]) + 
  geom_abline(intercept = fixef(m.2)[1] + ranef(m.2)[2,1], slope = fixef(m.2)[2] + ranef(m.2)[2,2], color=c.colors[2,1]) + 
  geom_abline(intercept = fixef(m.2)[1] + ranef(m.2)[3,1], slope = fixef(m.2)[2] + ranef(m.2)[3,2], color=c.colors[3,1]) + 
	geom_abline(intercept = fixef(m.2)[1], slope = fixef(m.2)[2], size = 2)

# we compare the models using liklihood ratio test
anova(m.0, m.1, m.2)
anova(m.0, m.2)

# we compare the explained variance of the various models
cor(df$y, predict(m.0))
cor(df$y, predict(m.1))
cor(df$y, predict(m.2))
# ... the school-type specific random slope increases model fit a lot

# we compare the fixed effects
fixef(m.0)
fixef(m.1)
fixef(m.2)
# ... in m.1 and m.2 we see the fixed parameters (intercept and slope)
# ... for model m.2 this is the fixed part of the parameters, random part has to be added to this

# we compare the random effects
ranef(m.0)
ranef(m.1)
ranef(m.2)
# ... for every school-type we get a specific intercept
# ... random intercepts of m.1 are almost the same when we add level 1 predictor x
# ... we can see the level 2 random slope component to be added to the fixed slope
# ... we can see the big differences in random slope, due to the school-type

# coefficients() works also with mle result objects
# it returns the aggregated coefficients
# fixed and random parts are already added
coefficients(m.0)
coefficients(m.1)
coefficients(m.2)


# a look at the variances
VarCorr(m.0)
VarCorr(m.1)
VarCorr(m.2)
# ... variances of residuals go down


# we look at ICC
VarCorr(m.0)
var.u0 = as.numeric(VarCorr(m.0)[1])
var.epsilon = as.numeric(VarCorr(m.0)[2])
ICC = var.u0 / (var.u0 + var.epsilon)
ICC

#require(psychometric)
#ICC1.lme(y, grp=schooltype, data=df)


```

Zur Interpretation der Ausgabe, speziell zu random StdDev 
[post](http://stats.stackexchange.com/questions/69882/why-computing-the-variance-of-the-extracted-random-effect-using-ranef-is-not-t)

sowie 

[post2](http://stats.stackexchange.com/questions/24452/how-to-interpret-variance-and-correlation-of-random-effects-in-a-mixed-effects-m)








## Example Kowles: Personality and School 

[example_personality_school](example_personality_school.html)

Verschiedene Herangehensweise an denselben Datensatz.





## Beispiel Schulen von Torres Reyna

[source](http://www.princeton.edu/~otorres/Multilevel101.pdf)

Variable | Beschreibung
--- | ---
school | Schule, 1...65
student | Schüler-Nummer, 1 ... 198, unter Schule genested, maximal 198 Schüler aus einer Schule
y | Score (of what?) vielleicht Abschlussnote
X1 | Reading test 
x2 | binär Female = 1
x3 | Type of school Schulform: 



```{r}
# the data file comes in Stata format *.dta
require(foreign)
dd <- read.dta("http://md.psych.bio.uni-goettingen.de/mv/data/div/schools.dta")
attributes(dd)$var.labels
# we rename vars according to labels
require(dplyr)
dd <- dplyr::rename(dd, score=y, reading.test=x1, gender.f=x2, schooltype=x3)
# and make school and student factors
dd$school <- factor(dd$school)
dd$student <- factor(dd$student)
dd$gender <- factor(dd$gender)

require(ggplot2)
# get a plot of variability
ggplot(data=dd, aes(x=school, y=score)) + geom_point() + geom_smooth()
# ... schools seem to be different in their score




# individual regressions
ggplot(data=dd, aes(x=reading.test, y=score, color=school)) + geom_point() + geom_smooth(method="lm", se=F)
# or a bit clearer
ggplot(data=dd, aes(x=reading.test, y=score, color=school))                + geom_smooth(method="lm", se=F)
# ... slopes seem to vary considerably


# global correlation
cor(dd$reading.test, dd$score)
# school specific correlation
require(dplyr)
cors <- dd %>% dplyr::group_by(school) %>% dplyr::summarize(COR=cor(score,reading.test))
# we sort correlations
cors <- dplyr::arrange(cors, COR)
head(cors)
# something is strange with school 48, we take a look at it
dd[dd$school == 48,]
# we exclude it, only two students took part
cors <- cors[2:length(cors$school),]
cors$o.nr <- 1:length(cors$school)
head(cors)
tail(cors)
# a graph of the correlations in schools
ggplot(data=cors, aes(x=o.nr, y=COR)) + geom_point() + geom_smooth(se=F)
# our correlations vary substancially




# we fit a simple anova
require(ez)
m.ez <- ezANOVA(data=dd, dv=score, wid=student, between=school, detailed=T, return_aov = T)
m.ez
summary(m.ez$aov)
# ... shows a significant factor school. And ...
# multiple correlation
cor(m.ez[["aov"]]$fitted.values, dd$score)
# squared multiple correlation
cor(m.ez[["aov"]]$fitted.values, dd$score)^2
# ... not too much explained variance in score
rm(m.ez)

require(nlme)
# zero multi level model
m.0 <- lme(fixed = score ~ 1, random = ~1 | school, data=dd, method='ML')
summary(m.0)
# get colors
# g <- ggplot_build(pplot)
# c.colors <- unique(g$data[[1]]["colour"])

# fitted value is the same for all subjects within a level2 group
head(predict(m.0), n=100)
unique(predict(m.0))
# but means vary
unique(predict(m.0))[order(unique(predict(m.0)))]
mean(dd$score[dd$school == 1])
# ... not surprising: each school has its own mean score which is the predicted value
ggplot(data=dd, aes(x=reading.test, y=predict(m.0), color=school)) + geom_point()

# level 1 predictor, random intercept
m.1 <- lme(fixed = score ~ 1 + reading.test, random = ~1 | school, data=dd, method='ML')
summary(m.1)
ggplot(data=dd, aes(x=reading.test, y=predict(m.1), color=school)) + geom_point()
ggplot(data=dd, aes(x=predict(m.1), y=score, color=school)) + geom_point()


# level 2 variable, random intercept and random slope
#m.2 <- lme(fixed = score ~ 1 + reading.test, random = ~ 1 + reading.test | school, data=dd, control=ctrl, method='ML')
m.2 <- lme(fixed = score ~ 1 + reading.test, random = ~ 1 + reading.test | school, data=dd, method='ML')
summary(m.2)

ggplot(data=dd, aes(x=reading.test, y=predict(m.2), color=school)) + geom_point()


# we compare the models using liklihood ratio test
anova(m.0, m.1, m.2)
anova(m.0, m.2)

# we compare the explained variance of the various models
cor(dd$score, predict(m.0))
cor(dd$score, predict(m.1))
cor(dd$score, predict(m.2))
# ... the school-type specific random slope increases model fit a lot

# we compare the fixed effects
fixef(m.0)
fixef(m.1)
fixef(m.2)
# ... in m.1 and m.2 we see the fixed parameters (intercept and slope)
# ... for model m.2 this is the fixed part of the parameters, random part has to be added to this

# we compare the random effects
head(ranef(m.0))
head(ranef(m.1))
head(ranef(m.2))
# ... for every school we get a specific intercept
# ... random intercepts of m.1 are almost the same when we add level 1 predictor x
# ... we can see the level 2 random slope component to be added to the fixed slope
# ... we can see the big differences in random slope, due to the school

# coefficients() works also with mle result objects
# it returns the aggregated coefficients
# fixed and random parts are already added
head(coefficients(m.0))
head(coefficients(m.1))
head(coefficients(m.2))


# a look at the variances
VarCorr(m.0)
VarCorr(m.1)
VarCorr(m.2)
# ... variances of residuals go down


# we look at ICC
VarCorr(m.0)
var.u0 = as.numeric(VarCorr(m.0)[1])
var.epsilon = as.numeric(VarCorr(m.0)[2])
ICC = var.u0 / (var.u0 + var.epsilon)
ICC

# ICC m.1
as.numeric(VarCorr(m.1)[1]) / (as.numeric(VarCorr(m.1)[1]) + as.numeric(VarCorr(m.1)[2]))
# ICC m.2
as.numeric(VarCorr(m.2)[1]) / (as.numeric(VarCorr(m.2)[1]) + as.numeric(VarCorr(m.2)[2]))

require(psychometric)
ICC1.lme(score, grp=school, data=dd)



```




## Beispiel Xavier Fernández-i-Marín (not available any more)

[source](http://xavier-fim.net/R/scripts/04-multilevel.R.html)

Leider keine Dokumentation zu den Daten.

Inhalt

Die SES erlaubt die Messung und differenzierte Beschreibung der subjektiv wahrgenommenen Schmerzen. Der Fragebogen besteht aus 24 Items, die sich 5 Teilskalen zuordnen lassen. Zwei Skalen beschreiben affektive Merkmale der Schmerzempfindung: «allgemeine affektive Schmerzangabe» und «Schmerzangabe der Hartnäckigkeit». Beide Merkmale aufsummiert bilden die Globalskala «SES-affektiv». Die drei weiteren Skalen beschreiben sensorische Aspekte der Schmerzempfindung: «sensorische Schmerzangabe der Rhythmik», «sensorische Schmerzangabe des lokalen Eindringens» und «sensorische Schmerzangabe der Temperatur». Die drei Merkmale aufsummiert bilden die Globalskala «SES-sensorisch». Um die Ökonomie der Durchführung bei umfassenden Studien zu erhöhen, sind drei Survey-Versionen, getrennt nach dem jeweiligen Beurteilungszeitraum, beziehbar. Instruktion und Items sind auf nur einem DIN-A4-Bogen abgebildet.

```{r eval = FALSE}
### !/usr/local/bin/R

# Xavier Fernández-i-Marín
# December 2004
# http://xavier-fim.net/

# load the example datasets from HLM/AppendixA/Examples
# Download them directly. Requires Internet connection
hsb1 <- read.table(file="http://xavier-fim.net/R/scripts/hsb1.dat")
hsb2 <- read.table(file="http://xavier-fim.net/R/scripts/hsb2.dat")

# get the variables and basic information of each dataset
names(hsb1)
length(hsb1$ID)
summary(as.data.frame(hsb1))

names(hsb2)
length(hsb2$ID)
summary(as.data.frame(hsb2))

# load the package to do Linear Mixed-Effects Models
library(nlme)

# create a new dataset with the merged data
ds <- merge(as.data.frame(hsb1), as.data.frame(hsb2))

### model specification

# 1
mod.1 <- lme(fixed = MATHACH ~ 1,
             random = ~ 1 | ID,
           data=ds, method="ML"
           )

# 2
mod.2 <- lme(fixed = MATHACH ~ SES,
            random = ~ 1 | ID,
             data=ds, method="ML"
             )

# 3
mod.3 <- lme(fixed = MATHACH ~ SES + SECTOR,
            random = ~ 1 | ID,
             data=ds, method="ML"
             )

# 4
mod.4 <- lme(fixed = MATHACH ~ SES + SECTOR,
             random = ~ 1 + SES | ID,
             data=ds, method="ML"
            )

# 5
# todo: convergence error
# mod.5 <- lme(fixed = MATHACH ~ SES * SECTOR,
#             random = ~ 1 + SES | ID,
#             data=ds, method="ML"
#              )


# Extract information about the models
print(mod.1)      # basic information
summary(mod.1)    # complete information
summary(mod.2)
summary(mod.3)
summary(mod.4)
# summary(mod.5)

# attributes of an 'lmeObject'
mod.1$dims     # dimensions 
mod.1$coef     # coefficients
mod.1$call     # the formula used to produce the object
head(mod.1$residuals)     # residuals of the model
head(mod.1$varFix)   # matrix with the fixed effects
mod.1$sigma    # estimated within-group error standard deviation

anova(mod.1)


# Plots
plot(mod.1)
    # by default, the plot of an lme object is 
    # form=resid(., type = "p") ~ fitted(.)
    # that corresponds to the standardized residuals vs predicted values

plot(mod.1, resid(.) ~ fitted(.) | SECTOR)
    # idem, but comparing the sectors 

qqnorm(mod.1)     # normality of the residuals

interaction.plot(ds$SES, ds$SECTOR, ds$MATHACH, col=c("blue", "red"))
     # show one of the plotting abilities of R

# standardized residuals vs predicted values, by sector
plot(mod.3, resid(., type = "p") ~ fitted(.) | SECTOR, abline = 0)

# boxplots with the residuals by school 
#  (it's just an example, there are too many objects in level 2)
plot(mod.4, ID ~ resid(.))

```


## Beispiel ...

[Quelle](https://stackoverflow.com/questions/1380694/how-to-fit-a-random-effects-model-with-subject-as-random-in-r)

```{r}
dd = structure(list(Score = c(1.84, 2.24, 3.8, 2.3, 3.8, 4.55, 1.13, 
2.49, 3.74, 2.84, 3.3, 4.82, 1.74, 2.89, 3.39, 2.08, 3.99, 4.07, 
1.93, 2.39, 3.63, 2.55, 3.09, 4.76), Subject = c(1L, 1L, 1L, 
2L, 2L, 2L, 3L, 3L, 3L, 4L, 4L, 4L, 5L, 5L, 5L, 6L, 6L, 6L, 7L, 
7L, 7L, 8L, 8L, 8L), Condition = c(0L, 0L, 0L, 1L, 1L, 1L, 0L, 
0L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 1L, 
1L), Time = c(1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 
1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L)), .Names = c("Score", 
"Subject", "Condition", "Time"), class = "data.frame", row.names = c(NA, 
-24L))

# or using data.frame()
dd <- data.frame(Score = c(1.84, 2.24, 3.8, 2.3, 3.8, 4.55, 1.13, 2.49, 3.74, 2.84, 3.3, 4.82, 1.74, 2.89, 3.39, 2.08, 3.99, 4.07, 1.93, 2.39, 3.63, 2.55, 3.09, 4.76), Subject = c(1L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 4L, 4L, 4L, 4L, 4L, 4L), Condition = c(0L, 0L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 1L, 1L, 0L, 0L, 0L, 1L, 1L, 1L), Time = c(1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L))


# or more realistic
dd = structure(list(Score = c(1.62, 2.18, 2.3, 3.46, 3.85, 4.7, 1.41, 
2.21, 3.32, 2.73, 3.34, 3.27, 2.14, 2.73, 2.74, 3.39, 3.59, 4.01, 
1.81, 1.83, 3.22, 3.64, 3.51, 4.26), Subject = structure(c(1L, 
1L, 1L, 2L, 2L, 2L, 3L, 3L, 3L, 4L, 4L, 4L, 5L, 5L, 5L, 6L, 6L, 
6L, 7L, 7L, 7L, 8L, 8L, 8L), .Label = c("A", "B", "C", "D", "E", 
"F", "G", "H"), class = "factor"), Condition = structure(c(1L, 
1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 2L, 1L, 1L, 1L, 2L, 2L, 
2L, 1L, 1L, 1L, 2L, 2L, 2L), .Label = c("No", "Yes"), class = "factor"), 
    Time = structure(c(1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 
    2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L, 1L, 2L, 3L), .Label = c("1PM", 
    "2PM", "3PM"), class = "factor")), .Names = c("Score", "Subject", 
"Condition", "Time"), class = "data.frame", row.names = c(NA, 
-24L))


# a first plot
library(ggplot2)
qplot(Time, Score, data = dd, geom = "line", group = Subject, colour = factor(Condition))

# using mle

library(nlme)
m1 <- lme(Score ~ Condition + Time + Condition*Time, data = dd, random = ~ 1 | Subject)
summary(m1)

# The intercept variance is basically 0, indicating no within subject effect, so this model is not capturing the between time relationship well. A random intercept model is rarely the type of model you want for a repeated measures design. A random intercept model assumes that the correlations between all time points are equal. i.e. the correlation between time 1 and time 2 is the same as between time 1 and time 3. Under normal circumstances (perhaps not those generating your fake data) we would expect the later to be less than the former. An auto regressive structure is usually a better way to go.

### m2<-gls(Score ~ Condition + Time + Condition*Time, data = dd, correlation = corAR1(form = ~ Time | Subject))
### summary(m2)

# Your data is showing a -.596 between time point correlation, which seems odd. normally there should, at a minimum be a positive correlation between time points. How was this data generated?
# addendum:
# With your new data we know that the data generating process is equivalent to a random intercept model (though that is not the most realistic for a longitudinal study. The visualization shows that the effect of time seems to be fairly linear, so we should feel comfortable treating it as a numeric variable.

library(nlme)
m1 <- lme(Score ~ Condition + as.numeric(Time) + Condition*as.numeric(Time), data = dd, random = ~ 1 | Subject)
summary(m1)

# We see a significant Condition effect, indicating that the 'yes' condition tends to have higher scores (by about 1.7), and a significant time effect, indicating that both groups go up over time. Supporting the plot, we find no differential effect of time between the two groups (the interaction). i.e. the slopes are the same.


# (using lme4 library) This fits your subject effect as random and also the variable that your random effects are grouped under. In this model the random effect is the intercept varying by subject.
require(lme4)
m <- lmer( Score ~ Condition + Time + (1|Subject), data=dd )
# To see the random effects you can just use
ranef(m)
# As Ian Fellows mentioned, your data may also have random Condition and Time components. You can test that with another model. In the one below Condition, Time, and the intercept are allowed to vary randomly by subject. It also evaluates their correlations.
### mi <- lmer( Score ~ Condition + Time + (Condition + Time|Subject), data=dd )
### summary(mi)
### ranef(mi)
# You could also test for this without correlations with the intercept, with interactions between Condition and Time, and numerous other models to see which best fits your data and / or your theory. Your question is a bit vague but these few commands should get you started.
# Note that Subject is your grouping factor so it's what you fit other effects as random under. It's not something you then explicitly fit as a predictor as well.


```



##  Beispiel

[Quelle](https://stats.stackexchange.com/questions/242724/multi-level-data-vs-multi-level-models)


```{r}
require(tidyverse)
require(foreign)
dd <- read.dta("https://stats.idre.ucla.edu/stat/stata/examples/mlm_ma_hox/popular.dta")
write_csv(dd, "popular.csv")
# ... see code and explanation in answers

```


<!--

### Studienfächer, ev. Standorte oder Länder und ...

Benotungskulturen
Reichtum
Autoritarismus


### Ländervergleich und ...

Religion bzw. religiöse Minderheiten
Gleichberechtigung


### Beispiel Split Plot

todo: think on split plot in biology: 

Anbauflächen in verschiedenen Gebieten, jeweils unterteilt in verschiedene Parzellen und in den Parzellen wird experimentell manipuliert, z. B. gedüngt o.ä.

Die Anbauflächen haben bestimmte Charakteristika, die die Wirkung der experimentellen Manipulation systematisch beeinflussen können: magerer Boden, sumpfig, sauer o. ä.

-->





## ICC in surgery example


```{r}
library(tidyverse)
library(nlme)
dd <- read_delim("https://md.psych.bio.uni-goettingen.de/mv/data/div/cosmetic_surgery.dat", "\t")
# or
dd <- readr::read_tsv("https://md.psych.bio.uni-goettingen.de/mv/data/div/cosmetic_surgery.dat")

dd$Gender <- factor(dd$Gender, levels = c(0, 1), labels = c("f", "m"))
dd$Surgery_Text <- factor(dd$Surgery_Text)
dd$Reason_Text <- factor(dd$Reason_Text)
dd$Gender_Text <- factor(dd$Gender_Text)

# model without interaction 
m.test <- nlme::lme(Post_QoL ~ Surgery,
                      data = dd,
                      random = ~Surgery | Clinic,
                      method = "ML")

# we look at ICC
VarCorr(m.test)
var.clinic = as.numeric(VarCorr(m.test)[1])
var.surgery = as.numeric(VarCorr(m.test)[2])
var.residual = as.numeric(VarCorr(m.test)[3])
cat("ICC clinic: ", var.clinic/(var.clinic + var.surgery + var.residual))
cat("ICC surgery: ", var.surgery/(var.clinic + var.surgery + var.residual))

# package psychometric offers a ICC calculation
require(psychometric)
cat("ICC clinic by package psychometric: ", ICC1.lme(Post_QoL, grp=Clinic, data=dd))



m.interaction <- update(m.test, .~. + Surgery:Reason)
m.interaction <- nlme::lme(Post_QoL ~ Surgery + Surgery:Reason,
                      data = dd,
                      random = ~Surgery | Clinic,
                      method = "ML")

# we look at ICC, they only refer to random effects ...
VarCorr(m.interaction)
var.clinic = as.numeric(VarCorr(m.interaction)[1])
var.surgery = as.numeric(VarCorr(m.interaction)[2])
var.residual = as.numeric(VarCorr(m.interaction)[3])
cat("ICC clinic: ", var.clinic/(var.clinic + var.surgery + var.residual))
cat("ICC surgery: ", var.surgery/(var.clinic + var.surgery + var.residual))

cat("Random effect surgery binds slightly less variance compared to m.test.\n")
cat("because the fixed interaction effect Surgery:Reason 'explains' part of the surgery variance.")

# the additional fixed effect 


fixef(m.test)

head(ranef(m.test))

```



## Example

Matthew Crump
PSYC 7709: Using R for Reproducible Research
[source](https://crumplab.github.io/psyc7709/book/docs/a-tutorial-for-using-the-lme-function-from-the-nlme-package-.html)

```{r}

```



 
# Links und Referenzen

Psychologische Rundschau (2006): [Mehrebenen.pdf](Mehrebenen.pdf)

Torres [paper](www.princeton.edu/~otorres/Multilevel101.pdf)

Zur Frage nach p-values und lme4 oder generell Signifikanztests bei Mixed-Designs:
[see](http://glmm.wikidot.com/faq)

[blog](http://www.r-bloggers.com/fitting-a-model-by-maximum-likelihood/) auf R-Bloggers.

[Erklärung](https://www.unc.edu/courses/2006spring/ecol/145/001/docs/lectures/lecture17.htm)
zu Akaike's weight bzw. zu Modellvergleichen via AIC.

Siehe hierzu auch den 
[Artikel](ejwagenmakers.com/2004/aic.pdf) 
von Wagenmakers & Farrel (2004) 


## Zum Multi Level Teil: 

Multi level online [course](http://www.bristol.ac.uk/cmm/learning/multilevel-models/)
of university of Bristol.


## formulae

[Vergleich](http://rpsychologist.com/r-guide-longitudinal-lme-lmer)
lme4 nlme formulae und Erklärungen dazu.
Leider ohne die Daten.

Rense Nieuwenhuis

[nlme](http://www.rensenieuwenhuis.nl/r-sessions-21-multilevel-model-specification-nlme/)

[lme4](http://www.rensenieuwenhuis.nl/r-sessions-16-multilevel-model-specification-lme4/)

[Einführung](https://rpubs.com/favstats/multilevel) von Fabio Votta

Good [post](https://stats.stackexchange.com/questions/24452/how-to-interpret-variance-and-correlation-of-random-effects-in-a-mixed-effects-m)
that explains varCorr and the output of lmer



<!--

# Plotting

## Plotting random effects in lme

with package sjPlot [source](https://mran.microsoft.com/snapshot/2017-04-27/web/packages/sjPlot/vignettes/sjplmer.html)


```{r eval = FALSE}
require(tidyverse)
require(nlme)
require(sjPlot)
require(sjmisc)
dd <- read_delim("https://md.psych.bio.uni-goettingen.de/mv/data/cosmetic_surgery.dat", "\t")
dd$Gender <- factor(dd$Gender, levels = c(0, 1), labels = c("f", "m"))
dd$Surgery_Text <- factor(dd$Surgery_Text)
dd$Reason_Text <- factor(dd$Reason_Text)
dd$Gender_Text <- factor(dd$Gender_Text)

m_random <- nlme::lme(Post_QoL ~ Base_QoL + Age + Gender + BDI + Surgery,
                      data = dd,
                      random = ~Surgery | Clinic,
                      method = "ML")






plot_model(m_random,               #
         facet.grid = F,           #
         sort.est = "sort.all",    #
         y.offset = 0.3,           #
         p.kr=F)       

plot_model(m_random, 
         facet.grid = FALSE, 
         type="est",
         show.legend = T)

plot_model(m_random, 
         facet.grid = TRUE, 
         type="eff",
         show.legend = T)

plot_model(m_random, 
         facet.grid = FALSE, 
         type="slope",
         show.legend = T)


```


## Test comparison of VarCorr in nlme and lme4

```{r eval = FALSE}
require(tidyverse)
require(nlme)
require(lme4)

dd <- read.delim("http://md.psych.bio.uni-goettingen.de/mv/data/div/oxboys.txt", fileEncoding = "UTF-8")
m.1 <- nlme::lme(height ~ occasion, data = dd, random = ~ 1 + occasion | subj, method = "ML")
# m.1 <- nlme::lme(height ~ occasion, data = dd, random = ~ 1 + occasion | subj)
VarCorr(m.1)
# VarCorr(m.1)$subj # not working in nlme::lme objects

ml.0 <- lmer(height ~ occasion + (1 | subj), data = dd)
ml.1 <- lmer(height ~ occasion + (1 + occasion | subj), data = dd)
anova(ml.0, ml.1)

summary(ml.1)

vc <- VarCorr(ml.1)
VarCorr(ml.1)$subj

# cov occasion 
sqrt(vc)
```

-->




# Check

- We present `library(nlme)` here, alternatively we might use `library(lme4)` in combination with `lmerTest()`

- `nlme::gls()`  vs. `nlme::lme()` 
  non random effect models vs. random effect models

- level specific variables: values are the same for all observations of a specific level

- `nlme::gls(..., method = "ML")` and `nlme:lme(..., method = "ML")`
  useful for comparison of non MLM and MLM via `anova()`

- Formulae
  2 parts: fixed and random effects have their own formula
  random: | separats random effects from the definition of the level structure
  
- 3 or more levels, separated by "/", the more right runs faster, 
  `lme(... | country/city )`
  
- fixef() und ranef()
  fixed effects: coefficients are equal for all levels (random) 
  random effects: coefficients for every level additional to the fixed effect (corrects the fixed effect) 

- modeling random effects binds (explains) variance in the outcome and thus reduces the residuals due to level specific coefficients,
  but we don't get level specific significance tests
  
- including a variable as fixed effect reveals in a significance test for the coefficient (for all dummy variables in categorial predictors)

- variables that define levels have to be categorial (usually of type factor)

- it does not make sense to use level defining variables as explanatory variables (or even reaction variables)

- we can model dependencies using the `corr` parameter
  this is especially relevant for time series (growth models)
  [examlple above](#corr)

- we can use `intervals(model)` to get informations on explanatory variables in form of confidence intervals,
  a variable might be relevant, if 0 is not included in the interval

- we can use ICC to get an idea whether a categorial variable might serve as level defining variable

- Vergleich Äquivalenz WW und random effects (passiert differenziert im sheet zu MLM)

- predict
  ohne die originalen Daten, die den (Random-)Koeffizienten zugrunde liegen, basieren die Vorhersagen bei MLM nur auf den fixed Effects
  

Variablen in Beispiel unten: #Beispiel_Schule,_Abschlussnote_und_Schuleingangstest-Niveau

```{r check}
library(tidyverse)
library(nlme)
# data come from library(mlmRev)
dd <- mlmRev::Exam
# normexam	Normalized exam score.
# standLRT	Standardised LR test score (London Reading Test - a reading test taken when they were 11 years old)

# categorial predictor school
mg.0 <- gls(normexam ~ 1,      data=dd, method="ML")
mg.1 <- gls(normexam ~ school, data=dd, method="ML")
car::Anova(mg.1, type=3)

# additional continuous predictor reading test: standLRT
mg.2 <- gls(normexam ~ school + standLRT, data=dd, method="ML", correlation=NULL)
car::Anova(mg.2, type=3)
# summary(mg.2)

# null model
m.0 <- lme(fixed = normexam ~ 1,
       random = ~ 1 | school,
       data=dd,
       method="ML")
summary(m.0)
nlme::intervals(m.0)

anova(mg.0, m.0)

# MLM 
m.1 <- lme(fixed = normexam ~ standLRT,
       random = ~ 1 | school,
       data=dd,
       method="ML")
m.2 <- lme(fixed = normexam ~ standLRT,
       random = ~ standLRT | school,
       data=dd,
       method="ML")

# visualization {}todo
dd.s <- dd %>% dplyr::filter(school %in% c(1,9))
p_m.1 <- predict(m.1)
p_m.2 <- predict(m.2)

dd.s %>% ggplot(aes(x=standLRT, y=normexam, color = school)) +
  geom_point() + 
  # geom_point(data = dd.s, aes(color=school))
  # geom_line(data=predict(m.1)) +
  # geom_line(data=predict(m.2))
  ggtitle("School exit")

# fixed effects coefficients
fixef(m.2)
# random effects coefficients, head, because we have them for every second level observation (school)
head(ranef(m.2))

# confidence intervals
nlme::intervals(m.2)

# log. likelhood deviance
m.2$logLik
-2*m.2$logLik
deviance(m.2) # -2ll   only available with method="ML" in lme()

# a look at the variances
VarCorr(m.0)
VarCorr(m.1)
VarCorr(m.2)
# ... variances of residuals go down



# using ICC to decide on levels
# we look at variance decoposition: proportion that variance of interest holds of total variance
# school
as.numeric(VarCorr(m.0)['(Intercept)',1]) / sum(as.numeric(VarCorr(m.0)[,1])) # ICC schools as level 2 on random intercepts
# schgend (mixed, boys, girls)
mm <- lme(fixed = normexam ~ 1, random = ~ 1 | schgend, data=dd, method="ML")
as.numeric(VarCorr(mm)['(Intercept)',1]) / sum(as.numeric(VarCorr(m.0)[,1])) # ICC schgend (boys, girls, mixed) as level 2 on random intercepts
# vr (verbal reasoning bottom 25%, mid 50%, and top 25%.)
mm <- lme(fixed = normexam ~ 1, random = ~ 1 | vr, data=dd, method="ML")
as.numeric(VarCorr(mm)['(Intercept)',1]) / sum(as.numeric(VarCorr(m.0)[,1])) # ICC vr (schools level in verbal reasoning at intake) as level 2 on random intercepts
# sex (F or M)
mm <- lme(fixed = normexam ~ 1, random = ~ 1 | sex, data=dd, method="ML")
as.numeric(VarCorr(mm)['(Intercept)',1]) / sum(as.numeric(VarCorr(m.0)[,1])) # ICC pupils gender as level 2 on random intercepts

# require(psychometric)
# school on level 2
psychometric::ICC1.lme(normexam, grp=school, data=dd)
# ... about the same as above

# predict is possible
head(predict(m.2))
# ... and therefore some manual check of "explained" variance
cor(predict(m.2), dd$normexam)
cor(predict(m.2), dd$normexam)^2

# using library(rsq) which takes into account the model used ...
rsq::rsq(m.2, adj=F) # same as rsq::rsq(m.2)
rsq::rsq(m.2, adj=T)



# some conceptual stuff
# including a categorial explanatory variable to have a continuous (normexam) and a categorial (vr) explanatory variable
# vr: Student level Verbal Reasoning (VR) score band at intake - a factor. Levels are bottom 25%, mid 50%, and top 25%.
m.f <- lme(fixed = normexam ~ standLRT + vr,
       random = ~ 1 | school,
       data=dd,
       method="ML")
m.fr <- lme(fixed = normexam ~ standLRT + vr,
       random = ~ standLRT + vr | school,
       data=dd,
       method="ML")
# ... in m.fr we have slope for continuous standLRT for every school which corrects fixed slope for standLRT
# ... we also have an intercept for categorial vr which corrects the fixed intercept for every school. As vr has 3 levels, we deal with two dummy variables.

# fixed effects coefficients (parameters) are the same for all observations (level 1)
fixef(m.f)
fixef(m.fr)

# random effects coefficients (parameters) are valid for a single unit of a level, here for every school
# they "correct" the fixed effects
# categorial random variables correct the fixed intercept, correct the level
# continuous random variables correct the fixed slope for their unit
head(ranef(m.f))
head(ranef(m.fr))

# we have as many random coefficients as units in level 2 (schools) 
nrow(ranef(m.fr))

```

## check parms

Unter 
[https://md.psych.bio.uni-goettingen.de/mv/unit/parms/parms.html#Some_points_to_remember_regarding_MLM_and_an_example](https://md.psych.bio.uni-goettingen.de/mv/unit/parms/parms.html#Some_points_to_remember_regarding_MLM_and_an_example)
findet sich eine Sammlung von beachtenswerten Punkten zu MLM.



<!--
```{r}

# ICCs
as.numeric(VarCorr(m.0)['(Intercept)',1]) / sum(as.numeric(VarCorr(m.0)[,1])) # ICC schools as level 2 on random intercepts
as.numeric(VarCorr(m.1)['(Intercept)',1]) / sum(as.numeric(VarCorr(m.0)[,1])) # ICC schools as level 2 on random intercepts with standLRT as fixed effect
as.numeric(VarCorr(m.2)['(Intercept)',1]) / sum(as.numeric(VarCorr(m.0)[,1])) # ICC schools as level 2 on random intercepts with standLRT as fixed effect and as random effect
as.numeric(VarCorr(m.2)['standLRT',1]) / sum(as.numeric(VarCorr(m.0)[,1])) # ICC schools as level 2 on random intercepts with standLRT as fixed effect and as random effect
sum(as.numeric(VarCorr(m.2)[c('(Intercept)', 'standLRT') ,1])) / sum(as.numeric(VarCorr(m.0)[,1])) # ICC schools as level 2 on random intercepts with standLRT as fixed effect and as random effect
# BUT '(Intercept)' and 'standLRT' are correlated ...
# unclear what happens here ...
```
-->




### Screencast 

- MLM - check - data structure and the concept of level variables
  [StudIP](https://studip.uni-goettingen.de/plugins.php/mediacastplugin/media/check/557da707f35cbedf27871677d1e8e982/576b6797fb8a1ed374fe1bf3c17dc2a4/126?v=mlm-levels-de-2020-07-08.mp4) - [ownCloud](https://owncloud.gwdg.de/index.php/s/FJeTd2yAx7V2RJ8)

- MLM - check - running lme() and looking at the output, also model tests
  [StudIP](https://studip.uni-goettingen.de/plugins.php/mediacastplugin/media/check/557da707f35cbedf27871677d1e8e982/ae2459569f82479b776b3d7637f352b7/126?v=mlm-modeltest-de-2020-07-08.mp4) - [ownCloud](https://owncloud.gwdg.de/index.php/s/bE5x3Ec6LfIdfxU)
  
- MLM - check - ICC, it helps to decide on which levels should be entered in the model
  [StudIP](https://studip.uni-goettingen.de/plugins.php/mediacastplugin/media/check/557da707f35cbedf27871677d1e8e982/4d808151952705961606604848b9c16b/126?v=mlm-icc-de-2020-07-08.mp4) - [ownCloud](https://owncloud.gwdg.de/index.php/s/ta7E7nAAewBxztx)
  
- MLM - check - something more conceptual, fixef and ranef in action
  [StudIP](https://studip.uni-goettingen.de/plugins.php/mediacastplugin/media/check/557da707f35cbedf27871677d1e8e982/713eb937d7e2909bc690c9f272da5c8a/126?v=mlm-fixef-ranef-de-2020-07-0.mp4) - [ownCloud](https://owncloud.gwdg.de/index.php/s/PqyVB9tygpHIwmG)


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