#### Claassen & Devine, The Politics of Imperial Nostalgia
#### Code to recode and scale attitudinal variables 

library(car)
library(dplyr)
library(tidyr)
library(forcats)
library(psych)
library(psychTools)
library(texreg)
library(arm)
library(lavaan)
library(semTable)
library(paran)
library(haven)
library(viridis)
library(RColorBrewer)
library(questionr)
library(list)
library(pewmethods)
library(glmnet)
library(vip)
library(party)
library(permimp)
library(pdp)
library(ranger)
library(foreach)
library(doParallel)
library(tth)
library(knitr)      
library(kableExtra) 

# folders
WD = rstudioapi::getActiveDocumentContext()$path 
setwd(dirname(WD))
print( getwd() )


## Read and merge data

# read data

datm = read.csv("rude_w1-3_merged_rep.csv")



### Party preferences----


## Wave 2

party_ques_w2 = c("Q2_Lab_w2", "Q2_Con_w2", "Q2_LibDem_w2", "Q2_Green_w2", "Q2_Reform_w2", 
                  "Q2_SNP_w2", "Q2_Plaid_w2")
sapply(datm[, party_ques_w2], table, useNA = "ifany")

# clean and recode NAs
party_w2 = c("sup_Lab_w2", "sup_Cons_w2", "sup_LibDem_w2", "sup_Green_w2", "sup_Reform_w2", 
             "sup_SNP_w2", "sup_Plaid_w2")
datm[, party_w2] = sapply(datm[, party_ques_w2], function(x) car::recode(x, "99=NA"))
datm[datm$Scotland == 0, "sup_SNP_w2"] = 0
datm[datm$Wales == 0, "sup_Plaid_w2"] = 0
sapply(datm[, party_w2], table, useNA = "ifany")

# find most-preferred party, breaking ties randomly
party_names = c("Lab", "Cons", "LibDem", "Green", "Reform", "SNP", "Plaid")
datm$pref_party_w2 = apply(datm[, party_w2], 1, function(row) {
  if (all(is.na(row))) {
    return(NA)
  } else {
    max_indices <- which(row == max(row, na.rm = TRUE))
    return(party_names[sample(max_indices, 1)])
  }
})
table(datm$pref_party_w2, useNA = "ifany")

# create left v right party preference indicator
datm$pref_party_di_w2 = car::recode(
  datm$pref_party_w2, "'Cons'='Right'; 'Reform'='Right'; NA=NA; else = 'Left'",
  as.factor = TRUE)


## Wave 3

party_ques_w3 = c("PartySupport_Lab_w3", "PartySupport_Con_w3", "PartySupport_LibDem_w3", 
                  "PartySupport_Green_w3", "PartySupport_Reform_w3", "PartySupport_SNP_w3", 
                  "PartySupport_Plaid_w3")
sapply(datm[, party_ques_w3], table, useNA = "ifany")

# clean and recode NAs
party_w3 = c("sup_Lab_w3", "sup_Cons_w3", "sup_LibDem_w3", "sup_Green_w3", 
             "sup_Reform_w3", "sup_SNP_w3", "sup_Plaid_w3")
datm[, party_w3] = sapply(datm[, party_ques_w3], function(x) car::recode(x, "99 = NA"))
datm[datm$Scotland == 0, "sup_SNP_w3"] = 0
datm[datm$Wales == 0, "sup_Plaid_w3"] = 0
sapply(datm[, party_w3], table, useNA = "ifany")

# most preferred party - ties broken randomly
party_names = c("Lab", "Cons", "LibDem", "Green", "Reform", "SNP", "Plaid")
datm$pref_party_w3 = apply(datm[, party_w3], 1, function(row) {
  if (all(is.na(row))) {
    return(NA)
  } else {
    max_indices <- which(row == max(row, na.rm = TRUE))
    return(party_names[sample(max_indices, 1)])
  }
})
table(datm$pref_party_w3, useNA = "ifany")

# create left v right party preference indicator
datm$pref_party_di_w3 = car::recode(
  datm$pref_party_w3, "'Cons'='Right'; 'Reform'='Right'; NA=NA; else = 'Left'",
  as.factor = TRUE)


### Vote intentions - w3

table(datm$voting_likelihood_w3, exclude = FALSE)
table(datm$voting_intention_w3, exclude = FALSE) # 26 NAs, 330 DKs, 43 others

print(sort(datm[datm$voting_intention_w3 == 9, "voting_intention_open_w3"]))

not_cit_resp = c("not a British citizen", "not vote", "not allowed to vote", 
                 "Nota", "I don’t qualify to vote as I don’t have a British passport")
datm$vote_intent_w3 = datm$voting_intention_w3
datm$vote_intent_w3 = 
  ifelse(datm$voting_likelihood_w3 == 1, 0,
         ifelse(datm$voting_intention_open_w3 %in% not_cit_resp, 0,
                ifelse(datm$voting_intention_w3 == 9 & 
                         datm$voting_intention_open_w3 == "__NA__", 99, 
                       datm$voting_intention_w3)))
table(datm$vote_intent_w3, exclude = FALSE) # 20 NAs, 317 DKs, 38 others
table(datm$voting_intention_w3, datm$vote_intent_w3, exclude = FALSE) # 26 NAs, 330 DKs, 43 others
datm$vote_intent_w3 = ifelse(datm$vote_intent_w3 == 99, NA, datm$vote_intent_w3)
datm$vote_intent_w3 = factor(
  datm$vote_intent_w3, levels = c(0, 1, 2, 3, 4, 5, 6, 7, 9), 
  labels = c("Not vote", "Cons", "Labour", "LibDem", "SNP", "Plaid", "Reform", "Green", "Other"))



### Imperial attitudes----


## wave 2

imp_att_ques_w2 = c("Q15_1_w2", "Q15_2_w2", "Q15_3_w2", "Q15_4_w2", "Q15_5_w2", "Q15_6_w2", 
                    "Q15_7_w2")
sapply(datm[, imp_att_ques_w2], table, useNA="ifany")

# recodes for plot
imp_att_qr_w2 = c("imp_nost_1_w2", "imp_nost_2_w2", "imp_nost_3_w2", "imp_nost_4_w2", 
                   "imp_nost_5_w2", "imp_nost_6_w2", "imp_nost_7_w2")
datm[, imp_att_qr_w2] = datm[, imp_att_ques_w2]
datm$imp_nost_6_w2 = ifelse(is.na(datm$Q15_6_w2), datm$Q1b_w2, datm$Q15_6_w2)
datm$imp_nost_7_w2 = ifelse(is.na(datm$Q15_7_w2), datm$Q1a_w2, datm$Q15_7_w2)
sapply(datm[, imp_att_qr_w2], table, useNA="ifany")

# items labels
imp_att_lab_w2 = c("The British Empire caused more\nharm than good to\ncolonised peoples",
                 "The British Empire had a great\ncivilising effect on the world",
                 "The British Empire advanced the\ninterests of humanity",
                 "The British Empire was\nresponsible for many atrocities",
                 "I wish Britain still had\nan empire",
                 "The British Empire was a golden\nage in our nation's history",
                 "The British Empire was a\nshameful period in our nation's\nhistory")

# colour palettes
lik5 = c(brewer.pal(5, "BrBG"), "#b4b1b1")
lik5_labs = c("strongly disagree", "somewhat disagree", "neither", "somewhat agree",
              "strongly agree", "don't know")
qual5 = c(brewer.pal(5, "Dark2"), "#b4b1b1")

# weighted proportions
wtd_props = list()
for(i in 1:length(imp_att_qr_w2)) {
  wtd_props[[i]] = prop.table(xtabs(datm$wt_w2 ~ datm[[imp_att_qr_w2[i]]], na.rm = TRUE)) * 100
}
wtd_props_matrix = do.call(cbind, lapply(wtd_props, function(x) as.numeric(x)))

# save summary table
grouped_props = as.data.frame(t(rbind(
  wtd_props_matrix[1,] + wtd_props_matrix[2,], # disagree
  wtd_props_matrix[3,] + wtd_props_matrix[6,], # unsure
  wtd_props_matrix[4,] + wtd_props_matrix[5,]  # agree
)))
colnames(grouped_props) = c("disagree", "unsure", "agree")
rownames(grouped_props) = c("harm_vs_good", "civilising_effect", "advanced_humanity", 
                            "responsible_atrocities", "wish_empire", "golden_age", "shameful_period")
grouped_props$net_pc_nost = NA
grouped_props$net_pc_nost[c(1, 4, 7)] = grouped_props$disagree[c(1, 4, 7)] - 
  grouped_props$agree[c(1, 4, 7)]
grouped_props$net_pc_nost[c(2, 3, 5, 6)] = grouped_props$agree[c(2, 3, 5, 6)] - 
  grouped_props$disagree[c(2, 3, 5, 6)]

write.csv(round(grouped_props, 1), "imperial_att_3cat_w2.csv", row.names = TRUE)

# Respondent pro vs anti summary
datm_recode = datm
for (item in imp_att_qr_w2) {
  if (item %in% imp_att_qr_w2[c(1, 4, 7)]) {
    datm_recode[[item]] = ifelse(datm[[item]] %in% 1:2, 1, ifelse(datm[[item]] %in% 4:5, -1, 0))
  } else {
    datm_recode[[item]] = ifelse(datm[[item]] %in% 1:2, -1, ifelse(datm[[item]] %in% 4:5, 1, 0))
  }
}
datm_recode$sum_recode = rowSums(datm_recode[imp_att_qr_w2], na.rm = TRUE)
sum(datm$wt_w2[datm_recode$sum_recode > 0], na.rm = TRUE) / 
  sum(datm$wt_w2, na.rm = TRUE) * 100 # pro-empire
sum(datm$wt_w2[datm_recode$sum_recode < 0], na.rm = TRUE) / 
  sum(datm$wt_w2, na.rm = TRUE) * 100 # anti-empire

# plot bar charts
pdf("imperial_att_dist_w2.pdf", height = 4, width = 6)
par(mfrow = c(2, 4), mar = c(4.5, 2, 0.5, 0.5), tcl = -0.2, cex = 0.75, las = 2, 
    mgp = c(1.8, 0.9, 0))
for(i in 1:length(imp_att_qr_w2)) {
  barplot(height = wtd_props[[i]], names.arg = lik5_labs, axes = FALSE, 
          col = lik5, cex.names = 0.6, mgp = c(1, 0.2, 0), ylim = c(0, 41))
  abline(h = seq(0, 40, by = 10), lty = 3, lwd = 0.75, col = rgb(0.5, 0.5, 0.5, 0.75))
  axis(side = 2, labels = TRUE, cex.axis = 0.6, mgp = c(1, 0.4, 0))
  barplot(height = wtd_props[[i]], col = lik5, axes = FALSE, names.arg = "", add = TRUE)
  text(x = 4, y = 36, imp_att_lab_w2[i], cex = 0.6, adj = 0.5)
}
dev.off()

# plot stacked barcharts
imp_att_lab_w2 = c("The British Empire caused\nmore harm than good\nto colonised peoples",
                 "The British Empire\nhad a great civilising\neffect on the world",
                 "The British Empire\nadvanced the interests\nof humanity",
                 "The British Empire\nwas responsible for\nmany atrocities",
                 "I wish Britain still\nhad an empire",
                 "The British Empire was\na golden age in\nour nation's history",
                 "The British Empire was\na shameful period in\nour nation's history")

lik5_labs = c("strongly\ndisagree", "somewhat\ndisagree", "neither", "somewhat\nagree",
              "strongly\nagree", "don't\nknow")

png("imperial_att_stack_w2.png", height = 1080, width = 1440)
par(mar = c(2, 10.5, 2, 1.5), tcl = -0.2, cex = 1.8, las = 0, mgp = c(3, 0.5, 0),
    xpd = TRUE)
barplot(height = wtd_props_matrix, names.arg = imp_att_lab_w2, col = lik5, 
        axes = FALSE, cex.names = 1, xlim = c(0, 100),
        beside = FALSE, las = 2, space = 0.3, horiz = TRUE)
segments(x0= seq(0, 100, by = 10), y0 = 0, y1 = 9, lty = 3, lwd = 1, 
         col = rgb(0, 0, 0, 0.75))
barplot(height = wtd_props_matrix, names.arg = NULL, col = lik5,
        axes = FALSE, beside = FALSE, add = TRUE, space = 0.3, horiz = TRUE)
axis(side = 1, at = seq(0, 100, by = 10), labels = TRUE, cex.axis = 1, mgp = c(3, 0.5, 0))
legend("top", legend = lik5_labs, fill = lik5, cex = 1, pt.cex = 1.2, 
       horiz = TRUE, bty= "n", inset = c(0, -0.05), x.intersp = 0.8)
dev.off()

# plot stacked barcharts for collapsed categories
lik3 = c(brewer.pal(5, "BrBG"))[c(1, 3, 5)]
lik3_labs = c("disagree", "unsure/dk", "agree")

png("imperial_att_collap_stack_w2.png", height = 1080, width = 1620)
layout(matrix(c(1, 2), ncol = 2), widths = c(4, 1.5))  # Adjust layout to have two panels

par(mar = c(2, 12, 1.5, 1.5), tcl = -0.2, cex = 1.9, las = 0, mgp = c(3, 0.5, 0),
    xpd = TRUE)
barplot(height = as.matrix(t(grouped_props[, 1:3])), names.arg = rep("", 7), 
        col = lik3, axes = FALSE, xlim = c(0, 100),
        beside = FALSE, las = 2, space = 0.3, horiz = TRUE)
segments(x0 = seq(0, 100, by = 10), y0 = 0, y1 = 9, lty = 3, 
         lwd = 1, col = rgb(0, 0, 0, 0.75))
barplot(height = as.matrix(t(grouped_props[, 1:3])), names.arg = imp_att_lab_w2, 
        col = lik3, axes = FALSE, cex.names = 1, las = 2, add = TRUE, xlim = c(0, 100),
        beside = FALSE, space = 0.3, horiz = TRUE)
