Extract coefficients of gamlss regressions

In this script I extract the coefficients for the time effects and intercepts and put them in a dataset together, to use them as input for the meta-regression.

rm(list=ls())
library(tidyverse)
library(dplyr)
library(gamlss)

ESS

# First load the empty lists
load("./data/final_data/regression_outputs/ess_list_empty_w.RData")

summary(ess_list_empty[[1]])

# Store the results in a new dataframe
ess_empty_df <- data.frame(dep_var = character())

        # Loop over the list of models
for (i in seq_along(ess_list_empty)) {

  # Extract the dependent variable so that we know for which dep var the effects are
  dep_var <- as.character(ess_list_empty[[i]]$mu.terms[[2]])
  
    # Summarize the model
  sum <- summary(ess_list_empty[[i]])
  
  # Extract the intercept and standard deviation for mu and sigma
  mu_int <- sum[[1]]
  mu_sd <- sum[[1,2]]
  sigma_int <- sum[[3]]
  sigma_sd <- sum[[3,2]]
  mu_time <- sum [[2]]
  mu_time_sd <- sum[[2,2]]
  sig_time <- sum[[4]]
  sig_time_sd <- sum[[4,2]]
  
  # Put the results 
  ess_empty_df <- rbind(ess_empty_df, data.frame(dep_var = dep_var, mu_intercept = mu_int, mu_sd = mu_sd, sigma_intercept = sigma_int, sigma_sd = sigma_sd, mu_time = mu_time, mu_time_sd = mu_time_sd, sig_time = sig_time, sig_time_sd = sig_time_sd))
}

# Of course this also has to be done for the other lists of the other models
load("./data/final_data/regression_outputs/ess_list_preds_w_new.RData")

# This does mean that the estimates of the sigma are in a different place
sum <- summary(ess_list_preds[[1]])

# Store the results in a new dataframe
ess_preds_df <- data.frame(dep_var = character())

# Loop over the list of models
for (i in seq_along(ess_list_preds)) {

  # Extract the dependent variable so that we know for which dep var the effects are
  dep_var <- as.character(ess_list_preds[[i]]$mu.terms[[2]])
  
    # Summarize the model
  sum <- summary(ess_list_preds[[i]])
  
  # Extract the intercept and standard deviation for mu and sigma
  mu_int <- sum[[1]]
  mu_sd <- sum[[1,2]]
  sigma_int <- sum[[7]]
  sigma_sd <- sum[[7,2]]
  mu_time <- sum [[2]]
  mu_time_sd <- sum[[2,2]]
  sig_time <- sum[[8]]
  sig_time_sd <- sum[[8,2]]
  
  # Put the results 
  ess_preds_df <- rbind(ess_preds_df, data.frame(dep_var = dep_var, mu_intercept_pred = mu_int, mu_sd_pred = mu_sd, sigma_intercept_pred = sigma_int, sigma_sd_pred = sigma_sd, mu_time_pred = mu_time, mu_time_sd_pred = mu_time_sd, sig_time_pred = sig_time, sig_time_sd_pred = sig_time_sd))
}

# And then for the interaction model
load("./data/final_data/regression_outputs/ess_list_interactions_w_new.RData")

# Store the results in a new dataframe
ess_interactions_df <- data.frame(dep_var = character())

# Loop over the number of dep vars in the list
for (i in seq_along(ess_list_interactions)) {

  # Extract the dependent variable so that we know for which dep var the effects are
  dep_var <- as.character(ess_list_interactions[[i]]$mu.terms[[2]])
  
    # Summarize the model
  sum <- summary(ess_list_interactions[[i]])
  
  # Extract the info I need from the summary
  mu_int <- sum[[1]]
  mu_sd <- sum[[1,2]]
  sigma_int <- sum[[11]]
  sigma_sd <- sum[[11,2]]
  mu_time <- sum [[2]]
  mu_time_sd <- sum[[2,2]]
  sig_time <- sum[[12]]
  sig_time_sd <- sum[[12,2]]
  
  # Put the results in a dataframe
  ess_interactions_df <- rbind(ess_interactions_df, data.frame(dep_var = dep_var, mu_intercept_int = mu_int, mu_sd_int = mu_sd, sigma_intercept_int = sigma_int, sigma_sd_int = sigma_sd, mu_time_int = mu_time, mu_time_sd_int = mu_time_sd, sig_time_int = sig_time, sig_time_sd_int = sig_time_sd))
}

ess_results_df <- left_join(ess_empty_df, ess_preds_df)
ess_results_df <- left_join(ess_results_df, ess_interactions_df)
ess_results_df$first_year <- 2016
ess_results_df$mean_year <- 2018
ess_results_df$data <- "ESS"

save(ess_results_df, file="./data/final_data/regression_outputs/ess_results_df_w_new.RData")

EVS

# Now the same for the EVS
rm(list=ls())
load("./data/final_data/regression_outputs/evs_list_empty_w.RData")

# Store the results in a new dataframe
evs_empty_df <- data.frame(dep_var = character())

# Loop over the list for all the dep vars
for (i in seq_along(evs_list_empty)) {

  # Extract the dependent variable so that we know for which dep var the effects are
  dep_var <- as.character(evs_list_empty[[i]]$mu.terms[[2]])
  
    # Summarize the model
  sum <- summary(evs_list_empty[[i]])
  
  # Extract the info
  mu_int <- sum[[1]]
  mu_sd <- sum[[1,2]]
  sigma_int <- sum[[3]]
  sigma_sd <- sum[[3,2]]
  mu_time <- sum [[2]]
  mu_time_sd <- sum[[2,2]]
  sig_time <- sum[[4]]
  sig_time_sd <- sum[[4,2]]
  
  # Put the results in a dataframe
  evs_empty_df <- rbind(evs_empty_df, data.frame(dep_var = dep_var, mu_intercept = mu_int, mu_sd = mu_sd, sigma_intercept = sigma_int, sigma_sd = sigma_sd, mu_time = mu_time, mu_time_sd = mu_time_sd, sig_time = sig_time, sig_time_sd = sig_time_sd))
}

# And the predictor model
load("./data/final_data/regression_outputs/evs_list_preds_w_new.RData")

# Store the results in a new dataframe
evs_preds_df <- data.frame(dep_var = character())

# Loop over the list of models
for (i in seq_along(evs_list_preds)) {

  # Extract the dependent variable so that we know for which dep var the effects are
  dep_var <- as.character(evs_list_preds[[i]]$mu.terms[[2]])
  
    # Summarize the model
  sum <- summary(evs_list_preds[[i]])
  
  # Extract the intercept and standard deviation for mu and sigma
  mu_int <- sum[[1]]
  mu_sd <- sum[[1,2]]
  sigma_int <- sum[[7]]
  sigma_sd <- sum[[7,2]]
  mu_time <- sum [[2]]
  mu_time_sd <- sum[[2,2]]
  sig_time <- sum[[8]]
  sig_time_sd <- sum[[8,2]]
  
  # Put the results 
  evs_preds_df <- rbind(evs_preds_df, data.frame(dep_var = dep_var, mu_intercept_pred = mu_int, mu_sd_pred = mu_sd, sigma_intercept_pred = sigma_int, sigma_sd_pred = sigma_sd, mu_time_pred = mu_time, mu_time_sd_pred = mu_time_sd, sig_time_pred = sig_time, sig_time_sd_pred = sig_time_sd))
}

# And lastly the interaction model
load("./data/final_data/regression_outputs/evs_list_interactions_w_new.RData")

# Store the results in a new dataframe
evs_interactions_df <- data.frame(dep_var = character())

