Gamlss regressions shorter time spans

In this script, I use loops to perform weighted gamlss regressions with shorter time spans, to use as a robustness check.

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

EVS

load(here("./data/final_data", "evssel.RData"))

# Select variables without missings 
evs_complete <- evssel %>% dplyr::select(surveyyear:marstat, weightvec, isced_cat)
evs_complete <- evs_complete %>% drop_na()
evs_complete$time <- evs_complete$surveyyear - min(evs_complete$surveyyear)

# Subset the data so that I only calculate the time effect between the first two waves
evs_1 <- subset(evs_complete, surveyyear == 1990 | surveyyear == 1999)

predictors <- c("time")
deps <- c("climate5")

evs_list_empty <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  evs_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = evs_1, weights = weightvec)
  
  evs_list_empty[[i]] <- evs_fit
} 

summary(evs_list_empty[[1]]) #Sig pos effect of time on intercept

# Save it, the 1 indicating that this is the first time effect
save(evs_list_empty, file="./data/final_data/regression_outputs/per_wave/evs_list_empty1.RData")

# Second model with predictors
predictors <- c("time", "sex", "isced_cat", "age")
deps <- c("climate5")

evs_list_preds <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  evs_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula,
                    family=NO(mu.link = "identity", sigma.link = "identity"), data = evs_1, weights = weightvec)
  
  evs_list_preds[[i]] <- evs_fit
} 

summary(evs_list_preds[[1]]) #Sig pos effect of time on intercept

save(evs_list_preds, file="./data/final_data/regression_outputs/per_wave/evs_list_preds1_new.RData")

# Third model with interaction effects as well
evs_1$age_time <- evs_1$age * evs_1$time
evs_1$sex_time <- as.numeric(evs_1$sex) * evs_1$time
evs_1 <- fastDummies::dummy_cols(evs_1, select_columns = "isced_cat")
evs_1$isced_basic_time <- evs_1$isced_cat_Basic * evs_1$time
evs_1$isced_intermediate_time <- evs_1$isced_cat_Intermediate * evs_1$time
evs_1$isced_advanced_time <- evs_1$isced_cat_Advanced * evs_1$time

interactions <- c("age_time", "sex_time", "isced_intermediate_time", "isced_advanced_time")

evs_list_interactions <- list()

 for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = " + ")))
  # Perform the linear regression
  evs_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = "+ "))), sigma.formula = sigma.formula," + ", family=NO(mu.link = "identity", sigma.link = "identity"), data = evs_1, weights = weightvec)
  
  evs_list_interactions[[i]] <- evs_fit
} 

summary(evs_list_interactions[[1]]) 

save(evs_list_interactions, file="./data/final_data/regression_outputs/per_wave/evs_list_interactions1_new.RData")

# Now the other waves
# Select the second and third wave to calculate the time effect between these two
evs_2 <- subset(evs_complete, surveyyear == 1999 | surveyyear == 2008)

predictors <- c("time")
deps <- c("climate5")

evs_list_empty <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regrevsion
  evs_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = evs_2, weights = weightvec)
  
  evs_list_empty[[i]] <- evs_fit
} 

summary(evs_list_empty[[1]]) #Sig pos effect of time on intercept

save(evs_list_empty, file="./data/final_data/regression_outputs/per_wave/evs_list_empty2.RData")

# Second model with predictors
predictors <- c("time", "sex", "isced_cat", "age")
deps <- c("climate5")

evs_list_preds <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  evs_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula,
                    family=NO(mu.link = "identity", sigma.link = "identity"), data = evs_2, weights = weightvec)
  
  evs_list_preds[[i]] <- evs_fit
} 

summary(evs_list_preds[[1]]) 

save(evs_list_preds, file="./data/final_data/regression_outputs/per_wave/evs_list_preds2_new.RData")

# Third model with interactions
evs_2$age_time <- evs_2$age * evs_2$time
evs_2$sex_time <- as.numeric(evs_2$sex) * evs_2$time
evs_2 <- fastDummies::dummy_cols(evs_2, select_columns = "isced_cat")
evs_2$isced_basic_time <- evs_2$isced_cat_Basic * evs_2$time
evs_2$isced_intermediate_time <- evs_2$isced_cat_Intermediate * evs_2$time
evs_2$isced_advanced_time <- evs_2$isced_cat_Advanced * evs_2$time

interactions <- c("age_time", "sex_time", "isced_intermediate_time", "isced_advanced_time")

evs_list_interactions <- list()

 for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = " + ")))
  # Perform the linear regression
  evs_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = "+ "))), sigma.formula = sigma.formula," + ", family=NO(mu.link = "identity", sigma.link = "identity"), data = evs_2, weights = weightvec)
  
  evs_list_interactions[[i]] <- evs_fit
} 


summary(evs_list_interactions[[1]]) 

save(evs_list_interactions, file="./data/final_data/regression_outputs/per_wave/evs_list_interactions2_new.RData")

I&O Research

# Repeat the same process for I&O research
load(here("./data/final_data", "io_total.RData"))
io_complete <- io_total %>% dplyr::select(worried:surveyyear, weightvec, isced_cat)
io_complete = subset(io_complete, select = -c(urban,IOINKOMEN, IOPOL2017) )
io_complete <- io_complete %>% drop_na()

io_complete$time <- io_complete$surveyyear - min(io_complete$surveyyear)

io_1 <- subset(io_complete, surveyyear == 2019 | surveyyear == 2020)

predictors <- c("time")
deps <- c("worried", "frontrunner", "worry_future", "resp_citiz", "dk_start", "do_gov", "buss_help", "min_contr", "human_resp")

io_list_empty <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  io_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = io_1, weights = weightvec)
  
  io_list_empty[[i]] <- io_fit
} 

save(io_list_empty, file="./data/final_data/regression_outputs/per_wave/io_list_empty1.RData")

# Model with predictors
predictors <- c("time", "sex", "isced_cat", "age")

io_list_preds <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  io_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula,
                   family=NO(mu.link = "identity", sigma.link = "identity"), data = io_1, weights = weightvec)
  
  io_list_preds[[i]] <- io_fit
} 

summary(io_list_preds[[1]]) 

save(io_list_preds, file="./data/final_data/regression_outputs/per_wave/io_list_preds1_new.RData")

# Model with interactions
io_1$age_time <- io_1$age * io_1$time
io_1$sex_time <- as.numeric(io_1$sex) * io_1$time
io_1 <- fastDummies::dummy_cols(io_1, select_columns = "isced_cat")
io_1$isced_basic_time <- io_1$isced_cat_Basic * io_1$time
io_1$isced_intermediate_time <- io_1$isced_cat_Intermediate * io_1$time
io_1$isced_advanced_time <- io_1$isced_cat_Advanced * io_1$time

interactions <- c("age_time", "sex_time", "isced_intermediate_time", "isced_advanced_time")

io_list_interactions <- list()

 for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = " + ")))
  # Perform the linear regression
  io_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = "+ "))), sigma.formula = sigma.formula," + ", family=NO(mu.link = "identity", sigma.link = "identity"), data = io_1, weights = weightvec)
  
  io_list_interactions[[i]] <- io_fit
} 

summary(io_list_interactions[[1]]) 

save(io_list_interactions, file="./data/final_data/regression_outputs/per_wave/io_list_interactions1_new.RData")

# The other waves
io_2 <- subset(io_complete, surveyyear == 2020 | surveyyear == 2022)

predictors <- c("time")
deps <- c("worried", "frontrunner", "worry_future", "resp_citiz", "dk_start", "do_gov", "buss_help", "min_contr", "human_resp")

io_list_empty <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  io_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = io_2, weights = weightvec)
  
  io_list_empty[[i]] <- io_fit
} 

save(io_list_empty, file="./data/final_data/regression_outputs/per_wave/io_list_empty2.RData")

# Model with predictors
predictors <- c("time", "sex", "isced_cat", "age")

io_list_preds <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  io_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula,
                   family=NO(mu.link = "identity", sigma.link = "identity"), data = io_2, weights = weightvec)
  
  io_list_preds[[i]] <- io_fit
} 

summary(io_list_preds[[1]]) 

save(io_list_preds, file="./data/final_data/regression_outputs/per_wave/io_list_preds2_new.RData")

# Model with interactions
io_2$age_time <- io_2$age * io_2$time
io_2$sex_time <- as.numeric(io_2$sex) * io_2$time
io_2 <- fastDummies::dummy_cols(io_2, select_columns = "isced_cat")
io_2$isced_basic_time <- io_2$isced_cat_Basic * io_2$time
io_2$isced_intermediate_time <- io_2$isced_cat_Intermediate * io_2$time
io_2$isced_advanced_time <- io_2$isced_cat_Advanced * io_2$time

interactions <- c("age_time", "sex_time", "isced_intermediate_time", "isced_advanced_time")

io_list_interactions <- list()

 for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = " + ")))
  # Perform the linear regression
  io_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = "+ "))), sigma.formula = sigma.formula," + ", family=NO(mu.link = "identity", sigma.link = "identity"), data = io_2, weights = weightvec)
  
  io_list_interactions[[i]] <- io_fit
} 


summary(io_list_interactions[[1]]) 

save(io_list_interactions, file="./data/final_data/regression_outputs/per_wave/io_list_interactions2_new.RData")

ISSP

# And the same procedure for ISSP
rm(list=ls())
load("./data/final_data/issptotal.Rdata")
issptotal$surveyyear[is.na(issptotal$surveyyear)] <- "2010"
issptotal$surveyyear <- as.numeric(issptotal$surveyyear)

# The issp has variables that are asked in all three waves, and variables that are asked only in 2 waves. So I have to run the code on subsets of the data. First I estimate the models for the statements that are asked throughout the three waves. The ones that are asked two waves don't have to be run again. 
issp_complete <- issptotal %>% dplyr::select(surveyyear:marstat, weightvec, isced_cat)
issp_complete <- issp_complete %>% drop_na()
issp_complete$time <- issp_complete$surveyyear - min(issp_complete$surveyyear)
issp_1 <- subset(issp_complete, surveyyear == 1993 | surveyyear == 2000)

predictors <- c("time")
deps <- c("worry", "lifeharm", "willing_price", "willing_tax", "willing_living", "do_right", "progharm", "econprotect", "people_decide")

issp_list_empty <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  issp_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula,
                     family=NO(mu.link = "identity", sigma.link = "identity"), data = issp_1, weights = weightvec)
  
  issp_list_empty[[i]] <- issp_fit
} 

save(issp_list_empty, file="./data/final_data/regression_outputs/per_wave/issp_list_empty1.RData")

# Now with predictors
predictors <- c("time", "sex", "isced_cat", "age")
deps <- c("worry", "lifeharm", "willing_price", "willing_tax", "willing_living", "do_right", "progharm", "econprotect", "people_decide" )

issp_list_preds <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  issp_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula,
                      family=NO(mu.link = "identity", sigma.link = "identity"), data = issp_1, weights = weightvec)
  
  issp_list_preds[[i]] <- issp_fit
} 

