# List of packages
packages <- c("tidyverse", "modelsummary", "fst", "viridis", "knitr", "rmarkdown", "questionr", "ggridges")

# Install packages if they aren't installed already
new_packages <- packages[!(packages %in% installed.packages()[,"Package"])]
if(length(new_packages)) install.packages(new_packages)

# Load the packages
lapply(packages, library, character.only = TRUE)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.4.4     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.0
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
## Loading required package: viridisLite
## [[1]]
##  [1] "lubridate" "forcats"   "stringr"   "dplyr"     "purrr"     "readr"    
##  [7] "tidyr"     "tibble"    "ggplot2"   "tidyverse" "stats"     "graphics" 
## [13] "grDevices" "utils"     "datasets"  "methods"   "base"     
## 
## [[2]]
##  [1] "modelsummary" "lubridate"    "forcats"      "stringr"      "dplyr"       
##  [6] "purrr"        "readr"        "tidyr"        "tibble"       "ggplot2"     
## [11] "tidyverse"    "stats"        "graphics"     "grDevices"    "utils"       
## [16] "datasets"     "methods"      "base"        
## 
## [[3]]
##  [1] "fst"          "modelsummary" "lubridate"    "forcats"      "stringr"     
##  [6] "dplyr"        "purrr"        "readr"        "tidyr"        "tibble"      
## [11] "ggplot2"      "tidyverse"    "stats"        "graphics"     "grDevices"   
## [16] "utils"        "datasets"     "methods"      "base"        
## 
## [[4]]
##  [1] "viridis"      "viridisLite"  "fst"          "modelsummary" "lubridate"   
##  [6] "forcats"      "stringr"      "dplyr"        "purrr"        "readr"       
## [11] "tidyr"        "tibble"       "ggplot2"      "tidyverse"    "stats"       
## [16] "graphics"     "grDevices"    "utils"        "datasets"     "methods"     
## [21] "base"        
## 
## [[5]]
##  [1] "knitr"        "viridis"      "viridisLite"  "fst"          "modelsummary"
##  [6] "lubridate"    "forcats"      "stringr"      "dplyr"        "purrr"       
## [11] "readr"        "tidyr"        "tibble"       "ggplot2"      "tidyverse"   
## [16] "stats"        "graphics"     "grDevices"    "utils"        "datasets"    
## [21] "methods"      "base"        
## 
## [[6]]
##  [1] "rmarkdown"    "knitr"        "viridis"      "viridisLite"  "fst"         
##  [6] "modelsummary" "lubridate"    "forcats"      "stringr"      "dplyr"       
## [11] "purrr"        "readr"        "tidyr"        "tibble"       "ggplot2"     
## [16] "tidyverse"    "stats"        "graphics"     "grDevices"    "utils"       
## [21] "datasets"     "methods"      "base"        
## 
## [[7]]
##  [1] "questionr"    "rmarkdown"    "knitr"        "viridis"      "viridisLite" 
##  [6] "fst"          "modelsummary" "lubridate"    "forcats"      "stringr"     
## [11] "dplyr"        "purrr"        "readr"        "tidyr"        "tibble"      
## [16] "ggplot2"      "tidyverse"    "stats"        "graphics"     "grDevices"   
## [21] "utils"        "datasets"     "methods"      "base"        
## 
## [[8]]
##  [1] "ggridges"     "questionr"    "rmarkdown"    "knitr"        "viridis"     
##  [6] "viridisLite"  "fst"          "modelsummary" "lubridate"    "forcats"     
## [11] "stringr"      "dplyr"        "purrr"        "readr"        "tidyr"       
## [16] "tibble"       "ggplot2"      "tidyverse"    "stats"        "graphics"    
## [21] "grDevices"    "utils"        "datasets"     "methods"      "base"
spain_data <- read.fst("spain_data.fst")
# not a bad practice to create a duplicate, if you ever wanted to trace back steps
df <- spain_data
 #First, we're creating new variables 'behave', 'secure', 'safety', 'tradition', and 'rules'
# These variables are derived from existing variables in the dataset (ipbhprp, impsafe, ipstrgv, imptrad, ipfrule)
df <- df %>%
  mutate(behave = ipbhprp,
         secure = impsafe,
         safety = ipstrgv,
         tradition = imptrad,
         rules = ipfrule) %>%
  
  # Next, we're applying 'na_if' to convert values 7, 8, and 9 to NA for all the new variables
  mutate(across(c("behave", "secure", "safety", "tradition", "rules"),
                ~ na_if(.x, 7) %>% na_if(8) %>% na_if(9))) %>%
  
  # Now, we're applying reverse coding for all the new variables by subtracting each value from 7
  # Important: this is because higher values on the scale, according to N+I's conceptualization, are 
  # indicative of being 'libertarian' and lower as 'authoritarian'
  mutate(across(c("behave", "secure", "safety", "tradition", "rules"), ~ 7 - .x ))

# Now we can calculate the 'auth' variable based on the sum of the reverse-coded variables
# We're using scales::rescale to rescale the sum to a range of 0 to 100, and we're removing NA values
df$auth <- scales::rescale(df$behave + 
                      df$secure + 
                      df$safety + 
                      df$tradition + 
                      df$rules, to=c(0,100), na.rm=TRUE)