# Loop over the number of dep vars in the list
for (i in seq_along(evs_list_interactions)) {

  # Extract the dependent variable so that we know for which dep var the effects are
  dep_var <- as.character(evs_list_interactions[[i]]$mu.terms[[2]])
  
    # Summarize the model
  sum <- summary(evs_list_interactions[[i]])
  
  # Extract the info I need from the summary
  mu_int <- sum[[1]]
  mu_sd <- sum[[1,2]]
  sigma_int <- sum[[11]]
  sigma_sd <- sum[[11,2]]
  mu_time <- sum [[2]]
  mu_time_sd <- sum[[2,2]]
  sig_time <- sum[[12]]
  sig_time_sd <- sum[[12,2]]
  
  # Put the results in a dataframe
  evs_interactions_df <- rbind(evs_interactions_df, data.frame(dep_var = dep_var, mu_intercept_int = mu_int, mu_sd_int = mu_sd, sigma_intercept_int = sigma_int, sigma_sd_int = sigma_sd, mu_time_int = mu_time, mu_time_sd_int = mu_time_sd, sig_time_int = sig_time, sig_time_sd_int = sig_time_sd))
}

# Merge the results
evs_results_df <- left_join(evs_empty_df, evs_preds_df)
evs_results_df <- left_join(evs_results_df, evs_interactions_df)
evs_results_df$first_year <- 1990
evs_results_df$mean_year <- 1999
evs_results_df$data <- "EVS"

save(evs_results_df, file="./data/final_data/regression_outputs/evs_results_df_w_new.RData")

I&O Research

# Now for I&O research
load("./data/final_data/regression_outputs/io_list_empty_w.RData")

# Store the results in a new dataframe
io_empty_df <- data.frame(dep_var = character())

for (i in seq_along(io_list_empty)) {

  # Extract the dependent variable so that we know for which dep var the effects are
  dep_var <- as.character(io_list_empty[[i]]$mu.terms[[2]])
  
    # Summarize the model
  sum <- summary(io_list_empty[[i]])
  
  # Extract the info
  mu_int <- sum[[1]]
  mu_sd <- sum[[1,2]]
  sigma_int <- sum[[3]]
  sigma_sd <- sum[[3,2]]
  mu_time <- sum [[2]]
  mu_time_sd <- sum[[2,2]]
  sig_time <- sum[[4]]
  sig_time_sd <- sum[[4,2]]
  
  # Put the results in a dataframe
  io_empty_df <- rbind(io_empty_df, data.frame(dep_var = dep_var, mu_intercept = mu_int, mu_sd = mu_sd, sigma_intercept = sigma_int, sigma_sd = sigma_sd, mu_time = mu_time, mu_time_sd = mu_time_sd, sig_time = sig_time, sig_time_sd = sig_time_sd))
}

# Predictor model
load("./data/final_data/regression_outputs/io_list_preds_w_new.RData")

# Store the results in a new dataframe
io_preds_df <- data.frame(dep_var = character())

# Loop over the list 
for (i in seq_along(io_list_preds)) {

  # Extract the dependent variable so that we know for which dep var the effects are
  dep_var <- as.character(io_list_preds[[i]]$mu.terms[[2]])
  
    # Summarize the model
  sum <- summary(io_list_preds[[i]])
  
  # Extract info
  mu_int <- sum[[1]]
  mu_sd <- sum[[1,2]]
  sigma_int <- sum[[7]]
  sigma_sd <- sum[[7,2]]
  mu_time <- sum [[2]]
  mu_time_sd <- sum[[2,2]]
  sig_time <- sum[[8]]
  sig_time_sd <- sum[[8,2]]
  
  # Put the results in a dataframe
  io_preds_df <- rbind(io_preds_df, data.frame(dep_var = dep_var, mu_intercept_pred = mu_int, mu_sd_pred = mu_sd, sigma_intercept_pred = sigma_int, sigma_sd_pred = sigma_sd, mu_time_pred = mu_time, mu_time_sd_pred = mu_time_sd, sig_time_pred = sig_time, sig_time_sd_pred = sig_time_sd))
}

# Interaction model
load("./data/final_data/regression_outputs/io_list_interactions_w_new.RData")

# Store the results in a new dataframe
io_interactions_df <- data.frame(dep_var = character())

# Loop over the number of dep vars in the list
for (i in seq_along(io_list_interactions)) {

  # Extract the dependent variable so that we know for which dep var the effects are
  dep_var <- as.character(io_list_interactions[[i]]$mu.terms[[2]])
  
    # Summarize the model
  sum <- summary(io_list_interactions[[i]])
  
  # Extract the info I need from the summary
  mu_int <- sum[[1]]
  mu_sd <- sum[[1,2]]
  sigma_int <- sum[[11]]
  sigma_sd <- sum[[11,2]]
  mu_time <- sum [[2]]
  mu_time_sd <- sum[[2,2]]
  sig_time <- sum[[12]]
  sig_time_sd <- sum[[12,2]]
  
  # Put the results in a dataframe
  io_interactions_df <- rbind(io_interactions_df, data.frame(dep_var = dep_var, mu_intercept_int = mu_int, mu_sd_int = mu_sd, sigma_intercept_int = sigma_int, sigma_sd_int = sigma_sd, mu_time_int = mu_time, mu_time_sd_int = mu_time_sd, sig_time_int = sig_time, sig_time_sd_int = sig_time_sd))
}

# Merge the results
io_results_df <- left_join(io_empty_df, io_preds_df)
io_results_df <- left_join(io_results_df, io_interactions_df)
io_results_df$first_year <- 2019
io_results_df$mean_year <- 2021
io_results_df$data <- "IO"

save(io_results_df, file="./data/final_data/regression_outputs/io_results_df_w_new.RData")

ISSP

# Now for ISSP, which has to happen twice (bc 2 substs based on dependent vars)
load("./data/final_data/regression_outputs/issp_list_empty_w.RData")

# Store the results in a new dataframe
issp_empty_df <- data.frame(dep_var = character())

# Loop over the list for all the dep vars
for (i in seq_along(issp_list_empty)) {

  # Extract the dependent variable so that we know for which dep var the effects are
  dep_var <- as.character(issp_list_empty[[i]]$mu.terms[[2]])
  
    # Summarize the model
  sum <- summary(issp_list_empty[[i]])
  
  # Extract the info
  mu_int <- sum[[1]]
  mu_sd <- sum[[1,2]]
  sigma_int <- sum[[3]]
  sigma_sd <- sum[[3,2]]
  mu_time <- sum [[2]]
  mu_time_sd <- sum[[2,2]]
  sig_time <- sum[[4]]
  sig_time_sd <- sum[[4,2]]
  
  # Put the results in a dataframe
  issp_empty_df <- rbind(issp_empty_df, data.frame(dep_var = dep_var, mu_intercept = mu_int, mu_sd = mu_sd, sigma_intercept = sigma_int, sigma_sd = sigma_sd, mu_time = mu_time, mu_time_sd = mu_time_sd, sig_time = sig_time, sig_time_sd = sig_time_sd))
}

# Predictor model 
load("./data/final_data/regression_outputs/issp_list_preds_w_new.RData")

# Store the results in a new dataframe
issp_preds_df <- data.frame(dep_var = character())