axis(side = 1, at = seq(0, 100, by = 10), labels = TRUE, cex.axis = 1, mgp = c(3, 0.5, 0))
legend("top", legend = lik3_labs, fill = lik3, cex = 1.1, pt.cex = 1.2, 
       horiz = TRUE, bty = "n", inset = c(0, -0.02), x.intersp = 0.8)

par(mar = c(2, 2, 1.5, 1.5), xpd = TRUE)
barplot(height = grouped_props$net_pc_nost, names.arg = NULL, horiz = TRUE, 
        col = "grey40", xlim = range(c(-50, 50)), space = 0.3)
segments(x0 = seq(-40, 40, by = 20), y0 = 0, y1 = 9, lty = 3, 
         lwd = 1, col = rgb(0, 0, 0, 0.75))
barplot(height = grouped_props$net_pc_nost, names.arg = NULL, horiz = TRUE, 
        col = "grey40", xlim = range(c(-50, 50)), space = 0.3, add = TRUE)
legend("top", legend = "Net % nostalgic", fill = "grey40", cex = 1.1, pt.cex = 1.2, 
       horiz = TRUE, bty = "n", inset = c(0, -0.02), x.intersp = 0.8)

dev.off()

# plot stacked barcharts with net % nostalgic
png("imperial_att_stack_netpc_w2.png", height = 1080, width = 1620)
layout(matrix(c(1, 2), ncol = 2), widths = c(4, 1.5))  

par(mar = c(2, 10.5, 1.5, 0.5), tcl = -0.2, cex = 1.9, las = 0, mgp = c(3, 0.5, 0),
    xpd = TRUE)
barplot(height = wtd_props_matrix, names.arg = imp_att_lab_w2, col = lik5, 
        axes = FALSE, cex.names = 1, xlim = c(0, 100),
        beside = FALSE, las = 2, space = 0.3, horiz = TRUE)
segments(x0= seq(0, 100, by = 10), y0 = 0, y1 = 9, lty = 3, lwd = 1, 
         col = rgb(0, 0, 0, 0.75))
barplot(height = wtd_props_matrix, names.arg = NULL, col = lik5,
        axes = FALSE, beside = FALSE, add = TRUE, space = 0.3, horiz = TRUE)
axis(side = 1, at = seq(0, 100, by = 10), labels = TRUE, cex.axis = 1, mgp = c(3, 0.5, 0))
legend("top", legend = lik5_labs, fill = lik5, cex = 1, pt.cex = 1.2, 
       horiz = TRUE, bty= "n", inset = c(0, -0.05), x.intersp = 0.8)

par(mar = c(2, 2, 1.5, 1.5), xpd = TRUE)
barplot(height = grouped_props$net_pc_nost, names.arg = NULL, horiz = TRUE, 
        col = "grey40", xlim = range(c(-50, 50)), space = 0.3)
segments(x0 = seq(-40, 40, by = 20), y0 = 0, y1 = 9, lty = 3, 
         lwd = 1, col = rgb(0, 0, 0, 0.75))
barplot(height = grouped_props$net_pc_nost, names.arg = NULL, horiz = TRUE, 
        col = "grey40", xlim = range(c(-50, 50)), space = 0.3, add = TRUE)
legend("top", legend = "Net % nostalgic", fill = "grey40", cex = 1.1, pt.cex = 1.2, 
       horiz = TRUE, bty = "n", inset = c(0, -0.02), x.intersp = 0.8)

dev.off()

# plot stacked barcharts for conjoint comparison
png("imperial_att_conjcomp_stack_w2.png", height = 450, width = 1440)
par(mar = c(2, 9, 3, 0.5), tcl = -0.2, cex = 1.8, las = 0, mgp = c(3, 0.5, 0),
    xpd = TRUE)
barplot(height = wtd_props_matrix[,c(2,4)], names.arg = imp_att_lab_w2[c(2,4)], 
        col = lik5, axes = FALSE, cex.names = 1, xlim = c(0, 100),
        beside = FALSE, las = 2, space = 0.3, horiz = TRUE)
segments(x0= seq(0, 100, by = 10), y0 = 0, y1 = 2, lty = 3, lwd = 1, 
         col = rgb(0, 0, 0, 0.75))
barplot(height = wtd_props_matrix[,c(2,4)], names.arg = NULL, col = lik5,
        axes = FALSE, beside = FALSE, add = TRUE, space = 0.3, horiz = TRUE)
axis(side = 1, at = seq(0, 100, by = 10), labels = TRUE, cex.axis = 1, mgp = c(3, 0.5, 0))
legend("top", legend = lik5_labs, fill = lik5, cex = 1, pt.cex = 1.2, 
       horiz = TRUE, bty= "n", inset = c(0, -0.25), x.intersp = 0.8)
dev.off()

# recodes for scaling
datm[, imp_att_qr_w2[c(1,4,7)]] = sapply(datm[, imp_att_ques_w2[c(1,4,7)]], 
                             function(x) car::recode(x, "1=5; 2=4; 3=3; 4=2; 5=1; 99=3"))
datm[, imp_att_qr_w2[c(2,3,5,6)]] = sapply(datm[, imp_att_ques_w2[c(2,3,5,6)]], 
                                      function(x) car::recode(x, "99=3"))
datm$imp_nost_6_w2 = ifelse(is.na(datm$Q15_6_w2), datm$Q1b_w2, datm$Q15_6_w2)
datm$imp_nost_7_w2 = ifelse(is.na(datm$Q15_7_w2), datm$Q1a_w2, datm$Q15_7_w2)
datm$imp_nost_6_w2 = car::recode(datm$imp_nost_6_w2, "99=3")
datm$imp_nost_7_w2 = car::recode(datm$imp_nost_7_w2, "1=5; 2=4; 3=3; 4=2; 5=1; 99=3")
sapply(datm[, imp_att_qr_w2], table, useNA = "ifany")

# plot pro- vs anti-imperial
datm$imp_att_add_w2 = rowMeans(datm[, imp_att_qr_w2], na.rm = TRUE)
imp_att_5c_w2 = cut(datm$imp_att_add_w2, breaks = c(0:5+0.5))
imp_att_3c_w2 = cut(datm$imp_att_add_w2, breaks = c(0.5, 2.5, 3.5, 5.5))
wtd_props = prop.table(xtabs(datm$wt_w2 ~ imp_att_3c_w2, na.rm = TRUE)) * 100

lik3 = c(brewer.pal(5, "BrBG"))[c(1, 3, 5)]
png("imp_att_3cat_w2.png", width = 600, height = 400)
par(mfrow = c(1, 1), mar = c(2, 3, 0.5, 0.5), las = 1, tcl = -0.2, cex = 1.7,
    mgp = c(1, 0.3, 0))
barplot(wtd_props, ylim = c(0, 50), axes = FALSE, names.arg = "", 
        xlab = "", ylab = "", border = NA, col = "white")
abline(h = seq(0, 100, by = 10), lty = 3, col = grey(0.8))
barplot(wtd_props, ylim = c(0, 60), add = TRUE, axes = FALSE, 
        names.arg = c("Anti-empire", "Neutral", "Pro-empire"),
        col = lik3)
axis(side = 2, at = seq(0, 100, by = 10), mgp = c(1, 0.4, 0), 
     cex.axis = 0.8)
mtext("%", side = 2, line = 1.7, cex = 1.6)
dev.off()

# plot pro- vs anti-imperial by party pref, left vs right
wtd_props_party = prop.table(xtabs(
  datm$wt_w2 ~ imp_att_3c_w2 + datm$pref_party_di_w2, na.rm = TRUE), margin = 2) * 100
lik3 = c(brewer.pal(5, "BrBG"))[c(1, 3, 5)]

png("imp_att_3cat_bypartyy_w2.png", width = 600, height = 400)

layout(matrix(c(1, 2), ncol = 2), widths = c(1, 0.85))  
par(mar = c(2, 2.5, 1.5, 1), las = 1, tcl = -0.2, cex = 1.7,
    mgp = c(1, 0.3, 0))
barplot(wtd_props_party[, 1], ylim = c(0, 50), axes = FALSE, names.arg = "", 
        xlab = "", ylab = "", border = NA, col = "white")
abline(h = seq(0, 100, by = 10), lty = 3, col = grey(0.8))
barplot(wtd_props_party[, 1], ylim = c(0, 60), add = TRUE, axes = FALSE, 
        names.arg = c("Anti-empire", "Neutral", "Pro-empire"),
        col = lik3)
axis(side = 2, at = seq(0, 100, by = 10), mgp = c(1, 0.4, 0), 
     cex.axis = 0.8)
mtext("%", side = 2, line = 1.7, cex = 1.6)
mtext("Prefers left party", side = 3, line = 0.3, cex = 1.6)

par(mar = c(2, 1, 1.5, 0.5), las = 1, tcl = -0.2, cex = 1.7,
    mgp = c(1, 0.3, 0))
barplot(wtd_props_party[, 2], ylim = c(0, 50), axes = FALSE, names.arg = "", 
        xlab = "", ylab = "", border = NA, col = "white")
abline(h = seq(0, 100, by = 10), lty = 3, col = grey(0.8))
barplot(wtd_props_party[, 2], ylim = c(0, 60), add = TRUE, axes = FALSE, 
        names.arg = c("Anti-empire", "Neutral", "Pro-empire"),
        col = lik3)
axis(side = 2, at = seq(0, 100, by = 10), mgp = c(1, 0.4, 0), 
     cex.axis = 0.8)
mtext("Prefers right party", side = 3, line = 0.3, cex = 1.6)

dev.off()

# reliability of scale
alpha(datm[, imp_att_qr_w2]) # alpha = 0.91; ave r = 0.60

# scree plot
png("imp_att_w2_scree.png", width = 720, height = 720)
in_cor_w2 = cor(datm[, imp_att_qr_w2], use = "pair")
par(mfrow = c(1, 1), mar = c(3, 3, 0.5, 0.5), las = 0, tcl = -0.2, cex = 1.6)
plot(eigen(in_cor_w2)$values, type = "n", axes = FALSE, xlab = "", ylab = "",
     xlim = c(1, 7), ylim = c(0, 5)) 
abline(h = 0:5, lty = 3, col = grey(0.8))
abline(v = 1:7, lty = 3, col = grey(0.8))
abline(h = 1, lty = 2, col = grey(0.5))
points(eigen(in_cor_w2)$values, type = "b", lwd = 1.5, pch = 19)
axis(side = 1, at = 1:7, mgp = c(1, 0.35, 0), cex.axis = 1)
axis(side = 2, at = 0:5, mgp = c(1, 0.45, 0), cex.axis = 1)
mtext("Component number", side = 1, line = 1.5, cex = 1.5)
mtext("Eigen values of components", side = 2, line = 1.7, cex = 1.5)
box()
dev.off()

# CFA
cfa_imp_att_w2 = '
imp_att_w2 =~ imp_nost_1_w2 + imp_nost_2_w2 + imp_nost_3_w2 + imp_nost_4_w2 
            + imp_nost_5_w2 + imp_nost_6_w2 + imp_nost_7_w2'
imp_att_cfa_w2 = cfa(cfa_imp_att_w2, data = datm,  estimator = "WLSMV", ordered = TRUE, 
                      std.lv = TRUE, missing = "pairwise")
semTable(imp_att_cfa_w2, paramSets = c("fits", "loadings", "latentvariances"),
         fits=c("chisq", "cfi", "rmsea", "srmr"), columns=c("est", "se", "p"),
         type="html", file="imp_att_cfa_w2_table.html")

imp_att_cfa_w2_fit = fitMeasures(
  imp_att_cfa_w2, output = "matrix", 
  fit.measures = c("cfi", "cfi.robust", "rmsea", "rmsea.robust", "srmr"))
write.csv(imp_att_cfa_w2_fit, file = "imp_att_cfa_w2_fit.csv", row.names = TRUE)

datm$imp_att_w2 = lavPredict(imp_att_cfa_w2)


## wave 3

imp_att_ques_w3 = c("empire_attitudes_1_w3", "empire_attitudes_2_w3", "empire_attitudes_3_w3", 
                    "empire_attitudes_4_w3", "empire_attitudes_5_w3", "empire_attitudes_6_w3", 
                    "empire_attitudes_7_w3")
sapply(datm[, imp_att_ques_w3], table, useNA="ifany")

# recodes for plot
imp_att_qr_w3 = c("imp_nost_1_w3", "imp_nost_2_w3", "imp_nost_3_w3", "imp_nost_4_w3", 
                "imp_nost_5_w3", "imp_nost_6_w3", "imp_nost_7_w3")
datm[, imp_att_qr_w3] = datm[, imp_att_ques_w3]
sapply(datm[, imp_att_qr_w3], table, useNA="ifany")

# items labels
imp_att_lab_w3 = c("The British Empire caused more\nharm than good to\ncolonised peoples",
                 "The British Empire had a great\ncivilising effect on the world",
                 "The British Empire advanced the\ninterests of humanity",
                 "The British Empire was\nresponsible for many atrocities",
                 "I wish Britain still had\nan empire",
                 "The British Empire was a golden\nage in our nation's history",
                 "The British Empire was a\nshameful period in our nation's\nhistory")

# colour palettes
lik5 = c(brewer.pal(5, "BrBG"), "#b4b1b1")
lik5_labs = c("strongly disagree", "somewhat disagree", "neither", "somewhat agree",
              "strongly agree", "don't know")
qual5 = c(brewer.pal(5, "Dark2"), "#b4b1b1")

# weighted proportions
wtd_props = list()
for(i in 1:length(imp_att_qr_w3)) {
  wtd_props[[i]] = prop.table(xtabs(datm$wt_w3 ~ datm[[imp_att_qr_w3[i]]], na.rm = TRUE)) * 100
}
wtd_props_matrix = do.call(cbind, lapply(wtd_props, function(x) as.numeric(x)))