summary(issp_list_preds[[2]]) 
summary(issp_list_preds[[3]]) 

save(issp_list_preds, file="./data/final_data/regression_outputs/per_wave/issp_list_preds1_new.RData")

# And now interactions
issp_1$age_time <- issp_1$age * issp_1$time
issp_1$sex_time <- as.numeric(issp_1$sex) * issp_1$time
issp_1 <- fastDummies::dummy_cols(issp_1, select_columns = "isced_cat")
issp_1$isced_basic_time <- issp_1$isced_cat_Basic * issp_1$time
issp_1$isced_intermediate_time <- issp_1$isced_cat_Intermediate * issp_1$time
issp_1$isced_advanced_time <- issp_1$isced_cat_Advanced * issp_1$time

interactions <- c("age_time", "sex_time", "isced_intermediate_time", "isced_advanced_time")

issp_list_interactions <- list()

 for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = " + ")))
  # Perform the linear regressisspn
  issp_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = "+ "))), sigma.formula = sigma.formula," + ", family=NO(mu.link = "identity", sigma.link = "identity"), data = issp_1, weights = weightvec)
  
  issp_list_interactions[[i]] <- issp_fit
} 

summary(issp_list_interactions[[1]]) 

save(issp_list_interactions, file="./data/final_data/regression_outputs/per_wave/issp_list_interactions1_new.RData")

# Second part
issp_2 <- subset(issp_complete, surveyyear == 2000 | surveyyear == 2010)

predictors <- c("time")
deps <- c("worry", "lifeharm", "willing_price", "willing_tax", "willing_living", "do_right", "progharm", "econprotect", "people_decide")

issp_list_empty <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  issp_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula,
                     family=NO(mu.link = "identity", sigma.link = "identity"), data = issp_2, weights = weightvec)
  
  issp_list_empty[[i]] <- issp_fit
} 

save(issp_list_empty, file="./data/final_data/regression_outputs/per_wave/issp_list_empty2.RData")

# Now with predictors
predictors <- c("time", "sex", "isced_cat", "age")
deps <- c("worry", "lifeharm", "willing_price", "willing_tax", "willing_living", "do_right", "progharm", "econprotect", "people_decide" )

issp_list_preds <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  issp_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula,
                      family=NO(mu.link = "identity", sigma.link = "identity"), data = issp_2, weights = weightvec)
  
  issp_list_preds[[i]] <- issp_fit
} 

save(issp_list_preds, file="./data/final_data/regression_outputs/per_wave/issp_list_preds2_new.RData")

# And now interactions
issp_2$age_time <- issp_2$age * issp_2$time
issp_2$sex_time <- as.numeric(issp_2$sex) * issp_2$time
issp_2 <- fastDummies::dummy_cols(issp_2, select_columns = "isced_cat")
issp_2$isced_basic_time <- issp_2$isced_cat_Basic * issp_2$time
issp_2$isced_intermediate_time <- issp_2$isced_cat_Intermediate * issp_2$time
issp_2$isced_advanced_time <- issp_2$isced_cat_Advanced * issp_2$time

interactions <- c("age_time", "sex_time", "isced_intermediate_time", "isced_advanced_time")

issp_list_interactions <- list()

 for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = " + ")))
  # Perform the linear regressisspn
  issp_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = "+ "))), sigma.formula = sigma.formula," + ", family=NO(mu.link = "identity", sigma.link = "identity"), data = issp_2, weights = weightvec)
  
  issp_list_interactions[[i]] <- issp_fit
} 

summary(issp_list_interactions[[1]]) 

save(issp_list_interactions, file="./data/final_data/regression_outputs/per_wave/issp_list_interactions2_new.RData")

EB

rm(list=ls())
# Now the most difficult one: Eurobarometer. No items are asked steadily throughout time, so have to work with a lot of subsets. 
load("./data/final_data/eb_tot.RData")

# first 1986 - 1995
eb_complete <- eb_tot %>% dplyr::select(surveyyear, env_ec_stat, env_prsimp, sex, isced: urban, weightvec, isced_cat)
# Religion is only asked in one wave
eb_complete <- eb_complete %>% drop_na()
eb_complete$time <- eb_complete$surveyyear - min(eb_complete$surveyyear)
eb_1 <- subset(eb_complete, surveyyear == 1986 | surveyyear == 1992)

# First empty model 
predictors <- c("time")
deps <- c("env_ec_stat", "env_prsimp")

eb_list_empty_1986 <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula,
                   family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_1, weights = weightvec)
  
  eb_list_empty_1986[[i]] <- eb_fit
} 

summary(eb_list_empty_1986[[1]]) 
summary(eb_list_empty_1986[[2]])

save(eb_list_empty_1986, file="./data/final_data/regression_outputs/per_wave/eb_list_empty_1986_1.RData")

# Model with predictors
predictors <- c("time", "sex", "isced_cat", "age")
deps <- c("env_ec_stat", "env_prsimp")

eb_list_preds_1986 <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula,
                   family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_1, weights = weightvec)
  
  eb_list_preds_1986[[i]] <- eb_fit
} 

summary(eb_list_preds_1986[[1]]) 
summary(eb_list_preds_1986[[2]])

save(eb_list_preds_1986, file="./data/final_data/regression_outputs/per_wave/eb_list_preds_1986_1_new.RData")

# And the model with interactions
eb_1$age_time <- eb_1$age * eb_1$time
eb_1$sex_time <- as.numeric(eb_1$sex) * eb_1$time
eb_1 <- fastDummies::dummy_cols(eb_1, select_columns = "isced_cat")
eb_1$isced_basic_time <- eb_1$isced_cat_Basic * eb_1$time
eb_1$isced_intermediate_time <- eb_1$isced_cat_Intermediate * eb_1$time
eb_1$isced_advanced_time <- eb_1$isced_cat_Advanced * eb_1$time

interactions <- c("age_time", "sex_time", "isced_intermediate_time", "isced_advanced_time")

eb_list_interactions_1986 <- list()

 for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = "+ "))), sigma.formula = sigma.formula," + ", family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_1, weights = weightvec)
  
  eb_list_interactions_1986[[i]] <- eb_fit
} 

summary(eb_list_interactions_1986[[1]]) 
summary(eb_list_interactions_1986[[2]])

save(eb_list_interactions_1986, file="./data/final_data/regression_outputs/per_wave/eb_list_interactions_1986_1.RData")

# Second part
eb_2 <- subset(eb_complete, surveyyear == 1992 | surveyyear == 1995)

predictors <- c("time")
deps <- c("env_ec_stat", "env_prsimp")

eb_list_empty_1986 <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula,
                   family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_2, weights = weightvec)
  
  eb_list_empty_1986[[i]] <- eb_fit
} 

summary(eb_list_empty_1986[[1]]) 
summary(eb_list_empty_1986[[2]])

save(eb_list_empty_1986, file="./data/final_data/regression_outputs/per_wave/eb_list_empty_1986_2.RData")

# Model with predictors
predictors <- c("time", "sex", "isced_cat", "age")
deps <- c("env_ec_stat", "env_prsimp")

eb_list_preds_1986 <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula,
                   family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_2, weights = weightvec)
  
  eb_list_preds_1986[[i]] <- eb_fit
} 

summary(eb_list_preds_1986[[1]]) 
summary(eb_list_preds_1986[[2]])

save(eb_list_preds_1986, file="./data/final_data/regression_outputs/per_wave/eb_list_preds_1986_2_new.RData")

# And the model with interactions
eb_2$age_time <- eb_2$age * eb_2$time
eb_2$sex_time <- as.numeric(eb_2$sex) * eb_2$time
eb_2 <- fastDummies::dummy_cols(eb_2, select_columns = "isced_cat")
eb_2$isced_basic_time <- eb_2$isced_cat_Basic * eb_2$time
eb_2$isced_intermediate_time <- eb_2$isced_cat_Intermediate * eb_2$time
eb_2$isced_advanced_time <- eb_2$isced_cat_Advanced * eb_2$time

interactions <- c("age_time", "sex_time", "isced_intermediate_time", "isced_advanced_time")


eb_list_interactions_1986 <- list()

 for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = "+ "))), sigma.formula = sigma.formula," + ", family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_2, weights = weightvec)
  
  eb_list_interactions_1986[[i]] <- eb_fit
} 


summary(eb_list_interactions_1986[[1]]) 
summary(eb_list_interactions_1986[[2]])

save(eb_list_interactions_1986, file="./data/final_data/regression_outputs/per_wave/eb_list_interactions_1986_2.RData")

# Now with the influence of environment on quality of life. This is asked 4 times, so I have to make 3 subsets. 
eb_complete <- eb_tot %>% dplyr::select(surveyyear, env_quallife, sex, isced, urban, age, weightvec, isced_cat)
# Lrplace and income are not asked in all waves

eb_complete <- eb_complete %>% drop_na()

eb_complete$time <- eb_complete$surveyyear - min(eb_complete$surveyyear)

eb_1 <- subset(eb_complete, surveyyear == 2004 | surveyyear == 2007)

# First empty model 
predictors <- c("time")
deps <- c("env_quallife")

eb_list_empty_quallife <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_1, weights = weightvec)
  
  eb_list_empty_quallife[[i]] <- eb_fit
} 

summary(eb_list_empty_quallife[[1]]) 

save(eb_list_empty_quallife, file="./data/final_data/regression_outputs/per_wave/eb_list_empty_quallife1.RData")

# Model with predictors
predictors <- c("time", "sex", "isced_cat", "age")
deps <- c("env_quallife")

eb_list_preds_quallife <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula,
                   family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_1, weights = weightvec)
  
  eb_list_preds_quallife[[i]] <- eb_fit
} 

summary(eb_list_preds_quallife[[1]]) 

save(eb_list_preds_quallife, file="./data/final_data/regression_outputs/per_wave/eb_list_preds_quallife1_new.RData")

# And the model with interactions
eb_1$age_time <- eb_1$age * eb_1$time
eb_1$sex_time <- as.numeric(eb_1$sex) * eb_1$time
eb_1 <- fastDummies::dummy_cols(eb_1, select_columns = "isced_cat")
eb_1$isced_basic_time <- eb_1$isced_cat_Basic * eb_1$time
eb_1$isced_intermediate_time <- eb_1$isced_cat_Intermediate * eb_1$time
eb_1$isced_advanced_time <- eb_1$isced_cat_Advanced * eb_1$time

interactions <- c("age_time", "sex_time", "isced_intermediate_time", "isced_advanced_time")

eb_list_interactions_quallife <- list()

 for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = "+ "))), sigma.formula = sigma.formula," + ", family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_1, weights = weightvec)
  
  eb_list_interactions_quallife[[i]] <- eb_fit
} 


summary(eb_list_interactions_quallife[[1]]) 

save(eb_list_interactions_quallife, file="./data/final_data/regression_outputs/per_wave/eb_list_interactions_quallife1_new.RData")

# Second part
eb_2 <- subset(eb_complete, surveyyear == 2007 | surveyyear == 2011)

