Meta analysis for independent variables

In this script, I perform meta-analyses for the models with the independent variables

rm(list=ls())
library(meta)
library(metafor)
library(here)
library(dplyr)
library(kableExtra)
set.seed(1)

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

# Effects of mu's 
# 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)

summary(gender_gen)

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

summary(age_gen)

# Educational level. 
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)

summary(edu_med_gen)

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)

summary(edu_high_gen)

# Now for variance
# 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)

summary(gender_var_gen)

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

summary(age_var_gen)

# 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)

summary(edu_var_med_gen)

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)

summary(edu_var_high_gen)

# Perform the meta-regression with meta-level variables on independent variables to see if there is something interesting

# For this dataset, the meta-level variables were not created yet, so that's needed first
total_indep_var_results$data[total_indep_var_results$dep_var=="climate5"] <- "EVS"
total_indep_var_results$data[total_indep_var_results$dep_var=="cause" | total_indep_var_results$dep_var == "cause" | total_indep_var_results$dep_var == "worry"] <- "ESS"
total_indep_var_results$data[total_indep_var_results$dep_var=="worried" | total_indep_var_results$dep_var == "frontrunner" | total_indep_var_results$dep_var == "worry_future" | total_indep_var_results$dep_var == "resp_citiz" | total_indep_var_results$dep_var == "dk_start" | total_indep_var_results$dep_var == "do_gov" | total_indep_var_results$dep_var == "buss_help" | total_indep_var_results$dep_var == "min_contr" | total_indep_var_results$dep_var == "human_resp"] <- "I&O"
total_indep_var_results$data[total_indep_var_results$dep_var=="worry" | total_indep_var_results$dep_var == "lifeharm" | total_indep_var_results$dep_var == "progharm" | total_indep_var_results$dep_var == "econprotect" | total_indep_var_results$dep_var == "growharm" | total_indep_var_results$dep_var == "willing_price" | total_indep_var_results$dep_var == "willing_tax" | total_indep_var_results$dep_var == "willing_living" | total_indep_var_results$dep_var == "dodiff" | total_indep_var_results$dep_var=="do_right" | total_indep_var_results$dep_var == "people_decide" | total_indep_var_results$dep_var == "moreimp" | total_indep_var_results$dep_var == "othersame" | total_indep_var_results$dep_var == "exag" | total_indep_var_results$dep_var == "country_effort"]  <- "ISSP"
total_indep_var_results$data[is.na(total_indep_var_results$data)] <- "EB"


total_indep_var_results$data[total_indep_var_results$dep_var=="cchange_mot" | total_indep_var_results$dep_var == "worried_mot" | total_indep_var_results$dep_var == "futuregen" | total_indep_var_results$dep_var == "nowor" | total_indep_var_results$dep_var == "ontime" | total_indep_var_results$dep_var == "gov" | total_indep_var_results$dep_var == "resp_gov" | total_indep_var_results$dep_var == "resp_comp" | total_indep_var_results$dep_var == "resp_mkb" | total_indep_var_results$dep_var=="resp_citiz_mot" | total_indep_var_results$dep_var == "resp_you" | total_indep_var_results$dep_var == "pers_resp_mot" | total_indep_var_results$dep_var == "sust_choice" | total_indep_var_results$dep_var == "contr" | total_indep_var_results$dep_var == "energy" | total_indep_var_results$dep_var == "energy" | total_indep_var_results$dep_var == "noidea" | total_indep_var_results$dep_var == "motiv"]  <- "MOT"
total_indep_var_results$data[total_indep_var_results$dep_var=="fut_gen_socon"] <- "SOCON"
total_indep_var_results$data[total_indep_var_results$dep_var=="lifestyle"] <- "LISS"

total_indep_var_results$dep_var[total_indep_var_results$dep_var == "worry" & total_indep_var_results$data == "ISSP"] <- "worry_issp"

