Descriptive statistics

In this script, I write the code to make some plots and tables that I use to in the data and analysis section. You will find the code for the plot with datasets per year, the Table with individual regression analysis for one variable of ISSP and the Table with all the pooled effects.

library(dplyr)
library(ggplot2)
library(kableExtra)
library(modelsummary)
library(tidyverse)
library(writexl)
library(readxl)
library(gamlss)
library(metafor)
library(meta)
library(ggrepel)
rm(list=ls())

Ggplot that shows the years that each dataset covers

datayear <- read_excel("./data/descriptives/dataset_year.xlsx")
datayear <- as.data.frame(datayear)

color_scale <- c ("#E02023", "#1676E3", "#20C66B", "#1ABC9C", "#85C1E9" , "#EB984E", "#AF7AC5" , "#F1C40F")
datayear$dataset <- factor(datayear$dataset, levels = c("MOT", "SOCON", "LISS", "I&O", "ESS", "ISSP", "EVS", "EB"))
datayear$color[datayear$dataset == "I&O"] <- "#1ABC9C"
datayear$color[datayear$dataset == "ESS"] <- "#85C1E9"
datayear$color[datayear$dataset == "ISSP"] <- "#EB984E"
datayear$color[datayear$dataset == "EVS"] <- "#AF7AC5"
datayear$color[datayear$dataset == "EB"] <- "#F1C40F"
datayear$color[datayear$dataset == "MOT"] <- "#E02023"
datayear$color[datayear$dataset == "SOCON"] <- "#1676E3"
datayear$color[datayear$dataset == "LISS"] <- "#20C66B"

first_last_year <- datayear %>%
  group_by(dataset) %>%
  summarize(first_year = min(year),
            last_year = max(year))

yearplot <- ggplot(data = datayear, aes(x = year, y = dataset, color = color)) +
  geom_point(size = 2) +
  geom_line(aes(group = dataset), size = 1) +
  geom_text(data = first_last_year, aes(x = first_year, y = dataset, label = as.character(first_year)),
            vjust = -1, color = "darkgrey", size = 3, show.legend = FALSE) + 
    geom_text(data = first_last_year, aes(x = last_year, y = dataset, label = as.character(last_year)),
            vjust = -1, color = "darkgrey", size = 3, show.legend = FALSE) +
  scale_color_manual(values = color_scale) +
  labs(x = "Year", y = "Dataset") +
  theme_classic() + 
  theme(legend.position = "none")
print(yearplot)

ggsave("yearplot.png", width = 7, height = 5, path = "./output/yearplot_DEC.png")

Table of the GAMLSS regression output

# Run this piece of code so that in the Table beneath it will not show NAs
options(knitr.kable.NA = '')

# Illustrate the gamlss regressions by showing one from ISSP with kable extra
load("./data/final_data/regression_outputs/issp_list_empty_w.Rdata")
load("./data/final_data/regression_outputs/issp_list_preds_w_new.Rdata")
load("./data/final_data/regression_outputs/issp_list_interactions_w_new.Rdata"))

# Use first variable, which is "worry"
model_time_only <- summary(issp_list_empty[[1]])
model_time_only_df <- as.data.frame(model_time_only)
model_time_only_df$mu_sigma <- c("Mu", "", "Sigma", "")
model_time_only_df$variables <- c("Intercept", "Time", "Intercept", "Time")
model_time_only_df <- model_time_only_df %>% dplyr:: select(mu_sigma, variables, everything())
rownames(model_time_only_df) <- NULL

model_indep_df <- as.data.frame(summary(issp_list_preds[[1]]))
model_indep_df$mu_sigma <- c("Mu", "", "", "", "", "","Sigma", "", "", "", "", "")
model_indep_df$variables <- c("Intercept", "Time", "Gender",  "Isced intermediate (ref=basic)", "Isced advanced (ref=basic)", "Age", "Intercept", "Time", "Gender",  "Isced intermediate (ref=basic)", "Isced advanced (ref=basic)", "Age")
model_indep_df <- model_indep_df %>% dplyr::select(mu_sigma, variables, everything())
rownames(model_indep_df) <- NULL

model_interact_df <- as.data.frame(summary(issp_list_interactions[[1]]))
model_interact_df$mu_sigma <- c("Mu", "", "", "", "", "", "", "", "", "", "Sigma", "", "", "", "", "", "", "", "", "")
model_interact_df$variables <- c("Intercept", "Time", "Gender",  "Isced intermediate (ref=basic)", "Isced advanced (ref=basic)", "Age", "Age x time", "Gender x time", "Isced intermediate x time", "Isced advanced x time", "Intercept", "Time", "Gender",  "Isced intermediate (ref=basic)", "Isced advanced (ref=basic)", "Age", "Age x time", "Gender x time", "Isced intermediate x time", "Isced advanced x time")
model_interact_df <- model_interact_df %>% dplyr::select(mu_sigma, variables, everything())
rownames(model_interact_df) <- NULL

# Difficult to columnbind them, because they're of different lengths (logically), so use a function
df_list <- list(model_time_only_df, model_indep_df, model_interact_df)

bind_cols_fill <- function(df_list) {

  max_rows <- map_int(df_list, nrow) %>% max()
  
  map(df_list, function(df) {
    if(nrow(df) == max_rows) return(df)
    first <- names(df)[1] %>% sym()
    df %>% add_row(!!first := rep(NA, max_rows - nrow(df)))
  }) %>% bind_cols()
}

model_total <- bind_cols_fill(df_list)

model_total[is.na(model_total)] <- ""

# The dataframe is binded, the order is not totally how i want it. The data are now all centered at the top, but for the 
# previous models, I need the data to be centered around the mu/intercept and sigma/intercept (with a lot of blank spaces). # Furthermore, some variables have a name that has special signs in them, which causes R not to recognize them. To make it #easy, I quickly fix these things in Excel. I don't make any changes to the data
#write_xlsx(model_total, "./data/descriptives/kable_model_total_gamlss.xlsx" )

model_new <- read_excel("./data/descriptives/kable_model_total_gamlss.xlsx" )

# I Want to round the numeric variables to 3 decimals
model_new <- dplyr::mutate(model_new, across(.cols = c(-mu_sigma, -variables), .fns= as.numeric))
numeric_cols <- sapply(model_new, is.numeric)
model_new[, numeric_cols] <- round(model_new[, numeric_cols], digits = 3)

# There are some annoying variable names that R can't recognize because it has spacing, so rename the columns
colnames(model_new) <- (c("mu_sigma", "Variables", "est_time", "std_err_time", "p_time", "est_indep", "std_err_indep", "p_indep", "est_inter", "std_error_inter", "p_inter", "p_text_time_only", "p_text_indep", "p_text_inter"))

# Create a new variable to use stars instead of the full p-values in the table
model_new$p_text_time_only <- case_when(
  model_new$p_time >= 0.05 ~ "",
model_new$p_time < 0.001 ~ "***",
model_new$p_time < 0.01 ~ "**",
model_new$p_time < 0.05 ~ "*"
)

