

# Some points to remember regarding MLM and an example

[Rmd](essence_mlm.Rmd)


- we have fixed and random effects
- fixed effects work equal (fix) for all observations, random effects work specific due to level-group membership
- the way effects work and can be interpreted depend on the scale-level of the variables
- we concentrate on categorial and interval scaled variables here
- in R, interval scaled variables are of type `numeric()`, categorial scaled variables are usually of type `factor()`
- the estimated parameters for interval scaled variables are slopes (that define a line in combination with an intercept)
- for categorial explanatory variables the estimated coefficients are differences that are added to every observation in a subgroup.
- consequently in a sub-category the estimated outcome is equal for all members of this sub-category
- it is the arithmetic mean of the outcome values in this sub-group
- therefore we talk about "mean differences" that are tested
- categorial explanatory variables with more than two sub-groups are transferred to dummy variables (number of sub-groups - 1)
- these dummy variables replace the original explanatory variable in the model
- sub-group membership is expressed by a combination of 1 (member) or 0 (no member) of all the dummy variables
- we have various types of dummy coding. A common one is reference coding. In the reference category we find 0 in all dummy variables. 
  This is default in R's `lm()` and in others.
- we have one effect and a specific coefficient for each dummy variable. 
  For reference coding it expresses the difference of means in comparison with the reference category.
- fixed effects result in a coefficient for every explanatory variable with a significance test included.
  In the case of categorial explanatory variables this is true for every dummy variable 
  that is used in the model in place of the categorial variable itself.
  
- random effects are effects that are valid in a sub-group of observations only
- these sub-groups are defined by one ore more leveling variable. It's levels define sublevel membership of our observations.
- with more than one leveling variable the levels are defined by a sub-level combination
- in MLM we get coefficients that are specific for every sub-level
- random effect coefficients have *no* significance test
- random effects "bind" or "explain" residual variability not explained by fixed effects 
- we get an impression of the importance of random effects via model tests or via confidence intervals `confint()`
- random coefficients can be seen as correctors of fixed effect coefficients, they are added to predict an observations outcome
- interpretation of random coefficients is the same like for fixed coefficients. They also depend on scale level.
- in `nlme::lme` we can use `fixef()` and `ranef()` to output fixed and random coefficients resp.
- the dummy variable decomposition of categorial variables holds true for random effects also. 
  The mean differences they express are valid for the sub-goup they refer to.
- categorial random effects generate (sub-level * dummy variables) coefficients. These can be a lot in comparison to the number of observations.
- it doesn't make sense to include a level defining variable as a fixed effect

- fixed coefficients can be used to predict outcome for new observations. We can use `predict()` to do so
- random coefficients can only be used for prediction, if the new observations hold the same leveling variables as were used for model estimation.

- in `nlme::lme()` sub-levels are defined after the `|` in in the random part of our model formula
- in front of the `|` we ask for the random effects we want integrate in the model
- the syntax for the fixed part in the model formula is the same like in `lm()`
- we also have a `summary()` for `nlme::lme()` result objects

## the above in German

- Es gibt fixed und random Effects
- Diese Effekte wirken entweder für alle Beobachtungen (Vpn) gleich, sind also "fix", oder sie wirken spezifisch für bestimmte Subgruppen, aus Sicht von Level 1 (Grundlevel) also "random".
- Effekte und Interpretation von Parametern hängen vom Skalenniveau der erkärenden Variablen ab. 
- Relevant sind in unserem Kontext kategoriale oder intervallskalierte (stetige) Variablen.
- Intervallskalierte Variablen haben in R den Datentyp numeric(), kategoriale Variablen den Datentyp factor().
- Bei intervallskalierten erklärenden Variablen ist der anzupassende Parameter eine Steigung über ihren Wertebereich (in Kombination mit einem Intercept, um die Gerade zu fixieren)
- Bei kategorialen erklärenden Variablen sind die Parameter Werte, die bei allen Beobachtungen einer Teilkategorie zum Schätzwert dazugezählt oder abgezogen werden.
  Pro Teilkategorie gibt es also einen konstanten Schätzwert für alle dieser Teilkategorie angehörenden Beobachtungen. 
  Dies ist auch der Mittelwert der vorhergesagten Werte in dieser Teilkategorie.
  Daher spricht man von Mittelwertsunterschieden.
- Bei kategorialen erklärenden Variablen mit mehr als zwei Teilkategorien werden neue erklärende Variablen generiert und benutzt and Stelle der tatsächlichen Werte der Variable. 
  Sie werden Dummy-Variablen genannt.
- Dummy-Variablen sorgen dafür, dass die Zugehörigkeit zu einer Teilkategorie durch 0 (nicht zugehörig) und 1 (zugehörig) in Kombination mit den anderen Dummy-Variablen ausgedrückt wird
  Es gibt `Teilkategorien - 1` Dummy-Variablen für eine kategoriale erklärende Variable.
  Meist wird eine der Teilkategorien als Referenz-Kategorie benutzt (Referenzkodierung). Beobachtungen aus dieser Referenzgruppe haben in allen  Dummy-Variablen eine 0
- Für jede Dummy-Variable wird ein eigener Effekt geschätzt. Dies ist bei Referenzkodierung die Abweichung der Teilkategorie vom Mittelwert der Referenzkategorie.
- Bei fixed Effects wird für jede erklärende Variable ein Koeffizient (Parameter) geschätzt und statistisch geprüft.
  Wenn kategoriale Variablen im Spiel sind, gilt das für jede Dummy-Variable.

- random Effects sind Effekte, die nur in speziellen Teilschichten wirken.
  Teilschichten sind Subgruppen von Beobachtungen, die nach einer bestimmten (oder mehreren) kategorialen Variable gebildet werden.
  in MLM sollen diese Teilschichten nicht als kategoriale erklärende Variable in den fixed Effects geprüft werden.
  Es werden pro Teilschicht ebenfalls Effekte geschätzt, die Parameter werden aber nicht einzeln statistisch geprüft.
- random Effects sind Effekte, die die zugehörigen fixed Effects für die jeweilige Teilschicht korrigieren.
- statistisch gehen die random Effects dadurch in die Testung ein, dass sie Varianz erklären, die sonst Residualvarianz wäre, also Varianz, die durch das Modell nicht erklärt wird.
- auch bei random Effekts ist die Interpretation der Koeffizienten vergleichbar mit der bei fixed Effects und ist abhängig vom Skalenniveau der erklärenden Variablen
- random Effects bzw. deren Koeffizienten für jede Teilschicht kann man sich bei nlme::lme() mit dem Befehl ranef(...) ausgeben lassen.
- bei intervallskaierten, erklärenden Variablen sind die ausgegebenen Koeffizienten Korrekturen der Steigung des korrespondierenden fixed Effects
- bei kategorialen erklärenden Variablen werden ebenfalls entsprechende Dummy-Effekte generiert und die Koeffizienten pro Teilschicht ermittelt
- auch bei kategorialen erklärenden Variablen können die random Koeffizienten als Abweichung der jeweiligen Teilschicht von den analogen fixed Effects interpretiert werden, die für alle Beobachtungen gelten.

- fixed Koeffizienten können für die Vorhersage von neuen Beobachtungen benutzt werden, beispielsweise mit dem Befehl `predict(...)`
- random Effects können nur dann für eine Vorhersage von neuen Beobachtungen benutzt werden, wenn dieselbe Teilschicht auch bei der Parameterschätzung aufgetreten ist
- wenn eine Versuchsperson die Teilschicht definiert (meist bei Messwiederholungs-Designs ), können also die random Koeffizienten nur für Vorhersagen bei derselben Versuchsperson benutzt werden
  