total_indep_var_results$attitude_cat[total_indep_var_results$dep_var == "dodiff" 
                               | total_indep_var_results$dep_var == "pers_resp"
                               | total_indep_var_results$dep_var == "worry"
                               | total_indep_var_results$dep_var == "worried"
                               | total_indep_var_results$dep_var == "worry_future"
                               | total_indep_var_results$dep_var == "frontrunner"
                               | total_indep_var_results$dep_var == "min_contr"
                               | total_indep_var_results$dep_var == "worried_mot"
                               | total_indep_var_results$dep_var == "futuregen"
                               | total_indep_var_results$dep_var == "nowor"
                               | total_indep_var_results$dep_var == "motivr"
                               | total_indep_var_results$dep_var == "fut_gen_socon"
                               | total_indep_var_results$dep_var == "pers_resp_mot"] <- "affective"

total_indep_var_results$attitude_cat[total_indep_var_results$dep_var == "willing_price" 
                               | total_indep_var_results$dep_var == "willing_tax"
                               | total_indep_var_results$dep_var == "willing_living"
                               | total_indep_var_results$dep_var == "do_right"
                               | total_indep_var_results$dep_var == "people_decide"
                               | total_indep_var_results$dep_var == "climate5"
                               | total_indep_var_results$dep_var == "prsaction"
                               | total_indep_var_results$dep_var == "cc_prsact"
                               | total_indep_var_results$dep_var == "buyprod"
                                | total_indep_var_results$dep_var == "sust_choice"
                                | total_indep_var_results$dep_var == "energy"
                                | total_indep_var_results$dep_var == "lifestyle"] <- "behavioral"

total_indep_var_results$attitude_cat[total_indep_var_results$dep_var == "worry_issp" 
                               | total_indep_var_results$dep_var == "lifeharm"
                               | total_indep_var_results$dep_var == "progharm"
                               | total_indep_var_results$dep_var == "econprotect"
                               | total_indep_var_results$dep_var == "growharm"
                               | total_indep_var_results$dep_var == "bus_decide"
                               | total_indep_var_results$dep_var == "moreimp"
                               | total_indep_var_results$dep_var == "othersame"
                               | total_indep_var_results$dep_var == "exag"
                               | total_indep_var_results$dep_var == "country_effort" 
                               | total_indep_var_results$dep_var == "cause"
                               | total_indep_var_results$dep_var == "resp_citiz"
                               | total_indep_var_results$dep_var == "dk_start"
                               | total_indep_var_results$dep_var == "do_gov"
                               | total_indep_var_results$dep_var == "buss_help"
                               | total_indep_var_results$dep_var == "human_resp"
                               | total_indep_var_results$dep_var == "env_ec_stat"
                               | total_indep_var_results$dep_var == "env_prsimp"
                                | total_indep_var_results$dep_var == "envp_eg"
                               | total_indep_var_results$dep_var == "effr_eg"
                               | total_indep_var_results$dep_var == "cchange"
                               | total_indep_var_results$dep_var == "cchange2"
                               | total_indep_var_results$dep_var == "cchangetot" 
                               | total_indep_var_results$dep_var == "ccpercept"
                               | total_indep_var_results$dep_var == "env_quallife"
                               | total_indep_var_results$dep_var == "doprot_comp"
                               | total_indep_var_results$dep_var == "doprot_region"
                               | total_indep_var_results$dep_var == "doprot_natgov"
                               | total_indep_var_results$dep_var == "doprot_city"
                               | total_indep_var_results$dep_var == "doprot_citiz"
                               | total_indep_var_results$dep_var == "doprot_eu"
                              | total_indep_var_results$dep_var == "cc_unstop"
                               | total_indep_var_results$dep_var == "cc_exag"
                               | total_indep_var_results$dep_var == "cc_poseu"
                               | total_indep_var_results$dep_var == "role_ind"
                               | total_indep_var_results$dep_var == "big_pol"
                               | total_indep_var_results$dep_var == "eff_daily"
                               | total_indep_var_results$dep_var == "pers_imp"
                              | total_indep_var_results$dep_var == "cchange_mot"
                              | total_indep_var_results$dep_var == "ontime"
                              | total_indep_var_results$dep_var == "gov"
                              | total_indep_var_results$dep_var == "resp_gov"
                              | total_indep_var_results$dep_var == "resp_comp"
                              | total_indep_var_results$dep_var == "resp_mkb"
                              | total_indep_var_results$dep_var == "resp_citiz_mot"
                              | total_indep_var_results$dep_var == "resp_you"
                              | total_indep_var_results$dep_var == "contr"
                              | total_indep_var_results$dep_var == "noidea"] <- "cognitive"

