# Laden der relevanten Pakete library(MASS) # Für Negativ-binomiale Regression library(marginaleffects) # Modellvorhersagen und Vergleiche library(performance) # Modellbeurteilung (hier vor allem R2 und Dispersion) 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(base_size = 12)) # Voreinstellung für Plots # Lesen und Aufbereiten des Datensatz von Van Erkel & Van Aelst d <- haven::read_stata("Vanerkel_Vanaelst_2021.dta") |> rename( Political_knowledge = PK, Personalized_news = personalized_news, Radio = News_channels_w4_1, Television = News_channels_w4_2, Newspapers = News_channels_w4_3, Online_news_sites = News_channels_w4_4, Twitter = News_channels_w4_5, Facebook = News_channels_w4_6, ) |> mutate( Gender = as_factor(Gender), Education = as_factor(Education), Age10 = Age / 10 ) # Beschreibung einzelen Fragen d |> select(starts_with("PK")) |> mutate(across(everything(), ~ factor(., labels = c("incorrect", "correct")))) |> report_sample() # Beschreibung Wissensindex d |> ggplot(aes(Political_knowledge)) + geom_bar() # Daten aus Faehnrich et al. laden d2 <- read_rds("Faehnrich_2020.rds") # Auszug aus dem Datensatz von Faehnrich et al. - Zaehlvariablen d2 |> select(id, likes_count, comments_count, shares_count) |> slice_sample(n = 10) # Plot Lineare Regression mit PK3 - Simuliertes Beispiel set.seed(1) d_fake <- d |> mutate(PK3_fake = rbinom(n = n(), size = 1, prob = plogis(-10 + 2 * Age10))) plt_bin_fake <- d_fake |> ggplot(aes(Age10, PK3_fake)) + geom_point(position = position_jitter(width = 0, height = 0.03, seed = 1), shape = 1) + labs(x = "Alter in Jahrzehnten (Age / 10)", y = "Frage 3 (0 = falsch, 1 = richtig)") + scale_y_continuous(limits = c(-0.2, 1.2), breaks = (0:4) / 4) plt_bin_fake + geom_smooth(method = "lm", linewidth = 2, se = F) # Plot Logistische Regression mit PK3 - Simuliertes Beispiel plt_bin_fake + geom_smooth( method = "glm", method.args = list(family = "binomial"), se = F, linewidth = 2 ) # Plot Lineare Regression mit PK3 plt_bin <- d |> ggplot(aes(Age10, PK3)) + geom_point(position = position_jitter(width = 0, height = 0.03, seed = 1), shape = 1) + labs(x = "Alter in Jahrzehnten (Age / 10)", y = "Frage 3 (0 = falsch, 1 = richtig)") + scale_y_continuous(limits = c(-0.2, 1.2), breaks = (0:4) / 4) plt_bin + geom_smooth(method = "lm", linewidth = 2, se = F) # Plot Logistische Regression mit PK3 plt_bin + geom_smooth( method = "glm", method.args = list(family = "binomial"), se = F, linewidth = 2 ) # Null-Modell Logistische Regression mit PK3 m0 <- glm( PK3 ~ 1, # Null-Modell nur mit Intercept (Konstante) family = binomial(link = "logit"), # aV ist binär, logit-Link = Logistische Regression data = d ) # Null-Modell Logistische Regression mit PK3 Tabelle m0 |> report_table(metrics = "R2", include_effectsize = FALSE) # Invers logit Formel 1 / (1 + exp(0.38)) # Invers logit Funktion plogis(-0.38) # PK 3 deskriptiv d |> select(PK3) |> mutate(across(everything(), ~ factor(., labels = c("incorrect", "correct")))) |> report_sample() # Bivariate Logistische Regression mit PK3 und Alter d <- d |> mutate( Age10_c = Age10 - mean(Age10) ) m1 <- glm( PK3 ~ Age10_c, # Regressionsgleichung family = binomial(link = "logit"), # aV ist binär, logit-Link = Logistische Regression data = d ) # Bivariate Logistische Regression mit PK3 und Alter Tabelle m1 |> report_table(metrics = "R2", include_effectsize = FALSE) # Bivariate Logistische Regression mit PK3 und Alter Baseline Wahrscheinlichkeit plogis(-0.39) # Bivariate Logistische Regression mit PK3 und Alter Average counterfactual comparison M+1 m1 |> avg_comparisons( variables = list(Age10_c = 1), # Welcher Vergleich? Hier: Mittelwert und M + 1 type = "response" # In welcher Metrik? Hier: Wahrscheinlichkeit ) # Bivariate Logistische Regression mit PK3 und Alter Average counterfactual comparison M-1 m1 |> avg_comparisons( variables = list(Age10_c = -1), # Welcher Vergleich? Hier: Mittelwert und M - 1 type = "response" # In welcher Metrik? Hier: Wahrscheinlichkeit ) # Bivariate Logistische Regression mit PK3 und Alter Average counterfactual comparison Junge m1 |> avg_comparisons( variables = list(Age10_c = c(min(d$Age10_c), min(d$Age10_c) + 1)), # Welcher Vergleich? Hier: minimum + 10 Jahre type = "response" # In welcher Metrik? Hier: Wahrscheinlichkeit ) # Bivariate Logistische Regression mit PK3 und Alter Average counterfactual comparison Alte m1 |> avg_comparisons( variables = list(Age10_c = c(max(d$Age10_c) - 1, max(d$Age10_c))), # Welcher Vergleich? Hier: maximum - 10 Jahre type = "response" # In welcher Metrik? Hier: Wahrscheinlichkeit ) # Bivariate Logistische Regression mit PK3 und Alter Odds Ratios m1 |> report_table(metrics = "R2", include_effectsize = FALSE, exponentiate = TRUE) # Multiple Logistische Regression mit PK3 d <- d |> mutate( Online_news_sites_c = Online_news_sites - mean(Online_news_sites), Twitter_c = Twitter - mean(Twitter), Facebook_c = Facebook - mean(Facebook), Age10_c = Age10 - mean(Age10) ) m2 <- glm( PK3 ~ Online_news_sites_c + Twitter_c + Facebook_c + Gender + Age10_c, # Regressionsgleichung family = binomial(link = "logit"), # aV ist binär, logit-Link = Logistische Regression data = d ) # Multiple Logistische Regression mit PK3 Tabelle m2 |> report_table(metrics = "R2", include_effectsize = FALSE) # Multiple Logistische Regression mit PK3 Baseline Wahrscheinlichkeit plogis(0.13) # Multiple Logistische Regression mit PK3 Average counterfactual comparison M+1 m2 |> avg_comparisons( variables = list( Online_news_sites_c = 1, Twitter_c = 1, Facebook_c = 1, Gender = "reference", Age10_c = 1 ), # Welcher Vergleich? Hier: Mittelwert und M + 1 und Faktor type = "response" # In welcher Metrik? Hier: Wahrscheinlichkeit ) # Multiple Logistische Regression mit PK3 Odds Ratios m2 |> report_table(metrics = "R2", include_effectsize = FALSE, exponentiate = TRUE) # Multiple Logistische Regression mit PK3 R2 Tabelle m2 |> report_table(metrics = "R2", include_effectsize = FALSE) # Multiple Logistische Regression mit PK3 R2 Text m2 |> report_performance() |> cat() # Varianten der Modellspezifikation für Test mit 5 Fragen m5q_v1 <- glm( cbind(Political_knowledge, 5 - Political_knowledge) ~ Age10_c, # Zahl der richtig und falsch beantworteten Fragen in cbind() family = binomial(link = "logit"), # eine Frage ist binär, logit-Link = Logistische Regression data = d ) m5q_v2 <- glm( Political_knowledge / 5 ~ Age10_c, # Anteil der richtig beantworteten Fragen als aV family = binomial(link = "logit"), # eine Frage ist binär, logit-Link = Logistische Regression weights = rep(5, nrow(d)), # Zahl der Fragen als Gewicht data = d ) # Ausgabe Logistische Regression mit 5 Fragen m5q_v1 |> report_table(metrics = "R2", include_effectsize = FALSE) # Intercept Logistische Regression mit 5 Fragen round(plogis(coef(m5q_v1))[1], 2) # Plot Lineare Regression mit Political_knowledge d |> ggplot(aes(Age10, Political_knowledge)) + geom_point(position = position_jitter(width = 0, height = 0.1, seed = 1), shape = 1) + labs(x = "Alter in Jahrzehnten (Age / 10)", y = "Richtig beantwortete Fragen") + geom_smooth(method = "lm", se = F, linewidth = 2) # Plot Logistische Regression mit Political_knowledge d |> ggplot(aes(Age10, Political_knowledge / 5, weight = 5)) + geom_point(position = position_jitter(width = 0, height = 0.02, seed = 1), shape = 1) + labs(x = "Alter in Jahrzehnten (Age / 10)", y = "Richtig beantwortete Fragen") + geom_smooth( method = "glm", method.args = list(family = "binomial"), se = F, linewidth = 2 ) # Mittelwert und Varianz der Beispiel-Variable d2_RU <- d2 |> mutate( ym = floor_date(ymd_hms(created_time), "month"), ym_int = as.integer(as.factor(ym)) - min(as.integer(as.factor(ym)), na.rm = TRUE) ) |> filter(uni == "UC San Francisco") |> na.omit() d2_RU |> summarise( M = mean(comments_count), Var = var(comments_count), n = n() ) # Plot Lineare Regression mit comments_count - Simuliertes Beispiel d_fake <- d2_RU |> mutate( lambda = exp(2 + 0.08 * ym_int + rnorm(n = n(), mean = 0, sd = 0.3)), comments_count_fake = rpois(n = n(), lambda = lambda) ) plt_count_fake <- d_fake |> ggplot(aes(ym, comments_count_fake)) + geom_point(position = position_jitter(width = 0, height = 0.1, seed = 1), shape = 1) + labs(x = "Zeit (in Monatsschritten)", y = "Zahl der Kommentare") plt_count_fake + geom_smooth(method = "lm", linewidth = 2, se = F) # Plot Poisson- und negativ-binomiale Regression mit comments_count - Simuliertes Beispiel plt_count_fake + geom_smooth( method = "glm", method.args = list(family = "poisson"), se = F, linewidth = 1, color = "green", linetype = 2 ) + geom_smooth( method = MASS::glm.nb, se = F, linewidth = 1 ) # Plot Lineare Regression mit comments_count plt_count <- d2_RU |> ggplot(aes(ym, comments_count)) + geom_point(position = position_jitter(width = 0, height = 0.1, seed = 1), shape = 1) + labs(x = "Zeit (in Monatsschritten)", y = "Zahl der Kommentare") plt_count + geom_smooth(method = "lm", linewidth = 2, se = F) # Plot Negativ-binomiale Regression mit comments_count plt_count + geom_smooth( method = MASS::glm.nb, se = F, linewidth = 2 ) # Null-Modell Negativ-binomiale Regression mit comments_count nbm0 <- glm.nb( # Negativ-binomiale Regression comments_count ~ 1, # Null-Modell nur mit Intercept (Konstante) data = d2_RU ) # Null-Modell Negativ-binomiale Regression mit comments_count Tabelle nbm0 |> report_table(metrics = "R2", include_effectsize = FALSE) # exp Funktion exp(coef(nbm0)) # comments_count deskriptiv mean(d2_RU$comments_count) # Bivariate negativ-binomiale Regression mit comments_count und Zeit d2_RU <- d2_RU |> mutate(ym_int_half = ym_int - (max(ym_int) + 1) / 2) nbm1 <- glm.nb( # Negativ-binomiale Regression comments_count ~ ym_int_half, # Regressionsgleichung data = d2_RU ) # Bivariate negativ-binomiale Regression mit comments_count und Zeit Tabelle nbm1 |> report_table(exponentiate = TRUE, metrics = "R2", include_effectsize = FALSE) # Bivariate negativ-binomiale Regression comments_count und Zeit Average counterfactual comparison M+1 nbm1 |> avg_comparisons( variables = list(ym_int_half = c(0, 1)), # Welcher Vergleich? Hier: Mitte der Zeit und M + 1 Monat type = "response" # In welcher Metrik? Hier: Erwartete Zahl der Kommentare ) # Bivariate negativ-binomiale Regression mit comments_count und Zeit Average counterfactual comparison M-1 nbm1 |> avg_comparisons( variables = list(ym_int_half = c(0, -1)), # Welcher Vergleich? Hier: Mitte der Zeit und M - 1 Monat type = "response" # In welcher Metrik? Hier: Erwartete Zahl der Kommentare ) # Bivariate negativ-binomiale Regression mit comments_count und Zeit Average counterfactual comparison Junge nbm1 |> avg_comparisons( variables = list(ym_int_half = c(min(d2_RU$ym_int_half), min(d2_RU$ym_int_half) + 1)), # Welcher Vergleich? Hier: minimum + 1 Monat type = "response" # In welcher Metrik? Hier: Erwartete Zahl der Kommentare ) # Bivariate negativ-binomiale Regression mit comments_count und Zeit Average counterfactual comparison Alte nbm1 |> avg_comparisons( variables = list(ym_int_half = c(max(d2_RU$ym_int_half) - 1, max(d2_RU$ym_int_half))), # Welcher Vergleich? Hier: maximum - 1 Monat type = "response" # In welcher Metrik? Hier: Erwartete Zahl der Kommentare ) # Multiple negativ-binomiale Regression mit comments_count d2_RU <- d2_RU |> mutate( ym_int_half = ym_int - (max(ym_int) + 1) / 2, word_count_100_c = (word_count - mean(word_count, rm.na = TRUE)) / 100 ) nbm2 <- glm.nb( comments_count ~ ym_int_half + topic_research + topic_teaching + word_count_100_c, # Regressionsgleichung data = d2_RU ) # Multiple negativ-binomiale Regression mit comments_count Tabelle nbm2 |> report_table(exponentiate = TRUE, metrics = "R2", include_effectsize = FALSE) # Multiple negativ-binomiale Regression mit comments_count Average counterfactual comparison nbm2 |> avg_comparisons( variables = list( ym_int_half = c(0, 1), topic_research = "reference", topic_teaching = "reference", word_count_100_c = 1 ), # Welche Vergleiche? type = "response" # In welcher Metrik? Hier: Erwartete Zahl der Likes ) # Multiple negativ-binomiale Regression mit comments_count R2 Tabelle nbm2 |> report_table(exponentiate = TRUE, metrics = "R2", include_effectsize = FALSE) # Multiple negativ-binomiale Regression mit comments_count R2 Text nbm2 |> report_performance() |> cat() # Poisson-Regression pm2 <- glm( comments_count ~ ym_int_half + topic_research + topic_teaching + word_count_100_c, # Regressionsgleichung family = poisson(link = "log"), # aV Zählvariable, log-Link = Poisson-Regression data = d2_RU ) # Modellvergleich Koeffizienten compare_models(Poisson = pm2, "Negativ-binomial" = nbm2, exponentiate = TRUE) # Test der einzelnen Modelle: Poisson-Regression check_overdispersion(pm2) # Test der einzelnen Modelle: Negativ-binomiale Regression check_overdispersion(nbm2) # Modellvergleich Informationkriterien compare_performance(pm2, nbm2, metrics = c("AIC", "BIC")) # Modellvergleich Test library(lmtest) lrtest(pm2, nbm2)