- bei nlme::lme() erfolgt die Definition der Teilschichten hinter dem | im random-Teil der Formula
- vor dem | im random Teil stehen die random Effects, die für die Teilschichten ermittelt werden sollen
- der fixed Teil entspricht den Formulae bei den lm(), glm() etc. Befehlen

- summary(...) funktioniert auch für lme() Ergebnisobjekte
- die genauen Parameter können wir mit `fixef(...)` und `ranef(...)` anfordern


## The same effekt fixed or random, categorial explanatory variables

This is shown using a four level categorial explanatory variable `grp`. 
The same grouping variable is used as fixed effect in one model and as random effect in another model.
To count for covariation effects we include interaction in the fixed effect model.

```{r}
library(tidyverse)
dd <- readRDS(gzcon(url("http://md.psych.bio.uni-goettingen.de/mv/data/div/parms_dat.rds")))
#dd <- readRDS("parms_dat.rds")

# we begin simplifying things by reading in pre-post data only and consider it as beeing a between subjects design
dd <- dd %>% dplyr::filter(time %in% c("t0", "t1")) %>% dplyr::filter(grp %in% c("0cg", "0i1", "ecg", "ei1"))
dd$grp <- factor(dd$grp)
dd$time <- factor(dd$time)
dd$gender <- factor(dd$gender)

#psych::describe.by()
Rmisc::summarySE(dd, "rv", c("grp", "time"))

library(nlme)
# seeing it as pure fixed effect model, no repeated measure design
m.a0 <- nlme::gls(rv ~ 1, data=dd, method='ML')
m.a1 <- nlme::gls(rv ~        grp, data=dd, method='ML')
m.a2 <- nlme::gls(rv ~ time*grp, data=dd, method='ML')
anova(m.a0, m.a1, m.a2)
# we might be interested in the 
# summary(m.a1)
# but now more in the coefficients
m.a1$coefficients
# so we see: alphabetic order applies to grp levels, our reference subgroup is 0cg, the coefficients are relative to 0cg
# coefficients can be interpreted as differences between group means
m.a2$coefficients

# visualization of m.a1
cc <- scales::hue_pal()(length(unique(dd$grp)))
displ <- 0.2
dd %>% 
  ggplot(aes(x=grp, y=rv, group=grp, color=grp)) +
  #geom_point(position="jitter") + 
  geom_jitter(width = 0.2, alpha=0.3) +
  geom_smooth(method=lm, se=FALSE, size=0.2) +
  stat_summary(fun = mean, geom="point", size=3) +
  stat_summary(fun.data = mean_se, geom = "errorbar", width=0.2) + 
  # we add the effects based on coefficients
  geom_hline(aes(yintercept=m.a1$coefficients['(Intercept)']), color = cc[1], linetype="dashed") +
  geom_hline(aes(yintercept=m.a1$coefficients['(Intercept)'] + m.a1$coefficients['grp0i1']), color = cc[2], linetype="dashed") +
  geom_hline(aes(yintercept=m.a1$coefficients['(Intercept)'] + m.a1$coefficients['grpecg']), color = cc[3], linetype="dashed") +
  geom_hline(aes(yintercept=m.a1$coefficients['(Intercept)'] + m.a1$coefficients['grpei1']), color = cc[4], linetype="dashed") +
  geom_segment(aes(x = 1.8, y = m.a1$coefficients['(Intercept)'], xend = 2-displ, yend = m.a1$coefficients['(Intercept)'] + m.a1$coefficients['grp0i1']), color = cc[2], arrow = arrow(length = unit(0.3, "cm")), show.legend=F) +
  geom_segment(aes(x = 2.8, y = m.a1$coefficients['(Intercept)'], xend = 3-displ, yend = m.a1$coefficients['(Intercept)'] + m.a1$coefficients['grpecg']), color = cc[3], arrow = arrow(length = unit(0.3, "cm")), show.legend=F) +
  geom_segment(aes(x = 3.8, y = m.a1$coefficients['(Intercept)'], xend = 4-displ, yend = m.a1$coefficients['(Intercept)'] + m.a1$coefficients['grpei1']), color = cc[4], arrow = arrow(length = unit(0.3, "cm")), show.legend=F) +
  ggtitle('Visualization m.a1')

# visualization of m.a2
m.a2$coefficients
# cc <- scales::hue_pal()(length(m.a2$coefficients))
cc <- scales::hue_pal()(length(unique(dd$grp)))
dd %>% 
  ggplot(aes(x=time, y=rv, group=grp, color=grp)) +
  geom_jitter(width = 0.2, alpha=0.3) +
  geom_smooth(method=lm, se=FALSE, size=0.2) +
  stat_summary(fun = mean, geom="point", size=3) +
  stat_summary(fun.data = mean_se, geom = "errorbar", width=0.1) +
  stat_summary(fun = mean, geom="line") +
  geom_hline(aes(yintercept=m.a2$coefficients['(Intercept)']), color = cc[1], linetype="dashed") +
  geom_segment(aes(x = 1.8, y = m.a2$coefficients['(Intercept)'] + m.a2$coefficients['timet1'], xend = 2.2, yend = m.a2$coefficients['(Intercept)'] + m.a2$coefficients['timet1']), color = cc[1], linetype="dotted") +
  geom_segment(aes(x = 1.8, y = m.a2$coefficients['(Intercept)'] + m.a2$coefficients['timet1'] + m.a2$coefficients['grp0i1'] + m.a2$coefficients['timet1:grp0i1'], xend = 2.2, yend = m.a2$coefficients['(Intercept)'] + m.a2$coefficients['timet1'] +  m.a2$coefficients['grp0i1'] + m.a2$coefficients['timet1:grp0i1']), color = cc[2], linetype="dotted") +
  # geom_hline(aes(yintercept=m.a2$coefficients[1] + m.a2$coefficients[2]), color = cc[2], linetype="dashed") +
  #geom_hline(aes(yintercept=m.a2$coefficients[1] + m.a2$coefficients[3]), color = cc[3], linetype="dashed") +
  #geom_hline(aes(yintercept=m.a2$coefficients[1] + m.a2$coefficients[4]), color = cc[4], linetype="dashed") +
  #geom_hline(aes(yintercept=m.a2$coefficients[1] + m.a2$coefficients[5]), color = cc[5], linetype="dashed") +
  #
  # 0i1
  geom_segment(aes(x = 0.7, y = m.a2$coefficients['(Intercept)'] + m.a2$coefficients['grp0i1'], xend = 2.3, yend = m.a2$coefficients['(Intercept)'] + m.a2$coefficients['grp0i1']), color = cc[2], linetype="dotted") +
  geom_segment(aes(x = 0.9, y = m.a2$coefficients['(Intercept)'], xend = 0.9, yend = m.a2$coefficients['(Intercept)'] + m.a2$coefficients['grp0i1']), color = cc[2], arrow = arrow(length = unit(0.3, "cm")), show.legend=F) +
  geom_segment(aes(x = 2.1, y = m.a2$coefficients['(Intercept)'] + m.a2$coefficients['grp0i1'], xend = 2.1, yend = m.a2$coefficients['(Intercept)'] + m.a2$coefficients['grp0i1'] + m.a2$coefficients['timet1']), color = cc[2], arrow = arrow(length = unit(0.3, "cm")), show.legend=F) +
  geom_segment(aes(x = 2.2, y = m.a2$coefficients['(Intercept)'] + m.a2$coefficients['grp0i1'] + m.a2$coefficients['timet1'], xend = 2.2, yend = m.a2$coefficients['(Intercept)'] + m.a2$coefficients['grp0i1'] + m.a2$coefficients['timet1'] + m.a2$coefficients['timet1:grp0i1']), color = cc[2], arrow = arrow(length = unit(0.3, "cm")), show.legend=F) +
  # we show details for grp ei1 time t1
  geom_segment(aes(x = 0.7, y = m.a2$coefficients['(Intercept)'] + m.a2$coefficients['grpei1'], xend = 2.3, yend = m.a2$coefficients['(Intercept)'] + m.a2$coefficients['grpei1']), color = cc[4], linetype="dotted") +
  geom_segment(aes(x = 0.8, y = m.a2$coefficients['(Intercept)'], xend = 0.8, yend = m.a2$coefficients['(Intercept)'] + m.a2$coefficients['grpei1']), color = cc[4], arrow = arrow(length = unit(0.3, "cm")), show.legend=F) +
  # annotate(geom="text", x=0.7, y=0, label="grpei1", color=cc[4]) +
  geom_text(label="grpei1", x=0.9, y=0, angle=270, size=3, color=cc[4]) +
  geom_segment(aes(x = 2.1, y = m.a2$coefficients['(Intercept)'] + m.a2$coefficients['grpei1'], xend = 2.1, yend = m.a2$coefficients['(Intercept)'] + m.a2$coefficients['grpei1'] + m.a2$coefficients['timet1']), color = cc[4], arrow = arrow(length = unit(0.3, "cm")), show.legend=F) +
  geom_text(label="timet1", x=2, y=-0.30, angle=270, size=3, color=cc[4]) +
    geom_segment(aes(x = 2.2, y = m.a2$coefficients['(Intercept)'] + m.a2$coefficients['grpei1'] + m.a2$coefficients['timet1'], xend = 2.2, yend = m.a2$coefficients['(Intercept)'] + m.a2$coefficients['grpei1'] + m.a2$coefficients['timet1'] + m.a2$coefficients['timet1:grpei1']), color = cc[4], arrow = arrow(length = unit(0.3, "cm")), show.legend=F) +
  geom_text(label="timet1:grpei1", x=2.3, y=-0.20, angle=90, size=3, color=cc[4]) +
    #geom_segment(aes(x = 2.8, y = m.a2$coefficients[1], xend = 3-displ, yend = m.a2$coefficients[1] + m.a2$coefficients[3]), color = cc[3], arrow = arrow(length = unit(0.3, "cm")), show.legend=F) +
  #geom_segment(aes(x = 3.8, y = m.a2$coefficients[1], xend = 4-displ, yend = m.a2$coefficients[1] + m.a2$coefficients[4]), color = cc[4], arrow = arrow(length = unit(0.3, "cm")), show.legend=F) +
  ggtitle('Visualization m.a2')
m.a2$coefficients

# cc <- scales::hue_pal()(length(m1$coefficients) - 1)
# dd %>% 
#   ggplot(aes(x=time, y=rv, group=gender, color=gender)) +
#   geom_jitter(width = 0.2, alpha=0.3) +
#   geom_smooth(method=lm, se=FALSE, size=0.2) +
#   stat_summary(fun = mean, geom="point", size=3) +
#   stat_summary(fun.data = mean_se, geom = "errorbar", width=0.2) +
#   stat_summary(fun = mean, geom="line") +
#   geom_hline(aes(yintercept=m1$coefficients[1]), color = cc[1], linetype="dashed") +
#   geom_hline(aes(yintercept=m1$coefficients[1] + m1$coefficients[2]), color = cc[2], linetype="dashed") +
#   geom_hline(aes(yintercept=m1$coefficients[1] + m1$coefficients[3]), color = cc[3], linetype="dashed") +
#   ggtitle('Visualization m1')
# 


# how about considering grp as random effect?
m.r0 <- lme(fixed = rv ~ 1, random = ~1 | grp, data=dd, method='ML')
# m.0 already has a fixed effect Intercept, which can be seen as the grand mean of all observations
# m.0 random effects are per subgroup variations of the fixed effect
m.r0$coefficients
# m.2 includes time (pre post) as fixed effect
m.r1 <- lme(fixed = rv ~ time, random = ~1  | grp, data=dd, method='ML')  
m.r2 <- lme(fixed = rv ~ time, random = ~1 + time | grp, data=dd, method='ML')
m.r2$coefficients
anova(m.r0, m.r1, m.r2)

# visualization of m.r2
# we find coefficients for fixed effects in m.r2$coefficients[[1]]
# and coefficients for random effects in m.r2$coefficients[[2]]
#cc <- scales::hue_pal()(length(m.r2$coefficients))
cc <- scales::hue_pal()(length(unique(dd$grp)))
dd %>% 
  ggplot(aes(x=time, y=rv, group=grp, color=grp)) +
  geom_jitter(width = 0.2, alpha=0.3) +
  geom_smooth(method=lm, se=FALSE, size=0.2) +
  stat_summary(fun = mean, geom="point", size=3) +
  stat_summary(fun.data = mean_se, geom = "errorbar", width=0.1) +
  stat_summary(fun = mean, geom="line") +
  geom_hline(aes(yintercept=m.r2$coefficients[['fixed']]['(Intercept)']), color = "darkblue") +
  geom_text(label="fixed (Intercept)", x=0.7, y=m.r2$coefficients[['fixed']]['(Intercept)'] + 0.02, color="darkblue") +
  geom_hline(aes(yintercept=m.r2$coefficients[['fixed']]['(Intercept)'] + m.r2$coefficients[['fixed']]['timet1']), color = "darkblue", linetype="dashed") +
  geom_segment(aes(x = 1.9, y = m.r2$coefficients[['fixed']]['(Intercept)'], xend = 1.9, yend = m.r2$coefficients[['fixed']]['(Intercept)'] + m.r2$coefficients[['fixed']]['timet1']), color = "darkblue", arrow = arrow(length = unit(0.3, "cm")), size=1.5, show.legend=F) +
  geom_text(label="fixed timet1", x=1.9, y=m.r2$coefficients[['fixed']]['(Intercept)'] + m.r2$coefficients[['fixed']]['timet1'] + 0.1, angle=90, color="darkblue") +
  geom_segment(aes(x = 1.1, y = m.r2$coefficients[['fixed']]['(Intercept)'], xend = 1.1, yend = m.r2$coefficients[['fixed']]['(Intercept)'] + m.r2$coefficients[['random']][[1]]['0cg','(Intercept)']), color = cc[1], arrow = arrow(length = unit(0.1, "cm")), show.legend=F) +
  geom_segment(aes(x = 1.2, y = m.r2$coefficients[['fixed']]['(Intercept)'], xend = 1.2, yend = m.r2$coefficients[['fixed']]['(Intercept)'] + m.r2$coefficients[['random']][[1]]['0i1','(Intercept)']), color = cc[2], arrow = arrow(length = unit(0.1, "cm")), show.legend=F) +
  geom_segment(aes(x = 1.3, y = m.r2$coefficients[['fixed']]['(Intercept)'], xend = 1.3, yend = m.r2$coefficients[['fixed']]['(Intercept)'] + m.r2$coefficients[['random']][[1]]['ecg','(Intercept)']), color = cc[3], arrow = arrow(length = unit(0.1, "cm")), show.legend=F) +
  geom_segment(aes(x = 1.4, y = m.r2$coefficients[['fixed']]['(Intercept)'], xend = 1.4, yend = m.r2$coefficients[['fixed']]['(Intercept)'] + m.r2$coefficients[['random']][[1]]['ei1','(Intercept)']), color = cc[4], arrow = arrow(length = unit(0.1, "cm")), show.legend=F) +
  geom_segment(aes(x = 1, y = m.r2$coefficients[['fixed']]['(Intercept)'] + m.r2$coefficients[['random']][[1]]['ei1','(Intercept)'], 
                   xend = 1.4, yend = m.r2$coefficients[['fixed']]['(Intercept)'] + m.r2$coefficients[['random']][[1]]['ei1','(Intercept)']), color = cc[4], linetype="dotted", show.legend=F) +
  geom_text(label="random intercept ei1", x=1.4, y=-0.35, angle=270, size=3, color=cc[4]) +
  geom_segment(aes(x = 2.05, y = m.r2$coefficients[['fixed']]['(Intercept)'] + m.r2$coefficients[['fixed']]['timet1'], 
                   xend = 2.05, yend = m.r2$coefficients[['fixed']]['(Intercept)'] + m.r2$coefficients[['fixed']]['timet1']  + m.r2$coefficients[['random']][[1]]['0cg','(Intercept)']), 
                     color = cc[1], arrow = arrow(length = unit(0.1, "cm")), show.legend=F) +
  geom_segment(aes(x = 2.1, y = m.r2$coefficients[['fixed']]['(Intercept)'] + m.r2$coefficients[['fixed']]['timet1']+ m.r2$coefficients[['random']][[1]]['0cg','(Intercept)'], 
                   xend = 2.1, yend = m.r2$coefficients[['fixed']]['(Intercept)'] + m.r2$coefficients[['fixed']]['timet1']  + m.r2$coefficients[['random']][[1]]['0cg','(Intercept)'] + m.r2$coefficients[['random']][[1]]['0cg','timet1']), 
                     color = cc[1], arrow = arrow(length = unit(0.1, "cm")), show.legend=F) +
  # 0i1
  geom_segment(aes(x = 2.2, y = m.r2$coefficients[['fixed']]['(Intercept)'] + m.r2$coefficients[['fixed']]['timet1'], 
                   xend = 2.2, yend = m.r2$coefficients[['fixed']]['(Intercept)'] + m.r2$coefficients[['fixed']]['timet1']  + m.r2$coefficients[['random']][[1]]['0i1','(Intercept)']), 
                     color = cc[2], arrow = arrow(length = unit(0.1, "cm")), show.legend=F) +
  geom_text(label="random intercept 0i1", x=2.2, y=0.55, angle=90, size=3, color=cc[2]) +
  geom_segment(aes(x = 2.25, y = m.r2$coefficients[['fixed']]['(Intercept)'] + m.r2$coefficients[['fixed']]['timet1']+ m.r2$coefficients[['random']][[1]]['0i1','(Intercept)'], 
                   xend = 2.25, yend = m.r2$coefficients[['fixed']]['(Intercept)'] + m.r2$coefficients[['fixed']]['timet1']  + m.r2$coefficients[['random']][[1]]['0i1','(Intercept)'] + m.r2$coefficients[['random']][[1]]['0i1','timet1']), 
                     color = cc[2], arrow = arrow(length = unit(0.1, "cm")), show.legend=F) +
  geom_segment(aes(x = 2, y = m.r2$coefficients[['fixed']]['(Intercept)'] + m.r2$coefficients[['fixed']]['timet1']+ m.r2$coefficients[['random']][[1]]['0i1','(Intercept)'] + m.r2$coefficients[['random']][[1]]['0i1','timet1'], 
                   xend = 2.4, yend = m.r2$coefficients[['fixed']]['(Intercept)'] + m.r2$coefficients[['fixed']]['timet1']  + m.r2$coefficients[['random']][[1]]['0i1','(Intercept)'] + m.r2$coefficients[['random']][[1]]['0i1','timet1']), 
                     color = cc[2], linetype="dotted", show.legend=F) +  
  geom_text(label="random timet1 0i1", x=2.25, y=0.82, angle=90, size=3, color=cc[2]) +
  # ei1
  geom_segment(aes(x = 2.35, y = m.r2$coefficients[['fixed']]['(Intercept)'] + m.r2$coefficients[['fixed']]['timet1'], 
                   xend = 2.35, yend = m.r2$coefficients[['fixed']]['(Intercept)'] + m.r2$coefficients[['fixed']]['timet1']  + m.r2$coefficients[['random']][[1]]['ei1','(Intercept)']), 
                     color = cc[4], arrow = arrow(length = unit(0.1, "cm")), show.legend=F) +
  geom_text(label="random intercept ei1", x=2.35, y=-0.2, angle=270, size=3, color=cc[4]) +
  geom_segment(aes(x = 2.4, y = m.r2$coefficients[['fixed']]['(Intercept)'] + m.r2$coefficients[['fixed']]['timet1']+ m.r2$coefficients[['random']][[1]]['ei1','(Intercept)'], 
                   xend = 2.4, yend = m.r2$coefficients[['fixed']]['(Intercept)'] + m.r2$coefficients[['fixed']]['timet1']  + m.r2$coefficients[['random']][[1]]['ei1','(Intercept)'] + m.r2$coefficients[['random']][[1]]['ei1','timet1']), 
                     color = cc[4], arrow = arrow(length = unit(0.1, "cm")), show.legend=F) +
  geom_segment(aes(x = 2, y = m.r2$coefficients[['fixed']]['(Intercept)'] + m.r2$coefficients[['fixed']]['timet1']+ m.r2$coefficients[['random']][[1]]['ei1','(Intercept)'] + m.r2$coefficients[['random']][[1]]['ei1','timet1'], 
                   xend = 2.4, yend = m.r2$coefficients[['fixed']]['(Intercept)'] + m.r2$coefficients[['fixed']]['timet1']  + m.r2$coefficients[['random']][[1]]['ei1','(Intercept)'] + m.r2$coefficients[['random']][[1]]['ei1','timet1']), 
                     color = cc[4], linetype="dotted", show.legend=F) +
  geom_text(label="random timet1 ei1", x=2.4, y=0.1, angle=90, size=3, color=cc[4]) +
  ggtitle('Visualization m.r2')  

m.r2$coefficients


```