# Loop over the list 
for (i in seq_along(issp_list_preds)) {

  # Extract the dependent variable so that we know for which dep var the effects are
  dep_var <- as.character(issp_list_preds[[i]]$mu.terms[[2]])
  
    # Summarize the model
  sum <- summary(issp_list_preds[[i]])
  
  # Extract info
  mu_int <- sum[[1]]
  mu_sd <- sum[[1,2]]
  sigma_int <- sum[[7]]
  sigma_sd <- sum[[7,2]]
  mu_time <- sum [[2]]
  mu_time_sd <- sum[[2,2]]
  sig_time <- sum[[8]]
  sig_time_sd <- sum[[8,2]]
  
  # Put the results in a dataframe
  issp_preds_df <- rbind(issp_preds_df, data.frame(dep_var = dep_var, mu_intercept_pred = mu_int, mu_sd_pred = mu_sd, sigma_intercept_pred = sigma_int, sigma_sd_pred = sigma_sd, mu_time_pred = mu_time, mu_time_sd_pred = mu_time_sd, sig_time_pred = sig_time, sig_time_sd_pred = sig_time_sd))
}

# Interaction model 
load("./data/final_data/regression_outputs/issp_list_interactions_w_new.RData")

# Store the results in a new dataframe
issp_interactions_df <- data.frame(dep_var = character())

# Loop over the number of dep vars in the list
for (i in seq_along(issp_list_interactions)) {

  # Extract the dependent variable so that we know for which dep var the effects are
  dep_var <- as.character(issp_list_interactions[[i]]$mu.terms[[2]])
  
    # Summarize the model
  sum <- summary(issp_list_interactions[[i]])
  
  # Extract the info I need from the summary
  mu_int <- sum[[1]]
  mu_sd <- sum[[1,2]]
  sigma_int <- sum[[11]]
  sigma_sd <- sum[[11,2]]
  mu_time <- sum [[2]]
  mu_time_sd <- sum[[2,2]]
  sig_time <- sum[[12]]
  sig_time_sd <- sum[[12,2]]
  
  # Put the results in a dataframe
  issp_interactions_df <- rbind(issp_interactions_df, data.frame(dep_var = dep_var, mu_intercept_int = mu_int, mu_sd_int = mu_sd, sigma_intercept_int = sigma_int, sigma_sd_int = sigma_sd, mu_time_int = mu_time, mu_time_sd_int = mu_time_sd, sig_time_int = sig_time, sig_time_sd_int = sig_time_sd))
}

# Merge the results
issp_results_df <- left_join(issp_empty_df, issp_preds_df)
issp_results_df <- left_join(issp_results_df, issp_interactions_df)
issp_results_df$first_year <- 1993
issp_results_df$mean_year <- 2002
issp_results_df$data <- "ISSP"

save(issp_results_df, file="./data/final_data/regression_outputs/issp_results_df_w_new.RData")

# And with the second model
load("./data/final_data/regression_outputs/issp_list_empty_2_waves_w.RData")

# Store the results in a new dataframe
issp_empty_2_df <- data.frame(dep_var = character())

# Loop over the list for all the dep vars
for (i in seq_along(issp_list_empty_2_waves)) {

  # Extract the dependent variable so that we know for which dep var the effects are
  dep_var <- as.character(issp_list_empty_2_waves[[i]]$mu.terms[[2]])
  
    # Summarize the model
  sum <- summary(issp_list_empty_2_waves[[i]])
  
  # Extract the info
  mu_int <- sum[[1]]
  mu_sd <- sum[[1,2]]
  sigma_int <- sum[[3]]
  sigma_sd <- sum[[3,2]]
  mu_time <- sum [[2]]
  mu_time_sd <- sum[[2,2]]
  sig_time <- sum[[4]]
  sig_time_sd <- sum[[4,2]]
  
  # Put the results in a dataframe
  issp_empty_2_df <- rbind(issp_empty_2_df, data.frame(dep_var = dep_var, mu_intercept = mu_int, mu_sd = mu_sd, sigma_intercept = sigma_int, sigma_sd = sigma_sd, mu_time = mu_time, mu_time_sd = mu_time_sd, sig_time = sig_time, sig_time_sd = sig_time_sd))
}

# Predictor model 
load("./data/final_data/regression_outputs/issp_list_preds_2_waves_w_new.RData")

# Store the results in a new dataframe
issp_preds_2_df <- data.frame(dep_var = character())

# Loop over the list 
for (i in seq_along(issp_list_preds_2_waves)) {

  # Extract the dependent variable so that we know for which dep var the effects are
  dep_var <- as.character(issp_list_preds_2_waves[[i]]$mu.terms[[2]])
  
    # Summarize the model
  sum <- summary(issp_list_preds_2_waves[[i]])
  
  # Extract info
  mu_int <- sum[[1]]
  mu_sd <- sum[[1,2]]
  sigma_int <- sum[[7]]
  sigma_sd <- sum[[7,2]]
  mu_time <- sum [[2]]
  mu_time_sd <- sum[[2,2]]
  sig_time <- sum[[8]]
  sig_time_sd <- sum[[8,2]]
  
  # Put the results in a dataframe
  issp_preds_2_df <- rbind(issp_preds_2_df, data.frame(dep_var = dep_var, mu_intercept_pred = mu_int, mu_sd_pred = mu_sd, sigma_intercept_pred = sigma_int, sigma_sd_pred = sigma_sd, mu_time_pred = mu_time, mu_time_sd_pred = mu_time_sd, sig_time_pred = sig_time, sig_time_sd_pred = sig_time_sd))
}

# Interaction model 
load("./data/final_data/regression_outputs/issp_list_interactions_2_waves_w.RData")

# Store the results in a new dataframe
issp_interactions_2_df <- data.frame(dep_var = character())

# Loop over the number of dep vars in the list
for (i in seq_along(issp_list_interactions_2_waves)) {

  # Extract the dependent variable so that we know for which dep var the effects are
  dep_var <- as.character(issp_list_interactions_2_waves[[i]]$mu.terms[[2]])
  
    # Summarize the model
  sum <- summary(issp_list_interactions_2_waves[[i]])
  
  # Extract the info I need from the summary
  mu_int <- sum[[1]]
  mu_sd <- sum[[1,2]]
  sigma_int <- sum[[11]]
  sigma_sd <- sum[[11,2]]
  mu_time <- sum [[2]]
  mu_time_sd <- sum[[2,2]]
  sig_time <- sum[[12]]
  sig_time_sd <- sum[[12,2]]
  
  # Put the results in a dataframe
  issp_interactions_2_df <- rbind(issp_interactions_2_df, data.frame(dep_var = dep_var, mu_intercept_int = mu_int, mu_sd_int = mu_sd, sigma_intercept_int = sigma_int, sigma_sd_int = sigma_sd, mu_time_int = mu_time, mu_time_sd_int = mu_time_sd, sig_time_int = sig_time, sig_time_sd_int = sig_time_sd))
}

# Merge the results
issp_results_2_df <- left_join(issp_empty_2_df, issp_preds_2_df)
issp_results_2_df <- left_join(issp_results_2_df, issp_interactions_2_df)
issp_results_2_df$first_year <- 2000
issp_results_2_df$mean_year <- 2005
issp_results_2_df$data <- "ISSP"

save(issp_results_2_df, file="./data/final_data/regression_outputs/issp_results_2_df_new.RData")

EB

# For all the other eurobarometer waves, exactly the same needs to be done each time. Therefore I want to loop over these different lists. 