#Create a variable that indicates whether the dependent variable can be interpreted in 2 ways (ambiguous)

total_indep_var_results$ambiguous[total_indep_var_results$dep_var != "dodiff" 
                               | total_indep_var_results$dep_var != "frontrunner"
                               | total_indep_var_results$dep_var != "min_contr"
                               | total_indep_var_results$dep_var != "people_decide"
                               | total_indep_var_results$dep_var != "econprotect"
                               | total_indep_var_results$dep_var != "growharm"
                               | total_indep_var_results$dep_var != "bus_decide"
                            | total_indep_var_results$dep_var != "othersame"
                            | total_indep_var_results$dep_var != "resp_citiz"
                            | total_indep_var_results$dep_var != "dk_start"
                            | total_indep_var_results$dep_var != "buss_help"
                            | total_indep_var_results$dep_var != "envp_eg"
                            | total_indep_var_results$dep_var != "effr_eg"
                            | total_indep_var_results$dep_var != "env_quallife"
                            | total_indep_var_results$dep_var != "cc_unstop"
                            | total_indep_var_results$dep_var != "cc_poseu"
                            | total_indep_var_results$dep_var != "big_pol"
                            | total_indep_var_results$dep_var == "noidea"] <- "No"
total_indep_var_results$ambiguous[total_indep_var_results$dep_var == "dodiff" 
                               | total_indep_var_results$dep_var == "frontrunner"
                               | total_indep_var_results$dep_var == "min_contr"
                               | total_indep_var_results$dep_var == "people_decide"
                               | total_indep_var_results$dep_var == "econprotect"
                               | total_indep_var_results$dep_var == "growharm"
                               | total_indep_var_results$dep_var == "bus_decide"
                            | total_indep_var_results$dep_var == "othersame"
                            | total_indep_var_results$dep_var == "resp_citiz"
                            | total_indep_var_results$dep_var == "dk_start"
                            | total_indep_var_results$dep_var == "buss_help"
                            | total_indep_var_results$dep_var == "envp_eg"
                            | total_indep_var_results$dep_var == "effr_eg"
                            | total_indep_var_results$dep_var == "env_quallife"
                            | total_indep_var_results$dep_var == "cc_unstop"
                            | total_indep_var_results$dep_var == "cc_poseu"
                            | total_indep_var_results$dep_var == "big_pol"
                            | total_indep_var_results$dep_var == "noidea"] <- "Yes"



# This time, year is a bit more difficult, so I have to assign that per dep_var
# And the average year 
total_indep_var_results$mean_year[total_indep_var_results$dep_var == "role_ind"] <- 2012
total_indep_var_results$mean_year[total_indep_var_results$dep_var == "big_pol"] <- 2012
total_indep_var_results$mean_year[total_indep_var_results$dep_var == "cc_unstop"] <- 2008
total_indep_var_results$mean_year[total_indep_var_results$dep_var == "cc_exag"] <- 2008
total_indep_var_results$mean_year[total_indep_var_results$dep_var == "cc_poseu"] <- 2008
total_indep_var_results$mean_year[total_indep_var_results$dep_var == "cc_prsact"] <- 2008
total_indep_var_results$mean_year[total_indep_var_results$dep_var == "cc_poseu"] <- 2008
total_indep_var_results$mean_year[total_indep_var_results$dep_var == "ccpercept"] <- 2014
total_indep_var_results$mean_year[total_indep_var_results$dep_var == "cchange"] <- 2015
total_indep_var_results$mean_year[total_indep_var_results$dep_var == "cchange2"] <- 2015
total_indep_var_results$mean_year[total_indep_var_results$dep_var == "cchangetot"] <- 2015
total_indep_var_results$mean_year[total_indep_var_results$dep_var == "envp_eg"] <- 2013
total_indep_var_results$mean_year[total_indep_var_results$dep_var == "effr_eg"] <- 2013
total_indep_var_results$mean_year[total_indep_var_results$dep_var == "buyprod"] <- 2010
total_indep_var_results$mean_year[total_indep_var_results$dep_var == "doprot_natgov"] <- 2013
total_indep_var_results$mean_year[total_indep_var_results$dep_var == "doprot_eu"] <- 2013
total_indep_var_results$mean_year[total_indep_var_results$dep_var == "doprot_region"] <- 2013
total_indep_var_results$mean_year[total_indep_var_results$dep_var == "doprot_comp"] <- 2013
total_indep_var_results$mean_year[total_indep_var_results$dep_var == "doprot_citiz"] <- 2013
total_indep_var_results$mean_year[total_indep_var_results$dep_var == "doprot_city"] <- 2015
total_indep_var_results$mean_year[total_indep_var_results$dep_var == "eff_daily"] <- 2012
total_indep_var_results$mean_year[total_indep_var_results$dep_var == "pers_imp"] <- 2009
total_indep_var_results$mean_year[total_indep_var_results$dep_var == "env_quallife"] <- 2009
total_indep_var_results$mean_year[total_indep_var_results$dep_var == "env_ec_stat"] <- 1991
total_indep_var_results$mean_year[total_indep_var_results$dep_var == "env_prsimp"] <- 1991
total_indep_var_results$mean_year[total_indep_var_results$dep_var == "prsaction"] <- 2014