We can interprete the parameters  `m.a2$coefficients` of model `m.a2` as differences of group means.
We start at the mean of the reference group and add the effect (the coefficients) that are caused by membership of the respective sub-category.

With MLM and in model `m.r2` the fixed intercept represents the mean at t0 for all sub-groups.
Each sub-group (sub-level) has a random intercept which expresses the difference of the sub-group from the fixed intercept.
`t1` values are the result of fixed intercept + fixed intercept t1.
We then add the random effects for the respective sub-group, i. e. random intercept plus random t1.

<!--
Modell `m.a2` werden die Parameter `m.a2$coefficients` im Rahmen des ALM interpretiert als Abstände zwischen Mittelwerten. 
Wir gehen aus von der Referenzgruppe und addieren die Effekte, also die Unterschiede, die durch die Zugehörigkeit zur jeweiligen Subgruppen-Kombination entstehen.

Im MLM Modell `m.r2` entspricht der fixed Intercept dem Mittelwert zu t0 über alle Teilgruppen. 
Jede Teilgruppe hat dann einen random Intercept, der die Abweichung der jeweiligen Teilgruppe vom fixed Intercept ausdrückt.
Die Werte zu t1 ergeben sich aus dem fixed Intercept + fixed effect t1. 
Hierzu werden die beiden random effects für die jeweilige Teilgruppe addiert, also random Intercept und random t1 der jeweiligen Teilgruppe.
-->