# First load all the eurobarometer lists with empty models in my environment
rm(list=ls())
load("./data/final_data/regression_outputs/issp_list_interactions_2_waves_w.RData")
load("./data/final_data/regression_outputs/eb_list_empty_1986_w.RData")
load("./data/final_data/regression_outputs/eb_list_empty_2007_w.RData")
load("./data/final_data/regression_outputs/eb_list_empty_2008_w.RData")
load("./data/final_data/regression_outputs/eb_list_empty_2009_w.RData")
load("./data/final_data/regression_outputs/eb_list_empty_2011_w.RData")
load("./data/final_data/regression_outputs/eb_list_empty_buyprod_w.RData")
load("./data/final_data/regression_outputs/eb_list_empty_cchange2_w.RData")
load("./data/final_data/regression_outputs/eb_list_empty_doprot_w.RData")
load("./data/final_data/regression_outputs/eb_list_empty_eff_daily_w.RData")
load("./data/final_data/regression_outputs/eb_list_empty_pers_imp_w.RData")
load("./data/final_data/regression_outputs/eb_list_empty_prsaction_w.RData")
load("./data/final_data/regression_outputs/eb_list_empty_quallife_w.RData")
load("./data/final_data/regression_outputs/eb_list_empty_doprot_city_w.RData")


# Store the results in a new dataframe
eb_results_df <- data.frame(dep_var = character())

# Get the lists from the environment
lists <- ls()

for (j in lists) { 
  if(is.list(get(j))) { 
    my_list <- get(j)

# Loop over the list for all the dep vars
for (i in seq_along(my_list)) {

  # Extract the dependent variable so that we know for which dep var the effects are
  dep_var <- as.character(my_list[[i]]$mu.terms[[2]])
  
    # Summarize the model
  sum <- summary(my_list[[i]])
  
  # Extract the info
  mu_int <- sum[[1]]
  mu_sd <- sum[[1,2]]
  sigma_int <- sum[[3]]
  sigma_sd <- sum[[3,2]]
  mu_time <- sum [[2]]
  mu_time_sd <- sum[[2,2]]
  sig_time <- sum[[4]]
  sig_time_sd <- sum[[4,2]]
  
  # Put the results in a dataframe
  eb_results_df <- rbind(eb_results_df, data.frame(dep_var = dep_var, mu_intercept = mu_int, mu_sd = mu_sd, sigma_intercept = sigma_int, sigma_sd = sigma_sd, mu_time = mu_time, mu_time_sd = mu_time_sd, sig_time = sig_time, sig_time_sd = sig_time_sd))
    }
  } 
} 

# It produces an error, but actually everything has worked as i wish?
eb_results_empty_df <- eb_results_df

# let's save it 
save(eb_results_empty_df, file= "./data/final_data/regression_outputs/eb_results_empty_df_w_new.RData" )

# And the numbers are indeed extracted correctly if i check this randomly
sum <- summary(eb_list_empty_quallife[[1]])

# Now i want the same for the models with the predictors
# First clear the environment because I don't want the empty lists iterated in the predictor loop
rm(list=ls())
load("./data/final_data/regression_outputs/eb_list_preds_1986_w_new.RData")
load("./data/final_data/regression_outputs/eb_list_preds_2007_w_new.RData")
load("./data/final_data/regression_outputs/eb_list_preds_2008_w_new.RData")
load("./data/final_data/regression_outputs/eb_list_preds_2009_w_new.RData")
load("./data/final_data/regression_outputs/eb_list_preds_2011_w_new.RData")
load("./data/final_data/regression_outputs/eb_list_preds_buyprod_w_new.RData")
load("./data/final_data/regression_outputs/eb_list_preds_cchange2_w_new.RData")
load("./data/final_data/regression_outputs/eb_list_preds_doprot_w_new.RData")
load("./data/final_data/regression_outputs/eb_list_preds_eff_daily_w_new.RData")
load("./data/final_data/regression_outputs/eb_list_preds_pers_imp_w_new.RData")
load("./data/final_data/regression_outputs/eb_list_preds_prsaction_w_new.RData")
load("./data/final_data/regression_outputs/eb_list_preds_doprot_city_w_new.RData")


# Store the results in a new dataframe
eb_results_preds_df <- data.frame(dep_var = character())

# Get the lists from the environment
lists <- ls()

for (j in lists) { 
  if(is.list(get(j))) { 
    my_list <- get(j)

# Loop over the list 
for (i in seq_along(my_list)) {

  # Extract the dependent variable so that we know for which dep var the effects are
  dep_var <- as.character(my_list[[i]]$mu.terms[[2]])
  
    # Summarize the model
  sum <- summary(my_list[[i]])
  
  # Extract info
  mu_int <- sum[[1]]
  mu_sd <- sum[[1,2]]
  sigma_int <- sum[[7]]
  sigma_sd <- sum[[7,2]]
  mu_time <- sum [[2]]
  mu_time_sd <- sum[[2,2]]
  sig_time <- sum[[8]]
  sig_time_sd <- sum[[8,2]]
  
  # Put the results in a dataframe
  eb_results_preds_df <- rbind(eb_results_preds_df, data.frame(dep_var = dep_var, mu_intercept_pred = mu_int, mu_sd_pred = mu_sd, sigma_intercept_pred = sigma_int, sigma_sd_pred = sigma_sd, mu_time_pred = mu_time, mu_time_sd_pred = mu_time_sd, sig_time_pred = sig_time, sig_time_sd_pred = sig_time_sd))
    }
  }
} 

# Again, an error, but the output is correct
save(eb_results_preds_df, file= "./data/final_data/regression_outputs/eb_results_preds_df_w_new.RData" )

# Lastly for the interaction models
rm(list=ls())
load("./data/final_data/regression_outputs/eb_list_interactions_1986_w_new.RData")
load("./data/final_data/regression_outputs/eb_list_interactions_2007_w_new.RData")
load("./data/final_data/regression_outputs/eb_list_interactions_2008_w_new.RData")
load("./data/final_data/regression_outputs/eb_list_interactions_2009_w_new.RData")
load("./data/final_data/regression_outputs/eb_list_interactions_2011_w_new.RData")
load("./data/final_data/regression_outputs/eb_list_interactions_buyprod_w_new.RData")
load("./data/final_data/regression_outputs/eb_list_interactions_cchange2_w_new.RData")
load("./data/final_data/regression_outputs/eb_list_interactions_doprot_w_new.RData")
load("./data/final_data/regression_outputs/eb_list_interactions_eff_daily_w_new.RData")
load("./data/final_data/regression_outputs/eb_list_interactions_pers_imp_w_new.RData")
load("./data/final_data/regression_outputs/eb_list_interactions_prsaction_w_new.RData")
load("./data/final_data/regression_outputs/eb_list_interactions_quallife_w_new.RData")
load("./data/final_data/regression_outputs/eb_list_interactions_doprot_city_w_new.RData")


# Store the results in a new dataframe
eb_results_interactions_df <- data.frame(dep_var = character())

# Get the lists from the environment
lists <- ls()

for (j in lists) { 
  if(is.list(get(j))) { 
    my_list <- get(j)

# Loop over the number of dep vars in the list
for (i in seq_along(my_list)) {

  # Extract the dependent variable so that we know for which dep var the effects are
  dep_var <- as.character(my_list[[i]]$mu.terms[[2]])
  
    # Summarize the model
  sum <- summary(my_list[[i]])
  
  # Extract the info I need from the summary
  mu_int <- sum[[1]]
  mu_sd <- sum[[1,2]]
  sigma_int <- sum[[11]]
  sigma_sd <- sum[[11,2]]
  mu_time <- sum [[2]]
  mu_time_sd <- sum[[2,2]]
  sig_time <- sum[[12]]
  sig_time_sd <- sum[[12,2]]
  
  # Put the results in a dataframe
  eb_results_interactions_df <- rbind(eb_results_interactions_df, data.frame(dep_var = dep_var, mu_intercept_int = mu_int, mu_sd_int = mu_sd, sigma_intercept_int = sigma_int, sigma_sd_int = sigma_sd, mu_time_int = mu_time, mu_time_sd_int = mu_time_sd, sig_time_int = sig_time, sig_time_sd_int = sig_time_sd))
    }
  } 
} 


