# Laden der relevanten Pakete
library(lme4) # Mehrebenenmodelle
library(lmerTest) # Signifikanztests für Mehrebenenmodelle
library(marginaleffects) # Modellvorhersagen und Vergleiche
library(performance) # Modellbeurteilung (hier vor allem ICC)
library(parameters) # Modellschätzungen extrahieren und darstellen
library(report) # Einfaches Erstellen von statistischen Berichten
library(tidyverse) # Datenmanagement und Visualisierung: https://www.tidyverse.org/
theme_set(theme_classic()) # Voreinstellung für Plots


# Daten laden
d <- read_rds("Faehnrich_2020.rds")


# Struktur der Stichprobe
d |>
  summarise(n_posts = n(), .by = uni)


# Stichprobe: Unabhängige Variablen (L1)
sample_l1 <- d |>
  select(timeofday, type, year, starts_with("topic"), word_count, day_type) |>
  report_sample()
sample_l1 |>
  _[1:12, ]


# Stichprobe: Unabhängige Variablen (L1) - Fortsetzung
sample_l1 |>
  _[13:nrow(sample_l1), ]


# Stichprobe: Abhängige Variablen (L1)
d |>
  select(likes_count, comments_count, shares_count) |>
  report_table()


# Stichprobe: Unabhängige Variablen (L2)
d |>
  distinct(uni, uni_fans, uni_us) |>
  select(-uni) |>
  report_sample()


# Stichprobe: log-Transformation - Code
# d <- d |>
#   mutate(
#     word_count_log = log1p(word_count),
#     uni_fans_log = log(uni_fans),
#     likes_count_log = log1p(likes_count),
#     comments_count_log = log1p(comments_count),
#     shares_count_log = log1p(shares_count)
#   )


# Stichprobe: log-Transformation - Grafik
d |>
  select(id, uni_fans, ends_with("count"), ends_with("log")) |>
  gather(key, value, -id) |>
  mutate(
    type = if_else(str_detect(key, "log"), "log", "original"),
    key = str_remove_all(key, "_log")
  ) |>
  spread(type, value) |>
  ggplot(aes(original, log)) +
  geom_point() +
  facet_wrap(vars(key), scales = "free") +
  scale_x_continuous(labels = scales::label_comma())


# Stichprobe: Datensatz
d |>
  select(comments_count_log, topic_research, uni, uni_fans_log) |>
  slice_head(n = 4, by = uni) |>
  _[1:10, ]


# Null-Modell schaetzen
m0 <- lmer( # Funktion für lineares Mehrebenenmodell
  comments_count_log ~ # abhängige Variable
    1 + # Gewichteter Gesamtmittelwert
    (1 | uni), # Struktur des Modells
  data = d
) # Daten


# Null-Modell Tabelle
m0 |>
  report_table(metrics = "R2", include_effectsize = FALSE)


# Null-Modell Vorhersage
plot_predictions(m0, condition = "uni") +
  coord_flip() +
  scale_x_discrete(limits = rev)


# ICC
icc(m0, by_group = TRUE)




# Random-Intercept-Modell mit L1-Praediktoren schaetzen
m1 <- lmer( # Funktion für lineares Mehrebenenmodell
  comments_count_log ~ # abhängige Variable
    topic_interact + topic_research + topic_teaching + # Prädiktoren
    (1 | uni), # Struktur des Modells
  data = d
)


# Random-Intercept-Modell mit L1-Praediktoren Tabelle
m1 |>
  report_table(metrics = "R2", include_effectsize = FALSE)


# Random-Intercept-Modell mit L1-Praediktoren Plot
m1 |>
  parameters(effects = "fixed") |>
  plot()


# Random-Intercept-Modell mit L2-Praediktoren schaetzen
m2 <- lmer(comments_count_log ~ topic_interact + topic_research +
  topic_teaching + uni_fans_log + (1 | uni), data = d)


# Random-Intercept-Modell mit L2-Praediktoren Tabelle
m2 |>
  report_table(metrics = "R2", include_effectsize = FALSE)


# Random-Intercept-Modell mit L2-Praediktoren Plot
m2 |>
  parameters(effects = "fixed") |>
  plot()


# Random-Slope-Modell schaetzen
m3 <- lmer(
  comments_count_log ~ topic_interact + topic_research +
    topic_teaching + uni_fans_log +
    (topic_research | uni), # Koeffizient von topic_research darf zwischen Uni-Seiten variieren
  data = d
)


# Random-Slope-Modell Modellvergleich
anova(m2, m3, model.names = c("Random Intercept", "Random Slope")) |>
  rownames_to_column(var = "Modell")


# Random-Slope-Modell Tabelle der Fixed Effects
m3 |>
  report_table(metrics = "R2", include_effectsize = FALSE) |>
  filter(Effects == "fixed") |>
  select(-Effects, -Group, -Fit)


# Random-Slope-Modell Tabelle der Random Effects
m3 |>
  parameters(effects = "random", ci_random = FALSE, group_level = FALSE)


# Random-Slope-Modell Plot der Random Effects
m3 |>
  plot_slopes(variables = "topic_research", condition = "uni") +
  coord_flip() +
  scale_x_discrete(limits = rev) +
  geom_hline(yintercept = 0, linetype = 2)


# Random-Slope-Modell mit Cross-Level-Interaktion schaetzen
m4 <- lmer(comments_count_log ~ topic_interact +
  topic_research * uni_fans_log + # Cross-Level-Interaktion
  topic_teaching + (topic_research | uni), data = d)


# Random-Slope-Modell mit Cross-Level-Interaktion Tabelle der Fixed Effects
m4 |>
  report_table(metrics = "R2", include_effectsize = FALSE) |>
  filter(Effects == "fixed") |>
  select(-Effects, -Group, -Fit)


# Random-Slope-Modell mit Cross-Level-Interaktion Plot der Interaktion
plot_slopes(m4, variables = "topic_research", condition = "uni_fans_log") +
  geom_hline(yintercept = 0, linetype = 2)