### Conclusion

We finally get at the same means in model `m.a2` and `m.r2`.
So we express the same effects as fixed effects with interaction or as random effects.
With fixed effects, every coefficient has it`s significance test.


## MLM continuous explanatory variable and levels

A fixed effect for a continuous explanatory variable can have two parameters: an intercept and an inclination.
These parameters are the same across all levels.
If we put the same explanatory variable as random effect, we get per level corrections of the fixed effects that we have to add.


### One continuous xplanatory variable and its fixed and random effects with a categorial level variable

- mv - reaction variable
- ntime - explanatory variable, type numeric, fixed
- ntime - explanatory variable, type numeric, random
- grp - level-variable, four sublevels, type categorial

We visualize the parameters and how to use or interprete them.

```{r}
# dd.w <- read.delim("http://md.psych.bio.uni-goettingen.de/data/virt/v_stait_exam_wide.txt")
dd <- readRDS(gzcon(url("http://md.psych.bio.uni-goettingen.de/mv/data/div/parms_dat.rds")))
#dd <- readRDS("parms_dat.rds")

# we begin simplifying things by reading in pre-post data only and consider it as beeing a between subjects design
#dd <- dd %>% dplyr::filter(time %in% c("t0", "t1")) %>% dplyr::filter(grp %in% c("0cg", "0i1", "ecg", "ei1"))
dd <- dd %>% dplyr::filter(time %in% c("t1", "tf1", "tf2")) %>% dplyr::filter(grp %in% c("0cg", "0i1", "ecg", "ei1"))
dd$grp <- factor(dd$grp)
# dd$time <- factor(dd$time) # we consider time to be constant, this would imply at least equal measure intervals ...
dd$time <- factor(dd$time, levels=c("0", "t1", "tf1", "tf2"))
dd$ntime <- as.numeric(factor(dd$time))
dd$gender <- factor(dd$gender)