predictors <- c("time")
deps <- c("env_quallife")

eb_list_empty_quallife <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_2, weights = weightvec)
  
  eb_list_empty_quallife[[i]] <- eb_fit
} 

summary(eb_list_empty_quallife[[1]]) 

save(eb_list_empty_quallife, file="./data/final_data/regression_outputs/per_wave/eb_list_empty_quallife2.RData")

# Model with predictors
predictors <- c("time", "sex", "isced_cat", "age")
deps <- c("env_quallife")

eb_list_preds_quallife <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula,
                   family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_2, weights = weightvec)
  
  eb_list_preds_quallife[[i]] <- eb_fit
} 

summary(eb_list_preds_quallife[[1]]) 

save(eb_list_preds_quallife, file="./data/final_data/regression_outputs/per_wave/eb_list_preds_quallife2_new.RData")

# And the model with interactions
eb_2$age_time <- eb_2$age * eb_2$time
eb_2$sex_time <- as.numeric(eb_2$sex) * eb_2$time
eb_2 <- fastDummies::dummy_cols(eb_2, select_columns = "isced_cat")
eb_2$isced_basic_time <- eb_2$isced_cat_Basic * eb_2$time
eb_2$isced_intermediate_time <- eb_2$isced_cat_Intermediate * eb_2$time
eb_2$isced_advanced_time <- eb_2$isced_cat_Advanced * eb_2$time

interactions <- c("age_time", "sex_time", "isced_intermediate_time", "isced_advanced_time")

eb_list_interactions_quallife <- list()

 for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = " + ")))
  # Perform the linear regressebn
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = "+ "))), sigma.formula = sigma.formula," + ", family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_2, weights = weightvec)
  
  eb_list_interactions_quallife[[i]] <- eb_fit
} 

summary(eb_list_interactions_quallife[[1]]) 

save(eb_list_interactions_quallife, file="./data/final_data/regression_outputs/per_wave/eb_list_interactions_quallife2_new.RData")

# And the third set
eb_3 <- subset(eb_complete, surveyyear == 2011 | surveyyear == 2014)

predictors <- c("time")
deps <- c("env_quallife")

eb_list_empty_quallife <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_3, weights = weightvec)
  
  eb_list_empty_quallife[[i]] <- eb_fit
} 

summary(eb_list_empty_quallife[[1]]) 

save(eb_list_empty_quallife, file="./data/final_data/regression_outputs/per_wave/eb_list_empty_quallife3.RData")

# Model with predictors
predictors <- c("time", "sex", "isced_cat", "age")
deps <- c("env_quallife")

eb_list_preds_quallife <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula,
                   family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_3, weights = weightvec)
  
  eb_list_preds_quallife[[i]] <- eb_fit
} 

summary(eb_list_preds_quallife[[1]]) 

save(eb_list_preds_quallife, file="./data/final_data/regression_outputs/per_wave/eb_list_preds_quallife3_new.RData")

# And the model with interactions
eb_3$age_time <- eb_3$age * eb_3$time
eb_3$sex_time <- as.numeric(eb_3$sex) * eb_3$time
eb_3 <- fastDummies::dummy_cols(eb_3, select_columns = "isced_cat")
eb_3$isced_basic_time <- eb_3$isced_cat_Basic * eb_3$time
eb_3$isced_intermediate_time <- eb_3$isced_cat_Intermediate * eb_3$time
eb_3$isced_advanced_time <- eb_3$isced_cat_Advanced * eb_3$time

interactions <- c("age_time", "sex_time", "isced_intermediate_time", "isced_advanced_time")

eb_list_interactions_quallife <- list()

 for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = " + ")))
  # Perform the linear regressebn
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = "+ "))), sigma.formula = sigma.formula," + ", family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_3, weights = weightvec)
  
  eb_list_interactions_quallife[[i]] <- eb_fit
} 

summary(eb_list_interactions_quallife[[1]]) 

save(eb_list_interactions_quallife, file="./data/final_data/regression_outputs/per_wave/eb_list_interactions_quallife3_new.RData")

# Now 2 statements that are asked 4 times
eb_complete <- eb_tot %>% dplyr::select(surveyyear, role_ind, big_pol, sex, isced, urban, age, weightvec, isced_cat)
eb_complete <- eb_complete %>% drop_na()
eb_complete$time <- eb_complete$surveyyear - min(eb_complete$surveyyear)
eb_1 <- subset(eb_complete, surveyyear == 2007 | surveyyear == 2011)

# First empty model 
predictors <- c("time")
deps <- c("role_ind", "big_pol")

eb_list_empty_2007 <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_1, weights = weightvec)
  
  eb_list_empty_2007[[i]] <- eb_fit
} 

summary(eb_list_empty_2007[[2]]) 

save(eb_list_empty_2007, file="./data/final_data/regression_outputs/eb_list_empty_2007_1.RData")

# Model with predictors
predictors <- c("time", "sex", "isced_cat", "age")
deps <- c("role_ind", "big_pol")

eb_list_preds_2007 <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_1, weights = weightvec)
  
  eb_list_preds_2007[[i]] <- eb_fit
} 

summary(eb_list_preds_2007[[1]]) 

save(eb_list_preds_2007, file="./data/final_data/regression_outputs/eb_list_preds_2007_1_new.RData")

# And the model with interactions
eb_1$age_time <- eb_1$age * eb_1$time
eb_1$sex_time <- as.numeric(eb_1$sex) * eb_1$time
eb_1 <- fastDummies::dummy_cols(eb_1, select_columns = "isced_cat")
eb_1$isced_basic_time <- eb_1$isced_cat_Basic * eb_1$time
eb_1$isced_intermediate_time <- eb_1$isced_cat_Intermediate * eb_1$time
eb_1$isced_advanced_time <- eb_1$isced_cat_Advanced * eb_1$time

interactions <- c("age_time", "sex_time", "isced_intermediate_time", "isced_advanced_time")

eb_list_interactions_2007 <- list()

 for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = " + ")))
  # Perform the linear regressebn
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = "+ "))), sigma.formula = sigma.formula," + ", family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_1, weights = weightvec)
  
  eb_list_interactions_2007[[i]] <- eb_fit
} 

summary(eb_list_interactions_2007[[1]]) 

save(eb_list_interactions_2007, file="./data/final_data/regression_outputs/eb_list_interactions_2007_1_new.RData")

# Second set
eb_2 <- subset(eb_complete, surveyyear == 2011 | surveyyear == 2014)

predictors <- c("time")
deps <- c("role_ind", "big_pol")

eb_list_empty_2007 <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_2, weights = weightvec)
  
  eb_list_empty_2007[[i]] <- eb_fit
} 

summary(eb_list_empty_2007[[2]]) 

save(eb_list_empty_2007, file="./data/final_data/regression_outputs/eb_list_empty_2007_2.RData")

# Model with predictors
predictors <- c("time", "sex", "isced_cat", "age")
deps <- c("role_ind", "big_pol")

eb_list_preds_2007 <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_2, weights = weightvec)
  
  eb_list_preds_2007[[i]] <- eb_fit
} 

summary(eb_list_preds_2007[[1]]) 

save(eb_list_preds_2007, file="./data/final_data/regression_outputs/eb_list_preds_2007_2_new.RData")

# And the model with interactions
eb_2$age_time <- eb_2$age * eb_2$time
eb_2$sex_time <- as.numeric(eb_2$sex) * eb_2$time
eb_2 <- fastDummies::dummy_cols(eb_2, select_columns = "isced_cat")
eb_2$isced_basic_time <- eb_2$isced_cat_Basic * eb_2$time
eb_2$isced_intermediate_time <- eb_2$isced_cat_Intermediate * eb_2$time
eb_2$isced_advanced_time <- eb_2$isced_cat_Advanced * eb_2$time

interactions <- c("age_time", "sex_time", "isced_intermediate_time", "isced_advanced_time")

eb_list_interactions_2007 <- list()

 for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = " + ")))
  # Perform the linear regressebn
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = "+ "))), sigma.formula = sigma.formula," + ", family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_2, weights = weightvec)
  
  eb_list_interactions_2007[[i]] <- eb_fit
} 

summary(eb_list_interactions_2007[[1]]) 

save(eb_list_interactions_2007, file="./data/final_data/regression_outputs/eb_list_interactions_2007_2_new.RData")

# Third set
eb_3 <- subset(eb_complete, surveyyear == 2014 | surveyyear == 2017)

predictors <- c("time")
deps <- c("role_ind", "big_pol")

eb_list_empty_2007 <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_3, weights = weightvec)
  
  eb_list_empty_2007[[i]] <- eb_fit
} 

summary(eb_list_empty_2007[[2]]) 

save(eb_list_empty_2007, file="./data/final_data/regression_outputs/eb_list_empty_2007_3.RData")

# Model with predictors
predictors <- c("time", "sex", "isced_cat", "age")
deps <- c("role_ind", "big_pol")

eb_list_preds_2007 <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_3, weights = weightvec)
  
  eb_list_preds_2007[[i]] <- eb_fit
} 

summary(eb_list_preds_2007[[1]]) 

save(eb_list_preds_2007, file="./data/final_data/regression_outputs/eb_list_preds_2007_3_new.RData")

# And the model with interactions
eb_3$age_time <- eb_3$age * eb_3$time
eb_3$sex_time <- as.numeric(eb_3$sex) * eb_3$time
eb_3 <- fastDummies::dummy_cols(eb_3, select_columns = "isced_cat")
eb_3$isced_basic_time <- eb_3$isced_cat_Basic * eb_3$time
eb_3$isced_intermediate_time <- eb_3$isced_cat_Intermediate * eb_3$time
eb_3$isced_advanced_time <- eb_3$isced_cat_Advanced * eb_3$time

interactions <- c("age_time", "sex_time", "isced_intermediate_time", "isced_advanced_time")

eb_list_interactions_2007 <- list()

 for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = " + ")))
  # Perform the linear regressebn
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = "+ "))), sigma.formula = sigma.formula," + ", family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_3, weights = weightvec)
  
  eb_list_interactions_2007[[i]] <- eb_fit
} 

summary(eb_list_interactions_2007[[1]]) 

save(eb_list_interactions_2007, file="./data/final_data/regression_outputs/eb_list_interactions_2007_3_new.RData")

# Statement about willingness to buy more expensive products that are environmentally friendly
eb_complete <- eb_tot %>% dplyr::select(surveyyear, buyprod, sex, isced, urban, age, weightvec, isced_cat)
eb_complete <- eb_complete %>% drop_na()
eb_complete$time <- eb_complete$surveyyear - min(eb_complete$surveyyear)
eb_1 <- subset(eb_complete, surveyyear == 2007| surveyyear == 2011)

# First empty model 
predictors <- c("time")
deps <- c("buyprod")

eb_list_empty_buyprod <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_1, weights = weightvec)
  
  eb_list_empty_buyprod[[i]] <- eb_fit
} 

summary(eb_list_empty_buyprod[[1]]) 