# save summary table
grouped_props = as.data.frame(t(rbind(
  wtd_props_matrix[1,] + wtd_props_matrix[2,], # disagree
  wtd_props_matrix[3,] + wtd_props_matrix[6,], # unsure
  wtd_props_matrix[4,] + wtd_props_matrix[5,]  # agree
)))
colnames(grouped_props) = c("disagree", "unsure", "agree")
rownames(grouped_props) = c("harm_vs_good", "civilising_effect", "advanced_humanity", 
                            "responsible_atrocities", "wish_empire", "golden_age", "shameful_period")
grouped_props$net_pc_nost = NA
grouped_props$net_pc_nost[c(1, 4, 7)] = grouped_props$disagree[c(1, 4, 7)] - 
  grouped_props$agree[c(1, 4, 7)]
grouped_props$net_pc_nost[c(2, 3, 5, 6)] = grouped_props$agree[c(2, 3, 5, 6)] - 
  grouped_props$disagree[c(2, 3, 5, 6)]
write.csv(round(grouped_props, 1), "imperial_att_3cat_w3.csv", row.names = TRUE)

# Respondent pro vs anti summary
datm_recode = datm
for (item in imp_att_qr_w3) {
  if (item %in% imp_att_qr_w3[c(1, 4, 7)]) {
    datm_recode[[item]] = ifelse(datm[[item]] %in% 1:2, 1, ifelse(datm[[item]] %in% 4:5, -1, 0))
  } else {
    datm_recode[[item]] = ifelse(datm[[item]] %in% 1:2, -1, ifelse(datm[[item]] %in% 4:5, 1, 0))
  }
}
datm_recode$sum_recode = rowSums(datm_recode[imp_att_qr_w3], na.rm = TRUE)
sum(datm$wt_w3[datm_recode$sum_recode > 0], na.rm = TRUE) / 
  sum(datm$wt_w3, na.rm = TRUE) * 100 # pro-empire
sum(datm$wt_w3[datm_recode$sum_recode < 0], na.rm = TRUE) / 
  sum(datm$wt_w3, na.rm = TRUE) * 100 # anti-empire

# plot
pdf("imperial_att_dist_w3.pdf", height = 4, width = 6)
par(mfrow = c(2, 4), mar = c(4.5, 2, 0.5, 0.5), tcl = -0.2, cex = 0.75, las = 2, 
    mgp = c(1.8, 0.9, 0))
for(i in 1:length(imp_att_qr_w3)) {
  barplot(height = wtd_props[[i]], names.arg = lik5_labs, axes = FALSE, 
          col = lik5, cex.names = 0.6, mgp = c(1, 0.2, 0), ylim = c(0, 41))
  abline(h = seq(0, 40, by = 10), lty = 3, lwd = 0.75, col = rgb(0.5, 0.5, 0.5, 0.75))
  axis(side = 2, labels = TRUE, cex.axis = 0.6, mgp = c(1, 0.4, 0))
  barplot(height = wtd_props[[i]], col = lik5, axes = FALSE, names.arg = "", add = TRUE)
  text(x = 4, y = 36, imp_att_lab_w3[i], cex = 0.6, adj = 0.5)
}
dev.off()

# plot stacked barcharts
imp_att_lab_w3 = c("The British Empire caused\nmore harm than good\nto colonised peoples",
                 "The British Empire\nhad a great civilising\neffect on the world",
                 "The British Empire\nadvanced the interests\nof humanity",
                 "The British Empire\nwas responsible for\nmany atrocities",
                 "I wish Britain still\nhad an empire",
                 "The British Empire was\na golden age in\nour nation's history",
                 "The British Empire was\na shameful period in\nour nation's history")

lik5_labs = c("strongly\ndisagree", "somewhat\ndisagree", "neither", "somewhat\nagree",
              "strongly\nagree", "don't\nknow")

png("imperial_att_stack_w3.png", height = 1080, width = 1440)
par(mar = c(2, 10.5, 2, 1.5), tcl = -0.2, cex = 1.8, las = 0, mgp = c(3, 0.5, 0))
barplot(height = wtd_props_matrix, names.arg = imp_att_lab_w3, col = lik5, 
        axes = FALSE, cex.names = 1, xlim = c(0, 100),
        beside = FALSE, las = 2, space = 0.3, horiz = TRUE)
segments(x0= seq(0, 100, by = 10), y0 = 0, y1 = 9, lty = 3, lwd = 1, 
         col = rgb(0, 0, 0, 0.75))
barplot(height = wtd_props_matrix, names.arg = NULL, col = lik5,
        axes = FALSE, beside = FALSE, add = TRUE, space = 0.3, horiz = TRUE)
axis(side = 1, at = seq(0, 100, by = 10), labels = TRUE, cex.axis = 1, mgp = c(3, 0.5, 0))
par(xpd = TRUE)
legend("top", legend = lik5_labs, fill = lik5, cex = 1, pt.cex = 1.2, 
       horiz = TRUE, bty= "n", inset = c(0, -0.05), x.intersp = 0.8)
dev.off()

# plot stacked barcharts for collapsed categories
lik3 = c(brewer.pal(5, "BrBG"))[c(1, 3, 5)]
lik3_labs = c("Disagree", "Unsure/DK", "Agree")

png("imperial_att_collap_stack_w3.png", height = 1080, width = 1620)
layout(matrix(c(1, 2), ncol = 2), widths = c(4, 1.5))  # Adjust layout to have two panels

par(mar = c(2, 12, 1.5, 1.5), tcl = -0.2, cex = 1.9, las = 0, mgp = c(3, 0.5, 0),
    xpd = TRUE)
barplot(height = as.matrix(t(grouped_props[, 1:3])), names.arg = rep("", 7), 
        col = lik3, axes = FALSE, xlim = c(0, 100),
        beside = FALSE, las = 2, space = 0.3, horiz = TRUE)
segments(x0 = seq(0, 100, by = 10), y0 = 0, y1 = 9, lty = 3, 
         lwd = 1, col = rgb(0, 0, 0, 0.75))
barplot(height = as.matrix(t(grouped_props[, 1:3])), names.arg = imp_att_lab_w2, 
        col = lik3, axes = FALSE, cex.names = 1, las = 2, add = TRUE, xlim = c(0, 100),
        beside = FALSE, space = 0.3, horiz = TRUE)
axis(side = 1, at = seq(0, 100, by = 10), labels = TRUE, cex.axis = 1, mgp = c(3, 0.5, 0))
legend("top", legend = lik3_labs, fill = lik3, cex = 1.1, pt.cex = 1.2, 
       horiz = TRUE, bty = "n", inset = c(0, -0.02), x.intersp = 0.8)

par(mar = c(2, 2, 1.5, 1.5), xpd = TRUE)
barplot(height = grouped_props$net_pc_nost, names.arg = NULL, horiz = TRUE, 
        col = "grey40", xlim = range(c(-50, 50)), space = 0.3)
segments(x0 = seq(-40, 40, by = 20), y0 = 0, y1 = 9, lty = 3, 
         lwd = 1, col = rgb(0, 0, 0, 0.75))
barplot(height = grouped_props$net_pc_nost, names.arg = NULL, horiz = TRUE, 
        col = "grey40", xlim = range(c(-50, 50)), space = 0.3, add = TRUE)
legend("top", legend = "Net % nostalgic", fill = "grey40", cex = 1.1, pt.cex = 1.2, 
       horiz = TRUE, bty = "n", inset = c(0, -0.02), x.intersp = 0.8)

dev.off()

# plot stacked barcharts with net % nostalgic
png("imperial_att_stack_netpc_w3.png", height = 1080, width = 1620)
layout(matrix(c(1, 2), ncol = 2), widths = c(4, 1.5))  

par(mar = c(2, 10.5, 1.5, 0.5), tcl = -0.2, cex = 1.9, las = 0, mgp = c(3, 0.5, 0),
    xpd = TRUE)
barplot(height = wtd_props_matrix, names.arg = imp_att_lab_w3, col = lik5, 
        axes = FALSE, cex.names = 1, xlim = c(0, 100),
        beside = FALSE, las = 2, space = 0.3, horiz = TRUE)
segments(x0= seq(0, 100, by = 10), y0 = 0, y1 = 9, lty = 3, lwd = 1, 
         col = rgb(0, 0, 0, 0.75))
barplot(height = wtd_props_matrix, names.arg = NULL, col = lik5,
        axes = FALSE, beside = FALSE, add = TRUE, space = 0.3, horiz = TRUE)
axis(side = 1, at = seq(0, 100, by = 10), labels = TRUE, cex.axis = 1, mgp = c(3, 0.5, 0))
legend("top", legend = lik5_labs, fill = lik5, cex = 1, pt.cex = 1.2, 
       horiz = TRUE, bty= "n", inset = c(0, -0.05), x.intersp = 0.8)

par(mar = c(2, 2, 1.5, 1.5), xpd = TRUE)
barplot(height = grouped_props$net_pc_nost, names.arg = NULL, horiz = TRUE, 
        col = "grey40", xlim = range(c(-50, 50)), space = 0.3)
segments(x0 = seq(-40, 40, by = 20), y0 = 0, y1 = 9, lty = 3, 
         lwd = 1, col = rgb(0, 0, 0, 0.75))
barplot(height = grouped_props$net_pc_nost, names.arg = NULL, horiz = TRUE, 
        col = "grey40", xlim = range(c(-50, 50)), space = 0.3, add = TRUE)
legend("top", legend = "Net % nostalgic", fill = "grey40", cex = 1.1, pt.cex = 1.2, 
       horiz = TRUE, bty = "n", inset = c(0, -0.02), x.intersp = 0.8)
dev.off()

# plot stacked barcharts for conjoint comparison
png("imperial_att_conjcomp_stack_w3.png", height = 450, width = 1440)
par(mar = c(2, 9, 3, 0.5), tcl = -0.2, cex = 1.8, las = 0, mgp = c(3, 0.5, 0),
    xpd = TRUE)
barplot(height = wtd_props_matrix[,c(2,4)], names.arg = imp_att_lab_w2[c(2,4)], 
        col = lik5, axes = FALSE, cex.names = 1, xlim = c(0, 100),
        beside = FALSE, las = 2, space = 0.3, horiz = TRUE)
segments(x0= seq(0, 100, by = 10), y0 = 0, y1 = 2, lty = 3, lwd = 1, 
         col = rgb(0, 0, 0, 0.75))
barplot(height = wtd_props_matrix[,c(2,4)], names.arg = NULL, col = lik5,
        axes = FALSE, beside = FALSE, add = TRUE, space = 0.3, horiz = TRUE)
axis(side = 1, at = seq(0, 100, by = 10), labels = TRUE, cex.axis = 1, mgp = c(3, 0.5, 0))
legend("top", legend = lik5_labs, fill = lik5, cex = 1, pt.cex = 1.2, 
       horiz = TRUE, bty= "n", inset = c(0, -0.25), x.intersp = 0.8)
dev.off()

# recodes for scales
datm[, imp_att_qr_w3[c(1,4,7)]] = sapply(datm[, imp_att_ques_w3[c(1,4,7)]], 
                                       function(x) car::recode(x, "1=5; 2=4; 3=3; 4=2; 5=1; 99=3"))
datm[, imp_att_qr_w3[c(2,3,5,6)]] = sapply(datm[, imp_att_ques_w3[c(2,3,5,6)]], 
                                         function(x) car::recode(x, "99=3"))
sapply(datm[, imp_att_qr_w3], table, useNA = "ifany")

# plot pro- vs anti-imperial
datm$imp_att_add_w3 = rowMeans(datm[, imp_att_qr_w3], na.rm = TRUE)
imp_att_5c_w3 = cut(datm$imp_att_add_w3, breaks = c(0:5+0.5))
imp_att_3c_w3 = cut(datm$imp_att_add_w3, breaks = c(0.5, 2.5, 3.5, 5.5))
wtd_props = prop.table(xtabs(datm$wt_w3 ~ imp_att_3c_w3, na.rm = TRUE)) * 100

lik3 = c(brewer.pal(5, "BrBG"))[c(1, 3, 5)]
png("imp_att_3cat_w3.png", width = 600, height = 400)
par(mfrow = c(1, 1), mar = c(2, 3, 0.5, 0.5), las = 1, tcl = -0.2, cex = 1.7,
    mgp = c(1, 0.3, 0))
barplot(wtd_props, ylim = c(0, 50), axes = FALSE, names.arg = "", 
        xlab = "", ylab = "", border = NA, col = "white")
abline(h = seq(0, 100, by = 10), lty = 3, col = grey(0.8))
barplot(wtd_props, ylim = c(0, 60), add = TRUE, axes = FALSE, 
        names.arg = c("Anti-empire", "Neutral", "Pro-empire"),
        col = lik3)
axis(side = 2, at = seq(0, 100, by = 10), mgp = c(1, 0.4, 0), 
     cex.axis = 0.8)
mtext("%", side = 2, line = 1.7, cex = 1.6)
dev.off()

# reliability scale
alpha(datm[, imp_att_qr_w3]) # alpha = 0.92; ave r = 0.63

# scree plot
png("imp_att_w3_scree.png", width = 720, height = 720)
in_cor_w3 = cor(datm[, imp_att_qr_w3], use = "pair")
par(mfrow = c(1, 1), mar = c(3, 3, 0.5, 0.5), las = 0, tcl = -0.2, cex = 1.6)
plot(eigen(in_cor_w3)$values, type = "n", axes = FALSE, xlab = "", ylab = "",
     xlim = c(1, 7), ylim = c(0, 5)) 
abline(h = 0:5, lty = 3, col = grey(0.8))
abline(v = 1:7, lty = 3, col = grey(0.8))
abline(h = 1, lty = 2, col = grey(0.5))
points(eigen(in_cor_w3)$values, type = "b", lwd = 1.5, pch = 19)
axis(side = 1, at = 1:7, mgp = c(1, 0.35, 0), cex.axis = 1)
axis(side = 2, at = 0:5, mgp = c(1, 0.45, 0), cex.axis = 1)
mtext("Component number", side = 1, line = 1.5, cex = 1.5)
mtext("Eigen values of components", side = 2, line = 1.7, cex = 1.5)
box()
dev.off()