#psych::describe.by()
Rmisc::summarySE(dd, "rv", c("grp", "time"))

m.c0  <- lme(fixed = rv ~ 1, random = ~1 | grp, data=dd, method='ML')
m.c0a <- lme(fixed = rv ~ 1, random = ~1 | subj, data=dd, method='ML')

# m.0 already has a fixed effect Intercept, which can be seen as the grand mean of all observations
# m.0 random effects are per subgroup variations of the fixed effect
m.c0$coefficients
# m.2 includes time (pre post) as fixed effect
m.c1a  <- lme(fixed = rv ~ ntime, random = ~1  | grp, data=dd, method='ML')
ctrl <- lmeControl(opt='optim');
m.c1 <- lme(fixed = rv ~ ntime, random = ~1 + ntime | grp,      data=dd, control=ctrl, method='ML')
cc <- scales::hue_pal()(length(unique(dd$grp)))

# one fixed continuous explanatory variable, the same variable as random effect, one level
mm  <- lme(fixed = rv ~ ntime, random = ~1 + ntime | grp, data=dd, control=ctrl, method='ML')
mm$coefficients
dd %>% 
  ggplot(aes(x=time, y=rv, group=grp, color=grp)) + 
  # scale_x_discrete(limits=c('0','t1', 'tf1', 'tf2')) +
  #ggplot(aes(x=time, y=rv, group=grp, color=grp), scale_x_discrete(limits=c('0','t1', 'tf1', 'tf2'))) +  #xlim('0','t1', 'tf1', 'tf2')
  geom_jitter(width = 0.2, alpha=0.3) +
  geom_smooth(method=lm, se=FALSE, size=0.2) +
  stat_summary(fun = mean, geom="point", size=3, alpha=0.3) +
  stat_summary(fun.data = mean_se, geom = "errorbar", width=0.1, alpha=0.3) +
  stat_summary(fun = mean, geom="line", alpha=0.3) +
  # geom_abline(aes(intercept=a, slope=b), data=df3) +
  geom_segment(aes(x = 0, y = 0, xend = 0, yend = mm$coefficients$fixed['(Intercept)']), color = "darkblue", arrow = arrow(length = unit(0.3, "cm")), size=1.5, show.legend=F) +
  geom_abline(aes(x=ntime, intercept=mm$coefficients$fixed['(Intercept)'], slope=mm$coefficients$fixed['ntime']), color = "darkblue", size=1.5) +
  geom_text(label="f (Intercept), f slope ntime", x=0.8, y=mm$coefficients$fixed['(Intercept)'] + 0.02, angle=-5, color="darkblue") +
  geom_text(label='0,0', x=-0, y=0, color="black", show.legend=F) +
  geom_segment(aes(x=-0.3, y=0, xend=0, yend=0), alpha=0, show.legend=F) +
  # we visualize random effect of grp 0i1
  geom_abline(aes(x=ntime, intercept=mm$coefficients$fixed['(Intercept)'] + mm$coefficients$random$grp['0i1','(Intercept)'], slope=mm$coefficients$fixed['ntime'] + mm$coefficients$random$grp['0i1','ntime']), color = cc[2], size=1.5, linetype="dotted") +
  geom_segment(aes(x = 0, y = mm$coefficients$fixed['(Intercept)'], xend = 0, yend = mm$coefficients$fixed['(Intercept)'] + mm$coefficients$random$grp['0i1','(Intercept)']), color = cc[2], arrow = arrow(length = unit(0.3, "cm")), size=1.5, show.legend=F) +
  ggtitle('Visualization of, Parameters of mm')  
  

m.c2 <- mm
m.c2$coefficients

```


Blue line

- the blue line visualizes the fixed effect of `ntime`
- it is valid for all four sublevels
- it is composed on the fixed intercept `r mm$coefficients$fixed['(Intercept)']` and the fixed slope `r mm$coefficients$fixed['ntime']`
- this line with a slope of  `r mm$coefficients$fixed['ntime']` is fixed at x=0 and y=`r mm$coefficients$fixed['(Intercept)']`

Green arrow and green dotted line

- the vislualizes the effect of `ntime` for sublevel `grp == '0i1'`
- it is composed on base of fixed _and_ random parameters
- the slope is the sum of fixed effect slope  and random effect slope

- we can see random parameters as a way to correct the fixed effect to a specific sublevel (subgroup)

  

### One continuous xplanatory variable and its fixed and random effects with two nested categorial level variables

- mv - reaction variable
- ntime - explanatory variable, type numeric, fixed
- ntime - explanatory variable, type numeric, random
- grp - level-variable, four sublevels, type categorial, level 3
- subj - level-variable that's nested under grp, level 2

We visualize the parameters and show how to use or interprete them.

We show only one random group effect and the effects of one subject in this group, just to show the idea.

```{r eval=FALSE}
dd <- readRDS(gzcon(url("http://md.psych.bio.uni-goettingen.de/mv/data/div/parms_dat.rds")))