save(eb_results_interactions_df, file= "./data/final_data/regression_outputs/eb_results_interactions_df_w_new.RData" )

# Make one dataframe of the three loops
load("./data/final_data/regression_outputs/eb_results_empty_df_w_new.RData")
load("./data/final_data/regression_outputs/eb_results_preds_df_w_new.RData")

# Merge the results
eb_results_df <- left_join(eb_results_empty_df, eb_results_preds_df)
eb_results_df <- left_join(eb_results_df, eb_results_interactions_df)
eb_results_df$data <- "EB"

# This time, year is a bit more difficult, so I have to assign that per dep_var
eb_results_df$first_year[eb_results_df$dep_var == "role_ind"] <- 2007
eb_results_df$first_year[eb_results_df$dep_var == "big_pol"] <- 2007
eb_results_df$first_year[eb_results_df$dep_var == "cc_unstop"] <- 2008
eb_results_df$first_year[eb_results_df$dep_var == "cc_exag"] <- 2008
eb_results_df$first_year[eb_results_df$dep_var == "cc_poseu"] <- 2008
eb_results_df$first_year[eb_results_df$dep_var == "cc_prsact"] <- 2008
eb_results_df$first_year[eb_results_df$dep_var == "cc_poseu"] <- 2008
eb_results_df$first_year[eb_results_df$dep_var == "ccpercept"] <- 2009
eb_results_df$first_year[eb_results_df$dep_var == "cchange"] <- 2009
eb_results_df$first_year[eb_results_df$dep_var == "cchange2"] <- 2011
eb_results_df$first_year[eb_results_df$dep_var == "cchangetot"] <- 2011
eb_results_df$first_year[eb_results_df$dep_var == "envp_eg"] <- 2011
eb_results_df$first_year[eb_results_df$dep_var == "effr_eg"] <- 2011
eb_results_df$first_year[eb_results_df$dep_var == "buyprod"] <- 2007
eb_results_df$first_year[eb_results_df$dep_var == "doprot_natgov"] <- 2009
eb_results_df$first_year[eb_results_df$dep_var == "doprot_eu"] <- 2009
eb_results_df$first_year[eb_results_df$dep_var == "doprot_region"] <- 2009
eb_results_df$first_year[eb_results_df$dep_var == "doprot_comp"] <- 2009
eb_results_df$first_year[eb_results_df$dep_var == "doprot_citiz"] <- 2009
eb_results_df$first_year[eb_results_df$dep_var == "doprot_city"] <- 2014
eb_results_df$first_year[eb_results_df$dep_var == "eff_daily"] <- 2007
eb_results_df$first_year[eb_results_df$dep_var == "pers_imp"] <- 2009
eb_results_df$first_year[eb_results_df$dep_var == "env_quallife"] <- 2004
eb_results_df$first_year[eb_results_df$dep_var == "env_ec_stat"] <- 1986
eb_results_df$first_year[eb_results_df$dep_var == "env_prsimp"] <- 1986
eb_results_df$first_year[eb_results_df$dep_var == "prsaction"] <- 2011

# And the average year 
eb_results_df$mean_year[eb_results_df$dep_var == "role_ind"] <- 2012
eb_results_df$mean_year[eb_results_df$dep_var == "big_pol"] <- 2012
eb_results_df$mean_year[eb_results_df$dep_var == "cc_unstop"] <- 2008
eb_results_df$mean_year[eb_results_df$dep_var == "cc_exag"] <- 2008
eb_results_df$mean_year[eb_results_df$dep_var == "cc_poseu"] <- 2008
eb_results_df$mean_year[eb_results_df$dep_var == "cc_prsact"] <- 2008
eb_results_df$mean_year[eb_results_df$dep_var == "cc_poseu"] <- 2008
eb_results_df$mean_year[eb_results_df$dep_var == "ccpercept"] <- 2014
eb_results_df$mean_year[eb_results_df$dep_var == "cchange"] <- 2015
eb_results_df$mean_year[eb_results_df$dep_var == "cchange2"] <- 2015
eb_results_df$mean_year[eb_results_df$dep_var == "cchangetot"] <- 2015
eb_results_df$mean_year[eb_results_df$dep_var == "envp_eg"] <- 2013
eb_results_df$mean_year[eb_results_df$dep_var == "effr_eg"] <- 2013
eb_results_df$mean_year[eb_results_df$dep_var == "buyprod"] <- 2010
eb_results_df$mean_year[eb_results_df$dep_var == "doprot_natgov"] <- 2013
eb_results_df$mean_year[eb_results_df$dep_var == "doprot_eu"] <- 2013
eb_results_df$mean_year[eb_results_df$dep_var == "doprot_region"] <- 2013
eb_results_df$mean_year[eb_results_df$dep_var == "doprot_comp"] <- 2013
eb_results_df$mean_year[eb_results_df$dep_var == "doprot_citiz"] <- 2013
eb_results_df$mean_year[eb_results_df$dep_var == "doprot_city"] <- 2015
eb_results_df$mean_year[eb_results_df$dep_var == "eff_daily"] <- 2012
eb_results_df$mean_year[eb_results_df$dep_var == "pers_imp"] <- 2009
eb_results_df$mean_year[eb_results_df$dep_var == "env_quallife"] <- 2009
eb_results_df$mean_year[eb_results_df$dep_var == "env_ec_stat"] <- 1991
eb_results_df$mean_year[eb_results_df$dep_var == "env_prsimp"] <- 1991
eb_results_df$mean_year[eb_results_df$dep_var == "prsaction"] <- 2014

save(eb_results_df, file= "./data/final_data/regression_outputs/eb_results_df_w_new.RData" )

Motivaction

rm(list=ls())
load("./data/final_data/regression_outputs/mot_list_empty_w.RData")

mot_empty_df <- data.frame(dep_var = character())


# Loop over the list of models
for (i in seq_along(mot_list_empty)) {

  # Extract the dependent variable so that we know for which dep var the effects are
  dep_var <- as.character(mot_list_empty[[i]]$mu.terms[[2]])
  
    # Summarize the model
  sum <- summary(mot_list_empty[[i]])
  
  # Extract the intercept and standard deviation for mu and sigma
  mu_int <- sum[[1]]
  mu_sd <- sum[[1,2]]
  sigma_int <- sum[[3]]
  sigma_sd <- sum[[3,2]]
  mu_time <- sum [[2]]
  mu_time_sd <- sum[[2,2]]
  sig_time <- sum[[4]]
  sig_time_sd <- sum[[4,2]]
  
  # Put the results 
  mot_empty_df <- rbind(mot_empty_df, data.frame(dep_var = dep_var, mu_intercept = mu_int, mu_sd = mu_sd, sigma_intercept = sigma_int, sigma_sd = sigma_sd, mu_time = mu_time, mu_time_sd = mu_time_sd, sig_time = sig_time, sig_time_sd = sig_time_sd))
}

load("./data/final_data/regression_outputs/mot_list_preds_w_new.RData")

# This does mean that the estimates of the sigma are in a different place
sum <- summary(mot_list_preds[[1]])

# Store the results in a new dataframe
mot_preds_df <- data.frame(dep_var = character())