# CFA 
cfa_imp_att_w3 = '
imp_att_w3 =~ imp_nost_1_w3 + imp_nost_2_w3 + imp_nost_3_w3 + imp_nost_4_w3 
            + imp_nost_5_w3 + imp_nost_6_w3 + imp_nost_7_w3'
imp_att_cfa_w3 = cfa(cfa_imp_att_w3, data = datm,  estimator = "WLSMV", ordered = TRUE, 
                      std.lv = TRUE, missing = "pairwise")
semTable(imp_att_cfa_w3, paramSets = c("fits", "loadings", "latentvariances"),
         fits=c("chisq", "cfi", "rmsea", "srmr"), columns=c("est", "se", "p"),
         type="html", file="imp_att_cfa_w3_table.html")
imp_att_cfa_w3_fit = fitMeasures(
  imp_att_cfa_w3, output = "matrix", 
  fit.measures = c("cfi", "cfi.robust", "rmsea", "rmsea.robust", "srmr"))
write.csv(imp_att_cfa_w3_fit, file = "imp_att_cfa_w3_fit.csv", row.names = TRUE)

datm$imp_att_w3 = lavPredict(imp_att_cfa_w3)

cor(datm$imp_att_w3, datm$imp_att_add_w3, use = "pair") # correlation = .994
cor(datm$imp_att_w2, datm$imp_att_w3, use = "pair") # test-retest correlation = .868



### Imperial emotions----


## Wave 2

imp_emot_ques_w2 = c("Q16_1_w2", "Q16_2_w2", "Q16_3_w2", "Q16_4_w2", "Q16_5_w2", 
                     "Q16_6_w2", "Q16_7_w2", "Q16_8_w2", "Q16_9_w2", "Q16_10_w2")
datm[, imp_emot_ques_w2] = sapply(datm[, imp_emot_ques_w2], function(x) car::recode(x, "99 = -1"))
sapply(datm[, imp_emot_ques_w2], table, useNA = "ifany")

# items
emot_items_w2 = c("embarrassment", "shame", "guilt", "anger", "pride", "nostalgia", "sadness", 
                  "disgust", "happiness", "satisfaction")

# colour palettes
col4 = c("#b4b1b1", brewer.pal(7, "BrBG")[4:7])
extent4_labs = c("DK", "Not at all", "a small extent", "a moderate extent", "a great extent")

# weighted proportions
wtd_props = list()
for(i in 1:length(imp_emot_ques_w2)) {
  wtd_props[[i]] = prop.table(xtabs(datm$wt_w2 ~ datm[[imp_emot_ques_w2[i]]], na.rm = TRUE)) * 100
}
wtd_props_matrix = do.call(cbind, lapply(wtd_props, function(x) as.numeric(x)))

# Respondent pro vs anti summary
datm_recode = datm
for (item in imp_emot_ques_w2) {
  if (item %in% imp_emot_ques_w2[c(5, 6, 9, 10)]) {
    datm_recode[[item]] = ifelse(datm[[item]] %in% 2:4, 1, 0)
  } else {
    datm_recode[[item]] = ifelse(datm[[item]] %in% 2:4, -1, 0)
  }
}
datm_recode$sum_recode = rowSums(datm_recode[imp_emot_ques_w2], na.rm = TRUE)
sum(datm$wt_w2[datm_recode$sum_recode > 0], na.rm = TRUE) / 
  sum(datm$wt_w2, na.rm = TRUE) * 100 # pro-empire
sum(datm$wt_w2[datm_recode$sum_recode < 0], na.rm = TRUE) / 
  sum(datm$wt_w2, na.rm = TRUE) * 100 # anti-empire

# re-order negative then positive  
emot_ord = c(1:4, 8, 7, 6, 5, 10, 9)

# plot
pdf("imperial_emot_dist_w2.pdf", height = 4, width = 6)
par(mfrow = c(2, 5), mar = c(4.5, 2, 0.5, 0.5), tcl = -0.2, cex = 0.75, las = 2, 
    mgp = c(1.8, 0.9, 0))
for(i in 1:length(imp_emot_ques_w2)) {
  barplot(height = wtd_props[[ emot_ord[i] ]], names.arg = extent4_labs, axes = FALSE, 
          col = col4, cex.names = 0.6, mgp = c(1, 0.2, 0), ylim = c(0, 50))
  abline(h = seq(0, 50, by = 10), lty = 3, lwd = 0.75, col = rgb(0.5, 0.5, 0.5, 0.75))
  axis(side = 2, labels = TRUE, cex.axis = 0.6, mgp = c(1, 0.4, 0))
  barplot(height = wtd_props[[ emot_ord[i] ]], col = col4, axes = FALSE, 
          names.arg = "", add = TRUE)
  text(x = 3, y = 48, emot_items_w2[emot_ord[i]], cex = 0.7, adj = 0.5)
}
dev.off()

# plot stacked barcharts
extent4_labs = c("Don't\nknow", "Not at\nall", "A small\nextent", 
                 "A moderate\nextent", "A great\nextent")

png("imperial_emot_stack_w2.png", height = 1080, width = 1440)
par(mar = c(2, 7.5, 2, 1.5), tcl = -0.2, cex = 1.8, las = 0, mgp = c(3, 0.5, 0))
barplot(height = wtd_props_matrix[, emot_ord], names.arg = emot_items_w2[emot_ord], 
        col = col4, axes = FALSE, cex.names = 1, xlim = c(0, 100),
        beside = FALSE, las = 2, space = 0.4, horiz = TRUE)
segments(x0= seq(0, 100, by = 10), y0 = 0, y1 = 13, lty = 3, lwd = 1, 
         col = rgb(0, 0, 0, 0.75))
barplot(height = wtd_props_matrix[, emot_ord], names.arg = NULL, col = col4,
        axes = FALSE, beside = FALSE, add = TRUE, space = 0.4, horiz = TRUE)
axis(side = 1, at = seq(0, 100, by = 10), labels = TRUE, cex.axis = 1, mgp = c(3, 0.5, 0))
par(xpd = TRUE)
legend("top", legend = extent4_labs, fill = col4, cex = 1, pt.cex = 1.2, 
       horiz = TRUE, bty= "n", inset = c(0, -0.05), x.intersp = 0.8)
dev.off()

# overall % nostalgic vs aversive
grouped_props = wtd_props_matrix[3:4, ] %>%
  colSums() %>%
  `*`(c(-1, -1, -1, -1, 1, 1, -1, -1, 1, 1))

# plot collapsed categories with valences
png("imperial_emot_collap_stack_w2.png", height = 1080, width = 1620)
layout(matrix(c(1, 2), ncol = 2), widths = c(5, 1.5)) 

par(mar = c(2, 7.5, 1.5, 1.5), tcl = -0.2, cex = 1.9, las = 0, mgp = c(3, 0.5, 0))
barplot(height = wtd_props_matrix[, emot_ord], names.arg = emot_items_w2[emot_ord], 
        col = col4, axes = FALSE, cex.names = 1, xlim = c(0, 100),
        beside = FALSE, las = 2, space = 0.4, horiz = TRUE)
segments(x0= seq(0, 100, by = 10), y0 = 0, y1 = 13, lty = 3, lwd = 1, 
         col = rgb(0, 0, 0, 0.75))
barplot(height = wtd_props_matrix[, emot_ord], names.arg = NULL, col = col4,
        axes = FALSE, beside = FALSE, add = TRUE, space = 0.4, horiz = TRUE)
axis(side = 1, at = seq(0, 100, by = 10), labels = TRUE, cex.axis = 1, mgp = c(3, 0.5, 0))
par(xpd = TRUE)
legend("top", legend = extent4_labs, fill = col4, cex = 1, pt.cex = 1.2, 
       horiz = TRUE, bty= "n", inset = c(0, -0.05), x.intersp = 0.8)

par(mar = c(2, 2, 1.5, 1.5), xpd = TRUE)
barplot(height = grouped_props[emot_ord], names.arg = NULL, horiz = TRUE, 
        col = "grey40", xlim = range(c(-30, 30)), space = 0.4)
segments(x0 = c(-20, 0, 20), y0 = 0, y1 = 14, lty = 3, 
         lwd = 1, col = rgb(0, 0, 0, 0.75))
barplot(height = grouped_props[emot_ord], names.arg = NULL, horiz = TRUE, 
        col = "grey40", xlim = range(c(-30, 30)), space = 0.4, add = TRUE)
legend("top", legend = "% nostalgic\n(moderate + great extent)", fill = "grey40", 
       cex = 1.1, pt.cex = 1.2, horiz = TRUE, bty = "n", 
       inset = c(0, -0.06), x.intersp = 0.8)

dev.off()

# plot by party preferences
wtd_props_party = list()
for(i in 1:length(imp_emot_ques_w2)) {
  wtd_props_party[[i]] = prop.table(xtabs(
    datm$wt_w2 ~ datm[[imp_emot_ques_w2[i]]] + datm$pref_party_di_w2, na.rm = TRUE),
    margin = 2) * 100
}

imp_emot_left = lapply(wtd_props_party, function(x) x[, 1])
imp_emot_left_matrix = do.call(cbind, lapply(imp_emot_left, function(x) as.numeric(x)))
imp_emot_right = lapply(wtd_props_party, function(x) x[, 2])
imp_emot_right_matrix = do.call(cbind, lapply(imp_emot_right, function(x) as.numeric(x)))

# plot collapsed categories with valences
png("imperial_emot_stack_byparty_w2.png", height = 1080, width = 1620)
layout(matrix(c(1, 3, 2, 3), ncol = 2), heights = c(7.5, 1), widths = c(1, 0.78))  

par(mar = c(1, 7.5, 1.2, 1), tcl = -0.2, cex = 1.9, las = 0, mgp = c(3, 0.5, 0))
barplot(height = imp_emot_left_matrix[,emot_ord], names.arg = emot_items_w2[emot_ord], 
        col = col4, axes = FALSE, cex.names = 1, xlim = c(0, 100),
        beside = FALSE, las = 2, space = 0.4, horiz = TRUE)
segments(x0= seq(0, 100, by = 25), y0 = 0, y1 = 13, lty = 3, lwd = 1, 
         col = rgb(0, 0, 0, 0.75))
barplot(height = imp_emot_left_matrix[,emot_ord], names.arg = NULL, col = col4,
        axes = FALSE, beside = FALSE, add = TRUE, space = 0.4, horiz = TRUE)
axis(side = 1, at = seq(0, 100, by = 25), labels = TRUE, cex.axis = 0.9, mgp = c(3, 0.3, 0))
mtext("Supports left party", side = 3, line = -0.2, cex = 2.2)
mtext("%", side = 1, line = 1.5, cex = 1.9)

par(mar = c(1, 1, 1.2, 1), tcl = -0.2, cex = 1.9, las = 0, mgp = c(3, 0.5, 0))
barplot(height = imp_emot_right_matrix[,emot_ord], names.arg = NULL, 
        col = col4, axes = FALSE, cex.names = 1, xlim = c(0, 100),
        beside = FALSE, las = 2, space = 0.4, horiz = TRUE)
segments(x0= seq(0, 100, by = 25), y0 = 0, y1 = 13, lty = 3, lwd = 1, 
         col = rgb(0, 0, 0, 0.75))
barplot(height = imp_emot_right_matrix[,emot_ord], names.arg = NULL, col = col4,
        axes = FALSE, beside = FALSE, add = TRUE, space = 0.4, horiz = TRUE)
axis(side = 1, at = seq(0, 100, by = 25), labels = TRUE, cex.axis = 0.9, mgp = c(3, 0.3, 0))
mtext("Supports right party", side = 3, line = -0.2, cex = 2.2)
mtext("%", side = 1, line = 1.5, cex = 1.9)

par(mar = c(0, 0, 0.5, 0), xpd = TRUE)  
plot.new()  
legend("center", legend = extent4_labs, fill = col4, cex = 1, horiz = TRUE, 
       bty = "n", x.intersp = 0.8)

dev.off()

# recodes for scales
imp_emot_qr_w2 = c("imp_emot_1_w2", "imp_emot_2_w2", "imp_emot_3_w2", "imp_emot_4_w2", "imp_emot_5_w2", 
                   "imp_emot_6_w2", "imp_emot_7_w2", "imp_emot_8_w2", "imp_emot_9_w2", "imp_emot_10_w2")
datm[, imp_emot_qr_w2] = sapply(datm[, imp_emot_ques_w2], function(x) car::recode(x, "-1 = NA"))
sapply(datm[, imp_emot_qr_w2], table, useNA = "ifany")

# plot pro- vs anti-imperial
datm$imp_emot_add_w2 = rowMeans(datm[, imp_emot_qr_w2[c(5, 6, 9, 10)]], na.rm = TRUE) -
  rowMeans(datm[, imp_emot_qr_w2[c(1, 2, 3, 4, 7, 8)]], na.rm = TRUE)
imp_emot_3c_w2 = cut(datm$imp_emot_add_w2, breaks = c(-3.5, -0.5, 0.5, 3.5))
wtd_props = prop.table(xtabs(datm$wt_w2 ~ imp_emot_3c_w2, na.rm = TRUE)) * 100

lik3 = c(brewer.pal(5, "BrBG"))[c(1, 3, 5)]
png("imp_emot_3cat_w2.png", width = 600, height = 400)
par(mfrow = c(1, 1), mar = c(2, 3, 0.5, 0.5), las = 1, tcl = -0.2, cex = 1.7,
    mgp = c(1, 0.3, 0))
barplot(wtd_props, ylim = c(0, 50), axes = FALSE, names.arg = "", 
        xlab = "", ylab = "", border = NA, col = "white")
abline(h = seq(0, 100, by = 10), lty = 3, col = grey(0.8))
barplot(wtd_props, ylim = c(0, 60), add = TRUE, axes = FALSE, 
        names.arg = c("Anti-empire", "Neutral", "Pro-empire"),
        col = lik3)
