Mixed Effects model visualisation

######## Miminal working example


# Bin nicht sicher, ob Du alle packages hier brauchst.
pacman::p_load(tidyverse, nlme, ggplot2, ggthemes)


therapie <- read_csv("https://raw.githubusercontent.com/methodenlehre/data/master/therapie.csv") %>% 
  mutate(id = as.factor(id),
         bedingung = as.factor(bedingung))

therapie_gr1 <- therapie %>% 
  filter(bedingung == "Therapiegruppe")

therapie05_gr1 <- therapie_gr1 %>% 
  filter(woche <= 5) %>% 
  mutate(woche_f = as.factor(woche))
opt <- lmeControl(opt = "optim")

# takes a bit to fit
random.coefficients.05.c <- lme(wohlbefinden ~ woche_f, 
                                random = ~woche_f | id, 
                                data = therapie05_gr1, 
                                method = "ML", control = opt)



# if random.coefficients.05.c would be a MerMod object (from lme4), I would simply do this:
#therapie05_gr1$woche.wohlbefinden <- predict(random.coefficients.05.c, allow.new.levels =  TRUE, newdata = {therapie05_gr1 %>% mutate(id = "new")})
# as we used nlme, we have to do this by hand:
cfs <- fixef(random.coefficients.05.c)
cfs[2:length(cfs)] <- cfs[2:length(cfs)]+cfs[1]

therapie05_gr1$woche.wohlbefinden <- cfs[therapie05_gr1$woche_f]


ggplot(aes(x = woche, y = wohlbefinden, group = id),
       data = therapie05_gr1 %>% filter(woche <= 5)) +
  geom_line(aes(color = id), size = .2) +
  geom_jitter(aes(y = wohlbefinden, color = id),
              alpha = .4, width = 0.1, height = 0.1) +
  geom_line(aes(y = woche.wohlbefinden, group = 1), color = "black") +
  ggtitle("Kategoriales Random-Coefficients-Modell (nur Therapiegruppe)") +
  guides(color = FALSE) +  #Legends removed
  xlab("Woche") +
  xlab("Woche") +
  ylab("Wohlbefinden") +
  theme_tufte()