model_new$p_text_indep <- case_when(
  model_new$p_indep >= 0.05 ~ "",
model_new$p_indep < 0.001 ~ "***",
model_new$p_indep < 0.01 ~ "**",
model_new$p_indep < 0.05 ~ "*"
)

model_new$p_text_inter <- case_when(
  model_new$p_inter >= 0.05 ~ "",
model_new$p_inter < 0.001 ~ "***",
model_new$p_inter < 0.01 ~ "**",
model_new$p_inter < 0.05 ~ "*"
)

# Now only select the variables that I want to show in the table
model_new <- dplyr::select(model_new, mu_sigma, Variables, est_time, std_err_time, p_text_time_only, est_indep, std_err_indep, p_text_indep, est_inter, std_error_inter, p_text_inter)

# Make and save the table 
kable(model_new, digits = 3, caption = "Table 2. GAMLSS model on worry ISSP (1993 - 2010)", 
      col.names = c( "", "Variables", "Estimate", "Std. Error", "", "Estimate", "Std. Error", "", "Estimate", "Std. Error", ""), align = "l") %>%
  kable_classic_2(html_font = "Times", fixed_thead = T, full_width = F) %>%
  add_header_above(c(" " = 2, "Time-only model" = 3, "Indepent var model" = 3, "Interaction model" = 3), align = "l") %>%
   footnote(
    general_title = "Note.",
    general = "*** = p < 0.001, ** = p < 0.01, * = p < 0.05",
    threeparttable = TRUE,
    footnote_as_chunk = TRUE
    ) 

Table 2. GAMLSS model on worry ISSP (1993 - 2010)

		Time-only model			Indepent var model			Interaction model
	Variables	Estimate	Std. Error		Estimate	Std. Error		Estimate	Std. Error
Mu	Intercept	3.358	0.029	***	3.453	0.086	***	3.602	0.120	***
	Time	-0.005	0.003		-0.011	0.003	***	-0.045	0.015	**
	Gender				0.183	0.038	***	0.173	0.055	**
	Isced intermediate (ref=basic)				0.468	0.040	***	0.397	0.055	***
	Isced advanced (ref=basic)				0.002	0.031		0.004	0.045
	Age				-0.008	0.001	***	-0.011	0.002	***
	Age x time							0.000	0.000
	Gender x time							0.002	0.005
	Isced intermediate x time							0.009	0.008
	Isced advanced x time							0.016	0.008
Sigma	Intercept	1.065	0.020	***	1.109	0.061	***	1.038	0.085	***
	Time	0.001	0.002		0.003	0.002		0.008	0.011
	Gender				-0.087	0.027	**	-0.070	0.039
	Isced intermediate (ref=basic)				-0.093	0.028	**	-0.122	0.039	**
	Isced advanced (ref=basic)				-0.028	0.022		-0.023	0.032
	Age				0.000	0.001		0.002	0.001
	Age x time							0.000	0.000
	Gender x time							-0.002	0.004
	Isced intermediate x time							0.006	0.006
	Isced advanced x time							0.006	0.006
Note. *** = p < 0.001, ** = p < 0.01, * = p < 0.05

#%>%
    #save_kable("/Users/anuschka/Documents/climatechange/climatechange/output/gamlss_table.pdf")

Table all pooled effects

rm(list=ls())
# I want to repeat the above Table, but instead of regressions I now want to put all the pooled effects in there
# Starting with the time only model
# You may notice that I use the rma function and metagen. The estimate TE.random of metagen (pooled effect) is the same as the beta in the rma model
set.seed(1) #(as I am performing analyses with a random aspect)

load("./data/meta_analysis/total_reg_results_gam_w_new.RData")

# Extract the pooled effects, first mean attitudes
model_step_1_att <- rma(yi = mu_time,
              sei = mu_time_sd,
              data = total_reg_results,
              method = "ML",
              test = "knha")

one <- summary(model_step_1_att)
one_df <- data.frame(matrix(ncol = 5, nrow = 1))
colnames(one_df) <- c("mu_sigma", "Variable", "Estimate", "SE", "p")

one_df$mu_sigma <- c("Mean attitudes")
one_df$Variable <- "Time"
one_df$Estimate <- as.numeric(one$b)
one_df$SE <- one$se
one_df$p <- one$pval

# And now variance
model_step_1_var <- rma(yi = sig_time,
              sei = sig_time_sd,
              data = total_reg_results,
              method = "ML",
              test = "knha")

two <- summary(model_step_1_var)
two_df <- data.frame(matrix(ncol = 5, nrow = 1))
colnames(two_df) <- c("mu_sigma", "Variable", "Estimate", "SE", "p")

two_df$mu_sigma <- c("Polarization attitudes")
two_df$Variable <- "Time"
two_df$Estimate <- as.numeric(two$b)
two_df$SE <- two$se
two_df$p <- two$pval

# Rowbind them, these are all the pooled effects of the time-only model
empty_df <- rbind(one_df, two_df)

# Now continue with the independent variable model, first the mean attitudes 
load("./data/meta_analysis/total_indep_var_results_new.RData")

m2_att <- rma(yi = mu_time_pred,
              sei = mu_time_sd_pred,
              data = total_reg_results,
              method = "ML",
              test = "knha")

# Gender 
gender_gen <- metagen(TE = mu_sex_est,
                 seTE = mu_sex_sd,
                 data = total_indep_var_results,
                 fixed = FALSE,
                 random = TRUE,
                 method.tau = "ML",
                 hakn = TRUE)

# Age 
age_gen <- metagen(TE = mu_age_est,
                 seTE = mu_age_sd,
                 data = total_indep_var_results,
                 fixed = FALSE,
                 random = TRUE,
                 method.tau = "ML",
                 hakn = TRUE)

# Educational level. This is now a dummy
edu_med_gen <- metagen(TE = mu_isced_med_est,
                 seTE = mu_isced_med_sd,
                 data = total_indep_var_results,
                 fixed = FALSE,
                 random = TRUE,
                 method.tau = "ML",
                 hakn = TRUE)

edu_high_gen <- metagen(TE = mu_isced_high_est,
                 seTE = mu_isced_high_sd,
                 data = total_indep_var_results,
                 fixed = FALSE,
                 random = TRUE,
                 method.tau = "ML",
                 hakn = TRUE)

# Extract all the data from the models and put them in a dataframe that will later be the input of the kable table
three_df <- data.frame(matrix(ncol = 5, nrow = 5))
colnames(three_df) <- c("mu_sigma_1", "Variable_1", "Estimate_1", "SE_1", "p_1")
three_df$mu_sigma_1 <- c("Mean attitudes", "", "", "", "")
three_df$Variable_1 <- c("Time", "Gender", "Age", "Isced intermediate (ref=basic)", "Isced advanced (ref=basic)")
three_df$Estimate_1 <- c(as.numeric(m2_att$b), gender_gen$TE.random, age_gen$TE.random, edu_med_gen$TE.random, edu_high_gen$TE.random)
three_df$SE_1 <- c(m2_att$se, gender_gen$seTE.random, age_gen$seTE.random, edu_med_gen$seTE.random, edu_high_gen$seTE.random)
three_df$p_1 <- c(m2_att$pval, gender_gen$pval.random, age_gen$pval.random, edu_med_gen$pval.random, edu_high_gen$pval.random)                       
                       