axis(side = 2, at = seq(0, 100, by = 10), mgp = c(1, 0.4, 0), 
     cex.axis = 0.8)
mtext("%", side = 2, line = 1.7, cex = 1.6)
dev.off()

# plot pro- vs anti-imperial by party pref, left vs right
wtd_props_party = prop.table(xtabs(
  datm$wt_w2 ~ imp_emot_3c_w2 + datm$pref_party_di_w2, na.rm = TRUE), margin = 2) * 100
lik3 = c(brewer.pal(5, "BrBG"))[c(1, 3, 5)]

png("imp_emot_3cat_bypartyy_w2.png", width = 800, height = 600)

layout(matrix(c(1, 2), ncol = 2), widths = c(1, 0.92))  
par(mar = c(3, 2.5, 1.5, 1), las = 1, tcl = -0.2, cex = 1.5,
    mgp = c(1, 0.5, 0))
bp1 = barplot(wtd_props_party[, 1], ylim = c(0, 50), axes = FALSE, names.arg = "", 
              xlab = "", ylab = "", border = NA, col = "white")
abline(h = seq(0, 100, by = 10), lty = 3, col = grey(0.8))
barplot(wtd_props_party[, 1], ylim = c(0, 60), add = TRUE, axes = FALSE, 
        names.arg = "", col = lik3)
axis(side = 2, at = seq(0, 100, by = 10), mgp = c(1, 0.4, 0), 
     cex.axis = 1)
mtext("%", side = 2, line = 1.7, cex = 1.6)
mtext("Prefers left party", side = 3, line = 0.3, cex = 1.6)
mtext(c("Anti-\nempire", "Neutral/\nmixed", "Pro-\nempire"), side = 1, at = bp1, 
      line = 1, cex = 1.5)

par(mar = c(3, 1, 1.5, 0.5), las = 1, tcl = -0.2, cex = 1.6,
    mgp = c(1, 0.5, 0))
bp2 = barplot(wtd_props_party[, 2], ylim = c(0, 50), axes = FALSE, names.arg = "",
              xlab = "", ylab = "", border = NA, col = "white")
abline(h = seq(0, 100, by = 10), lty = 3, col = grey(0.8))
barplot(wtd_props_party[, 2], ylim = c(0, 60), add = TRUE, axes = FALSE, 
        names.arg = "", col = lik3)
axis(side = 2, at = seq(0, 100, by = 10), mgp = c(1, 0.4, 0), 
     cex.axis = 1)
mtext("Prefers right party", side = 3, line = 0.3, cex = 1.6)
mtext(c("Anti-\nempire", "Neutral/\nmixed", "Pro-\nempire"), side = 1, at = bp2, 
      line = 1, cex = 1.5)

dev.off()

# dimensionality
png("imp_emot_w2_scree.png", width = 720, height = 720)
ie_cor_w2 = cor(datm[, imp_emot_qr_w2], use = "pair")
par(mfrow = c(1, 1), mar = c(3, 3, 0.5, 0.5), las = 0, tcl = -0.2, cex = 1.6)
plot(eigen(ie_cor_w2)$values, type = "n", axes = FALSE, xlab = "", ylab = "",
     xlim = c(1, 10), ylim = c(0, 6)) 
abline(h = 0:6, lty = 3, col = grey(0.8))
abline(v = 1:10, lty = 3, col = grey(0.8))
abline(h = 1, lty = 2, col = grey(0.5))
points(eigen(ie_cor_w2)$values, type = "b", lwd = 1.5, pch = 19)
axis(side = 1, at = 1:10, mgp = c(1, 0.35, 0), cex.axis = 1)
axis(side = 2, at = 0:6, mgp = c(1, 0.45, 0), cex.axis = 1)
mtext("Component number", side = 1, line = 1.5, cex = 1.5)
mtext("Eigen values of components", side = 2, line = 1.7, cex = 1.5)
box()
dev.off()

# EFA
imp_emot_w2_efa = fa(datm[, imp_emot_qr_w2], nfactors = 2) # clear 2D pattern 
writeLines(tth(fa2latex(imp_emot_w2_efa)), "imp_emot_w2_efa.html")

# CFA 1d
cfa_imp_emot_1d_w2 = '
  imp_emot_w2 =~ imp_emot_5_w2 + imp_emot_6_w2 + imp_emot_9_w2 + imp_emot_10_w2 + imp_emot_1_w2 + imp_emot_2_w2 + imp_emot_3_w2 + imp_emot_4_w2 
                     + imp_emot_7_w2 + imp_emot_8_w2 '
imp_emot_cfa_1d_w2 = cfa(cfa_imp_emot_1d_w2, data = datm,  estimator = "WLSMV", ordered = TRUE, 
                      std.lv = TRUE, missing = "pairwise")
semTable(imp_emot_cfa_1d_w2, paramSets = c("fits", "loadings", "latentvariances", "latentcovariances"),
         fits=c("chisq", "cfi", "rmsea", "srmr"), columns=c("est", "se", "p"),
         type="html", file="imp_emot_cfa_w2_table.html")
imp_emot_cfa_1d_w2_fit = fitMeasures(
  imp_emot_cfa_1d_w2, output = "matrix", 
  fit.measures = c("cfi", "cfi.robust", "rmsea", "rmsea.robust", "srmr"))
write.csv(imp_emot_cfa_1d_w2_fit, file = "imp_emot_cfa_1d_w2_fit.csv", row.names = TRUE)

# CFA 2d
cfa_imp_emot_w2 = '
  imp_emot_pos_w2 =~ imp_emot_5_w2 + imp_emot_6_w2 + imp_emot_9_w2 + imp_emot_10_w2
  imp_emot_neg_w2 =~ imp_emot_1_w2 + imp_emot_2_w2 + imp_emot_3_w2 + imp_emot_4_w2 
                     + imp_emot_7_w2 + imp_emot_8_w2 '
imp_emot_cfa_w2 = cfa(cfa_imp_emot_w2, data = datm,  estimator = "WLSMV", ordered = TRUE, 
                      std.lv = TRUE, missing = "pairwise")

semTable(imp_emot_cfa_w2, paramSets = c("fits", "loadings", "latentvariances", "latentcovariances"),
         fits=c("chisq", "cfi", "rmsea", "srmr"), columns=c("est", "se", "p"),
         type="html", file="imp_emot_cfa_w2_table.html")
imp_emot_cfa_w2_fit = fitMeasures(
  imp_emot_cfa_w2, output = "matrix", 
  fit.measures = c("cfi", "cfi.robust", "rmsea", "rmsea.robust", "srmr"))
write.csv(imp_emot_cfa_w2_fit, file = "imp_emot_cfa_w2_fit.csv", row.names = TRUE)

# save scores
ids = lavInspect(imp_emot_cfa_w2, "case.idx")
imp_emot_pred_w2 = lavPredict(imp_emot_cfa_w2)
for (i in colnames(imp_emot_pred_w2)) { datm[ids, i] = imp_emot_pred_w2[ , i] }
datm$imp_emot_diff_w2 = datm$imp_emot_pos_w2 - datm$imp_emot_neg_w2

cor(datm$imp_emot_diff_w2, datm$imp_emot_add_w2, use = "pair") # .989

# reliability scales
alpha(datm[, c("imp_emot_5_w2", "imp_emot_6_w2", "imp_emot_9_w2", "imp_emot_10_w2")]) # pos alpha = 0.92
alpha(datm[, c("imp_emot_1_w2", "imp_emot_2_w2", "imp_emot_3_w2", "imp_emot_4_w2", 
               "imp_emot_7_w2", "imp_emot_8_w2")]) # neg alpha = 0.95


# Imperial attitudes and emotions CFA 2d
cfa_imp_nost_2d_w2 = '
  imp_nost_pos_w2 =~ imp_emot_5_w2 + imp_emot_6_w2 + imp_emot_9_w2 + imp_emot_10_w2 + 
                     imp_nost_2_w2 + imp_nost_3_w2 + imp_nost_5_w2 + imp_nost_6_w2 
  imp_nost_neg_w2 =~ imp_emot_1_w2 + imp_emot_2_w2 + imp_emot_3_w2 + imp_emot_4_w2 +
                     imp_emot_7_w2 + imp_emot_8_w2 + imp_nost_1_w2 + imp_nost_4_w2 +
                     imp_nost_7_w2'
imp_nost_2d_cfa_w2 = cfa(cfa_imp_nost_2d_w2, data = datm, estimator = "WLSMV", 
                         ordered = TRUE, std.lv = TRUE, missing = "pairwise")

semTable(imp_nost_2d_cfa_w2, 
         paramSets = c("fits", "loadings", "latentvariances", "latentcovariances"),
         fits = c("chisq", "cfi", "rmsea", "srmr"), columns = c("est", "se", "p"),
         type = "html", file="imp_nost_2d_cfa_w2_table.html")
imp_nost_2d_cfa_w2_fit = fitMeasures(
  imp_nost_2d_cfa_w2, output = "matrix", 
  fit.measures = c("cfi", "cfi.robust", "rmsea", "rmsea.robust", "srmr"))
write.csv(imp_nost_2d_cfa_w2_fit, file = "imp_nost_2d_cfa_w2_fit.csv", row.names = TRUE)



## Wave 3

imp_emot_ques_w3 = c("empire_emotions_1_w3", "empire_emotions_2_w3", "empire_emotions_3_w3", 
                     "empire_emotions_4_w3", "empire_emotions_5_w3", "empire_emotions_6_w3", 
                     "empire_emotions_7_w3", "empire_emotions_8_w3", "empire_emotions_9_w3", 
                     "empire_emotions_10_w3")
datm[, imp_emot_ques_w3] = sapply(datm[, imp_emot_ques_w3], function(x) car::recode(x, "99 = -1"))
sapply(datm[, imp_emot_ques_w3], table, useNA = "ifany")

# items
emot_items_w3 = c("embarrassment", "shame", "guilt", "anger", "pride", "nostalgia", "sadness", 
                  "disgust", "happiness", "satisfaction")

# colour palettes
col4 = c("#b4b1b1", brewer.pal(7, "BrBG")[4:7])
extent4_labs = c("DK", "Not at all", "a small extent", "a moderate extent", "a great extent")

# re-order negative then positive  
emot_ord = c(1:4, 8, 7, 6, 5, 10, 9)

# weighted proportions
wtd_props = list()
for(i in 1:length(imp_emot_ques_w3)) {
  wtd_props[[i]] = prop.table(xtabs(datm$wt_w3 ~ datm[[imp_emot_ques_w3[i]]], na.rm = TRUE)) * 100
}
wtd_props_matrix = do.call(cbind, lapply(wtd_props, function(x) as.numeric(x)))

# Respondent pro vs anti summary
datm_recode = datm
for (item in imp_emot_ques_w3) {
  if (item %in% imp_emot_ques_w3[c(5, 6, 9, 10)]) {
    datm_recode[[item]] = ifelse(datm[[item]] %in% 2:4, 1, 0)
  } else {
    datm_recode[[item]] = ifelse(datm[[item]] %in% 2:4, -1, 0)
  }
}
datm_recode$sum_recode = rowSums(datm_recode[imp_emot_ques_w3], na.rm = TRUE)
sum(datm$wt_w3[datm_recode$sum_recode > 0], na.rm = TRUE) / 
  sum(datm$wt_w3, na.rm = TRUE) * 100 # pro-empire
sum(datm$wt_w3[datm_recode$sum_recode < 0], na.rm = TRUE) / 
  sum(datm$wt_w3, na.rm = TRUE) * 100 # anti-empire

# plot
pdf("imperial_emot_dist_w3.pdf", height = 4, width = 6)
par(mfrow = c(2, 5), mar = c(4.5, 2, 0.5, 0.5), tcl = -0.2, cex = 0.75, las = 2, 
    mgp = c(1.8, 0.9, 0))
for(i in 1:length(imp_emot_ques_w3)) {
  barplot(height = wtd_props[[ emot_ord[i] ]], names.arg = extent4_labs, axes = FALSE, 
          col = col4, cex.names = 0.6, mgp = c(1, 0.2, 0), ylim = c(0, 50))
  abline(h = seq(0, 50, by = 10), lty = 3, lwd = 0.75, col = rgb(0.5, 0.5, 0.5, 0.75))
  axis(side = 2, labels = TRUE, cex.axis = 0.6, mgp = c(1, 0.4, 0))
  barplot(height = wtd_props[[ emot_ord[i] ]], col = col4, axes = FALSE, names.arg = "", add = TRUE)
  text(x = 3, y = 48, emot_items_w2[emot_ord[i]], cex = 0.7, adj = 0.5)
}
dev.off()

# plot stacked barcharts
extent4_labs = c("Don't\nknow", "Not at\nall", "A small\nextent", 
                 "A moderate\nextent", "A great\nextent")

png("imperial_emot_stack_w3.png", height = 1080, width = 1440)
par(mar = c(2, 7.5, 2, 1.5), tcl = -0.2, cex = 1.8, las = 0, mgp = c(3, 0.5, 0))
barplot(height = wtd_props_matrix[, emot_ord], names.arg = emot_items_w3[emot_ord], 
        axes = FALSE, cex.names = 1, xlim = c(0, 100),
        beside = FALSE, las = 2, space = 0.4, horiz = TRUE)
segments(x0= seq(0, 100, by = 10), y0 = 0, y1 = 13, lty = 3, lwd = 1, 
         col = rgb(0, 0, 0, 0.75))
barplot(height = wtd_props_matrix[, emot_ord], names.arg = NULL, col = col4,
        axes = FALSE, beside = FALSE, add = TRUE, space = 0.4, horiz = TRUE)
axis(side = 1, at = seq(0, 100, by = 10), labels = TRUE, cex.axis = 1, mgp = c(3, 0.5, 0))
par(xpd = TRUE)
legend("top", legend = extent4_labs, fill = col4, cex = 1, pt.cex = 1.2, 
       horiz = TRUE, bty= "n", inset = c(0, -0.05), x.intersp = 0.8)
dev.off()