# Finally, we're filtering out rows where 'auth' is NA to retain only valid data
df <- df %>% filter(!is.na(auth))
table(df$auth)
## 
##    0    4    8   12   16   20   24   28   32   36   40   44   48   52   56   60 
##    7    3   11   10   32   45   63   92  144  189  268  397  462  573  807  985 
##   64   68   72   76   80   84   88   92   96  100 
## 1172 1532 1554 1670 2017 1396 1077  869  617  488
ggplot(df, aes(x = auth)) +
  geom_histogram(bins = 30, fill = "blue", color = "black") +
  theme_minimal() +
  labs(title = "Histogram of Authoritarian Scale",
       x = "Authoritarian Scale",
       y = "Count")

ggplot(df, aes(x = rules)) +
  geom_histogram(bins = 30, fill = "blue", color = "black") +
  theme_minimal() +
  labs(title = "Histogram of 'Rules' variable",
       x = "Scale",
       y = "Count")

df$year <- NA
replacements <- c(2002, 2004, 2006, 2008, 2010, 2012, 2014, 2016, 2018, 2020)
for(i in 1:10){
  df$year[df$essround == i] <- replacements[i]
}
auth_avg <- df %>%
  group_by(year) %>%
  summarize(auth_avg = mean(auth, na.rm = TRUE))

plot_zoom_auth <- ggplot(auth_avg, aes(x = year, y = auth_avg)) +
  geom_point(aes(size = auth_avg, color = auth_avg), alpha = 0.7) + 
  geom_line(aes(group = 1), color = "blue", linetype = "dashed") +  
  labs(title = "Authoritarian Scale Average by Year (Spain)",
       x = "Survey Year",
       y = "Authoritarian Scale Average") +
  theme_minimal() +
  theme(legend.position = "none")

print(plot_zoom_auth)

#Now note that we have zoomed in a lot to give us a sense of year-to-year fluctations. However, the scale can be misleading. Consider:
plot_full_scale_auth <- ggplot(auth_avg, aes(x = year, y = auth_avg)) +
  geom_point(aes(size = auth_avg, color = auth_avg), alpha = 0.7) + 
  geom_line(aes(group = 1), color = "blue", linetype = "dashed") +  
  labs(title = "Authoritarian Scale Average by Year (Spain)",
       x = "Survey Year",
       y = "Authoritarian Scale Average") +
  theme_minimal() +
  theme(legend.position = "none") +
  scale_y_continuous(limits = c(0, 100))

print(plot_full_scale_auth)

#Now to measure happiness

unique(spain_data$cntry)
## [1] "ES"
spain_happy <- spain_data %>%
  filter(cntry == "ES") %>%
  rename(happiness = happy)

spain_happy$y <- spain_happy$happiness


# Remove values 77, 88, 99 or else will alter results
spain_happy$y[spain_happy$y %in% 77:99] <- NA
# Calculating Mean and Median
mean_y <- mean(spain_happy$y, na.rm = TRUE)
cat("Mean of 'y' is:", mean_y, "\n")
## Mean of 'y' is: 7.496362
median_y <- median(spain_happy$y, na.rm = TRUE)
cat("Median of 'y' is:", median_y, "\n")
## Median of 'y' is: 8
# Finding the Mode
mode_y <- spain_happy %>%
  filter(!is.na(y)) %>%
  count(y) %>%
  arrange(desc(n)) %>%
  slice(1) %>%
  pull(y)

cat("\nMode of Y:", mode_y, "\n")
## 
## Mode of Y: 8
# Standard Deviation
sd_y <- sd(spain_happy$y, na.rm = TRUE)
cat("Standard Deviation of 'y':", sd_y, "\n")
## Standard Deviation of 'y': 1.793277
# Visualize the distribution of happiness
ggplot(spain_happy, aes(x = y)) +
  geom_histogram(binwidth = 1, fill = "skyblue", color = "black") +
  labs(
    title = "Distribution of Happiness in Spain",
    x = "Happiness Score",
    y = "Frequency"
  ) +
  theme_minimal()
## Warning: Removed 71 rows containing non-finite values (`stat_bin()`).

# Fit the regression model
model <- lm(lrscale ~ clsprty + yrbrn, data = spain_data)

# Get the coefficient estimates
coefficients <- coef(summary(model))

# Print the coefficients
print(coefficients)
##                 Estimate Std. Error   t value      Pr(>|t|)
## (Intercept) 2.6276182117  1.1587857  2.267562  2.336685e-02
## clsprty     6.5715146574  0.1849979 35.522113 7.684579e-268
## yrbrn       0.0008293065  0.0005622  1.475109  1.401993e-01
# Load the required libraries
library(ggplot2)
# Create the interaction effect plot using domicil instead of geo
interaction_plot <- ggplot(spain_data, aes(x = lrscale, y = happy, color = as.factor(domicil))) +
  geom_smooth(method = "lm", se = FALSE) +  # Add smoothed lines for each level of 'domicil'
  facet_wrap(~as.factor(domicil)) +  # Separate plots for each level of 'domicil'
  labs(x = "Average Authoritarianism", y = "Happiness", color = "Residence Type") +
  theme_minimal()

# Print the plot
print(interaction_plot)
## `geom_smooth()` using formula = 'y ~ x'