# Loop over the list of models
for (i in seq_along(mot_list_preds)) {

  # Extract the dependent variable so that we know for which dep var the effects are
  dep_var <- as.character(mot_list_preds[[i]]$mu.terms[[2]])
  
    # Summarize the model
  sum <- summary(mot_list_preds[[i]])
  
  # Extract the intercept and standard deviation for mu and sigma
  mu_int <- sum[[1]]
  mu_sd <- sum[[1,2]]
  sigma_int <- sum[[7]]
  sigma_sd <- sum[[7,2]]
  mu_time <- sum [[2]]
  mu_time_sd <- sum[[2,2]]
  sig_time <- sum[[8]]
  sig_time_sd <- sum[[8,2]]
  
  # Put the results 
  mot_preds_df <- rbind(mot_preds_df, data.frame(dep_var = dep_var, mu_intercept_pred = mu_int, mu_sd_pred = mu_sd, sigma_intercept_pred = sigma_int, sigma_sd_pred = sigma_sd, mu_time_pred = mu_time, mu_time_sd_pred = mu_time_sd, sig_time_pred = sig_time, sig_time_sd_pred = sig_time_sd))
}


# And then for the interaction model
load("./data/final_data/regression_outputs/mot_list_interactions_w_new.RData")

# Store the results in a new dataframe
mot_interactions_df <- data.frame(dep_var = character())

# Loop over the number of dep vars in the list
for (i in seq_along(mot_list_interactions)) {

  # Extract the dependent variable so that we know for which dep var the effects are
  dep_var <- as.character(mot_list_interactions[[i]]$mu.terms[[2]])
  
    # Summarize the model
  sum <- summary(mot_list_interactions[[i]])
  
  # Extract the info I need from the summary
  mu_int <- sum[[1]]
  mu_sd <- sum[[1,2]]
  sigma_int <- sum[[11]]
  sigma_sd <- sum[[11,2]]
  mu_time <- sum [[2]]
  mu_time_sd <- sum[[2,2]]
  sig_time <- sum[[12]]
  sig_time_sd <- sum[[12,2]]
  
  # Put the results in a dataframe
  mot_interactions_df <- rbind(mot_interactions_df, data.frame(dep_var = dep_var, mu_intercept_int = mu_int, mu_sd_int = mu_sd, sigma_intercept_int = sigma_int, sigma_sd_int = sigma_sd, mu_time_int = mu_time, mu_time_sd_int = mu_time_sd, sig_time_int = sig_time, sig_time_sd_int = sig_time_sd))
}


#For some reason left_join produces an error, but merge gives the same result
#mot_results_df <- left_join(mot_empty_df, mot_preds_df)
#mot_results_df <- left_join(mot_results_df, mot_interactions_df)

mot_results_df <- merge(mot_empty_df, mot_preds_df)
mot_results_df <- merge(mot_results_df, mot_interactions_df)
mot_results_df$first_year <- 2019
mot_results_df$mean_year <- 2020
mot_results_df$data <- "MOT"

save(mot_results_df, file="./data/final_data/regression_outputs/mot_results_df_w_new.RData")

LISS

rm(list = ls())
load("./data/final_data/regression_outputs/liss_list_empty_w.RData")

liss_empty_df <- data.frame(dep_var = character())


# Loop over the list of models
for (i in seq_along(liss_list_empty)) {

  # Extract the dependent variable so that we know for which dep var the effects are
  dep_var <- as.character(liss_list_empty[[i]]$mu.terms[[2]])
  
    # Summarize the model
  sum <- summary(liss_list_empty[[i]])
  
  # Extract the intercept and standard deviation for mu and sigma
  mu_int <- sum[[1]]
  mu_sd <- sum[[1,2]]
  sigma_int <- sum[[3]]
  sigma_sd <- sum[[3,2]]
  mu_time <- sum [[2]]
  mu_time_sd <- sum[[2,2]]
  sig_time <- sum[[4]]
  sig_time_sd <- sum[[4,2]]
  
  # Put the results 
  liss_empty_df <- rbind(liss_empty_df, data.frame(dep_var = dep_var, mu_intercept = mu_int, mu_sd = mu_sd, sigma_intercept = sigma_int, sigma_sd = sigma_sd, mu_time = mu_time, mu_time_sd = mu_time_sd, sig_time = sig_time, sig_time_sd = sig_time_sd))
}

load("./data/final_data/regression_outputs/liss_list_preds_w_new.RData")

# This does mean that the estimates of the sigma are in a different place
sum <- summary(liss_list_preds[[1]])

# Store the results in a new dataframe
liss_preds_df <- data.frame(dep_var = character())

# Loop over the list of models
for (i in seq_along(liss_list_preds)) {

  # Extract the dependent variable so that we know for which dep var the effects are
  dep_var <- as.character(liss_list_preds[[i]]$mu.terms[[2]])
  
    # Summarize the model
  sum <- summary(liss_list_preds[[i]])
  
  # Extract the intercept and standard deviation for mu and sigma
  mu_int <- sum[[1]]
  mu_sd <- sum[[1,2]]
  sigma_int <- sum[[7]]
  sigma_sd <- sum[[7,2]]
  mu_time <- sum [[2]]
  mu_time_sd <- sum[[2,2]]
  sig_time <- sum[[8]]
  sig_time_sd <- sum[[8,2]]
  
  # Put the results 
  liss_preds_df <- rbind(liss_preds_df, data.frame(dep_var = dep_var, mu_intercept_pred = mu_int, mu_sd_pred = mu_sd, sigma_intercept_pred = sigma_int, sigma_sd_pred = sigma_sd, mu_time_pred = mu_time, mu_time_sd_pred = mu_time_sd, sig_time_pred = sig_time, sig_time_sd_pred = sig_time_sd))
}


# And then for the interaction model
load("./data/final_data/regression_outputs/liss_list_interactions_w_new.RData")

# Store the results in a new dataframe
liss_interactions_df <- data.frame(dep_var = character())

# Loop over the number of dep vars in the list
for (i in seq_along(liss_list_interactions)) {

  # Extract the dependent variable so that we know for which dep var the effects are
  dep_var <- as.character(liss_list_interactions[[i]]$mu.terms[[2]])
  
    # Summarize the model
  sum <- summary(liss_list_interactions[[i]])
  
  # Extract the info I need from the summary
  mu_int <- sum[[1]]
  mu_sd <- sum[[1,2]]
  sigma_int <- sum[[11]]
  sigma_sd <- sum[[11,2]]
  mu_time <- sum [[2]]
  mu_time_sd <- sum[[2,2]]
  sig_time <- sum[[12]]
  sig_time_sd <- sum[[12,2]]
  
  # Put the results in a dataframe
  liss_interactions_df <- rbind(liss_interactions_df, data.frame(dep_var = dep_var, mu_intercept_int = mu_int, mu_sd_int = mu_sd, sigma_intercept_int = sigma_int, sigma_sd_int = sigma_sd, mu_time_int = mu_time, mu_time_sd_int = mu_time_sd, sig_time_int = sig_time, sig_time_sd_int = sig_time_sd))
}


#For some reason left_join produces an error, but merge gives the same result
#liss_results_df <- left_join(liss_empty_df, liss_preds_df)
#liss_results_df <- left_join(liss_results_df, liss_interactions_df)

liss_results_df <- merge(liss_empty_df, liss_preds_df)
liss_results_df <- merge(liss_results_df, liss_interactions_df)

liss_results_df$first_year <- 2019
liss_results_df$mean_year <- 2020
liss_results_df$data <- "LISS"

save(liss_results_df, file="./data/final_data/regression_outputs/liss_results_df_w_new.RData")

SOCON

rm(list=ls())
load("./data/final_data/regression_outputs/socon_list_empty_w.RData")

socon_empty_df <- data.frame(dep_var = character())