# The same for variance
m2_var <- rma(yi = sig_time_pred,
              sei = sig_time_sd_pred,
              data = total_reg_results,
              method = "ML",
              test = "knha")

# Gender
gender_var_gen <- metagen(TE = sig_sex_est,
                 seTE = sig_sex_sd,
                 data = total_indep_var_results,
                 fixed = FALSE,
                 random = TRUE,
                 method.tau = "ML",
                 hakn = TRUE)

# Age
age_var_gen <- metagen(TE = sig_age_est,
                 seTE = sig_age_sd,
                 data = total_indep_var_results,
                 fixed = FALSE,
                 random = TRUE,
                 method.tau = "ML",
                 hakn = TRUE)

# Education
edu_var_med_gen <- metagen(TE = sig_isced_med_est,
                 seTE = sig_isced_med_sd,
                 data = total_indep_var_results,
                 fixed = FALSE,
                 random = TRUE,
                 method.tau = "ML",
                 hakn = TRUE)

edu_var_high_gen <- metagen(TE = sig_isced_high_est,
                 seTE = sig_isced_high_sd,
                 data = total_indep_var_results,
                 fixed = FALSE,
                 random = TRUE,
                 method.tau = "ML",
                 hakn = TRUE)

four_df <- data.frame(matrix(ncol = 5, nrow = 5))
colnames(four_df) <- c("mu_sigma_1", "Variable_1", "Estimate_1", "SE_1", "p_1")
four_df$mu_sigma_1 <- c("Polarization attitudes", "", "", "", "")
four_df$Variable_1 <- c("Time", "Gender", "Age", "Isced intermediate (ref=basic)", "Isced advanced (ref=basic)")
four_df$Estimate_1 <- c(as.numeric(m2_var$b), gender_var_gen$TE.random, age_var_gen$TE.random, edu_var_med_gen$TE.random, edu_var_high_gen$TE.random)
four_df$SE_1 <- c(m2_att$se, gender_var_gen$seTE.random, age_var_gen$seTE.random, edu_var_med_gen$seTE.random, edu_var_high_gen$seTE.random)
four_df$p_1 <- c(m2_att$pval, gender_var_gen$pval.random, age_var_gen$pval.random, edu_var_med_gen$pval.random, edu_var_high_gen$pval.random)


# Bind rows of these two, these are all the pooled effects of the independent variables model
indep_df <- rbind(three_df, four_df)

#Now i want the interaction model as well
load("./data/meta_analysis/total_inter_results_new.RData")

# First mean attitudes
m3_att <- rma(yi = mu_time_int,
              sei = mu_time_sd_int,
              data = total_reg_results,
              method = "ML",
              test = "knha")

# Gender 
gender_gen_int <- metagen(TE = mu_sex_est,
                 seTE = mu_sex_sd,
                 data = total_inter_results,
                 fixed = FALSE,
                 random = TRUE,
                 method.tau = "ML",
                 hakn = TRUE)

# Age 
age_gen_int <- metagen(TE = mu_age_est,
                 seTE = mu_age_sd,
                 data = total_inter_results,
                 fixed = FALSE,
                 random = TRUE,
                 method.tau = "ML",
                 hakn = TRUE)

# Educational level
edu_med_gen_int <- metagen(TE = mu_isced_med_est,
                 seTE = mu_isced_med_sd,
                 data = total_inter_results,
                 fixed = FALSE,
                 random = TRUE,
                 method.tau = "ML",
                 hakn = TRUE)

edu_high_gen_int <- metagen(TE = mu_isced_high_est,
                 seTE = mu_isced_high_sd,
                 data = total_inter_results,
                 fixed = FALSE,
                 random = TRUE,
                 method.tau = "ML",
                 hakn = TRUE)

# And the interaction effects 
# Gender 
gender_int_eff <- metagen(TE = mu_sex_est_inter,
                 seTE = mu_sex_sd_inter,
                 data = total_inter_results,
                 fixed = FALSE,
                 random = TRUE,
                 method.tau = "ML",
                 hakn = TRUE)

# Age 
age_int_eff <- metagen(TE = mu_age_est_inter,
                 seTE = mu_age_sd_inter,
                 data = total_inter_results,
                 fixed = FALSE,
                 random = TRUE,
                 method.tau = "ML",
                 hakn = TRUE)

# Educational level
edu_med_int_eff <- metagen(TE = mu_isced_med_est_inter,
                 seTE = mu_isced_med_sd_inter,
                 data = total_inter_results,
                 fixed = FALSE,
                 random = TRUE,
                 method.tau = "ML",
                 hakn = TRUE)

edu_high_int_eff <- metagen(TE = mu_isced_high_est_inter,
                 seTE = mu_isced_high_sd_inter,
                 data = total_inter_results,
                 fixed = FALSE,
                 random = TRUE,
                 method.tau = "ML",
                 hakn = TRUE)

# Now extract the data for a table
five_df <- data.frame(matrix(ncol = 5, nrow = 9))
colnames(five_df) <- c("mu_sigma_2", "Variable_2", "Estimate_2", "SE_2", "p_2")
five_df$mu_sigma_2 <- c("Mean attitudes", "", "", "", "", "", "", "", "")
five_df$Variable_2 <- c("Time", "Gender", "Age", "Isced intermediate (ref=basic)", "Isced advanced (ref=basic)", "Gender x time", "Age x time", "Isced intermediate x time", "Isced advanced x time")
five_df$Estimate_2 <- c(as.numeric(m3_att$b), gender_gen_int$TE.random, age_gen_int$TE.random, edu_med_gen_int$TE.random, edu_high_gen_int$TE.random, gender_int_eff$TE.random, age_int_eff$TE.random, edu_med_int_eff$TE.random, edu_high_int_eff$TE.random)
five_df$SE_2 <- c(m3_att$se, gender_gen_int$seTE.random, age_gen_int$seTE.random, edu_med_gen_int$seTE.random, edu_high_gen_int$seTE.random, gender_int_eff$seTE.random, age_int_eff$seTE.random, edu_med_int_eff$seTE.random, edu_high_int_eff$seTE.random)
five_df$p_2 <- c(m3_att$pval, gender_gen_int$pval.random, age_gen_int$pval.random, edu_med_gen_int$pval.random, edu_high_gen_int$pval.random, gender_int_eff$pval.random, age_int_eff$pval.random, edu_med_int_eff$pval.random, edu_high_int_eff$pval.random)

# Now the same procedure for the variance
m3_var <- rma(yi = sig_time_int,
              sei = sig_time_sd_int,
              data = total_reg_results,
              method = "ML",
              test = "knha")