total_indep_var_results$mean_year[total_indep_var_results$data == "I&O"] <- 2020
total_indep_var_results$mean_year[total_indep_var_results$data == "EVS"] <- 1999
total_indep_var_results$mean_year[total_indep_var_results$data == "ESS"] <- 2018
total_indep_var_results$mean_year[total_indep_var_results$data == "MOT"] <- 2020
total_indep_var_results$mean_year[total_indep_var_results$data == "SOCON"] <- 2021
total_indep_var_results$mean_year[total_indep_var_results$data == "LISS"] <- 2020
# For the ISSP, have to specify it per dependent variable
total_indep_var_results$mean_year[total_indep_var_results$dep_var == "exag"] <- 2005
total_indep_var_results$mean_year[total_indep_var_results$dep_var == "moreimp"] <- 2005
total_indep_var_results$mean_year[total_indep_var_results$dep_var == "othersame"] <- 2005
total_indep_var_results$mean_year[total_indep_var_results$dep_var == "country_effort"] <- 2005
total_indep_var_results$mean_year[total_indep_var_results$dep_var == "growharm"] <- 2005
total_indep_var_results$mean_year[is.na(total_indep_var_results$mean_year)] <- 2002

total_indep_var_results$mean_year_centered <- total_indep_var_results$mean_year - mean(total_indep_var_results$mean_year)

save(total_indep_var_results, file= "./data/meta_analysis/total_indep_var_results_new.RData" )

m1_gender <- rma(yi = mu_sex_est,
              sei = mu_sex_sd,
              data = total_indep_var_results,
              mods = ~ mean_year_centered + attitude_cat + ambiguous,
              method = "REML",
              test = "knha")

m1_gender 

m1_age <- rma(yi = mu_age_est,
              sei = mu_age_sd,
              data = total_indep_var_results,
              mods = ~ mean_year_centered + attitude_cat + ambiguous,
              method = "REML",
              test = "knha")

m1_age

m1_isced_med <- rma(yi = mu_isced_med_est,
              sei = mu_isced_med_sd,
              data = total_indep_var_results,
              mods = ~ mean_year_centered + attitude_cat + ambiguous,
              method = "REML",
              test = "knha")

m1_isced_med

m1_isced_high <- rma(yi = mu_isced_high_est,
              sei = mu_isced_high_sd,
              data = total_indep_var_results,
              mods = ~ mean_year_centered + attitude_cat + ambiguous,
              method = "REML",
              test = "knha")

m1_isced_high

# Same for variance/polarization
m1_gender_var <- rma(yi = sig_sex_est,
              sei = sig_sex_sd,
              data = total_indep_var_results,
              mods = ~ mean_year_centered + attitude_cat + ambiguous,
              method = "REML",
              test = "knha")

m1_gender_var 

m1_age_var <- rma(yi = sig_age_est,
              sei = sig_age_sd,
              data = total_indep_var_results,
              mods = ~ mean_year_centered + attitude_cat + ambiguous,
              method = "REML",
              test = "knha")

m1_age_var