# Loop over the list of models
for (i in seq_along(socon_list_empty)) {

  # Extract the dependent variable so that we know for which dep var the effects are
  dep_var <- as.character(socon_list_empty[[i]]$mu.terms[[2]])
  
    # Summarize the model
  sum <- summary(socon_list_empty[[i]])
  
  # Extract the intercept and standard deviation for mu and sigma
  mu_int <- sum[[1]]
  mu_sd <- sum[[1,2]]
  sigma_int <- sum[[3]]
  sigma_sd <- sum[[3,2]]
  mu_time <- sum [[2]]
  mu_time_sd <- sum[[2,2]]
  sig_time <- sum[[4]]
  sig_time_sd <- sum[[4,2]]
  
  # Put the results 
  socon_empty_df <- rbind(socon_empty_df, data.frame(dep_var = dep_var, mu_intercept = mu_int, mu_sd = mu_sd, sigma_intercept = sigma_int, sigma_sd = sigma_sd, mu_time = mu_time, mu_time_sd = mu_time_sd, sig_time = sig_time, sig_time_sd = sig_time_sd))
}

load("./data/final_data/regression_outputs/socon_list_preds_w_new.RData")

# This does mean that the estimates of the sigma are in a different place
sum <- summary(socon_list_preds[[1]])

# Store the results in a new dataframe
socon_preds_df <- data.frame(dep_var = character())

# Loop over the list of models
for (i in seq_along(socon_list_preds)) {

  # Extract the dependent variable so that we know for which dep var the effects are
  dep_var <- as.character(socon_list_preds[[i]]$mu.terms[[2]])
  
    # Summarize the model
  sum <- summary(socon_list_preds[[i]])
  
  # Extract the intercept and standard deviation for mu and sigma
  mu_int <- sum[[1]]
  mu_sd <- sum[[1,2]]
  sigma_int <- sum[[7]]
  sigma_sd <- sum[[7,2]]
  mu_time <- sum [[2]]
  mu_time_sd <- sum[[2,2]]
  sig_time <- sum[[8]]
  sig_time_sd <- sum[[8,2]]
  
  # Put the results 
  socon_preds_df <- rbind(socon_preds_df, data.frame(dep_var = dep_var, mu_intercept_pred = mu_int, mu_sd_pred = mu_sd, sigma_intercept_pred = sigma_int, sigma_sd_pred = sigma_sd, mu_time_pred = mu_time, mu_time_sd_pred = mu_time_sd, sig_time_pred = sig_time, sig_time_sd_pred = sig_time_sd))
}


# And then for the interaction model
load("./data/final_data/regression_outputs/socon_list_interactions_w_new.RData")

# Store the results in a new dataframe
socon_interactions_df <- data.frame(dep_var = character())

# Loop over the number of dep vars in the list
for (i in seq_along(socon_list_interactions)) {

  # Extract the dependent variable so that we know for which dep var the effects are
  dep_var <- as.character(socon_list_interactions[[i]]$mu.terms[[2]])
  
    # Summarize the model
  sum <- summary(socon_list_interactions[[i]])
  
  # Extract the info I need from the summary
  mu_int <- sum[[1]]
  mu_sd <- sum[[1,2]]
  sigma_int <- sum[[11]]
  sigma_sd <- sum[[11,2]]
  mu_time <- sum [[2]]
  mu_time_sd <- sum[[2,2]]
  sig_time <- sum[[12]]
  sig_time_sd <- sum[[12,2]]
  
  # Put the results in a dataframe
  socon_interactions_df <- rbind(socon_interactions_df, data.frame(dep_var = dep_var, mu_intercept_int = mu_int, mu_sd_int = mu_sd, sigma_intercept_int = sigma_int, sigma_sd_int = sigma_sd, mu_time_int = mu_time, mu_time_sd_int = mu_time_sd, sig_time_int = sig_time, sig_time_sd_int = sig_time_sd))
}


#For some reason left_join produces an error, but merge gives the same result
#socon_results_df <- left_join(socon_empty_df, socon_preds_df)
#socon_results_df <- left_join(socon_results_df, socon_interactions_df)

socon_results_df <- merge(socon_empty_df, socon_preds_df)
socon_results_df <- merge(socon_results_df, socon_interactions_df)
utils::View(socon_results_df)

socon_results_df$first_year <- 2020
socon_results_df$mean_year <- 2021
socon_results_df$data <- "SOCON"

save(socon_results_df, file="./data/final_data/regression_outputs/socon_results_df_w_new.RData")

Merge everything in one dataset

# The final step is to merge all results into one dataframe
rm(list=ls())
load("./data/final_data/regression_outputs/ess_results_df_w_new.RData")
load("./data/final_data/regression_outputs/evs_results_df_w_new.RData")
load("./data/final_data/regression_outputs/io_results_df_w_new.RData")
load("./data/final_data/regression_outputs/issp_results_df_w_new.RData")
load("./data/final_data/regression_outputs/issp_results_2_df_new.RData")
load("./data/final_data/regression_outputs/eb_results_df_w_new.RData")
load("./data/final_data/regression_outputs/mot_results_df_w_new.RData")
load("./data/final_data/regression_outputs/liss_results_df_w_new.RData")
load("./data/final_data/regression_outputs/socon_results_df_w_new.RData")

total_reg_results <- rbind(eb_results_df, ess_results_df, evs_results_df, io_results_df, issp_results_df, issp_results_2_df, mot_results_df, liss_results_df, socon_results_df)

save(total_reg_results, file= "./data/meta_analysis/total_reg_results_gam_w_new.RData" )

rm(list=ls())
load("./data/meta_analysis/total_reg_results_gam_w_new.RData")

# I also want to add a variable that describes whether the dependent variable is mostly of affective, behavioral or cognitive nature 

# 2 variables that are called worry, not handy 
total_reg_results$dep_var[total_reg_results$dep_var == "worry" & total_reg_results$data == "ISSP"] <- "worry_issp"

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

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

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

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

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

table(total_reg_results$ambiguous, useNA = "always")

# I also want an indicator for whether it is a national or international dataset. 
total_reg_results$national[total_reg_results$data == "IO"|
                             total_reg_results$data == "MOT" |
                             total_reg_results$data == "LISS" |
                             total_reg_results$data == "SOCON"] <- 1
total_reg_results$national[total_reg_results$data == "EB"|
                             total_reg_results$data == "ESS" |
                             total_reg_results$data == "EVS" |
                             total_reg_results$data == "ISSP"] <- 0