# Gender 
gender_gen_int_var <- metagen(TE = sig_sex_est,
                 seTE = sig_sex_sd,
                 data = total_inter_results,
                 fixed = FALSE,
                 random = TRUE,
                 method.tau = "ML",
                 hakn = TRUE)

# Age 
age_gen_int_var <- metagen(TE = sig_age_est,
                 seTE = sig_age_sd,
                 data = total_inter_results,
                 fixed = FALSE,
                 random = TRUE,
                 method.tau = "ML",
                 hakn = TRUE)

# Educational level
edu_med_gen_int_var <- metagen(TE = sig_isced_med_est,
                 seTE = sig_isced_med_sd,
                 data = total_inter_results,
                 fixed = FALSE,
                 random = TRUE,
                 method.tau = "ML",
                 hakn = TRUE)

edu_high_gen_int_var <- metagen(TE = sig_isced_high_est,
                 seTE = sig_isced_high_sd,
                 data = total_inter_results,
                 fixed = FALSE,
                 random = TRUE,
                 method.tau = "ML",
                 hakn = TRUE)

# And the interaction effects 
# Gender 
gender_int_eff_var <- metagen(TE = sig_sex_est_inter,
                 seTE = sig_sex_sd_inter,
                 data = total_inter_results,
                 fixed = FALSE,
                 random = TRUE,
                 method.tau = "ML",
                 hakn = TRUE)

# Age 
age_int_eff_var <- metagen(TE = sig_age_est_inter,
                 seTE = sig_age_sd_inter,
                 data = total_inter_results,
                 fixed = FALSE,
                 random = TRUE,
                 method.tau = "ML",
                 hakn = TRUE)

# Educational level
edu_med_int_eff_var <- metagen(TE = sig_isced_med_est_inter,
                 seTE = sig_isced_med_sd_inter,
                 data = total_inter_results,
                 fixed = FALSE,
                 random = TRUE,
                 method.tau = "ML",
                 hakn = TRUE)

edu_high_int_eff_var <- metagen(TE = sig_isced_high_est_inter,
                 seTE = sig_isced_high_sd_inter,
                 data = total_inter_results,
                 fixed = FALSE,
                 random = TRUE,
                 method.tau = "ML",
                 hakn = TRUE)

six_df <- data.frame(matrix(ncol = 5, nrow = 9))
colnames(six_df) <- c("mu_sigma_2", "Variable_2", "Estimate_2", "SE_2", "p_2")
six_df$mu_sigma_2 <- c("Polarization attitudes", "", "", "", "", "", "", "", "")
six_df$Variable_2 <- c("Time", "Gender", "Age", "Isced intermediate (ref=basic)", "Isced advanced (ref=basic)", "Gender x time", "Age x time", "Isced intermediate x time", "Isced advanced x time")
six_df$Estimate_2 <- c(as.numeric(m3_var$b), gender_gen_int_var$TE.random, age_gen_int_var$TE.random, edu_med_gen_int_var$TE.random, edu_high_gen_int_var$TE.random, gender_int_eff_var$TE.random, age_int_eff_var$TE.random, edu_med_int_eff_var$TE.random, edu_high_int_eff_var$TE.random)
six_df$SE_2 <- c(m3_var$se, gender_gen_int_var$seTE.random, age_gen_int_var$seTE.random, edu_med_gen_int_var$seTE.random, edu_high_gen_int_var$seTE.random, gender_int_eff_var$seTE.random, age_int_eff_var$seTE.random, edu_med_int_eff_var$seTE.random, edu_high_int_eff_var$seTE.random)
six_df$p_2 <- c(m3_att$pval, gender_gen_int_var$pval.random, age_gen_int_var$pval.random, edu_med_gen_int_var$pval.random, edu_high_gen_int_var$pval.random, gender_int_eff_var$pval.random, age_int_eff_var$pval.random, edu_med_int_eff_var$pval.random, edu_high_int_eff_var$pval.random)

#Bind these two, this now has all the pooled effects of the interaction effects
inter_df <- rbind(five_df, six_df)

#Columnbind the three datasets
# Difficult to columnbind them, because they're of different lengths (logically)
df_list <- list(empty_df, indep_df, inter_df)

bind_cols_fill <- function(df_list) {

  max_rows <- map_int(df_list, nrow) %>% max()
  
  map(df_list, function(df) {
    if(nrow(df) == max_rows) return(df)
    first <- names(df)[1] %>% sym()
    df %>% add_row(!!first := rep(NA, max_rows - nrow(df)))
  }) %>% bind_cols()
}

total_df <- bind_cols_fill(df_list)

#Change the order a bit, just like in the previous table 
write_xlsx(total_df,  "./data/descriptives/kable_model_all_pooled_effects.xlsx" )

model_new <- read_excel("./data/descriptives/kable_model_all_pooled_effects.xlsx" )

# I Want to round the numeric variables to 3 decimals
model_new <- dplyr::mutate(model_new, across(.cols = c(-mu_sigma_2, -Variable_2), .fns = as.numeric))
numeric_cols <- sapply(model_new, is.numeric)
model_new[, numeric_cols] <- round(model_new[, numeric_cols], digits = 3)

# Create a new variable to use stars instead of the full p-values in the table
model_new$p_text <- case_when(
  model_new$p >= 0.05 ~ "",
model_new$p < 0.001 ~ "***",
model_new$p < 0.01 ~ "**",
model_new$p < 0.05 ~ "*"
)

model_new$p_text_1 <- case_when(
  model_new$p_1 >= 0.05 ~ "",
model_new$p_1 < 0.001 ~ "***",
model_new$p_1 < 0.01 ~ "**",
model_new$p_1 < 0.05 ~ "*"
)

model_new$p_text_2 <- case_when(
  model_new$p_2 >= 0.05 ~ "",
model_new$p_2 < 0.001 ~ "***",
model_new$p_2 < 0.01 ~ "**",
model_new$p_2 < 0.05 ~ "*"
)

# Now only select the variables that I want to show in the table
model_new <- dplyr::select(model_new, mu_sigma_2, Variable_2, Estimate, SE, p_text, Estimate_1, SE_1, p_text_1, Estimate_2, SE_2, p_text_2)

# Make and save the table as html file
kable(model_new, digits = 3, caption = "Table 3. All pooled effects mean attitudes and polarization", 
      col.names = c( "", "Variables", "Estimate", "Std. Error", "", "Estimate", "Std. Error", "", "Estimate", "Std. Error", "")) %>%
  kable_classic_2(html_font = "Times", fixed_thead = T, full_width = F) %>%
  add_header_above(c(" " = 2, "Time-only model" = 3, "Indepent var model" = 3, "Interaction model" = 3)) %>%
   footnote(
    general_title = "Note.",
    general = "*** = p < 0.001, ** = p < 0.01, * = p < 0.05",
    threeparttable = TRUE,
    footnote_as_chunk = TRUE
    ) #%>%