# we use only the data after the intervention, that is post treatment (t1) and two follow ups
# we also consider it a repeated measures design because we take care of the common influence in the data of one person
dd <- dd %>% dplyr::filter(time %in% c("t1", "tf1", "tf2")) %>% dplyr::filter(grp %in% c("0cg", "0i1", "ecg", "ei1"))
dd$grp <- factor(dd$grp)
# dd$time <- factor(dd$time, levels=c("0", "t1", "tf1", "tf2"))
dd$time <- factor(dd$time)
dd$ntime <- as.numeric(factor(dd$time))
dd$gender <- factor(dd$gender)
# we prepare colors
cc <- scales::hue_pal()(length(unique(dd$grp)))
# and the data of one subject
#ss <- '0i1_1'
ss <- '0i1_9'
gs <- '0i1/0i1_9'
s.data <- c(dd$rv[dd$subj==ss & dd$time=='t1'],
            dd$rv[dd$subj==ss & dd$time=='tf1'],
            dd$rv[dd$subj==ss & dd$time=='tf2']
            )
# Rmisc::summarySE(dd, "rv", c("grp", "time"))

# one fixed continuous explanatory variable, the same variable as random effect, one level
mm  <- lme(fixed = rv ~ ntime, random = ~1 + ntime | grp/subj, data=dd, control=ctrl, method='ML')
mm$coefficients$fixed
mm$coefficients$random$grp
head(mm$coefficients$random$subj)
dd %>% 
  ggplot(aes(x=time, y=rv, group=grp, color=grp)) + 
  # scale_x_discrete(limits=c('0','t1', 'tf1', 'tf2')) +
  #ggplot(aes(x=time, y=rv, group=grp, color=grp), scale_x_discrete(limits=c('0','t1', 'tf1', 'tf2'))) +  #xlim('0','t1', 'tf1', 'tf2')
  geom_jitter(width = 0.2, alpha=0.3) +
  geom_smooth(method=lm, se=FALSE, size=0.2) +
  stat_summary(fun = mean, geom="point", size=3, alpha=0.3) +
  stat_summary(fun.data = mean_se, geom = "errorbar", width=0.1, alpha=0.3) +
  stat_summary(fun = mean, geom="line", alpha=0.3) +
  # geom_abline(aes(intercept=a, slope=b), data=df3) +
  geom_segment(aes(x = 0, y = 0, xend = 0, yend = mm$coefficients$fixed['(Intercept)']), color = "darkblue", arrow = arrow(length = unit(0.3, "cm")), size=1.5, show.legend=F) +
  geom_abline(aes(x=ntime, intercept=mm$coefficients$fixed['(Intercept)'], slope=mm$coefficients$fixed['ntime']), color = "darkblue", size=1.5) +
  geom_text(label="f (Intercept), f slope ntime", x=0.8, y=mm$coefficients$fixed['(Intercept)'] + 0.02, angle=-5, color="darkblue") +
  geom_text(label='0,0', x=-0, y=0, color="black", show.legend=F) +
  geom_segment(aes(x=-0.3, y=0, xend=0, yend=0), alpha=0, show.legend=F) +
  # we visualize random effect of grp 0i1
  geom_abline(aes(x=ntime, intercept=mm$coefficients$fixed['(Intercept)'] + mm$coefficients$random$grp['0i1','(Intercept)'], slope=mm$coefficients$fixed['ntime'] + mm$coefficients$random$grp['0i1','ntime']), color = cc[2], size=1.5, linetype="dashed") +
  geom_segment(aes(x = 0, y = mm$coefficients$fixed['(Intercept)'], xend = 0, yend = mm$coefficients$fixed['(Intercept)'] + mm$coefficients$random$grp['0i1','(Intercept)']), color = cc[2], arrow = arrow(length = unit(0.3, "cm")), size=1.5, show.legend=F) +
  #geom_point(data=dd[subj==0i1_1], size=3, color=cc[2])
  # geom_point(stat = "identity", aes(y=s.data, x=c("t1", "tf1", "tf2")), size=3, color=cc[2]) +
  geom_point(stat = "identity", aes(y=s.data[1], x="t1"), size=3, color="lightgreen") +
  geom_point(stat = "identity", aes(y=s.data[2], x="tf1"), size=3, color="lightgreen") +
  geom_point(stat = "identity", aes(y=s.data[3], x="tf2"), size=3, color="lightgreen") +
  # geom_line(aes())
  # intercept for one subject
  geom_segment(aes(x = 0, 
                   y = mm$coefficients$fixed['(Intercept)'] + 
                     mm$coefficients$random$grp['0i1', '(Intercept)'],
                   xend = 0, 
                   yend = mm$coefficients$fixed['(Intercept)'] + 
                     mm$coefficients$random$grp['0i1','(Intercept)'] + 
                     mm$coefficients$random$subj[gs,'(Intercept)']), 
                   #color = cc[2], 
                color="lightgreen",
                arrow = arrow(length = unit(0.3, "cm")), 
                show.legend=F) +
  # regression line with slope for the same subject
  geom_abline(aes(x=ntime, 
                  intercept=mm$coefficients$fixed['(Intercept)'] + 
                    mm$coefficients$random$grp['0i1','(Intercept)'] + 
                    mm$coefficients$random$subj[gs,'(Intercept)'],
                  slope=mm$coefficients$fixed['ntime'] + 
                    mm$coefficients$random$grp['0i1','ntime'] +
                    mm$coefficients$random$subj[gs,'ntime']), 
                  # color = cc[2], 
                color="lightgreen",
                size=1.2,
                linetype="dashed") +
    ggtitle('Visualization of, Parameters of mm')  
  

m.c2 <- mm
m.c2$coefficients




















# two nested levels
mm  <- lme(fixed = rv ~ ntime, random = ~1 + ntime | grp/subj, data=dd, control=ctrl, method='ML')


# examples to access fixed effects using names
m.c2$coefficients$fixed['ev1']
m.c2$coefficients$fixed['(Intercept)']  # caution: brackets are part of the name, R names are case sensitive (capital I in Intercept)

# examples to access random effects by name
m.c2$coefficients$random$grp['0i1', 'timetf2']             #random effect of timetf2 for grp 0i1 of model m.c2
m.c2$coefficients$random$subj['ei1/ei1_9', 'timetf1']  # random effect of timetf1 for subject ei1_9 in grp ei1 of model m.c2