save(eb_list_empty_buyprod, file="./data/final_data/regression_outputs/eb_list_empty_buyprod1.RData")

# Model with predictors
predictors <- c("time", "sex", "isced_cat", "age")
deps <- c("buyprod")

eb_list_preds_buyprod <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_1, weights = weightvec)
  
  eb_list_preds_buyprod[[i]] <- eb_fit
} 

summary(eb_list_preds_buyprod[[1]]) 

save(eb_list_preds_buyprod, file="./data/final_data/regression_outputs/eb_list_preds_buyprod1_new.RData")

# And the model with interactions
eb_1$age_time <- eb_1$age * eb_1$time
eb_1$sex_time <- as.numeric(eb_1$sex) * eb_1$time
eb_1 <- fastDummies::dummy_cols(eb_1, select_columns = "isced_cat")
eb_1$isced_basic_time <- eb_1$isced_cat_Basic * eb_1$time
eb_1$isced_intermediate_time <- eb_1$isced_cat_Intermediate * eb_1$time
eb_1$isced_advanced_time <- eb_1$isced_cat_Advanced * eb_1$time

interactions <- c("age_time", "sex_time", "isced_intermediate_time", "isced_advanced_time")

eb_list_interactions_buyprod <- list()

 for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = "+ "))), sigma.formula = sigma.formula," + ", family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_1, weights = weightvec)
  
  eb_list_interactions_buyprod[[i]] <- eb_fit
} 


summary(eb_list_interactions_buyprod[[1]]) 

save(eb_list_interactions_buyprod, file="./data/final_data/regression_outputs/eb_list_interactions_buyprod1_new.RData")

# Second set
eb_2 <- subset(eb_complete, surveyyear == 2011| surveyyear == 2014)

# First empty model 
predictors <- c("time")
deps <- c("buyprod")

eb_list_empty_buyprod <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_2, weights = weightvec)
  
  eb_list_empty_buyprod[[i]] <- eb_fit
} 

summary(eb_list_empty_buyprod[[1]]) 

save(eb_list_empty_buyprod, file="./data/final_data/regression_outputs/eb_list_empty_buyprod2.RData")

# Model with predictors
predictors <- c("time", "sex", "isced_cat", "age")
deps <- c("buyprod")

eb_list_preds_buyprod <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_2, weights = weightvec)
  
  eb_list_preds_buyprod[[i]] <- eb_fit
} 

summary(eb_list_preds_buyprod[[1]]) 

save(eb_list_preds_buyprod, file="./data/final_data/regression_outputs/eb_list_preds_buyprod2_new.RData")

# And the model with interactions
eb_2$age_time <- eb_2$age * eb_2$time
eb_2$sex_time <- as.numeric(eb_2$sex) * eb_2$time
eb_2 <- fastDummies::dummy_cols(eb_2, select_columns = "isced_cat")
eb_2$isced_basic_time <- eb_2$isced_cat_Basic * eb_2$time
eb_2$isced_intermediate_time <- eb_2$isced_cat_Intermediate * eb_2$time
eb_2$isced_advanced_time <- eb_2$isced_cat_Advanced * eb_2$time

interactions <- c("age_time", "sex_time", "isced_intermediate_time", "isced_advanced_time")

eb_list_interactions_buyprod <- list()

 for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = "+ "))), sigma.formula = sigma.formula," + ", family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_2, weights = weightvec)
  
  eb_list_interactions_buyprod[[i]] <- eb_fit
} 


summary(eb_list_interactions_buyprod[[1]]) 

save(eb_list_interactions_buyprod, file="./data/final_data/regression_outputs/eb_list_interactions_buyprod2_new.RData")

# Affected daily
eb_complete <- eb_tot %>% dplyr::select(surveyyear, eff_daily, sex, isced, urban, age, weightvec, isced_cat)
eb_complete <- eb_complete %>% drop_na()
eb_complete$time <- eb_complete$surveyyear - min(eb_complete$surveyyear)
eb_1 <- subset(eb_complete, surveyyear == 2007 | surveyyear == 2011)

# First empty model 
predictors <- c("time")
deps <- c("eff_daily")

eb_list_empty_eff_daily <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_1, weights = weightvec)
  
  eb_list_empty_eff_daily[[i]] <- eb_fit
} 

summary(eb_list_empty_eff_daily[[1]]) 

save(eb_list_empty_eff_daily, file="./data/final_data/regression_outputs/eb_list_empty_eff_daily1.RData")

# Model with predictors
predictors <- c("time", "sex", "isced_cat", "age")
deps <- c("eff_daily")

eb_list_preds_eff_daily <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_1, weights = weightvec)
  
  eb_list_preds_eff_daily[[i]] <- eb_fit
} 

summary(eb_list_preds_eff_daily[[1]]) 

save(eb_list_preds_eff_daily, file="./data/final_data/regression_outputs/eb_list_preds_eff_daily1_new.RData")

# And the model with interactions
eb_1$age_time <- eb_1$age * eb_1$time
eb_1$sex_time <- as.numeric(eb_1$sex) * eb_1$time
eb_1 <- fastDummies::dummy_cols(eb_1, select_columns = "isced_cat")
eb_1$isced_basic_time <- eb_1$isced_cat_Basic * eb_1$time
eb_1$isced_intermediate_time <- eb_1$isced_cat_Intermediate * eb_1$time
eb_1$isced_advanced_time <- eb_1$isced_cat_Advanced * eb_1$time

interactions <- c("age_time", "sex_time", "isced_intermediate_time", "isced_advanced_time")

eb_list_interactions_eff_daily <- list()

 for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = "+ "))), sigma.formula = sigma.formula," + ", family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_1, weights = weightvec)
  
  eb_list_interactions_eff_daily[[i]] <- eb_fit
} 

summary(eb_list_interactions_eff_daily[[1]]) 

save(eb_list_interactions_eff_daily, file="./data/final_data/regression_outputs/eb_list_interactions_eff_daily1_new.RData")

# Second set
eb_2 <- subset(eb_complete, surveyyear == 2011 | surveyyear == 2014)

# First empty model 
predictors <- c("time")
deps <- c("eff_daily")

eb_list_empty_eff_daily <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_2, weights = weightvec)
  
  eb_list_empty_eff_daily[[i]] <- eb_fit
} 

summary(eb_list_empty_eff_daily[[1]]) 

save(eb_list_empty_eff_daily, file="./data/final_data/regression_outputs/eb_list_empty_eff_daily2.RData")

# Model with predictors
predictors <- c("time", "sex", "isced_cat", "age")
deps <- c("eff_daily")

eb_list_preds_eff_daily <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_2, weights = weightvec)
  
  eb_list_preds_eff_daily[[i]] <- eb_fit
} 

summary(eb_list_preds_eff_daily[[1]]) 

save(eb_list_preds_eff_daily, file="./data/final_data/regression_outputs/eb_list_preds_eff_daily2_new.RData")

# And the model with interactions
eb_2$age_time <- eb_2$age * eb_2$time
eb_2$sex_time <- as.numeric(eb_2$sex) * eb_2$time
eb_2 <- fastDummies::dummy_cols(eb_2, select_columns = "isced_cat")
eb_2$isced_basic_time <- eb_2$isced_cat_Basic * eb_2$time
eb_2$isced_intermediate_time <- eb_2$isced_cat_Intermediate * eb_2$time
eb_2$isced_advanced_time <- eb_2$isced_cat_Advanced * eb_2$time

interactions <- c("age_time", "sex_time", "isced_intermediate_time", "isced_advanced_time")

eb_list_interactions_eff_daily <- list()

 for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = "+ "))), sigma.formula = sigma.formula," + ", family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_2, weights = weightvec)
  
  eb_list_interactions_eff_daily[[i]] <- eb_fit
} 


summary(eb_list_interactions_eff_daily[[1]]) 

save(eb_list_interactions_eff_daily, file="./data/final_data/regression_outputs/eb_list_interactions_eff_daily2_new.RData")

# Third set
eb_3 <- subset(eb_complete, surveyyear == 2014 | surveyyear == 2017)

# First empty model 
predictors <- c("time")
deps <- c("eff_daily")

eb_list_empty_eff_daily <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_3, weights = weightvec)
  
  eb_list_empty_eff_daily[[i]] <- eb_fit
} 

summary(eb_list_empty_eff_daily[[1]]) 

save(eb_list_empty_eff_daily, file="./data/final_data/regression_outputs/eb_list_empty_eff_daily3.RData")

# Model with predictors
predictors <- c("time", "sex", "isced_cat", "age")
deps <- c("eff_daily")

eb_list_preds_eff_daily <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_3, weights = weightvec)
  
  eb_list_preds_eff_daily[[i]] <- eb_fit
} 

summary(eb_list_preds_eff_daily[[1]]) 

save(eb_list_preds_eff_daily, file="./data/final_data/regression_outputs/eb_list_preds_eff_daily3_new.RData")

# And the model with interactions
eb_3$age_time <- eb_3$age * eb_3$time
eb_3$sex_time <- as.numeric(eb_3$sex) * eb_3$time
eb_3 <- fastDummies::dummy_cols(eb_3, select_columns = "isced_cat")
eb_3$isced_basic_time <- eb_3$isced_cat_Basic * eb_3$time
eb_3$isced_intermediate_time <- eb_3$isced_cat_Intermediate * eb_3$time
eb_3$isced_advanced_time <- eb_3$isced_cat_Advanced * eb_3$time

interactions <- c("age_time", "sex_time", "isced_intermediate_time", "isced_advanced_time")

eb_list_interactions_eff_daily <- list()

 for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = "+ "))), sigma.formula = sigma.formula," + ", family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_3, weights = weightvec)
  
  eb_list_interactions_eff_daily[[i]] <- eb_fit
} 


summary(eb_list_interactions_eff_daily[[1]]) 

save(eb_list_interactions_eff_daily, file="./data/final_data/regression_outputs/eb_list_interactions_eff_daily3_new.RData")

# Climate change perception asked in six waves, so have to use 5 subsets
eb_complete <- eb_tot %>% dplyr::select(surveyyear, ccpercept, cchange, sex, isced, urban, age, weightvec, isced_cat)
eb_complete <- eb_complete %>% drop_na()
eb_complete$time <- eb_complete$surveyyear - min(eb_complete$surveyyear)


eb_1 <- subset(eb_complete, surveyyear == 2009 | surveyyear == 2011)
# First empty model 
predictors <- c("time")
deps <- c("ccpercept", "cchange")

eb_list_empty_2009 <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_1, weights = weightvec)
  
  eb_list_empty_2009[[i]] <- eb_fit
} 

summary(eb_list_empty_2009[[1]]) 

save(eb_list_empty_2009, file="./data/final_data/regression_outputs/eb_list_empty_2009_1.RData")

# Model with predictors
predictors <- c("time", "sex", "isced_cat", "age")
deps <- c("ccpercept", "cchange")

eb_list_preds_2009 <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_1, weights = weightvec)
  
  eb_list_preds_2009[[i]] <- eb_fit
} 