Table 3. All pooled effects mean attitudes and polarization

		Time-only model			Indepent var model			Interaction model
	Variables	Estimate	Std. Error		Estimate	Std. Error		Estimate	Std. Error
Mean attitudes	Time	0.017	0.016		0.014	0.017		-0.025	0.019
	Gender				0.128	0.016	***	0.104	0.025	***
	Age				0.002	0.001	**	0.000	0.001
	Isced intermediate (ref=basic)				0.161	0.024	***	0.117	0.033	**
	Isced advanced (ref=basic)				0.035	0.009	***	0.029	0.013
	Gender x time							0.001	0.004
	Age x time							0.001	0.000
	Isced intermediate x time							0.008	0.002	**
	Isced advanced x time							0.018	0.009
Polarization attitudes	Time	-0.009	0.004		-0.011	0.017		0.001	0.005
	Gender				-0.092	0.009	***	0.008	0.002	**
	Age				0.000	0.000		0.000	0.005
	Isced intermediate (ref=basic)				-0.041	0.008	***	0.018	0.009
	Isced advanced (ref=basic)				-0.004	0.004		1.183	0.043	***
	Gender x time							-0.001	0.001
	Age x time							0.000	0.000
	Isced intermediate x time							0.001	0.002
	Isced advanced x time							-0.001	0.003
Note. *** = p < 0.001, ** = p < 0.01, * = p < 0.05

    #save_kable("./output/pooled_effects_table_DECEMBER.html")

Appendix D. Pooled effects age-only vs gender and education-only model

rm(list=ls())
set.seed(1) #(as I am performing analyses with a random aspect)
load("./data/meta_analysis/total_indep_var_results_new_edugen.RData")

# The models in which I only included education and gender 
m2_att <- rma(yi = mu_time,
              sei = mu_time_sd,
              data = total_indep_var_results,
              method = "REML",
              test = "knha")

# Gender 
gender_gen <- metagen(TE = mu_sex_est,
                 seTE = mu_sex_sd,
                 data = total_indep_var_results,
                 fixed = FALSE,
                 random = TRUE,
                 method.tau = "REML",
                 hakn = TRUE)

# Educational level. This is now a dummy
edu_med_gen <- metagen(TE = mu_isced_med_est,
                 seTE = mu_isced_med_sd,
                 data = total_indep_var_results,
                 fixed = FALSE,
                 random = TRUE,
                 method.tau = "REML",
                 hakn = TRUE)

edu_high_gen <- metagen(TE = mu_isced_high_est,
                 seTE = mu_isced_high_sd,
                 data = total_indep_var_results,
                 fixed = FALSE,
                 random = TRUE,
                 method.tau = "REML",
                 hakn = TRUE)

# Extract all the data from the models and put them in a dataframe that will later be the input of the kable table
three_df <- data.frame(matrix(ncol = 6, nrow = 4))
colnames(three_df) <- c("mu_sigma_1", "age_educ_1", "Variable_1", "Estimate_1", "SE_1", "p_1")
three_df$mu_sigma_1 <- c("Mean attitudes", "", "", "")
three_df$age_educ_1 <- c("Gender and education", "", "", "")
three_df$Variable_1 <- c("Time", "Gender", "Isced intermediate (ref=basic)", "Isced advanced (ref=basic)")
three_df$Estimate_1 <- c(as.numeric(m2_att$b), gender_gen$TE.random,edu_med_gen$TE.random, edu_high_gen$TE.random)
three_df$SE_1 <- c(m2_att$se, gender_gen$seTE.random, edu_med_gen$seTE.random, edu_high_gen$seTE.random)
three_df$p_1 <- c(m2_att$pval, gender_gen$pval.random, edu_med_gen$pval.random, edu_high_gen$pval.random)                       
                       
# The same for variance
m2_var <- rma(yi = sig_time,
              sei = sig_time_sd,
              data = total_indep_var_results,
              method = "REML",
              test = "knha")

# Gender
gender_var_gen <- metagen(TE = sig_sex_est,
                 seTE = sig_sex_sd,
                 data = total_indep_var_results,
                 fixed = FALSE,
                 random = TRUE,
                 method.tau = "REML",
                 hakn = TRUE)

# Education
edu_var_med_gen <- metagen(TE = sig_isced_med_est,
                 seTE = sig_isced_med_sd,
                 data = total_indep_var_results,
                 fixed = FALSE,
                 random = TRUE,
                 method.tau = "REML",
                 hakn = TRUE)

edu_var_high_gen <- metagen(TE = sig_isced_high_est,
                 seTE = sig_isced_high_sd,
                 data = total_indep_var_results,
                 fixed = FALSE,
                 random = TRUE,
                 method.tau = "REML",
                 hakn = TRUE)

four_df <- data.frame(matrix(ncol = 6, nrow = 4))
colnames(four_df) <- c("mu_sigma_1", "age_educ_1", "Variable_1", "Estimate_1", "SE_1", "p_1")
four_df$mu_sigma_1 <- c("Polarization attitudes", "", "", "")
four_df$age_educ_1 <- c("Gender and education", "", "", "")
four_df$Variable_1 <- c("Time", "Gender", "Isced intermediate (ref=basic)", "Isced advanced (ref=basic)")
four_df$Estimate_1 <- c(as.numeric(m2_var$b), gender_var_gen$TE.random, edu_var_med_gen$TE.random, edu_var_high_gen$TE.random)
four_df$SE_1 <- c(m2_att$se, gender_var_gen$seTE.random, edu_var_med_gen$seTE.random, edu_var_high_gen$seTE.random)
four_df$p_1 <- c(m2_att$pval, gender_var_gen$pval.random, edu_var_med_gen$pval.random, edu_var_high_gen$pval.random)

rm(total_indep_var_results)
# Do the same, but then for the models with only age included
load("./data/meta_analysis/total_indep_var_results_new_age.RData")
total_indep_var_results_age <- as.data.frame(total_indep_var_results)

# The models in which I only included education and gender 
m2_att_age <- rma(yi = mu_time,
              sei = mu_time_sd,
              data = total_indep_var_results_age,
              method = "REML",
              test = "knha")

# Age
age_gen <- metagen(TE = mu_age_est,
                 seTE = mu_age_sd,
                 data = total_indep_var_results_age,
                 fixed = FALSE,
                 random = TRUE,
                 method.tau = "REML",
                 hakn = TRUE)

# Extract all the data from the models and put them in a dataframe that will later be the input of the kable table
five_df <- data.frame(matrix(ncol = 6, nrow = 2))
colnames(five_df) <- c("mu_sigma_1", "age_educ_1", "Variable_1", "Estimate_1", "SE_1", "p_1")
five_df$mu_sigma_1 <- c("", "")
five_df$age_educ_1 <- c("Age", "")
five_df$Variable_1 <- c("Time", "Age")
five_df$Estimate_1 <- c(as.numeric(m2_att_age$b), age_gen$TE.random)
five_df$SE_1 <- c(m2_att_age$se, age_gen$seTE.random)
five_df$p_1 <- c(m2_att_age$pval, age_gen$pval.random)                       
                       