# overall % nostalgic vs aversive
grouped_props = wtd_props_matrix[3:4, ] %>%
  colSums() %>%
  `*`(c(-1, -1, -1, -1, 1, 1, -1, -1, 1, 1))

# plot collapsed categories with valences
png("imperial_emot_collap_stack_w3.png", height = 1080, width = 1620)
layout(matrix(c(1, 2), ncol = 2), widths = c(4, 1.5))  # Adjust layout to have two panels

par(mar = c(2, 7.5, 1.5, 1.5), tcl = -0.2, cex = 1.9, las = 0, mgp = c(3, 0.5, 0))
barplot(height = wtd_props_matrix[, emot_ord], names.arg = emot_items_w3[emot_ord], col = col4, 
        axes = FALSE, cex.names = 1, xlim = c(0, 100),
        beside = FALSE, las = 2, space = 0.4, horiz = TRUE)
segments(x0= seq(0, 100, by = 10), y0 = 0, y1 = 13, lty = 3, lwd = 1, 
         col = rgb(0, 0, 0, 0.75))
barplot(height = wtd_props_matrix[, emot_ord], names.arg = NULL, col = col4,
        axes = FALSE, beside = FALSE, add = TRUE, space = 0.4, horiz = TRUE)
axis(side = 1, at = seq(0, 100, by = 10), labels = TRUE, cex.axis = 1, mgp = c(3, 0.5, 0))
par(xpd = TRUE)
legend("top", legend = extent4_labs, fill = col4, cex = 1, pt.cex = 1.2, 
       horiz = TRUE, bty= "n", inset = c(0, -0.05), x.intersp = 0.8)

par(mar = c(2, 2, 1.5, 1.5), xpd = TRUE)
barplot(height = grouped_props[emot_ord], names.arg = NULL, horiz = TRUE, 
        col = "grey40", xlim = range(c(-30, 30)), space = 0.4)
segments(x0 = c(-20, 0, 20), y0 = 0, y1 = 14, lty = 3, 
         lwd = 1, col = rgb(0, 0, 0, 0.75))
barplot(height = grouped_props[emot_ord], names.arg = NULL, horiz = TRUE, 
        col = "grey40", xlim = range(c(-30, 30)), space = 0.4, add = TRUE)
legend("top", legend = "% nostalgic (moderate + great extent)", fill = "grey40", 
       cex = 1.1, pt.cex = 1.2, 
       horiz = TRUE, bty = "n", inset = c(0, -0.02), x.intersp = 0.8)

dev.off()

# recodes for scales
imp_emot_qr_w3 = c("imp_emot_1_w3", "imp_emot_2_w3", "imp_emot_3_w3", "imp_emot_4_w3", "imp_emot_5_w3", 
                   "imp_emot_6_w3", "imp_emot_7_w3", "imp_emot_8_w3", "imp_emot_9_w3", "imp_emot_10_w3")
datm[, imp_emot_qr_w3] = sapply(datm[, imp_emot_ques_w3], function(x) car::recode(x, "-1 = NA"))
sapply(datm[, imp_emot_qr_w3], table, useNA = "ifany")

# plot pro- vs anti-imperial
datm$imp_emot_add_w3 = rowMeans(datm[, imp_emot_qr_w3[c(5, 6, 9, 10)]], na.rm = TRUE) -
  rowMeans(datm[, imp_emot_qr_w3[c(1, 2, 3, 4, 7, 8)]], na.rm = TRUE)
imp_emot_3c_w3 = cut(datm$imp_emot_add_w3, breaks = c(-3.5, -0.5, 0.5, 3.5))
wtd_props = prop.table(xtabs(datm$wt_w3 ~ imp_emot_3c_w3, na.rm = TRUE)) * 100

lik3 = c(brewer.pal(5, "BrBG"))[c(1, 3, 5)]
png("imp_emot_3cat_w3.png", width = 600, height = 400)
par(mfrow = c(1, 1), mar = c(2, 3, 0.5, 0.5), las = 1, tcl = -0.2, cex = 1.7,
    mgp = c(1, 0.3, 0))
barplot(wtd_props, ylim = c(0, 50), axes = FALSE, names.arg = "", 
        xlab = "", ylab = "", border = NA, col = "white")
abline(h = seq(0, 100, by = 10), lty = 3, col = grey(0.8))
barplot(wtd_props, ylim = c(0, 60), add = TRUE, axes = FALSE, 
        names.arg = c("Anti-empire", "Neutral", "Pro-empire"),
        col = lik3)
axis(side = 2, at = seq(0, 100, by = 10), mgp = c(1, 0.4, 0), 
     cex.axis = 0.8)
mtext("%", side = 2, line = 1.7, cex = 1.6)
dev.off()

# scree plot
png("imp_emot_w3_scree.png", width = 720, height = 720)
ie_cor_w3 = cor(datm[, imp_emot_qr_w3], use = "pair")
par(mfrow = c(1, 1), mar = c(3, 3, 0.5, 0.5), las = 0, tcl = -0.2, cex = 1.6)
plot(eigen(ie_cor_w3)$values, type = "n", axes = FALSE, xlab = "", ylab = "",
     xlim = c(1, 10), ylim = c(0, 6)) 
abline(h = 0:6, lty = 3, col = grey(0.8))
abline(v = 1:10, lty = 3, col = grey(0.8))
abline(h = 1, lty = 2, col = grey(0.5))
points(eigen(ie_cor_w3)$values, type = "b", lwd = 1.5, pch = 19)
axis(side = 1, at = 1:10, mgp = c(1, 0.35, 0), cex.axis = 1)
axis(side = 2, at = 0:6, mgp = c(1, 0.45, 0), cex.axis = 1)
mtext("Component number", side = 1, line = 1.5, cex = 1.5)
mtext("Eigen values of components", side = 2, line = 1.7, cex = 1.5)
box()
dev.off()

# EFA
imp_emot_w3_efa = fa(datm[, imp_emot_qr_w3], nfactors = 2) # clear 2D pattern 
writeLines(tth(fa2latex(imp_emot_w3_efa)), "imp_emot_w3_efa.html")

# CFA 1d
cfa_imp_emot_1d_w3 = '
  imp_emot_w3 =~ imp_emot_5_w3 + imp_emot_6_w3 + imp_emot_9_w3 + imp_emot_10_w3 + imp_emot_1_w3 + imp_emot_2_w3 + imp_emot_3_w3 + imp_emot_4_w3 
                     + imp_emot_7_w3 + imp_emot_8_w3 '
imp_emot_cfa_1d_w3 = cfa(cfa_imp_emot_1d_w3, data = datm,  estimator = "WLSMV", ordered = TRUE, 
                         std.lv = TRUE, missing = "pairwise")
semTable(imp_emot_cfa_1d_w3, paramSets = c("fits", "loadings", "latentvariances", "latentcovariances"),
         fits=c("chisq", "cfi", "rmsea", "srmr"), columns=c("est", "se", "p"),
         type="html", file="imp_emot_cfa_w3_table.html")
imp_emot_cfa_1d_w3_fit = fitMeasures(
  imp_emot_cfa_1d_w3, output = "matrix", 
  fit.measures = c("cfi", "cfi.robust", "rmsea", "rmsea.robust", "srmr"))
write.csv(imp_emot_cfa_1d_w3_fit, file = "imp_emot_cfa_1d_w3_fit.csv", row.names = TRUE)

# CFA 2d
cfa_imp_emot_w3 = '
  imp_emot_pos_w3 =~ imp_emot_5_w3 + imp_emot_6_w3 + imp_emot_9_w3 + imp_emot_10_w3
  imp_emot_neg_w3 =~ imp_emot_1_w3 + imp_emot_2_w3 + imp_emot_3_w3 + imp_emot_4_w3 
                     + imp_emot_7_w3 + imp_emot_8_w3 '
imp_emot_cfa_w3 = cfa(cfa_imp_emot_w3, data = datm,  estimator = "WLSMV", ordered = TRUE, 
                      std.lv = TRUE, missing = "pairwise")

semTable(imp_emot_cfa_w3, paramSets = c("fits", "loadings", "latentvariances", "latentcovariances"),
         fits=c("chisq", "cfi", "rmsea", "srmr"), columns=c("est", "se", "p"),
         type="html", file="imp_emot_cfa_w3_table.html")
imp_emot_cfa_w3_fit = fitMeasures(
  imp_emot_cfa_w3, output = "matrix", 
  fit.measures = c("cfi", "cfi.robust", "rmsea", "rmsea.robust", "srmr"))
write.csv(imp_emot_cfa_w3_fit, file = "imp_emot_cfa_w3_fit.csv", row.names = TRUE)

# save scores
ids = lavInspect(imp_emot_cfa_w3, "case.idx")
imp_emot_pred_w3 = lavPredict(imp_emot_cfa_w3)
for (i in colnames(imp_emot_pred_w3)) { datm[ids, i] = imp_emot_pred_w3[ , i] }
datm$imp_emot_diff_w3 = datm$imp_emot_pos_w3 - datm$imp_emot_neg_w3

cor(datm$imp_emot_diff_w3, datm$imp_emot_add_w3, use = "pair") # .989

# reliability scales
alpha(datm[, c("imp_emot_5_w3", "imp_emot_6_w3", "imp_emot_9_w3", "imp_emot_10_w3")]) # pos alpha = 0.94
alpha(datm[, c("imp_emot_1_w3", "imp_emot_2_w3", "imp_emot_3_w3", "imp_emot_4_w3", 
               "imp_emot_7_w3", "imp_emot_8_w3")]) # neg alpha = 0.95

# correlations
write.csv(lowerCor(datm[, c("imp_att_w2", "imp_att_w3", "imp_emot_pos_w2", "imp_emot_pos_w3", 
                            "imp_emot_neg_w2", "imp_emot_neg_w3", "imp_emot_diff_w2", "imp_emot_diff_w3")]),
          file = "imp_nost_corr.csv", row.names = TRUE)



### Immigration----


## wave 2

immig_ques_w2 = c("Q9_scale_w2", "Q10_scale_w2")
sapply(datm[, immig_ques_w2], table, useNA = "ifany")

# recodes
immig_qr_w2 = c("immig_eco_w2", "immig_cult_w2")
datm[, immig_qr_w2] = sapply(datm[, immig_ques_w2], function(x) car::recode(x, "997=NA"))
table(datm$immig_eco, datm$immig_cult, useNA = 'ifany')

# test scale
alpha(datm[, immig_qr_w2], check.keys = TRUE) # a = 0.94, r = 0.89

# cfa
cfa_immig_supp = 'immig_supp_w2 =~ 1*immig_eco_w2 + 1*immig_cult_w2'
immig_supp_cfa_w2= cfa(cfa_immig_supp, data = datm, estimator = "WLSMV", std.lv = TRUE, 
                       missing = "pairwise")
summary(immig_supp_cfa_w2, fit.measures = TRUE)
ids = lavInspect(immig_supp_cfa_w2, "case.idx")
immig_supp_pred = lavPredict(immig_supp_cfa_w2)
datm[ids, "immig_supp_w2"] = immig_supp_pred

## wave 3

immig_ques_w3 = c("Immigration_w3", "Immigration_culture_w3")
sapply(datm[, immig_ques_w3], table, useNA = "ifany")

# recodes
immig_qr_w3 = c("immig_eco_w3", "immig_cult_w3")
datm[, immig_qr_w3] = sapply(datm[, immig_ques_w3], function(x) car::recode(x, "99=NA"))
table(datm$immig_eco, datm$immig_cult, useNA = 'ifany')

# test scale
alpha(datm[, immig_qr_w3], check.keys = TRUE) # a = 0.94, r = 0.89

# create variables
cfa_immig_supp = 'immig_supp_w3 =~ 1*immig_eco_w3 + 1*immig_cult_w3'
immig_supp_cfa_w3= cfa(cfa_immig_supp, data = datm, estimator = "WLSMV", std.lv = TRUE, 
                       missing = "pairwise")
summary(immig_supp_cfa_w3, fit.measures = TRUE)
ids = lavInspect(immig_supp_cfa_w3, "case.idx")
immig_supp_pred = lavPredict(immig_supp_cfa_w3)
datm[ids, "immig_supp_w3"] = immig_supp_pred



### Political values----

## Wave 2

val_ques_w2 = c("Q22_1_w2", "Q22_2_w2", "Q22_3_w2", "Q22_4_w2", "Q22_5_w2", "Q22_6_w2", 
             "Q22_7_w2", "Q22_8_w2")
sapply(datm[, val_ques_w2], table, useNA="ifany")

# recodes
val_qr_w2 = c("pol_val_1_w2", "pol_val_2_w2", "pol_val_3_w2", "pol_val_4_w2", "pol_val_5_w2", 
              "pol_val_6_w2", "pol_val_7_w2", "pol_val_8_w2")
datm[, val_qr_w2] = sapply(datm[, val_ques_w2], function(x) car::recode(x, "99=3"))
sapply(datm[, val_qr_w2], table, useNA="ifany")

# test scale
alpha(datm[, val_qr_w2[1:4]]) # alpha = 0.79
alpha(datm[, val_qr_w2[5:8]]) # alpha = 0.81 
plot(eigen(cor(datm[, val_qr_w2], use="pair"))$values, type="b") # 2 dims
fa(datm[, val_qr_w2], nfactors = 2) # 2 dim

# polarization on scale
eco_val_add_w2 = rowMeans(
  datm[, c("pol_val_1_w2", "pol_val_2_w2", "pol_val_3_w2", "pol_val_4_w2")], na.rm = TRUE)
hist(eco_val_add_w2, c(0.5, 1.5, 2.5, 3.5, 4.5, 5.5))

# cfa
cfa_mod_val_w2 = '
  left_ecoval_w2 =~ pol_val_1_w2 + pol_val_2_w2 + pol_val_3_w2 + pol_val_4_w2
  auth_val_w2 =~ pol_val_5_w2 + pol_val_6_w2 + pol_val_7_w2 + pol_val_8_w2 '