m.c2  <- lme(fixed = rv ~ time + ev1, random = ~1 + time | grp/subj, data=dd, control=ctrl, method='ML')
dd %>% 
  ggplot(aes(x=time, y=rv, group=grp, color=grp)) +
  geom_jitter(width = 0.2, alpha=0.3) +
  geom_smooth(method=lm, se=FALSE, size=0.2) +
  stat_summary(fun = mean, geom="point", size=3) +
  stat_summary(fun.data = mean_se, geom = "errorbar", width=0.1) +
  stat_summary(fun = mean, geom="line", alpha=0.3) +
  geom_hline(aes(yintercept=m.c2$coefficients$fixed['(Intercept)']), color = "darkblue") +
  geom_text(label="f (Intercept)", x=0.7, y=m.c2$coefficients$fixed['(Intercept)'] + 0.02, color="darkblue") +
  geom_segment(aes(x = 1.7, y=m.c2$coefficients$fixed['(Intercept)'] + m.c2$coefficients$fixed['timetf1'], xend=2.3, yend=m.c2$coefficients$fixed['(Intercept)'] + m.c2$coefficients$fixed['timetf1']), color = "darkblue", linetype="dashed") +
  geom_segment(aes(x = 1.9, y = m.c2$coefficients$fixed['(Intercept)'], xend = 1.9, yend = m.c2$coefficients$fixed['(Intercept)'] + m.c2$coefficients$fixed['timetf1']), color = "darkblue", arrow = arrow(length = unit(0.3, "cm")), size=1.5, show.legend=F) +
  geom_text(label="f timetf1", x=1.9, y=m.c2$coefficients$fixed['(Intercept)'] + m.c2$coefficients$fixed['timetf1'] + 0.1, angle=90, color="darkblue") +
  #geom_hline(aes(yintercept=m.c2$coefficients$fixed['(Intercept)'] + m.c2$coefficients$fixed['timetf2']), color = "darkblue", linetype="dashed") +
  geom_segment(aes(x = 2.7, y=m.c2$coefficients$fixed['(Intercept)'] + m.c2$coefficients$fixed['timetf2'], xend=3.3, yend=m.c2$coefficients$fixed['(Intercept)'] + m.c2$coefficients$fixed['timetf2']), color = "darkblue", linetype="dashed") +
  geom_segment(aes(x = 2.9, y = m.c2$coefficients$fixed['(Intercept)'], xend = 2.9, yend = m.c2$coefficients$fixed['(Intercept)'] + m.c2$coefficients$fixed['timetf2']), color = "darkblue", arrow = arrow(length = unit(0.3, "cm")), size=1.5, show.legend=F) +
  geom_text(label="f timetf2", x=2.9, y=m.c2$coefficients$fixed['(Intercept)'] + m.c2$coefficients$fixed['timetf2'] + 0.1, angle=90, color="darkblue") +
  ggtitle('Visualization m.c2')  
  
# m.c2$coefficients[['random']][['random$subj']]['ei1/ei1_9']
```

m.c2$coefficients

Only parameters are used to visualize the regression lines and effects ect. in the graph above.
This is the reason for the somewhat lengthy code.



## Visualization of the modelling of one single observation


```{r}
dd <- readRDS(gzcon(url("http://md.psych.bio.uni-goettingen.de/mv/data/div/parms_dat.rds")))
##readRDS("parms_dat.rds")
# we use only the data after the intervention, that is post treatment (t1) and two follow ups
# we also consider it a repeated measures design because we take care of the common influence in the data of one person
dd <- dd %>% dplyr::filter(time %in% c("t1", "tf1", "tf2", "tf3")) %>% dplyr::filter(grp %in% c("0cg", "0i1", "ecg", "ei1"))
dd$grp <- factor(dd$grp)
# dd$time <- factor(dd$time, levels=c("0", "t1", "tf1", "tf2"))
dd$time <- factor(dd$time)
dd$ntime <- as.numeric(factor(dd$time))
dd$gender <- factor(dd$gender)
# we prepare colors
cc <- scales::hue_pal()(length(unique(dd$grp)))
# and the data of one subject
#ss <- '0i1_1'
#ss <- '0i1_9'
ss <- '0i1_3'
gs <- '0i1/0i1_9'
s.data <- c(dd$rv[dd$subj==ss & dd$time=='t1'],
            dd$rv[dd$subj==ss & dd$time=='tf1'],
            dd$rv[dd$subj==ss & dd$time=='tf2'],
            dd$rv[dd$subj==ss & dd$time=='tf3']
            )
s.dd <- dd[dd$subj == ss,]
# Rmisc::summarySE(dd, "rv", c("grp", "time"))

ctrl <- lmeControl(opt='optim')
ctrl <- lmeControl(opt="nlminb")
ctrl <- lmeControl(maxIter = 500, msMaxIter = 500, tolerance = 1e-12, niterEM = 250,
          msMaxEval = 2000,
          msTol = 1e-14, msVerbose = FALSE,
          returnObject = FALSE, gradHess = TRUE, apVar = TRUE,
          .relStep = .Machine$double.eps^(1/3), minAbsParApVar = 0.005,
          #opt = "optim",
          optimMethod = "BFGS", natural = TRUE,
          sigma = NULL,
          allow.n.lt.q = FALSE
      )
##mm  <- lme(fixed = rv ~ ntime, random = ~1 + ntime | subj, data=dd, control=ctrl, method='ML')
##mm  <- lme(fixed = rv ~ ntime, random = ~1 | subj, data=dd, method='ML')

mm  <- lme(fixed = rv ~ ntime, random = ~1 + ntime | subj, control=list(opt="optim"), data=dd, method='ML')
mm$coefficients$fixed
##mm$coefficients$random$grp
##head(mm$coefficients$random$subj)
# random parameters for our selected subject
mm$coefficients$random$subj[ss,]
dd %>% 
  ggplot(aes(x=time, y=rv, group=grp, color=grp)) + 
  geom_jitter(width = 0.2, alpha=0.3) +
  geom_smooth(method=lm, se=FALSE, size=0.2) +
  stat_summary(fun = mean, geom="point", size=3, alpha=0.3) +
  stat_summary(fun.data = mean_se, geom = "errorbar", width=0.1, alpha=0.3) +
  stat_summary(fun = mean, geom="line", alpha=0.3) +
  geom_segment(aes(x = 0, y = 0, xend = 0, yend = mm$coefficients$fixed['(Intercept)']), color = "darkblue", arrow = arrow(length = unit(0.3, "cm")), size=1.5, show.legend=F) +
  geom_abline(aes(x=ntime, intercept=mm$coefficients$fixed['(Intercept)'], slope=mm$coefficients$fixed['ntime']), color = "darkblue", size=1.5) +
  geom_text(label='0,0', x=0, y=0, color="black", show.legend=F) +
  geom_segment(aes(x=-0.3, y=0, xend=0, yend=0), alpha=0, show.legend=F) +
  # we visualize random effect of grp 0i1
  #geom_abline(aes(x=ntime, intercept=mm$coefficients$fixed['(Intercept)'] + mm$coefficients$random$grp['0i1','(Intercept)'], slope=mm$coefficients$fixed['ntime'] + mm$coefficients$random$grp['0i1','ntime']), color = cc[2], size=1.5, linetype="dashed") +
  #geom_segment(aes(x = 0, y = mm$coefficients$fixed['(Intercept)'], xend = 0, yend = mm$coefficients$fixed['(Intercept)'] + mm$coefficients$random$grp['0i1','(Intercept)']), color = cc[2], arrow = arrow(length = unit(0.3, "cm")), size=1.5, show.legend=F) +
  #geom_point(data=dd[subj==0i1_1], size=3, color=cc[2])
  # geom_point(stat = "identity", aes(y=s.data, x=c("t1", "tf1", "tf2")), size=3, color=cc[2]) +
  geom_point(stat = "identity", aes(y=s.data[1], x="t1"), size=3, color="lightgreen") +
  geom_point(stat = "identity", aes(y=s.data[2], x="tf1"), size=3, color="lightgreen") +
  geom_point(stat = "identity", aes(y=s.data[3], x="tf2"), size=3, color="lightgreen") +
  geom_point(stat = "identity", aes(y=s.data[4], x="tf3"), size=3, color="lightgreen") +
  geom_point(data=s.dd, color="lightblue") +
  geom_smooth(data=s.dd, method=lm, se=FALSE, size=0.2, color="lightblue") +
  # geom_line(aes())
  # intercept for one subject
  geom_segment(aes(x = 0, 
                   y = mm$coefficients$fixed['(Intercept)'],
                     #+ mm$coefficients$random$grp['0i1', '(Intercept)'],
                   xend = 0, 
                   yend = mm$coefficients$fixed['(Intercept)'] + 
                     #mm$coefficients$random$grp['0i1','(Intercept)'] + 
                     mm$coefficients$random$subj[ss,'(Intercept)']), 
                   #color = cc[2], 
                color="lightgreen",
                arrow = arrow(length = unit(0.3, "cm")), 
                show.legend=F) +
  # regression line with slope for the same subject
  geom_abline(aes(x=ntime, 
                  intercept=mm$coefficients$fixed['(Intercept)'] + 
                    #mm$coefficients$random$grp['0i1','(Intercept)'] + 
                    mm$coefficients$random$subj[ss,'(Intercept)'],
                  slope=mm$coefficients$fixed['ntime'] + 
                    #mm$coefficients$random$grp['0i1','ntime'] +
                    mm$coefficients$random$subj[ss,'ntime']), 
                  # color = cc[2], 
                color="lightgreen",
                size=1.2,
                linetype="dashed") +
    ggtitle('Fixed time, random time, level subj')