# The same for variance
m2_var_age <- rma(yi = sig_time,
              sei = sig_time_sd,
              data = total_indep_var_results_age,
              method = "REML",
              test = "knha")

# Gender
age_var_gen <- metagen(TE = sig_age_est,
                 seTE = sig_age_sd,
                 data = total_indep_var_results_age,
                 fixed = FALSE,
                 random = TRUE,
                 method.tau = "REML",
                 hakn = TRUE)

six_df <- data.frame(matrix(ncol = 6, nrow = 2))
colnames(six_df) <- c("mu_sigma_1", "age_educ_1", "Variable_1", "Estimate_1", "SE_1", "p_1")
six_df$mu_sigma_1 <- c("", "")
six_df$age_educ_1 <- c("Age", "")
six_df$Variable_1 <- c("Time", "Age")
six_df$Estimate_1 <- c(as.numeric(m2_var_age$b), age_var_gen$TE.random)
six_df$SE_1 <- c(m2_var_age$se, age_var_gen$seTE.random)
six_df$p_1 <- c(m2_var_age$pval, age_var_gen$pval.random)

# The order I want in the tables is three -> five -> four -> six
total_df <- rbind(three_df, five_df, four_df, six_df)

# I Want to round the numeric variables to 3 decimals
total_df <- dplyr::mutate(total_df, across(.cols = c(-mu_sigma_1, -Variable_1, -age_educ_1), .fns = as.numeric))
numeric_cols <- sapply(total_df, is.numeric)
total_df[, numeric_cols] <- round(total_df[, numeric_cols], digits = 3)

# Create a new variable to use stars instead of the full p-values in the table
total_df$p_text_1 <- case_when(
  total_df$p_1 >= 0.05 ~ "",
total_df$p_1 < 0.001 ~ "***",
total_df$p_1 < 0.01 ~ "**",
total_df$p_1 < 0.05 ~ "*"
)

# Now only select the variables that I want to show in the table
total_df <- dplyr::select(total_df, mu_sigma_1, age_educ_1, Variable_1, Estimate_1, SE_1, p_text_1)

# Make and save the table as html file
kable(total_df, digits = 3, caption = "Appendix D. Pooled effects mean attitudes and polarization age, education and gender.", 
      col.names = c( "", "", "Variables", "Estimate", "Std. Error", "")) %>%
  kable_classic_2(html_font = "Times", fixed_thead = T, full_width = F) %>%
  add_header_above(c(" " = 3, "Indepent var model" = 3)) %>%
   footnote(
    general_title = "Note.",
    general = "*** = p < 0.001, ** = p < 0.01, * = p < 0.05",
    threeparttable = TRUE,
    footnote_as_chunk = TRUE
    )#%>%

Appendix D. Pooled effects mean attitudes and polarization age, education and gender.

			Indepent var model
		Variables	Estimate	Std. Error
Mean attitudes	Gender and education	Time	0.020	0.017
		Gender	0.134	0.016	***
		Isced intermediate (ref=basic)	0.167	0.025	***
		Isced advanced (ref=basic)	0.034	0.009	***
	Age	Time	0.016	0.016
		Age	0.000	0.001
Polarization attitudes	Gender and education	Time	-0.019	0.017
		Gender	-0.105	0.010	***
		Isced intermediate (ref=basic)	-0.071	0.009	***
		Isced advanced (ref=basic)	-0.003	0.004
	Age	Time	-0.009	0.004
		Age	0.001	0.000
Note. *** = p < 0.001, ** = p < 0.01, * = p < 0.05

 #save_kable("./output/appendix/pooled_effects_edugenage_DECEMBER.html")

0.1 Repeat main analyses without outlier(s)

# Load the data
load("./data/meta_analysis/total_reg_results_gam_w_new.RData")

# Run the first model and check for outliers
model_step_1_att <- rma(yi = mu_time,
              sei = mu_time_sd,
              data = total_reg_results,
              method = "ML",
              test = "knha")

model_step_1_att

## 
## Random-Effects Model (k = 71; tau^2 estimator: ML)
## 
## tau^2 (estimated amount of total heterogeneity): 0.0155 (SE = 0.0027)
## tau (square root of estimated tau^2 value):      0.1247
## I^2 (total heterogeneity / total variability):   99.82%
## H^2 (total variability / sampling variability):  557.16
## 
## Test for Heterogeneity:
## Q(df = 70) = 4216.1447, p-val < .0001
## 
## Model Results:
## 
## estimate      se    tval  df    pval    ci.lb   ci.ub    
##   0.0165  0.0160  1.0296  70  0.3067  -0.0155  0.0485    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

inf <- influence(model_step_1_att) # There is one outlier

# Delete outlier
total_reg_results <- subset(total_reg_results, subset = dep_var != "cc_unstop")

# Run the same analyses
model_step_1_att <- rma(yi = mu_time,
              sei = mu_time_sd,
              data = total_reg_results,
              method = "ML",
              test = "knha")

model_step_1_att

## 
## Random-Effects Model (k = 70; tau^2 estimator: ML)
## 
## tau^2 (estimated amount of total heterogeneity): 0.0044 (SE = 0.0008)
## tau (square root of estimated tau^2 value):      0.0660
## I^2 (total heterogeneity / total variability):   99.37%
## H^2 (total variability / sampling variability):  159.19
## 
## Test for Heterogeneity:
## Q(df = 69) = 3896.6240, p-val < .0001
## 
## Model Results:
## 
## estimate      se    tval  df    pval   ci.lb   ci.ub     
##   0.0265  0.0086  3.0758  69  0.0030  0.0093  0.0437  ** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

model_step_1_var <- rma(yi = sig_time,
              sei = sig_time_sd,
              data = total_reg_results,
              method = "ML",
              test = "knha")

model_step_1_var

## 
## Random-Effects Model (k = 70; tau^2 estimator: ML)
## 
## tau^2 (estimated amount of total heterogeneity): 0.0009 (SE = 0.0002)
## tau (square root of estimated tau^2 value):      0.0300
## I^2 (total heterogeneity / total variability):   98.50%
## H^2 (total variability / sampling variability):  66.55
## 
## Test for Heterogeneity:
## Q(df = 69) = 1444.9153, p-val < .0001
## 
## Model Results:
## 
## estimate      se     tval  df    pval    ci.lb    ci.ub    
##  -0.0096  0.0039  -2.4697  69  0.0160  -0.0173  -0.0018  * 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

model_step_2_att <- rma(yi = mu_time,
              sei = mu_time_sd,
              data = total_reg_results,
              mods = ~ attitude_cat,
              method = "ML",
              test = "knha")

summary(model_step_2_att)