val_cfa_w2 = cfa(cfa_mod_val_w2, data = datm,  estimator = "WLSMV", ordered = TRUE, 
                 std.lv = TRUE, missing = "pairwise")
summary(val_cfa_w2, fit.measures = TRUE)
ids = lavInspect(val_cfa_w2, "case.idx")
val_pred_w2 = lavPredict(val_cfa_w2)
for (i in colnames(val_pred_w2)) { datm[ids, i] = val_pred_w2[ , i] }


## Wave 3

val_ques_w3 = c("Political_values_1_w3", "Political_values_2_w3", "Political_values_3_w3", 
                "Political_values_4_w3", "Political_values_5_w3", "Political_values_6_w3", 
                "Political_values_7_w3", "Political_values_8_w3")
sapply(datm[, val_ques_w3], table, useNA="ifany")

# recodes
val_qr_w3 = c("pol_val_1_w3", "pol_val_2_w3", "pol_val_3_w3", "pol_val_4_w3", "pol_val_5_w3", 
              "pol_val_6_w3", "pol_val_7_w3", "pol_val_8_w3")
datm[, val_qr_w3] = sapply(datm[, val_ques_w3], function(x) car::recode(x, "99=3"))
sapply(datm[, val_qr_w3], table, useNA="ifany")

# test scale
alpha(datm[, val_qr_w3[1:4]]) # alpha = 0.83
alpha(datm[, val_qr_w3[5:8]]) # alpha = 0.81 
plot(eigen(cor(datm[, val_qr_w3], use="pair"))$values, type="b") # 2 dims
fa(datm[, val_qr_w3], nfactors = 2) # 2 dim

# cfa
cfa_mod_val_w3 = '
  left_ecoval_w3 =~ pol_val_1_w3 + pol_val_2_w3 + pol_val_3_w3 + pol_val_4_w3
  auth_val_w3 =~ pol_val_5_w3 + pol_val_6_w3 + pol_val_7_w3 + pol_val_8_w3 '
val_cfa_w3 = cfa(cfa_mod_val_w3, data = datm,  estimator = "WLSMV", ordered = TRUE, 
                 std.lv = TRUE, missing = "pairwise")
summary(val_cfa_w3, fit.measures = TRUE)
ids = lavInspect(val_cfa_w3, "case.idx")
val_pred_w3 = lavPredict(val_cfa_w3)
for (i in colnames(val_pred_w3)) { datm[ids, i] = val_pred_w3[ , i] }

# test-retest
lowerCor(datm[, c("left_ecoval_w2", "left_ecoval_w3", "auth_val_w2", "auth_val_w3")])



### Democratic support----

## wave 2

dem_supp_q_w2 = c("Q13_1_w2", "Q13_2_w2", "Q13_3_w2", "Q13_4_w2", "Q13_5_w2", "Q13_6_w2",
                  "Q13_7_w2", "Q13_8_w2")
sapply(datm[, dem_supp_q_w2], table, useNA="ifany")

# recodes
dem_supp_qr_w2 = c("democ_frexp2_w2", "democ_frassc1_w2", "democ_unisuff1a_w2", 
                   "democ_frelect2a_w2", "democ_decelec2a_w2", 
                   "democ_judcnstr2a_w2", "democ_legcnstr1a_w2", "democ_eqlaw1_w2")
datm[, dem_supp_qr_w2] = sapply(datm[, dem_supp_q_w2], 
                                function(x) car::recode(x, "1=5; 2=4; 3=3; 4=2; 5=1; 99=3"))
sapply(datm[, dem_supp_qr_w2], table, useNA="ifany")

# test scale
alpha(datm[, dem_supp_qr_w2], check.keys = TRUE) # alpha = 0.83
fa(datm[, dem_supp_qr_w2], nfactors = 1) 

# cfa
cfa_dem_supp_w2 = '
  dem_sup_w2 =~ democ_frexp2_w2 + democ_frassc1_w2 + democ_unisuff1a_w2 + democ_frelect2a_w2 
                    + democ_decelec2a_w2 + democ_judcnstr2a_w2 + democ_legcnstr1a_w2
                    + democ_eqlaw1_w2'
demsup_cfa_w2 = cfa(cfa_dem_supp_w2, data = datm,  estimator = "WLSMV", ordered = TRUE, 
                    std.lv = TRUE, missing = "pairwise")
semTable(demsup_cfa_w2, paramSets = c("fits", "loadings", "latentvariances"),
         fits=c("chisq", "cfi", "rmsea", "srmr"), columns=c("est", "se", "p"),
         type="html", file="demsup_cfa_w2_table.html")
demsup_cfa_w2_fit = fitMeasures(
  demsup_cfa_w2, output = "matrix", 
  fit.measures = c("cfi", "cfi.robust", "rmsea", "rmsea.robust", "srmr"))
write.csv(demsup_cfa_w2_fit, file = "demsup_cfa_w2_fit.csv", row.names = TRUE)

# save scores
ids = lavInspect(demsup_cfa_w2, "case.idx")
demsup_pred_w2 = lavPredict(demsup_cfa_w2)
datm[ids, "dem_sup_w2"] = demsup_pred_w2



### Hostile sexism----

## wave 2

sexism_qr_w2 = c("Q17_1_w2", "Q17_2_w2", "Q17_3_w2")
sapply(datm[, sexism_qr_w2], table, useNA="ifany")

# recodes
hos_sexism_qr = c("sexism_1_w2", "sexism_2_w2", "sexism_3_w2")
datm[, hos_sexism_qr] = sapply(datm[, sexism_qr_w2], function(x) car::recode(x, "99=3"))
sapply(datm[, hos_sexism_qr], table, useNA="ifany")

# test scale
alpha(datm[, hos_sexism_qr], check.keys = TRUE) # alpha = 0.85; ave r = 0.66
fa(datm[, hos_sexism_qr], nfactors = 1) 

# CFA
cfa_mod_sexism = 'sexism =~ sexism_1_w2 + sexism_2_w2 + sexism_3_w2'
sexism_cfa_w2 = cfa(cfa_mod_sexism, data = datm,  estimator = "WLSMV", ordered = TRUE, 
                    std.lv = TRUE, missing = "pairwise")
summary(sexism_cfa_w2, fit.measures = TRUE)
ids = lavInspect(sexism_cfa_w2, "case.idx")
sexism_pred = lavPredict(sexism_cfa_w2)
datm[ids, "hos_sexism_w2"] = sexism_pred



### Group status----

## wave 2

grp_stat_ques = c("Q5_scale_w2", "Q6_scale_w2", "Q7_scale_w2")
sapply(datm[, grp_stat_ques], table, useNA = "ifany")

# recodes
grp_stat_qr = c("grp_stat_1_w2", "grp_stat_2_w2", "grp_stat_3_w2")
datm[, grp_stat_qr] = sapply(datm[, grp_stat_ques], function(x) car::recode(x, "997=NA"))
sapply(datm[, grp_stat_qr], table, useNA = "ifany")

# test scale
lowerCor(datm[, grp_stat_qr])
alpha(datm[, grp_stat_qr], check.keys = TRUE) # 2 & 3 cohere, not 1
fa(datm[, grp_stat_qr], nfactors = 1) 

# create variables
datm$grp_polstat_w2 = rowMeans(datm[, grp_stat_qr[2:3]])



### Trust, satisfaction, ideological position----

## wave 1

sys_ques = c("core_plt_1_1_w1", "core_plt_1_2_w1", "core_plt_1_3_w1", "core_plt_1_4_w1", 
             "core_sat_1_w1", "core_pos_9_w1")
sapply(datm[, sys_ques], table, useNA="ifany")

# recodes
sys_qr = c("trust_gov_w1", "trust_parl_w1", "trust_legsys_w1", "trust_pols_w1", 
           "satis_dem_w1", "left_right_w1")
datm[, sys_qr] = sapply(datm[, sys_ques], function(x) car::recode(x, "12=NA"))
sapply(datm[, sys_qr], table, useNA="ifany")

# test trust scale
alpha(datm[, sys_qr[1:4]], check.keys = TRUE) # a = 0.87, r = 0.62
fa(datm[, sys_qr[1:4]], nfactors = 1) 

# create variables
datm$pol_trust_w1 = rowMeans(datm[, sys_qr[1:4]]) - 1


## wave 2

sys_ques = c("Q18scale_w2", "Q19scale_w2", "Q20scale_w2", "q20_scale_w2", "Q11scale_w2")
sapply(datm[, sys_ques], table, useNA="ifany")

# recodes
sys_qr = c("satis_dem_w2", "trust_pols_w2", "govt_appr_w2", "left_right_w2", "EU_indep_w2")
datm[, sys_qr] = sapply(datm[, sys_ques], function(x) car::recode(x, "997=NA"))
sapply(datm[, sys_qr], table, useNA="ifany")



### Populism----

## wave 2

popul_ques = c("Q14_1_w2", "Q14_2_w2", "Q14_3_w2", "Q14_4_w2")
sapply(datm[, popul_ques], table, useNA = "ifany")

# recodes
popul_qr = c("pop_1_w2", "pop_2_w2", "pop_3_w2", "pop_4_w2")
datm[, popul_qr] = sapply(datm[, popul_ques], function(x) car::recode(x, "99 = 3"))
sapply(datm[, popul_qr], table, useNA = "ifany")

# test scale
alpha(datm[, popul_qr], check.keys = TRUE) # a = 0.78, r = 0.50
fa(datm[, popul_qr], nfactors = 1) 

# CFA
cfa_mod_popul_att = 'popul_att_w2 =~ pop_1_w2 + pop_2_w2 + pop_3_w2 + pop_4_w2'
popul_att_cfa_1 = cfa(cfa_mod_popul_att, data = datm, estimator = "WLSMV", ordered = TRUE, 
                      std.lv = TRUE, missing = "pairwise")
summary(popul_att_cfa_1, fit.measures = TRUE)
ids = lavInspect(popul_att_cfa_1, "case.idx")
popul_att_pred = lavPredict(popul_att_cfa_1)
datm[ids, "popul_att_w2"] = popul_att_pred

# life better 50 years ago
datm$uk_ethant_7_w1 = car::recode(datm$uk_ethant_7_w1, "99 = 3; 999 = NA")


### National identities----

## wave 2

natid_ques = c("Q31_Scottish_w2", "Q31_English_w2", "Q31_Welsh_w2", "Q31_British_w2", 
               "Q31_Irish_w2", "Q31_NorthernIrish_w2", "Q31_European_w2")
sapply(datm[, natid_ques], table, useNA = "ifany")

# recodes
nat_id_qr = c("nat_id_scot_w2", "nat_id_eng_w2", "nat_id_wel_w2", "nat_id_brit_w2", 
             "nat_id_ir_w2", "nat_id_ni_w2", "nat_id_eur_w2")
datm[, nat_id_qr] = sapply(datm[, natid_ques], function(x) car::recode(x, "99 = NA"))
sapply(datm[, nat_id_qr], table, useNA = "ifany")
lowerCor(datm[, nat_id_qr])


## wave 3

nat_id_ques = c("national_identity_Scot_w3", "national_identity_Eng_w3", 
                "national_identity_Wel_w3", "national_identity_Brit_w3", 
                "national_identity_Irish_w3", "national_identity_NI_w3", 
                "national_identity_Euro_w3")
sapply(datm[, nat_id_ques], table, useNA = "ifany")

# recodes
nat_id_qr = c("nat_id_scot_w3", "nat_id_eng_w3", "nat_id_wel_w3", "nat_id_brit_w3", 
              "nat_id_ir_w3", "nat_id_ni_w3", "nat_id_eur_w3")
datm[, nat_id_qr] = sapply(datm[, nat_id_ques], function(x) car::recode(x, "99 = NA"))
sapply(datm[, nat_id_qr], table, useNA = "ifany")
lowerCor(datm[, nat_id_qr])



### National pride / chauvinism----

## wave 3

chauv_ques = c("chauvinism_1_w3", "chauvinism_2_w3", "chauvinism_3_w3", "chauvinism_4_w3", 
               "chauvinism_5_w3", "chauvinism_6_w3")
sapply(datm[, chauv_ques], table, useNA="ifany")

nat_chauv_qr = c("nat_chauv_1_w3", "nat_chauv_2_w3", "nat_chauv_3_w3", "nat_chauv_4_w3", 
                 "nat_chauv_5_w3", "nat_chauv_6_w3")
datm[, nat_chauv_qr] = sapply(datm[, chauv_ques], function(x) car::recode(x, "99 = 3"))
sapply(datm[, nat_chauv_qr], table, useNA="ifany")

# test scale - 2D pride/chauvinism apparent
fa(datm[, nat_chauv_qr], nfactors = 1)
fa(datm[, nat_chauv_qr], nfactors = 2)
psych::alpha(datm[, nat_chauv_qr])

# CFA
cfa_mod_nat_chauv = '
 nat_pride_w3 =~ nat_chauv_2_w3 + nat_chauv_5_w3 + nat_chauv_6_w3 
 nat_chauv_w3 =~ nat_chauv_1_w3 + nat_chauv_3_w3 + nat_chauv_4_w3 
'
nat_chauv_cfa_1 = cfa(cfa_mod_nat_chauv, data = datm, estimator = "WLSMV", ordered = TRUE, 
                      std.lv = TRUE, missing = "pairwise")
semTable(nat_chauv_cfa_1, paramSets = c("fits", "loadings", "latentvariances", "latentcovariances"),
         fits=c("chisq", "cfi", "rmsea", "srmr"), columns=c("est", "se", "p"),
         type="html", file="nat_chauv_cfa_w3_table.html")

# save scores
ids = lavInspect(nat_chauv_cfa_1, "case.idx")
nat_chauv_pred_w3 = lavPredict(nat_chauv_cfa_1)
for (i in colnames(nat_chauv_pred_w3)) { datm[ids, i] = nat_chauv_pred_w3[ , i] }



### Imperial associations of parties----

## wave 3