summary(eb_list_preds_2009[[1]]) 

save(eb_list_preds_2009, file="./data/final_data/regression_outputs/eb_list_preds_2009_1_new.RData")

# And the model with interactions
eb_1$age_time <- eb_1$age * eb_1$time
eb_1$sex_time <- as.numeric(eb_1$sex) * eb_1$time
eb_1 <- fastDummies::dummy_cols(eb_1, select_columns = "isced_cat")
eb_1$isced_basic_time <- eb_1$isced_cat_Basic * eb_1$time
eb_1$isced_intermediate_time <- eb_1$isced_cat_Intermediate * eb_1$time
eb_1$isced_advanced_time <- eb_1$isced_cat_Advanced * eb_1$time

interactions <- c("age_time", "sex_time", "isced_intermediate_time", "isced_advanced_time")

eb_list_interactions_2009 <- list()

 for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = "+ "))), sigma.formula = sigma.formula," + ", family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_1, weights = weightvec)
  
  eb_list_interactions_2009[[i]] <- eb_fit
} 

summary(eb_list_interactions_2009[[1]]) 

save(eb_list_interactions_2009, file="./data/final_data/regression_outputs/eb_list_interactions_2009_1_new.RData")

# Second set 
eb_2 <- subset(eb_complete, surveyyear == 2011 | surveyyear == 2013)
# First empty model 
predictors <- c("time")
deps <- c("ccpercept", "cchange")

eb_list_empty_2009 <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_2, weights = weightvec)
  
  eb_list_empty_2009[[i]] <- eb_fit
} 

summary(eb_list_empty_2009[[1]]) 

save(eb_list_empty_2009, file="./data/final_data/regression_outputs/eb_list_empty_2009_2.RData")

# Model with predictors
predictors <- c("time", "sex", "isced_cat", "age")
deps <- c("ccpercept", "cchange")

eb_list_preds_2009 <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_2, weights = weightvec)
  
  eb_list_preds_2009[[i]] <- eb_fit
} 

summary(eb_list_preds_2009[[1]]) 

save(eb_list_preds_2009, file="./data/final_data/regression_outputs/eb_list_preds_2009_2_new.RData")

# And the model with interactions
eb_2$age_time <- eb_2$age * eb_2$time
eb_2$sex_time <- as.numeric(eb_2$sex) * eb_2$time
eb_2 <- fastDummies::dummy_cols(eb_2, select_columns = "isced_cat")
eb_2$isced_basic_time <- eb_2$isced_cat_Basic * eb_2$time
eb_2$isced_intermediate_time <- eb_2$isced_cat_Intermediate * eb_2$time
eb_2$isced_advanced_time <- eb_2$isced_cat_Advanced * eb_2$time

interactions <- c("age_time", "sex_time", "isced_intermediate_time", "isced_advanced_time")


eb_list_interactions_2009 <- list()

 for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = "+ "))), sigma.formula = sigma.formula," + ", family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_2, weights = weightvec)
  
  eb_list_interactions_2009[[i]] <- eb_fit
} 

summary(eb_list_interactions_2009[[1]]) 

save(eb_list_interactions_2009, file="./data/final_data/regression_outputs/eb_list_interactions_2009_2_new.RData")

# Third set 
eb_3 <- subset(eb_complete, surveyyear == 2013 | surveyyear == 2015)
# First empty model 
predictors <- c("time")
deps <- c("ccpercept", "cchange")

eb_list_empty_2009 <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_3, weights = weightvec)
  
  eb_list_empty_2009[[i]] <- eb_fit
} 

summary(eb_list_empty_2009[[1]]) 

save(eb_list_empty_2009, file="./data/final_data/regression_outputs/eb_list_empty_2009_3.RData")

# Model with predictors
predictors <- c("time", "sex", "isced_cat", "age")
deps <- c("ccpercept", "cchange")

eb_list_preds_2009 <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_3, weights = weightvec)
  
  eb_list_preds_2009[[i]] <- eb_fit
} 

summary(eb_list_preds_2009[[1]]) 

save(eb_list_preds_2009, file="./data/final_data/regression_outputs/eb_list_preds_2009_3_new.RData")

# And the model with interactions
eb_3$age_time <- eb_3$age * eb_3$time
eb_3$sex_time <- as.numeric(eb_3$sex) * eb_3$time
eb_3 <- fastDummies::dummy_cols(eb_3, select_columns = "isced_cat")
eb_3$isced_basic_time <- eb_3$isced_cat_Basic * eb_3$time
eb_3$isced_intermediate_time <- eb_3$isced_cat_Intermediate * eb_3$time
eb_3$isced_advanced_time <- eb_3$isced_cat_Advanced * eb_3$time

interactions <- c("age_time", "sex_time", "isced_intermediate_time", "isced_advanced_time")

eb_list_interactions_2009 <- list()

 for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = "+ "))), sigma.formula = sigma.formula," + ", family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_3, weights = weightvec)
  
  eb_list_interactions_2009[[i]] <- eb_fit
} 

summary(eb_list_interactions_2009[[1]]) 

save(eb_list_interactions_2009, file="./data/final_data/regression_outputs/eb_list_interactions_2009_3_new.RData")

# Fourth set
eb_4 <- subset(eb_complete, surveyyear == 2015 | surveyyear == 2017)
# First empty model 
predictors <- c("time")
deps <- c("ccpercept", "cchange")

eb_list_empty_2009 <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_4, weights = weightvec)
  
  eb_list_empty_2009[[i]] <- eb_fit
} 

summary(eb_list_empty_2009[[1]]) 

save(eb_list_empty_2009, file="./data/final_data/regression_outputs/eb_list_empty_2009_4.RData")

# Model with predictors
predictors <- c("time", "sex", "isced_cat", "age")
deps <- c("ccpercept", "cchange")

eb_list_preds_2009 <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_4, weights = weightvec)
  
  eb_list_preds_2009[[i]] <- eb_fit
} 

summary(eb_list_preds_2009[[1]]) 

save(eb_list_preds_2009, file="./data/final_data/regression_outputs/eb_list_preds_2009_4_new.RData")

# And the model with interactions
eb_4$age_time <- eb_4$age * eb_4$time
eb_4$sex_time <- as.numeric(eb_4$sex) * eb_4$time
eb_4 <- fastDummies::dummy_cols(eb_4, select_columns = "isced_cat")
eb_4$isced_basic_time <- eb_4$isced_cat_Basic * eb_4$time
eb_4$isced_intermediate_time <- eb_4$isced_cat_Intermediate * eb_4$time
eb_4$isced_advanced_time <- eb_4$isced_cat_Advanced * eb_4$time

interactions <- c("age_time", "sex_time", "isced_intermediate_time", "isced_advanced_time")

eb_list_interactions_2009 <- list()

 for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = "+ "))), sigma.formula = sigma.formula," + ", family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_4, weights = weightvec)
  
  eb_list_interactions_2009[[i]] <- eb_fit
} 

summary(eb_list_interactions_2009[[1]]) 

save(eb_list_interactions_2009, file="./data/final_data/regression_outputs/eb_list_interactions_2009_4_new.RData")

# Fifth set
eb_5 <- subset(eb_complete, surveyyear == 2017 | surveyyear == 2021)
# First empty model 
predictors <- c("time")
deps <- c("ccpercept", "cchange")

eb_list_empty_2009 <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_5, weights = weightvec)
  
  eb_list_empty_2009[[i]] <- eb_fit
} 

summary(eb_list_empty_2009[[1]]) 

save(eb_list_empty_2009, file="./data/final_data/regression_outputs/eb_list_empty_2009_5.RData")

# Model with predictors
predictors <- c("time", "sex", "isced_cat", "age")
deps <- c("ccpercept", "cchange")

eb_list_preds_2009 <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_5, weights = weightvec)
  
  eb_list_preds_2009[[i]] <- eb_fit
} 

summary(eb_list_preds_2009[[1]]) 

save(eb_list_preds_2009, file="./data/final_data/regression_outputs/eb_list_preds_2009_5_new.RData")

# And the model with interactions
eb_5$age_time <- eb_5$age * eb_5$time
eb_5$sex_time <- as.numeric(eb_5$sex) * eb_5$time
eb_5 <- fastDummies::dummy_cols(eb_5, select_columns = "isced_cat")
eb_5$isced_basic_time <- eb_5$isced_cat_Basic * eb_5$time
eb_5$isced_intermediate_time <- eb_5$isced_cat_Intermediate * eb_5$time
eb_5$isced_advanced_time <- eb_5$isced_cat_Advanced * eb_5$time

interactions <- c("age_time", "sex_time", "isced_intermediate_time", "isced_advanced_time")

eb_list_interactions_2009 <- list()

 for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = "+ "))), sigma.formula = sigma.formula," + ", family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_5, weights = weightvec)
  
  eb_list_interactions_2009[[i]] <- eb_fit
} 

summary(eb_list_interactions_2009[[1]]) 

save(eb_list_interactions_2009, file="./data/final_data/regression_outputs/eb_list_interactions_2009_5_new.RData")

# Personal importance asked in 3 waves
eb_complete <- eb_tot %>% dplyr::select(surveyyear, pers_imp, sex, isced, urban, age, weightvec, isced_cat)
eb_complete <- eb_complete %>% drop_na()

eb_complete$time <- eb_complete$surveyyear - min(eb_complete$surveyyear)

eb_1 <- subset(eb_complete, surveyyear == 2009 | surveyyear == 2011)

# First empty model 
predictors <- c("time")
deps <- c("pers_imp")

eb_list_empty_pers_imp <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula,family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_1, weights = weightvec)
  
  eb_list_empty_pers_imp[[i]] <- eb_fit
} 

summary(eb_list_empty_pers_imp[[1]]) 

save(eb_list_empty_pers_imp, file="./data/final_data/regression_outputs/eb_list_empty_pers_imp1.RData")

# Model with predictors
predictors <- c("time", "sex", "isced_cat", "age")
deps <- c("pers_imp")

eb_list_preds_pers_imp <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula,
                   family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_1, weights = weightvec)
  
  eb_list_preds_pers_imp[[i]] <- eb_fit
} 

summary(eb_list_preds_pers_imp[[1]]) 

save(eb_list_preds_pers_imp, file="./data/final_data/regression_outputs/eb_list_preds_pers_imp1_new.RData")

# And the model with interactions
eb_1$age_time <- eb_1$age * eb_1$time
eb_1$sex_time <- as.numeric(eb_1$sex) * eb_1$time
eb_1 <- fastDummies::dummy_cols(eb_1, select_columns = "isced_cat")
eb_1$isced_basic_time <- eb_1$isced_cat_Basic * eb_1$time
eb_1$isced_intermediate_time <- eb_1$isced_cat_Intermediate * eb_1$time
eb_1$isced_advanced_time <- eb_1$isced_cat_Advanced * eb_1$time

interactions <- c("age_time", "sex_time", "isced_intermediate_time", "isced_advanced_time")