## 
## Mixed-Effects Model (k = 70; tau^2 estimator: ML)
## 
##    logLik   deviance        AIC        BIC       AICc   
##   85.6804   348.9783  -163.3608  -154.3668  -162.7454   
## 
## tau^2 (estimated amount of residual heterogeneity):     0.0042 (SE = 0.0008)
## tau (square root of estimated tau^2 value):             0.0649
## I^2 (residual heterogeneity / unaccounted variability): 99.33%
## H^2 (unaccounted variability / sampling variability):   149.30
## R^2 (amount of heterogeneity accounted for):            3.41%
## 
## Test for Residual Heterogeneity:
## QE(df = 67) = 2852.9376, p-val < .0001
## 
## Test of Moderators (coefficients 2:3):
## F(df1 = 2, df2 = 67) = 0.6171, p-val = 0.5425
## 
## Model Results:
## 
##                        estimate      se    tval  df    pval    ci.lb   ci.ub    
## intrcpt                  0.0153  0.0208  0.7334  67  0.4659  -0.0263  0.0568    
## attitude_cataffective    0.0312  0.0297  1.0534  67  0.2959  -0.0280  0.0904    
## attitude_catcognitive    0.0089  0.0234  0.3827  67  0.7032  -0.0377  0.0556    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

model_step_2_var <- rma(yi = sig_time,
              sei = sig_time_sd,
              data = total_reg_results,
              mods = ~ attitude_cat, # cogn as ref cat, sig
              method = "ML",
              test = "knha")

model_step_2_var

## 
## Mixed-Effects Model (k = 70; tau^2 estimator: ML)
## 
## tau^2 (estimated amount of residual heterogeneity):     0.0008 (SE = 0.0002)
## tau (square root of estimated tau^2 value):             0.0285
## I^2 (residual heterogeneity / unaccounted variability): 98.28%
## H^2 (unaccounted variability / sampling variability):   58.29
## R^2 (amount of heterogeneity accounted for):            9.95%
## 
## Test for Residual Heterogeneity:
## QE(df = 67) = 1247.9526, p-val < .0001
## 
## Test of Moderators (coefficients 2:3):
## F(df1 = 2, df2 = 67) = 3.4755, p-val = 0.0366
## 
## Model Results:
## 
##                        estimate      se     tval  df    pval    ci.lb    ci.ub    
## intrcpt                  0.0102  0.0089   1.1423  67  0.2574  -0.0076   0.0280    
## attitude_cataffective   -0.0156  0.0129  -1.2095  67  0.2307  -0.0414   0.0102    
## attitude_catcognitive   -0.0260  0.0100  -2.5912  67  0.0117  -0.0460  -0.0060  * 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

model_step_3_att <- rma(yi = mu_time,
              sei = mu_time_sd,
              data = total_reg_results,
              mods = ~ attitude_cat + mean_year_centered,
              method = "ML",
              test = "knha")

model_step_3_att

## 
## Mixed-Effects Model (k = 70; tau^2 estimator: ML)
## 
## tau^2 (estimated amount of residual heterogeneity):     0.0042 (SE = 0.0008)
## tau (square root of estimated tau^2 value):             0.0646
## I^2 (residual heterogeneity / unaccounted variability): 99.32%
## H^2 (unaccounted variability / sampling variability):   146.29
## R^2 (amount of heterogeneity accounted for):            4.19%
## 
## Test for Residual Heterogeneity:
## QE(df = 66) = 2660.4810, p-val < .0001
## 
## Test of Moderators (coefficients 2:4):
## F(df1 = 3, df2 = 66) = 0.4952, p-val = 0.6869
## 
## Model Results:
## 
##                        estimate      se    tval  df    pval    ci.lb   ci.ub    
## intrcpt                  0.0182  0.0217  0.8382  66  0.4049  -0.0252  0.0616    
## attitude_cataffective    0.0239  0.0332  0.7186  66  0.4749  -0.0424  0.0901    
## attitude_catcognitive    0.0064  0.0240  0.2678  66  0.7897  -0.0415  0.0544    
## mean_year_centered       0.0006  0.0012  0.5099  66  0.6118  -0.0018  0.0031    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

model_step_4_att <- rma(yi = mu_time,
              sei = mu_time_sd,
              data = total_reg_results,
              mods = ~ attitude_cat + mean_year_centered + ambiguous + national + or_scale + pec_miss + data + nr_waves,
              method = "ML",
              test = "knha")

## Warning: 9 studies with NAs omitted from model fitting.

## Warning: Redundant predictors dropped from the model.

model_step_4_att

## 
## Mixed-Effects Model (k = 61; tau^2 estimator: ML)
## 
## tau^2 (estimated amount of residual heterogeneity):     0.0032 (SE = 0.0006)
## tau (square root of estimated tau^2 value):             0.0565
## I^2 (residual heterogeneity / unaccounted variability): 99.13%
## H^2 (unaccounted variability / sampling variability):   114.47
## R^2 (amount of heterogeneity accounted for):            32.49%
## 
## Test for Residual Heterogeneity:
## QE(df = 47) = 1960.4389, p-val < .0001
## 
## Test of Moderators (coefficients 2:14):
## F(df1 = 13, df2 = 47) = 1.6737, p-val = 0.0988
## 
## Model Results:
## 
##                        estimate      se     tval  df    pval    ci.lb    ci.ub     
## intrcpt                  0.1583  0.0631   2.5093  47  0.0156   0.0314   0.2852   * 
## attitude_cataffective   -0.0115  0.0425  -0.2699  47  0.7884  -0.0970   0.0740     
## attitude_catcognitive   -0.0067  0.0256  -0.2622  47  0.7943  -0.0583   0.0449     
## mean_year_centered       0.0014  0.0024   0.5861  47  0.5606  -0.0035   0.0063     
## ambiguousYes            -0.0219  0.0253  -0.8629  47  0.3926  -0.0728   0.0291     
## national                -0.0130  0.0841  -0.1540  47  0.8783  -0.1822   0.1563     
## or_scale                 0.0175  0.0074   2.3749  47  0.0217   0.0027   0.0324   * 
## pec_miss                -0.8852  0.6657  -1.3298  47  0.1900  -2.2244   0.4540     
## dataESS                 -0.1184  0.0641  -1.8466  47  0.0711  -0.2474   0.0106   . 
## dataEVS                 -0.1154  0.0755  -1.5286  47  0.1331  -0.2673   0.0365     
## dataISSP                -0.0852  0.0304  -2.8001  47  0.0074  -0.1464  -0.0240  ** 
## dataLISS                -0.0739  0.1160  -0.6375  47  0.5269  -0.3072   0.1594     
## dataMOT                 -0.0889  0.0783  -1.1346  47  0.2623  -0.2465   0.0687     
## nr_waves                -0.0426  0.0127  -3.3547  47  0.0016  -0.0681  -0.0171  ** 
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

model_step_3_var <- rma(yi = sig_time,
              sei = sig_time_sd,
              data = total_reg_results,
              mods = ~  attitude_cat + mean_year_centered,
              method = "ML",
              test = "knha")

model_step_3_var