# I would like to change the names of the variables so that they are a bit more descriptive and look better
total_reg_results$dep_var_name[total_reg_results$dep_var == "dodiff"] <- "Do things differently" 
total_reg_results$dep_var_name[total_reg_results$dep_var == "frontrunner"] <- "NL frontrunner"
total_reg_results$dep_var_name[total_reg_results$dep_var == "min_contr"] <- "Minimal contribution"
total_reg_results$dep_var_name[total_reg_results$dep_var == "people_decide"] <- "Let people decide"
total_reg_results$dep_var_name[total_reg_results$dep_var == "econprotect"] <- "Need economic growth"
total_reg_results$dep_var_name[total_reg_results$dep_var == "growharm"] <- "Growth harms env"
total_reg_results$dep_var_name[total_reg_results$dep_var == "bus_decide"] <- "Let business decide"
total_reg_results$dep_var_name[total_reg_results$dep_var == "othersame"] <- "Others do same"
total_reg_results$dep_var_name[total_reg_results$dep_var == "resp_citiz"] <- "Responsibility citizen"
total_reg_results$dep_var_name[total_reg_results$dep_var == "dk_start"] <- "Don't know where start"
total_reg_results$dep_var_name[total_reg_results$dep_var == "buss_help"] <- "Large companies reduce"
total_reg_results$dep_var_name[total_reg_results$dep_var == "envp_eg"] <- "Protection boosts growth"
total_reg_results$dep_var_name[total_reg_results$dep_var == "effr_eg"] <- "Use natural resources"
total_reg_results$dep_var_name[total_reg_results$dep_var == "env_quallife"] <- "Env affects quality life"
total_reg_results$dep_var_name[total_reg_results$dep_var == "cc_unstop"] <- "Climate change unstoppable"
total_reg_results$dep_var_name[total_reg_results$dep_var == "cc_poseu"] <- "Fighting cc pos impact EU"
total_reg_results$dep_var_name[total_reg_results$dep_var == "big_pol"] <- "Big polluters responsible"
total_reg_results$dep_var_name[total_reg_results$dep_var == "pers_resp"] <- "Personal responsiblity" 
total_reg_results$dep_var_name[total_reg_results$dep_var == "worry"] <- "Worried climate change"
total_reg_results$dep_var_name[total_reg_results$dep_var == "worried"] <- "Worry CO2 and climate change"
total_reg_results$dep_var_name[total_reg_results$dep_var == "worry_future"] <- "Bleak future generations"
total_reg_results$dep_var_name[total_reg_results$dep_var == "willing_price"] <- "Willing pay higher prices"
total_reg_results$dep_var_name[total_reg_results$dep_var == "willing_tax"] <- "Willing pay higher taxes"
total_reg_results$dep_var_name[total_reg_results$dep_var == "willing_living"] <- "Willing cut std of living"
total_reg_results$dep_var_name[total_reg_results$dep_var == "do_right"] <- "Do what is right"
total_reg_results$dep_var_name[total_reg_results$dep_var == "climate5"] <- "Would give part income" 
total_reg_results$dep_var_name[total_reg_results$dep_var == "prsaction"] <- "Taken personal action"
total_reg_results$dep_var_name[total_reg_results$dep_var == "cc_prsact"] <- "Taken personal action"
total_reg_results$dep_var_name[total_reg_results$dep_var == "buyprod"] <- "Buy environmental products"
total_reg_results$dep_var_name[total_reg_results$dep_var == "worry_issp"] <- "Worry too much about future"
total_reg_results$dep_var_name[total_reg_results$dep_var == "lifeharm"] <- "Everything we do harms env"
total_reg_results$dep_var_name[total_reg_results$dep_var == "progharm"] <- "Worry too much about harm"
total_reg_results$dep_var_name[total_reg_results$dep_var == "moreimp"] <- "Other things more important"
total_reg_results$dep_var_name[total_reg_results$dep_var == "exag"] <- "Claims exaggerated" 
total_reg_results$dep_var_name[total_reg_results$dep_var == "country_effort"] <- "Countries doing enough"
total_reg_results$dep_var_name[total_reg_results$dep_var == "cause"] <- "Cause climate change"
total_reg_results$dep_var_name[total_reg_results$dep_var == "human_resp"] <- "Humans responsible"
total_reg_results$dep_var_name[total_reg_results$dep_var == "env_ec_stat"] <- "Env vs economy"
total_reg_results$dep_var_name[total_reg_results$dep_var == "env_prsimp"] <- "Urgent problem"
total_reg_results$dep_var_name[total_reg_results$dep_var == "cchange"] <- "Most important problem"
total_reg_results$dep_var_name[total_reg_results$dep_var == "cchange2"] <- "2nd most important problem"
total_reg_results$dep_var_name[total_reg_results$dep_var == "cchangetot"] <- "Mentioned as important problem"
total_reg_results$dep_var_name[total_reg_results$dep_var == "ccpercept"] <- "Serious climate change"
total_reg_results$dep_var_name[total_reg_results$dep_var == "doprot_comp"] <- "Companies do enough" 
total_reg_results$dep_var_name[total_reg_results$dep_var == "doprot_citiz"] <- "Citizens do enough"
total_reg_results$dep_var_name[total_reg_results$dep_var == "doprot_region"] <- "Region does enough"
total_reg_results$dep_var_name[total_reg_results$dep_var == "do_gov"] <- "Government does enough"
total_reg_results$dep_var_name[total_reg_results$dep_var == "doprot_natgov"] <- "Government does enough"
total_reg_results$dep_var_name[total_reg_results$dep_var == "doprot_eu"] <- "EU does enough"
total_reg_results$dep_var_name[total_reg_results$dep_var == "doprot_city"] <- "City does enough"
total_reg_results$dep_var_name[total_reg_results$dep_var == "cc_exag"] <- "Seriousness exaggerated"
total_reg_results$dep_var_name[total_reg_results$dep_var == "role_ind"] <- "Individuals can play role"
total_reg_results$dep_var_name[total_reg_results$dep_var == "eff_daily"] <- "Env issues affect daily life"
total_reg_results$dep_var_name[total_reg_results$dep_var == "pers_imp"] <- "Personal importance"
total_reg_results$dep_var_name[total_reg_results$dep_var == "cchange_mot"] <- "Think about climate change"
total_reg_results$dep_var_name[total_reg_results$dep_var == "worried_mot"] <- "Worry status climate"
total_reg_results$dep_var_name[total_reg_results$dep_var == "futuregen"] <- "Worry future generations"
total_reg_results$dep_var_name[total_reg_results$dep_var == "nowor"] <- "Don't worry global warming"
total_reg_results$dep_var_name[total_reg_results$dep_var == "ontime"] <- "Still on time to do sth"
total_reg_results$dep_var_name[total_reg_results$dep_var == "gov"] <- "Evaluation government action"
total_reg_results$dep_var_name[total_reg_results$dep_var == "resp_gov"] <- "Responsibility gov"
total_reg_results$dep_var_name[total_reg_results$dep_var == "resp_comp"] <- "Responsibility big comp"
total_reg_results$dep_var_name[total_reg_results$dep_var == "resp_mkb"] <- "Responsibility SME"
total_reg_results$dep_var_name[total_reg_results$dep_var == "resp_citiz_mot"] <- "Responsibility citizen"
total_reg_results$dep_var_name[total_reg_results$dep_var == "resp_you"] <- "Your responsibility"
total_reg_results$dep_var_name[total_reg_results$dep_var == "pers_resp_mot"] <- "Your responsibility CO2"
total_reg_results$dep_var_name[total_reg_results$dep_var == "sust_choice"] <- "Sustainability important"
total_reg_results$dep_var_name[total_reg_results$dep_var == "contr"] <- "Contribute to sustainable society"
total_reg_results$dep_var_name[total_reg_results$dep_var == "energy"] <- "Energy choices fight climate change"
total_reg_results$dep_var_name[total_reg_results$dep_var == "noidea"] <- "No idea what to do"
total_reg_results$dep_var_name[total_reg_results$dep_var == "motiv"] <- "Motivated to fight climate change"
total_reg_results$dep_var_name[total_reg_results$dep_var == "fut_gen_socon"] <- "Worry future generations"
total_reg_results$dep_var_name[total_reg_results$dep_var == "lifestyle"] <- "Importance changing lifestyle"


# I also calculated the percentage of missings per variable, in how many waves it was asked and how many categories the original scale had. I collected these data in Excel.
metavar <- readxl::read_excel("./data/meta_var.xlsx")
total_reg_results <- merge(total_reg_results, metavar)


table(total_reg_results$data)
save(total_reg_results, file= "./data/meta_analysis/total_reg_results_gam_w_new.RData" )