eb_list_interactions_pers_imp <- list()

 for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = "+ "))), sigma.formula = sigma.formula," + ", family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_1, weights = weightvec)
  
  eb_list_interactions_pers_imp[[i]] <- eb_fit
} 

summary(eb_list_interactions_pers_imp[[1]]) 

save(eb_list_interactions_pers_imp, file="./data/final_data/regression_outputs/eb_list_interactions_pers_imp1_new.RData")

# Second set
eb_2 <- subset(eb_complete, surveyyear == 2011 | surveyyear == 2014)

# First empty model 
predictors <- c("time")
deps <- c("pers_imp")

eb_list_empty_pers_imp <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula,family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_2, weights = weightvec)
  
  eb_list_empty_pers_imp[[i]] <- eb_fit
} 

summary(eb_list_empty_pers_imp[[1]]) 

save(eb_list_empty_pers_imp, file="./data/final_data/regression_outputs/eb_list_empty_pers_imp2.RData")

# Model with predictors
predictors <- c("time", "sex", "isced_cat", "age")
deps <- c("pers_imp")

eb_list_preds_pers_imp <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula,
                   family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_2, weights = weightvec)
  
  eb_list_preds_pers_imp[[i]] <- eb_fit
} 

summary(eb_list_preds_pers_imp[[1]]) 

save(eb_list_preds_pers_imp, file="./data/final_data/regression_outputs/eb_list_preds_pers_imp2_new.RData")

# And the model with interactions
eb_2$age_time <- eb_2$age * eb_2$time
eb_2$sex_time <- as.numeric(eb_2$sex) * eb_2$time
eb_2 <- fastDummies::dummy_cols(eb_2, select_columns = "isced_cat")
eb_2$isced_basic_time <- eb_2$isced_cat_Basic * eb_2$time
eb_2$isced_intermediate_time <- eb_2$isced_cat_Intermediate * eb_2$time
eb_2$isced_advanced_time <- eb_2$isced_cat_Advanced * eb_2$time

interactions <- c("age_time", "sex_time", "isced_intermediate_time", "isced_advanced_time")

eb_list_interactions_pers_imp <- list()

 for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = "+ "))), sigma.formula = sigma.formula," + ", family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_2, weights = weightvec)
  
  eb_list_interactions_pers_imp[[i]] <- eb_fit
} 

summary(eb_list_interactions_pers_imp[[1]]) 

save(eb_list_interactions_pers_imp, file="./data/final_data/regression_outputs/eb_list_interactions_pers_imp2_new.RData")

# Questions about who does enough asked in 3 waves
eb_complete <- eb_tot %>% dplyr::select(surveyyear, doprot_natgov:doprot_citiz, sex, isced, urban, age, weightvec, isced_cat)
eb_complete <- eb_complete %>% drop_na()
eb_complete$time <- eb_complete$surveyyear - min(eb_complete$surveyyear)
eb_1 <- subset(eb_complete, surveyyear == 2009 | surveyyear == 2014)

# First empty model 
predictors <- c("time")
deps <- c("doprot_natgov", "doprot_eu", "doprot_region", "doprot_comp", "doprot_citiz")

eb_list_empty_doprot <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_1, weights = weightvec)
  
  eb_list_empty_doprot[[i]] <- eb_fit
} 

summary(eb_list_empty_doprot[[1]]) 

save(eb_list_empty_doprot, file="./data/final_data/regression_outputs/eb_list_empty_doprot1.RData")

# Model with predictors
predictors <- c("time", "sex", "isced_cat", "age")
deps <- c("doprot_natgov", "doprot_eu", "doprot_region", "doprot_comp", "doprot_citiz")

eb_list_preds_doprot <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_1, weights = weightvec)
  
  eb_list_preds_doprot[[i]] <- eb_fit
} 

summary(eb_list_preds_doprot[[1]]) 

save(eb_list_preds_doprot, file="./data/final_data/regression_outputs/eb_list_preds_doprot1_new.RData")

# And the model with interactions
eb_1$age_time <- eb_1$age * eb_1$time
eb_1$sex_time <- as.numeric(eb_1$sex) * eb_1$time
eb_1 <- fastDummies::dummy_cols(eb_1, select_columns = "isced_cat")
eb_1$isced_basic_time <- eb_1$isced_cat_Basic * eb_1$time
eb_1$isced_intermediate_time <- eb_1$isced_cat_Intermediate * eb_1$time
eb_1$isced_advanced_time <- eb_1$isced_cat_Advanced * eb_1$time

interactions <- c("age_time", "sex_time", "isced_intermediate_time", "isced_advanced_time")

eb_list_interactions_doprot <- list()

 for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = "+ "))), sigma.formula = sigma.formula," + ", family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_1, weights = weightvec)
  
  eb_list_interactions_doprot[[i]] <- eb_fit
} 

summary(eb_list_interactions_doprot[[1]]) 

save(eb_list_interactions_doprot, file="./data/final_data/regression_outputs/eb_list_interactions_doprot1_new.RData")

# Second set
eb_2 <- subset(eb_complete, surveyyear == 2014 | surveyyear == 2017)

# First empty model 
predictors <- c("time")
deps <- c("doprot_natgov", "doprot_eu", "doprot_region", "doprot_comp", "doprot_citiz")

eb_list_empty_doprot <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_2, weights = weightvec)
  
  eb_list_empty_doprot[[i]] <- eb_fit
} 

summary(eb_list_empty_doprot[[1]]) 

save(eb_list_empty_doprot, file="./data/final_data/regression_outputs/eb_list_empty_doprot2.RData")

# Model with predictors
predictors <- c("time", "sex", "isced_cat", "age")
deps <- c("doprot_natgov", "doprot_eu", "doprot_region", "doprot_comp", "doprot_citiz")

eb_list_preds_doprot <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_2, weights = weightvec)
  
  eb_list_preds_doprot[[i]] <- eb_fit
} 

summary(eb_list_preds_doprot[[1]]) 

save(eb_list_preds_doprot, file="./data/final_data/regression_outputs/eb_list_preds_doprot2_new.RData")

# And the model with interactions
eb_2$age_time <- eb_2$age * eb_2$time
eb_2$sex_time <- as.numeric(eb_2$sex) * eb_2$time
eb_2 <- fastDummies::dummy_cols(eb_2, select_columns = "isced_cat")
eb_2$isced_basic_time <- eb_2$isced_cat_Basic * eb_2$time
eb_2$isced_intermediate_time <- eb_2$isced_cat_Intermediate * eb_2$time
eb_2$isced_advanced_time <- eb_2$isced_cat_Advanced * eb_2$time

interactions <- c("age_time", "sex_time", "isced_intermediate_time", "isced_advanced_time")

eb_list_interactions_doprot <- list()

 for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = "+ "))), sigma.formula = sigma.formula," + ", family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_2, weights = weightvec)
  
  eb_list_interactions_doprot[[i]] <- eb_fit
} 

summary(eb_list_interactions_doprot[[1]]) 

save(eb_list_interactions_doprot, file="./data/final_data/regression_outputs/eb_list_interactions_doprot2_new.RData")

# Questions about naming climate change a problem
eb_complete <- eb_tot %>% dplyr::select(surveyyear, cchange2, cchangetot, sex, isced, urban, age, weightvec, isced_cat)
eb_complete <- eb_complete %>% drop_na()
eb_complete$time <- eb_complete$surveyyear - min(eb_complete$surveyyear)

eb_1 <- subset(eb_complete, surveyyear == 2011 | surveyyear == 2013)

# First empty model 
predictors <- c("time")
deps <- c("cchange2", "cchangetot")

eb_list_empty_cchange2 <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_1, weights = weightvec)
  
  eb_list_empty_cchange2[[i]] <- eb_fit
} 

summary(eb_list_empty_cchange2[[1]]) 

save(eb_list_empty_cchange2, file="./data/final_data/regression_outputs/eb_list_empty_cchange2_1.RData")

# Model with predictors
predictors <- c("time", "sex", "isced_cat", "age")
deps <- c("cchange2", "cchangetot")

eb_list_preds_cchange2 <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula,
                   family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_1, weights = weightvec)
  
  eb_list_preds_cchange2[[i]] <- eb_fit
} 

summary(eb_list_preds_cchange2[[1]]) 

save(eb_list_preds_cchange2, file="./data/final_data/regression_outputs/eb_list_preds_cchange2_1_new.RData")

# And the model with interactions
eb_1$age_time <- eb_1$age * eb_1$time
eb_1$sex_time <- as.numeric(eb_1$sex) * eb_1$time
eb_1 <- fastDummies::dummy_cols(eb_1, select_columns = "isced_cat")
eb_1$isced_basic_time <- eb_1$isced_cat_Basic * eb_1$time
eb_1$isced_intermediate_time <- eb_1$isced_cat_Intermediate * eb_1$time
eb_1$isced_advanced_time <- eb_1$isced_cat_Advanced * eb_1$time

interactions <- c("age_time", "sex_time", "isced_intermediate_time", "isced_advanced_time")

eb_list_interactions_cchange2 <- list()

 for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = "+ "))), sigma.formula = sigma.formula," + ", family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_1, weights = weightvec)
  
  eb_list_interactions_cchange2[[i]] <- eb_fit
} 

summary(eb_list_interactions_cchange2[[1]]) 

save(eb_list_interactions_cchange2, file="./data/final_data/regression_outputs/eb_list_interactions_cchange2_1_new.RData")

# Second set
eb_2 <- subset(eb_complete, surveyyear == 2013 | surveyyear == 2015)

# First empty model 
predictors <- c("time")
deps <- c("cchange2", "cchangetot")

eb_list_empty_cchange2 <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_2, weights = weightvec)
  
  eb_list_empty_cchange2[[i]] <- eb_fit
} 

summary(eb_list_empty_cchange2[[1]]) 

save(eb_list_empty_cchange2, file="./data/final_data/regression_outputs/eb_list_empty_cchange2_2.RData")

# Model with predictors
predictors <- c("time", "sex", "isced_cat", "age")
deps <- c("cchange2", "cchangetot")

eb_list_preds_cchange2 <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula,
                   family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_2, weights = weightvec)
  
  eb_list_preds_cchange2[[i]] <- eb_fit
} 

summary(eb_list_preds_cchange2[[1]]) 

save(eb_list_preds_cchange2, file="./data/final_data/regression_outputs/eb_list_preds_cchange2_2_new.RData")

# And the model with interactions
eb_2$age_time <- eb_2$age * eb_2$time
eb_2$sex_time <- as.numeric(eb_2$sex) * eb_2$time
eb_2 <- fastDummies::dummy_cols(eb_2, select_columns = "isced_cat")
eb_2$isced_basic_time <- eb_2$isced_cat_Basic * eb_2$time
eb_2$isced_intermediate_time <- eb_2$isced_cat_Intermediate * eb_2$time
eb_2$isced_advanced_time <- eb_2$isced_cat_Advanced * eb_2$time

interactions <- c("age_time", "sex_time", "isced_intermediate_time", "isced_advanced_time")