# ... we see:
# blue arrow is fixed intercept
# blue line is fixed time slope
# light green arrow is random intercept for one chosen subject
# light green dashed line is random slope over time for selected subject 

# we add group membership as fixed effect

ss <- '0i1_4'
mm  <- lme(fixed = rv ~ ntime + grp, random = ~1 + ntime | subj, data=dd, control=list(opt="optim"), method='ML')
mm$coefficients$fixed
##mm$coefficients$random$grp
##head(mm$coefficients$random$subj)
# random parameters for our selected subject
mm$coefficients$random$subj[ss,]
dd %>% 
  ggplot(aes(x=time, y=rv, group=grp, color=grp)) + 
  geom_jitter(width = 0.2, alpha=0.3) +
  geom_smooth(method=lm, se=FALSE, size=0.2) +
  stat_summary(fun = mean, geom="point", size=3, alpha=0.3) +
  stat_summary(fun.data = mean_se, geom = "errorbar", width=0.1, alpha=0.3) +
  stat_summary(fun = mean, geom="line", alpha=0.3) +
  geom_segment(aes(x = 0, y = 0, xend = 0, yend = mm$coefficients$fixed['(Intercept)']), color = "darkblue", arrow = arrow(length = unit(0.3, "cm")), size=1.5, show.legend=F) +
  geom_abline(aes(x=ntime, intercept=mm$coefficients$fixed['(Intercept)'], slope=mm$coefficients$fixed['ntime']), color = "darkblue", size=1.5) +
  geom_text(label='0,0', x=0, y=0, color="black", show.legend=F) +
  geom_segment(aes(x=-0.3, y=0, xend=0, yend=0), alpha=0, show.legend=F) +
  # we visualize the fixed effect of grp 0i1
  geom_segment(aes(x = 0, y = mm$coefficients$fixed['(Intercept)'], xend = 0, yend = mm$coefficients$fixed['(Intercept)'] + mm$coefficients$fixed['grp0i1']), color = cc[2], arrow = arrow(length = unit(0.3, "cm")), size=1.5, show.legend=F) +
  geom_abline(aes(x=ntime, intercept=mm$coefficients$fixed['(Intercept)'] + mm$coefficients$fixed['grp0i1'], slope=mm$coefficients$fixed['ntime']), color = cc[2], size=1.5, linetype="dashed") +
  #geom_segment(aes(x = 0, y = mm$coefficients$fixed['(Intercept)'], xend = 0, yend = mm$coefficients$fixed['(Intercept)'] + mm$coefficients$random$grp['0i1','(Intercept)']), color = cc[2], arrow = arrow(length = unit(0.3, "cm")), size=1.5, show.legend=F) +
  #geom_point(data=dd[subj==0i1_1], size=3, color=cc[2])
  # geom_point(stat = "identity", aes(y=s.data, x=c("t1", "tf1", "tf2")), size=3, color=cc[2]) +
  geom_point(stat = "identity", aes(y=s.data[1], x="t1"), size=3, color="lightgreen") +
  geom_point(stat = "identity", aes(y=s.data[2], x="tf1"), size=3, color="lightgreen") +
  geom_point(stat = "identity", aes(y=s.data[3], x="tf2"), size=3, color="lightgreen") +
  geom_point(stat = "identity", aes(y=s.data[4], x="tf3"), size=3, color="lightgreen") +
  geom_point(data=s.dd, color="lightblue") +
  geom_smooth(data=s.dd, method=lm, se=FALSE, size=0.2, color="lightblue") +
  # geom_line(aes())
  # intercept for one subject
  geom_segment(aes(x = 0, 
                   y = mm$coefficients$fixed['(Intercept)'] + mm$coefficients$fixed['grp0i1'], 
                     #+ mm$coefficients$random$grp['0i1', '(Intercept)'],
                   xend = 0, 
                   yend = mm$coefficients$fixed['(Intercept)'] + mm$coefficients$fixed['grp0i1'] +
                     #mm$coefficients$random$grp['0i1','(Intercept)'] + 
                     mm$coefficients$random$subj[ss,'(Intercept)']), 
                   #color = cc[2], 
                color="lightgreen",
                arrow = arrow(length = unit(0.3, "cm")), 
                show.legend=F) +
  # regression line with slope for the same subject
  geom_abline(aes(x=ntime, 
                  intercept=mm$coefficients$fixed['(Intercept)'] + 
                    mm$coefficients$fixed['grp0i1'] +
                    #mm$coefficients$random$grp['0i1','(Intercept)'] + 
                    mm$coefficients$random$subj[ss,'(Intercept)'],
                  slope=mm$coefficients$fixed['ntime'] + 
                    #mm$coefficients$random$grp['0i1','ntime'] +
                    mm$coefficients$random$subj[ss,'ntime']), 
                  # color = cc[2], 
                color="lightgreen",
                size=1.2,
                linetype="dashed") +
    ggtitle('Fixed group, fixed time, random time, level subj')

# ... we see:
# blue arrow is fixed intercept
# blue line is fixed time slope
# olive arrow is fixed group intercept for selected group 0i1
# olive dashed line is fixed group slope for selected group 0i1
# light green arrow is random intercept for one chosen subject
# light green dashed line is random slope over time for selected subject 

# prediction
predict(mm, s.dd)
m.c2 <- mm
head(m.c2$coefficients)


```

The two graphs above want to show the composition of the prediction of a specific observation. 
The first graph is based on a model with fixed effect time and random effect time.
So the individual regression line (light green, dashed) is a "correction" of 
the dark blue regression line constructed on base of fixed parameters 
by the random parameters, so that finally results the light green, dashed regression line.
You can clearly see, that it is adapted to the emphasized individual data points of the selected observation (fat light green points).

The second graph shows an extension.
Here the resulting one observation regression line is composed of the fixed effect time,
of the fixed effect grp which results in an group specific line (olive green).
The selected observations regression line is build on base of the above mentioned two fixed effects plus it's own random coefficients.
The final individual line is the same, but built differently.

The `ggplot()` code is a bit lengthy to show, that all the visualizations can be built on base of coefficients, provided by our model `mm`.

