Relationship between engagement and news quality

results_model1 <- list()
results_model2 <- list()

for (platform in platforms) {
  data_subset <- ALL_Platforms_eng[ALL_Platforms_eng$platform == platform, ]
  
  # Model 1: pc1 only
  model1 <- feglm(scale(log10(engagement + 1)) ~ scale(pc1) | username,cluster = "username",data = data_subset)
  coef1 <- coef(model1)
  conf1 <- confint(model1)
  
  results_model1[[platform]] <- data.frame(
    platform = platform,
    variable = "scale(pc1)",
    estimate = coef1["scale(pc1)"],
    ci_lb = conf1["scale(pc1)", 1],
    ci_ub = conf1["scale(pc1)", 2],
    color = "black",
    shape = 15)
  
  # Model 2: pc1 + tox + certainty
  model2 <- feglm(scale(log10(engagement + 1)) ~ scale(pc1) + scale(toxicity_ntile) + scale(certainty_avg) | username, cluster = "username",data = data_subset)
  coef2 <- coef(model2)
  conf2 <- confint(model2)
  
  results_model2[[platform]] <- data.frame(
    platform = platform,
    variable = c("scale(pc1)", "scale(toxicity_ntile)", "scale(certainty_avg)"),
    estimate = coef2[c("scale(pc1)", "scale(toxicity_ntile)", "scale(certainty_avg)")],
    ci_lb = conf2[c("scale(pc1)", "scale(toxicity_ntile)", "scale(certainty_avg)"), 1],
    ci_ub = conf2[c("scale(pc1)", "scale(toxicity_ntile)", "scale(certainty_avg)"), 2],
    color = c("darkred", "darkgreen", "darkblue"),
    shape = 16
  )
}

# Combine results
results_df_model1 <- do.call(rbind, results_model1)
results_df_model2 <- do.call(rbind, results_model2)

# Variance for meta-analysis
results_df_model1$vi <- ((results_df_model1$ci_ub - results_df_model1$ci_lb) / (2 * 1.96))^2
results_df_model2$vi <- ((results_df_model2$ci_ub - results_df_model2$ci_lb) / (2 * 1.96))^2

# Meta-analysis for pc1 only
res_model1_pc1 <- rma(estimate, vi, data = results_df_model1)

# Meta-analysis for pc1 in extended model
res_model2_pc1 <- rma(estimate, vi, data = results_df_model2, subset = (variable == "scale(pc1)"))
res_model2_toxicity <- rma(estimate, vi, data = results_df_model2, subset = (variable == "scale(toxicity_ntile)"))
res_model2_certainty <- rma(estimate, vi, data = results_df_model2, subset = (variable == "scale(certainty_avg)"))

# Order platforms by pc1 effect size from Model 1
order_df <- results_df_model1[order(results_df_model1$estimate), ]
sorted_platforms <- order_df$platform
results_df_model1 <- results_df_model1[order(factor(results_df_model1$platform, levels = sorted_platforms)), ]
results_df_model2 <- results_df_model2[order(factor(results_df_model2$platform, levels = sorted_platforms)), ]

# Plot
par(mfrow = c(1, 2))

# Model 1 forest plot
forest(
  results_df_model1$estimate,
  ci.lb = results_df_model1$ci_lb,
  ci.ub = results_df_model1$ci_ub,
  slab = results_df_model1$platform,
  xlab = "Effect Size",
  refline = 0,
  shade = TRUE,
  col = results_df_model1$color,
  ylim = c(-3, nrow(results_df_model1) + 4),
  digits = 3
)
addpoly(res_model1_pc1, col = "black", cex = 1, border = "black", lwd = 2, row = -1, mlab = 'pc1')
title(main = "Model 1: Only pc1")

# Model 2 forest plot
forest(
  results_df_model2$estimate,
  ci.lb = results_df_model2$ci_lb,
  ci.ub = results_df_model2$ci_ub,
  slab = paste(results_df_model2$platform, results_df_model2$variable),
  xlab = "Effect Size",
  refline = 0,
  shade = TRUE,
  col = results_df_model2$color,
  ylim = c(-3, nrow(results_df_model2) + 4),
  digits = 3
)
addpoly(res_model2_pc1, col = "darkred", cex = 1, border = "darkred", lwd = 2, row = -1, mlab = 'pc1')
addpoly(res_model2_toxicity, col = "darkgreen", cex = 1, border = "darkgreen", lwd = 2, row = -2, mlab = 'toxicity')
addpoly(res_model2_certainty, col = "darkblue", cex = 1, border = "darkblue", lwd = 2, row = -3, mlab = 'certainty')
title(main = "Model 2: pc1 + toxicity + confidence ")