eb_list_interactions_cchange2 <- list()

 for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = "+ "))), sigma.formula = sigma.formula," + ", family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_2, weights = weightvec)
  
  eb_list_interactions_cchange2[[i]] <- eb_fit
} 

summary(eb_list_interactions_cchange2[[1]]) 

save(eb_list_interactions_cchange2, file="./data/final_data/regression_outputs/eb_list_interactions_cchange2_2_new.RData")

# Third set
eb_3 <- subset(eb_complete, surveyyear == 2015 | surveyyear == 2017)

# First empty model 
predictors <- c("time")
deps <- c("cchange2", "cchangetot")

eb_list_empty_cchange2 <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_3, weights = weightvec)
  
  eb_list_empty_cchange2[[i]] <- eb_fit
} 

summary(eb_list_empty_cchange2[[1]]) 

save(eb_list_empty_cchange2, file="./data/final_data/regression_outputs/eb_list_empty_cchange2_3.RData")

# Model with predictors
predictors <- c("time", "sex", "isced_cat", "age")
deps <- c("cchange2", "cchangetot")

eb_list_preds_cchange2 <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula,
                   family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_3, weights = weightvec)
  
  eb_list_preds_cchange2[[i]] <- eb_fit
} 

summary(eb_list_preds_cchange2[[1]]) 

save(eb_list_preds_cchange2, file="./data/final_data/regression_outputs/eb_list_preds_cchange2_3_new.RData")

# And the model with interactions
eb_3$age_time <- eb_3$age * eb_3$time
eb_3$sex_time <- as.numeric(eb_3$sex) * eb_3$time
eb_3 <- fastDummies::dummy_cols(eb_3, select_columns = "isced_cat")
eb_3$isced_basic_time <- eb_3$isced_cat_Basic * eb_3$time
eb_3$isced_intermediate_time <- eb_3$isced_cat_Intermediate * eb_3$time
eb_3$isced_advanced_time <- eb_3$isced_cat_Advanced * eb_3$time

interactions <- c("age_time", "sex_time", "isced_intermediate_time", "isced_advanced_time")

eb_list_interactions_cchange2 <- list()

 for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = "+ "))), sigma.formula = sigma.formula," + ", family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_3, weights = weightvec)
  
  eb_list_interactions_cchange2[[i]] <- eb_fit
} 

summary(eb_list_interactions_cchange2[[1]]) 

save(eb_list_interactions_cchange2, file="./data/final_data/regression_outputs/eb_list_interactions_cchange2_3_new.RData")

# Fourth set
eb_4 <- subset(eb_complete, surveyyear == 2017 | surveyyear == 2021)

# First empty model 
predictors <- c("time")
deps <- c("cchange2", "cchangetot")

eb_list_empty_cchange2 <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_4, weights = weightvec)
  
  eb_list_empty_cchange2[[i]] <- eb_fit
} 

summary(eb_list_empty_cchange2[[1]]) 

save(eb_list_empty_cchange2, file="./data/final_data/regression_outputs/eb_list_empty_cchange2_4.RData")

# Model with predictors
predictors <- c("time", "sex", "isced_cat", "age")
deps <- c("cchange2", "cchangetot")

eb_list_preds_cchange2 <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula,
                   family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_4, weights = weightvec)
  
  eb_list_preds_cchange2[[i]] <- eb_fit
} 

summary(eb_list_preds_cchange2[[1]]) 

save(eb_list_preds_cchange2, file="./data/final_data/regression_outputs/eb_list_preds_cchange2_4_new.RData")

# And the model with interactions
eb_4$age_time <- eb_4$age * eb_4$time
eb_4$sex_time <- as.numeric(eb_4$sex) * eb_4$time
eb_4 <- fastDummies::dummy_cols(eb_4, select_columns = "isced_cat")
eb_4$isced_basic_time <- eb_4$isced_cat_Basic * eb_4$time
eb_4$isced_intermediate_time <- eb_4$isced_cat_Intermediate * eb_4$time
eb_4$isced_advanced_time <- eb_4$isced_cat_Advanced * eb_4$time

interactions <- c("age_time", "sex_time", "isced_intermediate_time", "isced_advanced_time")

eb_list_interactions_cchange2 <- list()

 for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = "+ "))), sigma.formula = sigma.formula," + ", family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_4, weights = weightvec)
  
  eb_list_interactions_cchange2[[i]] <- eb_fit
} 

summary(eb_list_interactions_cchange2[[1]]) 

save(eb_list_interactions_cchange2, file="./data/final_data/regression_outputs/eb_list_interactions_cchange2_4_new.RData")

# Question about personal action taken
eb_complete <- eb_tot %>% dplyr::select(surveyyear, prsaction, sex, isced, urban, age, weightvec, isced_cat)
eb_complete <- eb_complete %>% drop_na()
eb_complete$time <- eb_complete$surveyyear - min(eb_complete$surveyyear)
eb_1 <- subset(eb_complete, surveyyear == 2011 | surveyyear == 2013)

# First empty model 
predictors <- c("time")
deps <- c("prsaction")

eb_list_empty_prsaction <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_1, weights = weightvec)
  
  eb_list_empty_prsaction[[i]] <- eb_fit
} 

summary(eb_list_empty_prsaction[[1]]) 

save(eb_list_empty_prsaction, file="./data/final_data/regression_outputs/eb_list_empty_prsaction1.RData")

# Model with predictors
predictors <- c("time", "sex", "isced_cat", "age")
deps <- c("prsaction")

eb_list_preds_prsaction <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_1, weights = weightvec)
  
  eb_list_preds_prsaction[[i]] <- eb_fit
} 

summary(eb_list_preds_prsaction[[1]]) 

save(eb_list_preds_prsaction, file="./data/final_data/regression_outputs/eb_list_preds_prsaction1_new.RData")

# And the model with interactions
eb_1$age_time <- eb_1$age * eb_1$time
eb_1$sex_time <- as.numeric(eb_1$sex) * eb_1$time
eb_1 <- fastDummies::dummy_cols(eb_1, select_columns = "isced_cat")
eb_1$isced_basic_time <- eb_1$isced_cat_Basic * eb_1$time
eb_1$isced_intermediate_time <- eb_1$isced_cat_Intermediate * eb_1$time
eb_1$isced_advanced_time <- eb_1$isced_cat_Advanced * eb_1$time

interactions <- c("age_time", "sex_time", "isced_intermediate_time", "isced_advanced_time")

eb_list_interactions_prsaction <- list()

 for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = "+ "))), sigma.formula = sigma.formula," + ", family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_1, weights = weightvec)
  
  eb_list_interactions_prsaction[[i]] <- eb_fit
} 

summary(eb_list_interactions_prsaction[[1]]) 

save(eb_list_interactions_prsaction, file="./data/final_data/regression_outputs/eb_list_interactions_prsaction1_new.RData")

# Second set
eb_2 <- subset(eb_complete, surveyyear == 2013 | surveyyear == 2015)

# First empty model 
predictors <- c("time")
deps <- c("prsaction")

eb_list_empty_prsaction <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_2, weights = weightvec)
  
  eb_list_empty_prsaction[[i]] <- eb_fit
} 

summary(eb_list_empty_prsaction[[1]]) 

save(eb_list_empty_prsaction, file="./data/final_data/regression_outputs/eb_list_empty_prsaction2.RData")

# Model with predictors
predictors <- c("time", "sex", "isced_cat", "age")
deps <- c("prsaction")

eb_list_preds_prsaction <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_2, weights = weightvec)
  
  eb_list_preds_prsaction[[i]] <- eb_fit
} 

summary(eb_list_preds_prsaction[[1]]) 

save(eb_list_preds_prsaction, file="./data/final_data/regression_outputs/eb_list_preds_prsaction2_new.RData")

# And the model with interactions
eb_2$age_time <- eb_2$age * eb_2$time
eb_2$sex_time <- as.numeric(eb_2$sex) * eb_2$time
eb_2 <- fastDummies::dummy_cols(eb_2, select_columns = "isced_cat")
eb_2$isced_basic_time <- eb_2$isced_cat_Basic * eb_2$time
eb_2$isced_intermediate_time <- eb_2$isced_cat_Intermediate * eb_2$time
eb_2$isced_advanced_time <- eb_2$isced_cat_Advanced * eb_2$time

interactions <- c("age_time", "sex_time", "isced_intermediate_time", "isced_advanced_time")

eb_list_interactions_prsaction <- list()

 for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = "+ "))), sigma.formula = sigma.formula," + ", family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_2, weights = weightvec)
  
  eb_list_interactions_prsaction[[i]] <- eb_fit
} 

summary(eb_list_interactions_prsaction[[1]]) 

save(eb_list_interactions_prsaction, file="./data/final_data/regression_outputs/eb_list_interactions_prsaction2_new.RData")

# Third set
eb_3 <- subset(eb_complete, surveyyear == 2015 | surveyyear == 2017)

# First empty model 
predictors <- c("time")
deps <- c("prsaction")

eb_list_empty_prsaction <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_3, weights = weightvec)
  
  eb_list_empty_prsaction[[i]] <- eb_fit
} 

summary(eb_list_empty_prsaction[[1]]) 

save(eb_list_empty_prsaction, file="./data/final_data/regression_outputs/eb_list_empty_prsaction3.RData")

# Model with predictors
predictors <- c("time", "sex", "isced_cat", "age")
deps <- c("prsaction")

eb_list_preds_prsaction <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_3, weights = weightvec)
  
  eb_list_preds_prsaction[[i]] <- eb_fit
} 

summary(eb_list_preds_prsaction[[1]]) 

save(eb_list_preds_prsaction, file="./data/final_data/regression_outputs/eb_list_preds_prsaction3_new.RData")

# And the model with interactions
eb_3$age_time <- eb_3$age * eb_3$time
eb_3$sex_time <- as.numeric(eb_3$sex) * eb_3$time
eb_3 <- fastDummies::dummy_cols(eb_3, select_columns = "isced_cat")
eb_3$isced_basic_time <- eb_3$isced_cat_Basic * eb_3$time
eb_3$isced_intermediate_time <- eb_3$isced_cat_Intermediate * eb_3$time
eb_3$isced_advanced_time <- eb_3$isced_cat_Advanced * eb_3$time

interactions <- c("age_time", "sex_time", "isced_intermediate_time", "isced_advanced_time")

eb_list_interactions_prsaction <- list()

 for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = " + ")))
  # Perform the linear regression
  eb_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = "+ "))), sigma.formula = sigma.formula," + ", family=NO(mu.link = "identity", sigma.link = "identity"), data = eb_3, weights = weightvec)
  
  eb_list_interactions_prsaction[[i]] <- eb_fit
} 

summary(eb_list_interactions_prsaction[[1]]) 

save(eb_list_interactions_prsaction, file="./data/final_data/regression_outputs/eb_list_interactions_prsaction3_new.RData")

SOCON

# Repeat the same process for SOCON
load("./data/final_data/socontotal.RData")

socon_complete <- socontotal %>% drop_na()