## 
## Mixed-Effects Model (k = 70; tau^2 estimator: ML)
## 
## tau^2 (estimated amount of residual heterogeneity):     0.0008 (SE = 0.0002)
## tau (square root of estimated tau^2 value):             0.0285
## I^2 (residual heterogeneity / unaccounted variability): 98.26%
## H^2 (unaccounted variability / sampling variability):   57.43
## R^2 (amount of heterogeneity accounted for):            10.19%
## 
## Test for Residual Heterogeneity:
## QE(df = 66) = 1233.5279, p-val < .0001
## 
## Test of Moderators (coefficients 2:4):
## F(df1 = 3, df2 = 66) = 2.3349, p-val = 0.0818
## 
## Model Results:
## 
##                        estimate      se     tval  df    pval    ci.lb    ci.ub    
## intrcpt                  0.0092  0.0094   0.9754  66  0.3329  -0.0096   0.0279    
## attitude_cataffective   -0.0132  0.0145  -0.9080  66  0.3672  -0.0422   0.0158    
## attitude_catcognitive   -0.0252  0.0103  -2.4363  66  0.0175  -0.0458  -0.0045  * 
## mean_year_centered      -0.0002  0.0005  -0.3786  66  0.7062  -0.0013   0.0009    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

model_step_4_var <- rma(yi = sig_time,
              sei = sig_time_sd,
              data = total_reg_results,
              mods = ~  attitude_cat +  mean_year_centered + ambiguous + national + or_scale + pec_miss + data + nr_waves,
              method = "ML",
              test = "knha")

## Warning: 9 studies with NAs omitted from model fitting.

## Warning: Redundant predictors dropped from the model.

model_step_4_var

## 
## Mixed-Effects Model (k = 61; tau^2 estimator: ML)
## 
## tau^2 (estimated amount of residual heterogeneity):     0.0006 (SE = 0.0001)
## tau (square root of estimated tau^2 value):             0.0236
## I^2 (residual heterogeneity / unaccounted variability): 97.53%
## H^2 (unaccounted variability / sampling variability):   40.44
## R^2 (amount of heterogeneity accounted for):            30.37%
## 
## Test for Residual Heterogeneity:
## QE(df = 47) = 792.3181, p-val < .0001
## 
## Test of Moderators (coefficients 2:14):
## F(df1 = 13, df2 = 47) = 1.2166, p-val = 0.2980
## 
## Model Results:
## 
##                        estimate      se     tval  df    pval    ci.lb   ci.ub    
## intrcpt                 -0.0179  0.0282  -0.6324  47  0.5302  -0.0747  0.0389    
## attitude_cataffective   -0.0050  0.0194  -0.2580  47  0.7975  -0.0440  0.0340    
## attitude_catcognitive   -0.0195  0.0111  -1.7538  47  0.0860  -0.0420  0.0029  . 
## mean_year_centered       0.0006  0.0010   0.5695  47  0.5718  -0.0015  0.0026    
## ambiguousYes             0.0100  0.0109   0.9240  47  0.3602  -0.0118  0.0319    
## national                -0.0487  0.0384  -1.2675  47  0.2112  -0.1260  0.0286    
## or_scale                -0.0023  0.0031  -0.7254  47  0.4718  -0.0085  0.0040    
## pec_miss                -0.0677  0.2920  -0.2320  47  0.8175  -0.6552  0.5197    
## dataESS                  0.0162  0.0281   0.5784  47  0.5658  -0.0402  0.0727    
## dataEVS                  0.0120  0.0316   0.3796  47  0.7059  -0.0516  0.0756    
## dataISSP                 0.0278  0.0131   2.1236  47  0.0390   0.0015  0.0540  * 
## dataLISS                 0.0658  0.0515   1.2790  47  0.2072  -0.0377  0.1694    
## dataMOT                  0.0630  0.0361   1.7458  47  0.0874  -0.0096  0.1356  . 
## nr_waves                 0.0080  0.0056   1.4323  47  0.1587  -0.0032  0.0193    
## 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Put them in a table
models <- list("M1" = model_step_1_att, "M2" = model_step_2_att, "M3" = model_step_3_att, "M4" = model_step_4_att, "M1" =  model_step_1_var, "M2"= model_step_2_var, "M3"=  model_step_3_var, "M4" = model_step_4_var)

modelsummary(models, output = "kableExtra", statistic = 'std.error', stars = TRUE, shape = term ~ model + statistic, title = "Appendix X. Meta-regression excluding outlier", fmt = fmt_statistic(estimate = 3, std.error =3), coef_rename = c("overall" ="Intercept", 
           "intercept" = "Intercept", 
           "attitude_cataffective" = "Affective attitude (ref = beh)",
           "attitude_catcognitive" = "Cognitive attitude", 
           "mean_year_centered" =   "Mean year centered", 
           "ambiguousYes" = "Ambiguous (ref = no)",
           "national" = "National (ref = no)", 
           "or_scale" = "Original scale", 
            "pec_miss" = "Perc. missings", 
           "dataESS"= "ESS (ref = EB)", 
           "dataEVS"= "EVS", 
           "dataISSP" = "ISSP", 
           "dataLISS"= "LISS", 
           "dataMOT" = "MOT", 
           "nr_waves" = "Nr. of waves"), gof_map = NA) %>%
  kable_classic_2(html_font = "Times", fixed_thead = T, full_width = F) %>%
   add_header_above(c(" " = 1, "Mu" = 8, "Sigma" = 8))#%>%

Appendix X. Meta-regression excluding outlier

	Mu								Sigma
	M1		M2		M3		M4		M1		M2		M3		M4
	Est.	S.E.	Est.	S.E.	Est.	S.E.	Est.	S.E.	Est.	S.E.	Est.	S.E.	Est.	S.E.	Est.	S.E.
Intercept	0.026**	0.009	0.015	0.021	0.018	0.022	0.158*	0.063	-0.010*	0.004	0.010	0.009	0.009	0.009	-0.018	0.028
Affective attitude (ref = beh)			0.031	0.030	0.024	0.033	-0.011	0.042			-0.016	0.013	-0.013	0.015	-0.005	0.019
Cognitive attitude			0.009	0.023	0.006	0.024	-0.007	0.026			-0.026*	0.010	-0.025*	0.010	-0.020+	0.011
Mean year centered					0.001	0.001	0.001	0.002					0.000	0.001	0.001	0.001
Ambiguous (ref = no)							-0.022	0.025							0.010	0.011
National (ref = no)							-0.013	0.084							-0.049	0.038
Original scale							0.018*	0.007							-0.002	0.003
Perc. missings							-0.885	0.666							-0.068	0.292
ESS (ref = EB)							-0.118+	0.064							0.016	0.028
EVS							-0.115	0.076							0.012	0.032
ISSP							-0.085**	0.030							0.028*	0.013
LISS							-0.074	0.116							0.066	0.051
MOT							-0.089	0.078							0.063+	0.036
Nr. of waves							-0.043**	0.013							0.008	0.006
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001

    #save_kable("./output/appendix/meta_regression_table_excluding_outlier.html")