m1_isced_med_var <- rma(yi = sig_isced_med_est,
              sei = sig_isced_med_sd,
              data = total_indep_var_results,
              mods = ~ mean_year_centered + attitude_cat + ambiguous,
              method = "REML",
              test = "knha")

m1_isced_med_var

m1_isced_high_var <- rma(yi = sig_isced_high_est,
              sei = sig_isced_high_sd,
              data = total_indep_var_results,
              mods = ~ mean_year_centered + attitude_cat + ambiguous,
              method = "REML",
              test = "knha")

m1_isced_high_var

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

# Make a table of the above regression output
indep_df <- data.frame(matrix(ncol = 14, nrow = 5))
colnames(indep_df) <- c("Mu/Sigma", "Variable", "Estimate", "SE", "p", "Estimate_1", "SE_1", "p_1", "Estimate_2", "SE_2", "p_2", "Estimate_3", "SE_3", "p_3")
indep_df$Variable <- c("Intercept", "Mean year (centered around mean)", "Attitude affective (ref = behavioral)", "Attitude cognitive", "Ambiguous (ref = no)")
indep_df$`Mu/Sigma` <- c("Mean attitudes", "", "", "", "")

indep_df$Estimate <- as.numeric(m1_gender$beta) #estimate
indep_df$SE <- m1_gender$se #se indeed but as a row
indep_df$p <- m1_gender$pval # same as for SE

indep_df$Estimate_1 <- as.numeric(m1_age$beta) #estimate
indep_df$SE_1 <- m1_age$se #se indeed but as a row
indep_df$p_1 <- m1_age$pval # same as for SE

indep_df$Estimate_2 <- as.numeric(m1_isced_med$beta) #estimate
indep_df$SE_2 <- m1_isced_med$se #se indeed but as a row
indep_df$p_2 <- m1_isced_med$pval # same as for SE

indep_df$Estimate_3 <- as.numeric(m1_isced_high$beta) #estimate
indep_df$SE_3 <- m1_isced_high$se #se indeed but as a row
indep_df$p_3 <- m1_isced_high$pval # same as for SE

indep_df$p_text <- case_when(
indep_df$p >= 0.05 ~ "",
indep_df$p < 0.001 ~ "***",
indep_df$p < 0.01 ~ "**",
indep_df$p < 0.05 ~ "*"
)

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

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

indep_df$p_text_3 <- case_when(
indep_df$p_3 >= 0.05 ~ "",
indep_df$p_3 < 0.001 ~ "***",
indep_df$p_3 < 0.01 ~ "**",
indep_df$p_3 < 0.05 ~ "*"
)

# I also want the R2, I2 and QE in the table
model_1_fit  <- data.frame(matrix(ncol = 14, nrow = 3))
colnames(model_1_fit) <- c("Mu/Sigma", "Variable", "Estimate", "SE", "p_text", "Estimate_1", "SE_1", "p_text_1", "Estimate_2", "SE_2", "p_text_2", "Estimate_3", "SE_3", "p_text_3")
model_1_fit$Variable <- c("R2", "I2", "QE")
model_1_fit$Estimate <- c(m1_gender$R2, m1_gender$I2, m1_gender$QE)
model_1_fit$Estimate_1 <- c(m1_age$R2, m1_age$I2, m1_age$QE)
model_1_fit$Estimate_2 <- c(m1_isced_med$R2, m1_isced_med$I2, m1_isced_med$QE)
model_1_fit$Estimate_3 <- c(m1_isced_high$R2, m1_isced_high$I2, m1_isced_high$QE)

indep_df <- select(indep_df, `Mu/Sigma`, Variable, Estimate, SE, p_text, Estimate_1, SE_1, p_text_1, Estimate_2, SE_2, p_text_2, Estimate_3, SE_3, p_text_3)
indep_df <- rbind(indep_df, model_1_fit)

# Variance
indep_var_df <- data.frame(matrix(ncol = 14, nrow = 5))
colnames(indep_var_df) <- c("Mu/Sigma", "Variable", "Estimate", "SE", "p", "Estimate_1", "SE_1", "p_1", "Estimate_2", "SE_2", "p_2", "Estimate_3", "SE_3", "p_3")
indep_var_df$Variable <- c("Intercept", "Mean year (centered around mean)", "Attitude affective (ref = behavioral)", "Attitude cognitive", "Ambiguous (ref = no)")
indep_var_df$`Mu/Sigma` <- c("Polarization in attitudes", "", "", "", "")