party_emp_ques = c("imperial_associations_Lab_w3", "imperial_associations_Con_w3", 
                   "imperial_associations_Lib_w3", "imperial_associations_Green_w3", 
                   "imperial_associations_Reform_w3", "imperial_associations_SNP_w3", 
                   "imperial_associations_Plaid_w3")
sapply(datm[, party_emp_ques], table, useNA = "ifany")

imp_ass_r = c("imp_ass_Lab_w3", "imp_ass_Cons_w3", "imp_ass_LibDem_w3", "imp_ass_Green_w3", 
              "imp_ass_Reform_w3", "imp_ass_SNP_w3", "imp_ass_Plaid_w3")
datm[, imp_ass_r] = sapply(datm[, party_emp_ques], function(x) car::recode(x, "99 = NA"))
sapply(datm[, imp_ass_r], table, useNA = "ifany")
sapply(datm[, imp_ass_r], mean, na.rm = TRUE)

# items labels
imp_ass_party = c("Labour", "Conservative", "Liberal\nDemocrat", "Green", "Reform", "SNP", "Plaid")

# colour palettes
col5 = c(brewer.pal(5, "BrBG"))

# weighted means
wtd_mean_mat = matrix(nrow = length(imp_ass_r), ncol = 2)
colnames(wtd_mean_mat) = c("coef", "SE")
rownames(wtd_mean_mat) = imp_ass_party
for (i in seq_along(imp_ass_r)) {
  model = lm(as.formula(paste(imp_ass_r[i], "~ 1")), weights = wt_w3, data = datm)
  wtd_mean_mat[i, "coef"] = summary(model)$coefficients[1, "Estimate"]
  wtd_mean_mat[i, "SE"] = summary(model)$coefficients[1, "Std. Error"]
  wtd_mean_mat[i, "SE"] = summary(model)$coefficients[1, "Std. Error"]
}
wtd_mean_mat = wtd_mean_mat[order(wtd_mean_mat[, "coef"]), ]

# plot dotplot
png("imperial_assoc_party_w3.png", height = 720, width = 1080)
par(mar = c(2, 6.2, 1, 1), tcl = -0.2, cex = 1.8, las = 1, mgp = c(3, 0.5, 0))
plot(wtd_mean_mat[, "coef"], 1:length(imp_ass_party), pch = 16, col = col5[5], 
     axes = FALSE, ylab = "", cex = 1.6, xlim = c(1, 9), las = 2)
abline(h = seq(1, 7, by = 1), lty = 1, lwd = 1, col = rgb(0.7, 0.7, 0.7, 0.75))
abline(v = seq(1, 9, by = 1), lty = 3, lwd = 1, col = rgb(0.7, 0.7, 0.7, 0.75))
points(wtd_mean_mat[, "coef"], 1:length(imp_ass_party), pch = 16, cex = 1.6, col = col5[5])
segments(x0 = wtd_mean_mat[, "coef"] + 1.96 * wtd_mean_mat[, "SE"], 
         x1 = wtd_mean_mat[, "coef"] - 1.96 * wtd_mean_mat[, "SE"],
         y0 = 1:7, lty = 1, lwd = 3, col = col5[5])
axis(side = 1, at = seq(1, 9, by = 1), labels = TRUE, cex.axis = 1, mgp = c(3, 0.5, 0))
axis(side = 2, at = 1:7, labels = row.names(wtd_mean_mat), cex.axis = 1.2, 
     mgp = c(1, 0.4, 0))
dev.off()



### Nostalgic deprivation----

## wave 1

nost_dep_ques = c("uk_poldep_n_w1", "uk_poldep_p_w1", "uk_poldep_f_w1", 
                  "uk_ecodep_n_w1", "uk_ecodep_p_w1", "uk_ecodep_f_w1")
sapply(datm[, nost_dep_ques], table, useNA = "ifany")

datm[, nost_dep_ques] = sapply(datm[, nost_dep_ques], 
                               function(x) car::recode(x, "12 = NA; 99 = NA"))
sapply(datm[, nost_dep_ques], table, useNA = "ifany")

# recode
datm$pol_stat_decl_w1 = datm$uk_poldep_n_w1 - datm$uk_poldep_p_w1
datm$eco_stat_decl_w1 = datm$uk_ecodep_n_w1 - datm$uk_ecodep_p_w1



### EFA / convergent & divergent validity----

vars2 = c(
  "imp_nost_1_w2", "imp_nost_2_w2", "imp_nost_3_w2", "imp_nost_4_w2", "imp_nost_5_w2", 
  "imp_nost_6_w2", "imp_nost_7_w2", "imp_emot_1_w2", "imp_emot_2_w2", "imp_emot_3_w2", 
  "imp_emot_4_w2", "imp_emot_5_w2", "imp_emot_6_w2", "imp_emot_7_w2", "imp_emot_8_w2", 
  "imp_emot_9_w2", "imp_emot_10_w2"
)

scree(datm[, vars2], factors = TRUE, pc = FALSE)
fa.parallel(datm[, vars2])
efa2_8 = fa(datm[, vars2], fm = "minres", cor = "mixed", nfactors = 2)
print(efa2_8, cut = 0.3)

## Wave 2

vars2 = c(
  "imp_nost_1_w2", "imp_nost_2_w2", "imp_nost_3_w2", "imp_nost_4_w2", "imp_nost_5_w2", 
  "imp_nost_6_w2", "imp_nost_7_w2", "immig_eco_w2", "immig_cult_w2", "pol_val_1_w2", 
  "pol_val_2_w2", "pol_val_3_w2", "pol_val_4_w2", "pol_val_5_w2", "pol_val_6_w2", 
  "pol_val_7_w2", "pol_val_8_w2", "sexism_1_w2", "sexism_2_w2", "sexism_3_w2", 
  "democ_frexp2_w2", "democ_frassc1_w2", "democ_unisuff1a_w2", 
  "democ_frelect2a_w2", "democ_decelec2a_w2", "democ_judcnstr2a_w2", 
  "democ_legcnstr1a_w2", "democ_eqlaw1_w2", "satis_dem_w2", "trust_pols_w2", 
  "govt_appr_w2", "left_right_w2", "EU_indep_w2", "pop_1_w2", "pop_2_w2", "pop_3_w2", 
  "pop_4_w2", "grp_stat_1_w2", "grp_stat_2_w2", "grp_stat_3_w2"
  )

scree(datm[, vars2], factors = TRUE, pc = FALSE)
fa.parallel(datm[, vars2], fa = "fa") # 8 factors
efa2_8 = fa(datm[, vars2], fm = "minres", max.iter = 200, nfactors = 8)
print(efa2_8, cut = 0.3)

# Extract loadings and edit
loadings2_8 = as.matrix(efa2_8$loadings)  
loadings2_8 = apply(loadings2_8, 2, round, 2)
loadings2_8 = apply(loadings2_8, 2, function(x) ifelse(abs(x) < 0.3, "", x))
loadings2_8 = as.data.frame(loadings2_8)

# Save EFA
efa2_8_table <- kable(loadings2_8, format = "latex", booktabs = TRUE, 
                caption = "EFA") %>%
  kable_styling(latex_options = c("hold_position"))
writeLines(efa2_8_table, "efa_table_all_w2.tex")


## Wave 3

vars3 = c(
  "imp_nost_1_w3", "imp_nost_2_w3", "imp_nost_3_w3", "imp_nost_4_w3", 
  "imp_nost_5_w3", "imp_nost_6_w3", "imp_nost_7_w3", "immig_eco_w3", "immig_cult_w3", 
  "pol_val_1_w3", "pol_val_2_w3", "pol_val_3_w3", "pol_val_4_w3", "pol_val_5_w3", 
  "pol_val_6_w3", "pol_val_7_w3", "pol_val_8_w3", "nat_id_eng_w3", "nat_id_brit_w3", 
  "nat_chauv_1_w3", "nat_chauv_2_w3", "nat_chauv_3_w3", "nat_chauv_4_w3", 
  "nat_chauv_5_w3", "nat_chauv_6_w3"
  )

scree(datm[, vars3], factors = TRUE, pc = FALSE)
fa.parallel(datm[, vars3], fa = "fa") # 6 factors
efa3_6 = fa(datm[, vars3], fm = "minres", max.iter = 200, nfactors = 5)
print(efa3_6, cut = 0.3)

# Extract loadings and edit
loadings3_6 = as.matrix(efa3_6$loadings)  
loadings3_6 = apply(loadings3_6, 2, round, 2)
loadings3_6 = apply(loadings3_6, 2, function(x) ifelse(abs(x) < 0.3, "", x))
loadings3_6 = as.data.frame(loadings3_6)

# Save EFA
efa3_6_table <- kable(loadings3_6, format = "latex", booktabs = TRUE, 
                      caption = "EFA") %>%
  kable_styling(latex_options = c("hold_position"))
writeLines(efa3_6_table, "efa_table_all_w3.tex")



## Both waves

vars_b = c(
  "imp_nost_1_w2", "imp_nost_2_w2", "imp_nost_3_w2", "imp_nost_4_w2", "imp_nost_5_w2", 
  "imp_nost_6_w2", "imp_nost_7_w2", "immig_eco_w2", "immig_cult_w2", "pol_val_1_w2", 
  "pol_val_2_w2", "pol_val_3_w2", "pol_val_4_w2", "pol_val_5_w2", "pol_val_6_w2", 
  "pol_val_7_w2", "pol_val_8_w2", "sexism_1_w2", "sexism_2_w2", "sexism_3_w2", 
  "democ_frexp2_w2", "democ_frassc1_w2", "democ_unisuff1a_w2", 
  "democ_frelect2a_w2", "democ_decelec2a_w2", "democ_judcnstr2a_w2", 
  "democ_legcnstr1a_w2", "democ_eqlaw1_w2", "satis_dem_w2", "trust_pols_w2", 
  "govt_appr_w2", "left_right_w2", "EU_indep_w2", "pop_1_w2", "pop_2_w2", "pop_3_w2", 
  "pop_4_w2", "nat_chauv_1_w3", "nat_chauv_2_w3", "nat_chauv_3_w3", "nat_chauv_4_w3", 
  "nat_chauv_5_w3", "nat_chauv_6_w3"
)

scree(datm[, vars_b], factors = TRUE, pc = FALSE)
fa.parallel(datm[, vars_b], fa = "fa", cor = "mixed", sim = FALSE) # 8 factors
efa_b = fa(datm[, vars_b], fm = "minres", cor = "mixed", rotate = "promax", 
           max.iter = 200, nfactors = 8)
print(efa_b, cut = 0.3)

# Extract loadings and edit
loadings_b = as.matrix(efa_b$loadings)  
loadings_b = apply(loadings_b, 2, round, 2)
loadings_b = apply(loadings_b, 2, function(x) ifelse(abs(x) < 0.3, "", x))
loadings_b = as.data.frame(loadings_b)

# Save EFA
efa_b_table <- kable(loadings_b, format = "latex", booktabs = TRUE, 
                      caption = "EFA") %>%
  kable_styling(latex_options = c("hold_position"))
writeLines(efa_b_table, "efa_table_all_both.tex")

# Extract inter-factor correlations
cors_b = efa_b$Phi
cors_b = as.data.frame(round(cors_b, 2))
efa_cors_b_table <- kable(cors_b, format = "latex", booktabs = TRUE, 
                     caption = "EFA inter-factor correlations") %>%
  kable_styling(latex_options = c("hold_position"))
writeLines(efa_cors_b_table, "efa_cors_all_both.tex")

# Heatmap of loadings

plot_cols <- viridis::viridis(100, option = "mako")[100:1]
loadings_b <- as.matrix(efa_b$loadings) 
loadings_abs <- abs(loadings_b)

item_labels <- c(
  "imperial nost. 1", "imperial nost. 2", "imperial nost. 3", "imperial nost. 4", "imperial nost. 5", 
  "imperial nost. 6", "imperial nost. 7", "immigr. eco.", "immigr. cult.", "eco. values 1", 
  "eco. values 2", "eco. values 3", "eco. values 4", "auth. values 1", "auth. values 2", 
  "auth. values 3", "auth. values 4", "sexism 1", "sexism 2", "sexism 3", "sup. democ. 1", 
  "sup. democ. 2", "sup. democ. 3", "sup. democ. 4", "sup. democ. 5", "sup. democ. 6", 
  "sup. democ. 7", "sup. democ. 8", "satis. democ.", "trust MPs", "gov. approval", 
  "left-right", "EU indepen.", "populism 1", "populism 2", "populism 3", "populism 4", 
  "chauv. nation. 1", "nation. pride 1", "chauv. nation. 2", "chauv. nation. 3", 
  "nation. pride 2", "nation. pride 3"
)
item_labs_rev <- rev(item_labels)
col_labs = c("F1", "F2", "F3", "F4", "F5", "F6", "F7", "F8")

png("EFA_loadings_heatmap.png", width = 720, height = 1080)
par(mar = c(0.5, 4.5, 1, 0.5), cex = 2.4, tcl = -0.1)
image(x = 1:ncol(loadings_abs), y = 1:nrow(loadings_abs), 
  z = t(loadings_abs)[, nrow(loadings_abs):1], 
  col = plot_cols, axes = FALSE, xlab = "", ylab = "", zlim = c(0, 1)
)
axis(3, at = 1:ncol(loadings_b), labels = col_labs, tcl = 0, las = 1, 
     mgp = c(1, 0.1, 0), cex.axis = 0.7)
axis(2, at = 1:nrow(loadings_abs), labels = item_labs_rev, las = 2, 
     mgp = c(1, 0.4, 0), cex.axis = 0.6)
for (i in 1:nrow(loadings_abs)) {
  for (j in 1:ncol(loadings_abs)) {
    value <- loadings_b[i, j]
    label <- ifelse(abs(value) >= 0.30, sprintf("%.2f", value), "") # Suppress loadings < |0.3|
    text(j, nrow(loadings_abs) - i + 1, label, cex = 0.6, col = "grey90")
  }
}

box()
dev.off()




## save combined

write.csv(datm, "rude_w1-3_merg_scaled_rep.csv", row.names = FALSE)