socon_complete$time <- socon_complete$surveyyear - min(socon_complete$surveyyear)

socon_1 <- subset(socon_complete, surveyyear == 2020 | surveyyear == 2021)

predictors <- c("time")
deps <- c("fut_gen_socon")

socon_list_empty <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regresssoconn
  socon_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = socon_1, weights = weightvec)
  
  socon_list_empty[[i]] <- socon_fit
} 

save(socon_list_empty, file="./data/final_data/regression_outputs/per_wave/socon_list_empty1.RData")

# Model with predictors
predictors <- c("time", "sex", "isced_cat", "age")

socon_list_preds <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regresssoconn
  socon_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula,
                   family=NO(mu.link = "identity", sigma.link = "identity"), data = socon_1, weights = weightvec)
  
  socon_list_preds[[i]] <- socon_fit
} 

summary(socon_list_preds[[1]]) 

save(socon_list_preds, file="./data/final_data/regression_outputs/per_wave/socon_list_preds1_new.RData")

# Model with interactions
socon_1$age_time <- socon_1$age * socon_1$time
socon_1$sex_time <- as.numeric(socon_1$sex) * socon_1$time
socon_1 <- fastDummies::dummy_cols(socon_1, select_columns = "isced_cat")
socon_1$isced_basic_time <- socon_1$isced_cat_Basic * socon_1$time
socon_1$isced_intermediate_time <- socon_1$isced_cat_Intermediate * socon_1$time
socon_1$isced_advanced_time <- socon_1$isced_cat_Advanced * socon_1$time

interactions <- c("age_time", "sex_time", "isced_intermediate_time", "isced_advanced_time")

socon_list_interactions <- list()

 for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = " + ")))
  # Perform the linear regresssoconn
  socon_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = "+ "))), sigma.formula = sigma.formula," + ", family=NO(mu.link = "identity", sigma.link = "identity"), data = socon_1, weights = weightvec)
  
  socon_list_interactions[[i]] <- socon_fit
} 

summary(socon_list_interactions[[1]]) 

save(socon_list_interactions, file="./data/final_data/regression_outputs/per_wave/socon_list_interactions1_new.RData")

# Repeat the same process for SOCON
socon_2 <- subset(socon_complete, surveyyear == 2021 | surveyyear == 2022)

predictors <- c("time")
deps <- c("fut_gen_socon")

socon_list_empty <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regresssoconn
  socon_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = socon_2, weights = weightvec)
  
  socon_list_empty[[i]] <- socon_fit
} 

save(socon_list_empty, file="./data/final_data/regression_outputs/per_wave/socon_list_empty2.RData")

# Model with predictors
predictors <- c("time", "sex", "isced_cat", "age")

socon_list_preds <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regresssoconn
  socon_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula,
                   family=NO(mu.link = "identity", sigma.link = "identity"), data = socon_2, weights = weightvec)
  
  socon_list_preds[[i]] <- socon_fit
} 

summary(socon_list_preds[[1]]) 

save(socon_list_preds, file="./data/final_data/regression_outputs/per_wave/socon_list_preds2_new.RData")

# Model with interactions
socon_2$age_time <- socon_2$age * socon_2$time
socon_2$sex_time <- as.numeric(socon_2$sex) * socon_2$time
socon_2 <- fastDummies::dummy_cols(socon_2, select_columns = "isced_cat")
socon_2$isced_basic_time <- socon_2$isced_cat_Basic * socon_2$time
socon_2$isced_intermediate_time <- socon_2$isced_cat_Intermediate * socon_2$time
socon_2$isced_advanced_time <- socon_2$isced_cat_Advanced * socon_2$time

interactions <- c("age_time", "sex_time", "isced_intermediate_time", "isced_advanced_time")

socon_list_interactions <- list()

 for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = " + ")))
  # Perform the linear regresssoconn
  socon_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = "+ "))), sigma.formula = sigma.formula," + ", family=NO(mu.link = "identity", sigma.link = "identity"), data = socon_2, weights = weightvec)
  
  socon_list_interactions[[i]] <- socon_fit
} 

summary(socon_list_interactions[[1]]) 

save(socon_list_interactions, file="./data/final_data/regression_outputs/per_wave/socon_list_interactions2_new.RData")

LISS

# Repeat the same process for LISS
rm(list=ls())
load("./data/final_data/lisstotal.RData")

liss_complete <- lisstotal %>% drop_na()

liss_complete$time <- liss_complete$surveyyear - min(liss_complete$surveyyear)
table(liss_complete$surveyyear)

liss_1 <- subset(liss_complete, surveyyear == 2019 | surveyyear == 2020.58)

predictors <- c("time")
deps <- c("lifestyle")

liss_list_empty <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regresslissn
  liss_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = liss_1, weights = weightvec)
  
  liss_list_empty[[i]] <- liss_fit
} 

save(liss_list_empty, file="./data/final_data/regression_outputs/per_wave/liss_list_empty1.RData")

# Model with predictors
predictors <- c("time", "sex", "isced_cat", "age")

liss_list_preds <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regresslissn
  liss_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula,
                   family=NO(mu.link = "identity", sigma.link = "identity"), data = liss_1, weights = weightvec)
  
  liss_list_preds[[i]] <- liss_fit
} 

summary(liss_list_preds[[1]]) 

save(liss_list_preds, file="./data/final_data/regression_outputs/per_wave/liss_list_preds1_new.RData")

# Model with interactions
liss_1$age_time <- liss_1$age * liss_1$time
liss_1$sex_time <- as.numeric(liss_1$sex) * liss_1$time
liss_1 <- fastDummies::dummy_cols(liss_1, select_columns = "isced_cat")
liss_1$isced_basic_time <- liss_1$isced_cat_Basic * liss_1$time
liss_1$isced_intermediate_time <- liss_1$isced_cat_Intermediate * liss_1$time
liss_1$isced_advanced_time <- liss_1$isced_cat_Advanced * liss_1$time

interactions <- c("age_time", "sex_time", "isced_intermediate_time", "isced_advanced_time")

liss_list_interactions <- list()

 for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = " + ")))
  # Perform the linear regresslissn
  liss_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = "+ "))), sigma.formula = sigma.formula," + ", family=NO(mu.link = "identity", sigma.link = "identity"), data = liss_1, weights = weightvec)
  
  liss_list_interactions[[i]] <- liss_fit
} 

summary(liss_list_interactions[[1]]) 

save(liss_list_interactions, file="./data/final_data/regression_outputs/per_wave/liss_list_interactions1_new.RData")

# Repeat the same process for LISS
table(liss_complete$surveyyear)

liss_2 <- subset(liss_complete, surveyyear == 2020.58 | surveyyear == 2020.83)

predictors <- c("time")
deps <- c("lifestyle")

liss_list_empty <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regresslissn
  liss_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = liss_2, weights = weightvec)
  
  liss_list_empty[[i]] <- liss_fit
} 

save(liss_list_empty, file="./data/final_data/regression_outputs/per_wave/liss_list_empty2.RData")

# Model with predictors
predictors <- c("time", "sex", "isced_cat", "age")

liss_list_preds <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regresslissn
  liss_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula,
                   family=NO(mu.link = "identity", sigma.link = "identity"), data = liss_2, weights = weightvec)
  
  liss_list_preds[[i]] <- liss_fit
} 

summary(liss_list_preds[[1]]) 

save(liss_list_preds, file="./data/final_data/regression_outputs/per_wave/liss_list_preds2_new.RData")

# Model with interactions
liss_2$age_time <- liss_2$age * liss_2$time
liss_2$sex_time <- as.numeric(liss_2$sex) * liss_2$time
liss_2 <- fastDummies::dummy_cols(liss_2, select_columns = "isced_cat")
liss_2$isced_basic_time <- liss_2$isced_cat_Basic * liss_2$time
liss_2$isced_intermediate_time <- liss_2$isced_cat_Intermediate * liss_2$time
liss_2$isced_advanced_time <- liss_2$isced_cat_Advanced * liss_2$time

interactions <- c("age_time", "sex_time", "isced_intermediate_time", "isced_advanced_time")

liss_list_interactions <- list()

 for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = " + ")))
  # Perform the linear regresslissn
  liss_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = "+ "))), sigma.formula = sigma.formula," + ", family=NO(mu.link = "identity", sigma.link = "identity"), data = liss_2, weights = weightvec)
  
  liss_list_interactions[[i]] <- liss_fit
} 

summary(liss_list_interactions[[1]]) 

save(liss_list_interactions, file="./data/final_data/regression_outputs/per_wave/liss_list_interactions2_new.RData")

# Repeat the same process for LISS
table(liss_complete$surveyyear)

liss_3 <- subset(liss_complete, surveyyear == 2020.83 | surveyyear == 2021)

predictors <- c("time")
deps <- c("lifestyle")

liss_list_empty <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regresslissn
  liss_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula, family=NO(mu.link = "identity", sigma.link = "identity"), data = liss_3, weights = weightvec)
  
  liss_list_empty[[i]] <- liss_fit
} 

save(liss_list_empty, file="./data/final_data/regression_outputs/per_wave/liss_list_empty3.RData")

# Model with predictors
predictors <- c("time", "sex", "isced_cat", "age")

liss_list_preds <- list()

# Loop through each dependent variable
  for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + ")))
  # Perform the linear regresslissn
  liss_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "))), sigma.formula = sigma.formula,
                   family=NO(mu.link = "identity", sigma.link = "identity"), data = liss_3, weights = weightvec)
  
  liss_list_preds[[i]] <- liss_fit
} 

summary(liss_list_preds[[1]]) 

save(liss_list_preds, file="./data/final_data/regression_outputs/per_wave/liss_list_preds3_new.RData")

# Model with interactions
liss_3$age_time <- liss_3$age * liss_3$time
liss_3$sex_time <- as.numeric(liss_3$sex) * liss_3$time
liss_3 <- fastDummies::dummy_cols(liss_3, select_columns = "isced_cat")
liss_3$isced_basic_time <- liss_3$isced_cat_Basic * liss_3$time
liss_3$isced_intermediate_time <- liss_3$isced_cat_Intermediate * liss_3$time
liss_3$isced_advanced_time <- liss_3$isced_cat_Advanced * liss_3$time

interactions <- c("age_time", "sex_time", "isced_intermediate_time", "isced_advanced_time")

liss_list_interactions <- list()

 for (i in 1:length(deps)) {
  # Specify the sigma formula
  sigma.formula <- as.formula(paste("~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = " + ")))
  # Perform the linear regresslissn
  liss_fit <- gamlss(as.formula(paste(deps[i], "~", paste(predictors, collapse = " + "), " + ", paste(interactions, collapse = "+ "))), sigma.formula = sigma.formula," + ", family=NO(mu.link = "identity", sigma.link = "identity"), data = liss_3, weights = weightvec)
  
  liss_list_interactions[[i]] <- liss_fit
} 

summary(liss_list_interactions[[1]]) 

save(liss_list_interactions, file="./data/final_data/regression_outputs/per_wave/liss_list_interactions3_new.RData")