indep_var_df$Estimate <- as.numeric(m1_gender_var$beta) #estimate
indep_var_df$SE <- m1_gender_var$se #se indeed but as a row
indep_var_df$p <- m1_gender_var$pval # same as for SE

indep_var_df$Estimate_1 <- as.numeric(m1_age_var$beta) #estimate
indep_var_df$SE_1 <- m1_age_var$se #se indeed but as a row
indep_var_df$p_1 <- m1_age_var$pval # same as for SE

indep_var_df$Estimate_2 <- as.numeric(m1_isced_med_var$beta) #estimate
indep_var_df$SE_2 <- m1_isced_med_var$se #se indeed but as a row
indep_var_df$p_2 <- m1_isced_med_var$pval # same as for SE

indep_var_df$Estimate_3 <- as.numeric(m1_isced_high_var$beta) #estimate
indep_var_df$SE_3 <- m1_isced_high_var$se #se indeed but as a row
indep_var_df$p_3 <- m1_isced_high_var$pval # same as for SE

indep_var_df$p_text <- case_when(
indep_var_df$p >= 0.05 ~ "",
indep_var_df$p < 0.001 ~ "***",
indep_var_df$p < 0.01 ~ "**",
indep_var_df$p < 0.05 ~ "*"
)

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

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

indep_var_df$p_text_3 <- case_when(
indep_var_df$p_3 >= 0.05 ~ "",
indep_var_df$p_3 < 0.001 ~ "***",
indep_var_df$p_3 < 0.01 ~ "**",
indep_var_df$p_3 < 0.05 ~ "*"
)

# I also want the R2, I2 and QE in the table
model_1_fit_var  <- data.frame(matrix(ncol = 14, nrow = 3))
colnames(model_1_fit_var) <- c("Mu/Sigma", "Variable", "Estimate", "SE", "p_text", "Estimate_1", "SE_1", "p_text_1", "Estimate_2", "SE_2", "p_text_2", "Estimate_3", "SE_3", "p_text_3")
model_1_fit_var$Variable <- c("R2", "I2", "QE")
model_1_fit_var$Estimate <- c(m1_gender_var$R2, m1_gender_var$I2, m1_gender_var$QE)
model_1_fit_var$Estimate_1 <- c(m1_age_var$R2, m1_age_var$I2, m1_age_var$QE)
model_1_fit_var$Estimate_2 <- c(m1_isced_med_var$R2, m1_isced_med_var$I2, m1_isced_med_var$QE)
model_1_fit_var$Estimate_3 <- c(m1_isced_high_var$R2, m1_isced_high_var$I2, m1_isced_high_var$QE)

indep_var_df <- select(indep_var_df, `Mu/Sigma`, Variable, Estimate, SE, p_text, Estimate_1, SE_1, p_text_1, Estimate_2, SE_2, p_text_2, Estimate_3, SE_3, p_text_3)
indep_var_df <- rbind(indep_var_df, model_1_fit_var)

#Now a total dataframe with both models together, and then make some changes before making the table 
indep_df_total <- rbind(indep_df, indep_var_df)

indep_df_total <- dplyr::mutate(indep_df_total, across(.cols =c(-`Mu/Sigma`, -Variable, -p_text, -p_text_1, -p_text_2, -p_text_3), .fns = as.numeric))

numeric_cols <- sapply(indep_df_total, is.numeric)
indep_df_total[, numeric_cols] <- round(indep_df_total[, numeric_cols], digits = 3)

kable(indep_df_total, digits = 3, caption = "Appendix E. Meta-regression on climate change attitudes and polarization", 
      col.names = c( "", "Variables", "Estimate", "Std. Error", "", "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, "Gender" = 3, "Age" = 3, "Isced intermediate (ref=basic)" = 3, "Isced advanced (ref=basic)" = 3)) %>%
   footnote(
    general_title = "Note.",
    general = "*** = p < 0.001, ** = p < 0.01, * = p < 0.05",
    threeparttable = TRUE,
    footnote_as_chunk = TRUE
    ) %>%
    save_kable("./output/meta_regression_indepvar_table_NOVEMBER.html")

`