### This project replicates the analyses reported in Toshkov and Krouwel (2022) 'Beyond the U-Curve: Citizen Preferences on European Integration in Multidimensional Political Space' in European Union Politics ### Part I: This script prepares and analyzes the data from the 2019 EP elections in The Netherlands ### script developed by Dimiter Toshkov # Clean up VAA data and make variables ------------------------------------------------------------------------- ## This part of the script shows the preparation of the datafile that is shared from the original datafile that is proprietary # Load the data nl.all <- read_sav('./data in/nl ep 2019/2019 VAA data AN.sav') # Recode variables and filter nl.all<- nl.all %>% mutate(lr = popup_answer_1040 - 5, proeu = popup_answer_1041 - 5, progr = popup_answer_1042 - 5, ep2014.vote = recode (popup_answer_1039, '0' = 'D66', '1' = 'CDA', '2' = 'PVV', '3' = 'VVD', '4' = 'SP', '5' = 'PvdA', '6' = 'CU', '7' = 'GL', '8' = 'PvdD', '9' = 'Other', '10' = 'Other', '11' = NA_character_, '12' = 'No vote', '13' = NA_character_), ep2019.vote.intent = recode (popup_answer_1038, '0' = '50+', '1' = 'CDA', '2' = 'CU', '3' = 'D66', '4' = 'Other', '5' = 'DENK', '6' = 'FvD', '7' = 'GL', '8' = 'Other', '9' = 'PvdA', '10' = 'PvdD', '11' = 'PVV', '12' = 'SP', '13' = 'Other', '14' = 'Other', '15' = 'VVD', '16' = NA_character_, '17' = 'No vote', '18' = NA_character_), tk2017.vote = recode (popup_answer_1036, '0' = 'VVD', '1' = 'PVV', '2' = 'CDA', '3' = 'D66', '4' = 'GL', '5' = 'SP', '6' = 'PvdA', '7' = 'CU', '8' = 'PvdD', '9' = '50+', '10' = 'SGP', '11' = 'DENK', '12' = 'FvD', '13' = 'Other', '14' = NA_character_, '15' = 'No vote', '16' = NA_character_, '17' = NA_character_), vote2017 = recode (popup_answer_1036, '0' = 'VVD', '1' = 'PVV', '2' = 'CDA', '3' = 'D66', '4' = 'GL', '5' = 'SP', '6' = 'PvdA', '7' = 'CU', '8' = 'PvdD', '9' = '50+', '10' = 'SGP', '11' = 'DENK', '12' = 'FvD', '13' = NA_character_, '14' = NA_character_, '15' = NA_character_, '16' = NA_character_, '17' = NA_character_), ep.2019.vote.sure = recode (popup_answer_1037, '0' = 'Sure', '1' = 'Doubt', '2' = 'Dont_know', '3' = 'No vote'), age = as.numeric(as.character(ifelse (bg_question_445==94, NA_character_, 2019 - (2003 - bg_question_445)))), age_cat3 = ifelse(age<35, 'young', ifelse(age<65, 'midage', 'old')), sex = recode (bg_question_416, '0' = 'man', '1' = 'woman', '2' = NA_character_, '3' = NA_character_), edu = recode (bg_question_447, '0' = '4.high', '1' = '4.high', '2' = '3.midhigh', '3' = '2.midlow', '4' = '2.midlow', '5' = '1.lower', '6' = '1.lower', '7' = NA_character_), edu_cat3 = ifelse(edu=='4.high', 'high_edu', ifelse(edu=='3.midhigh', 'mid_edu', 'low_edu')), themes = case_when ((relevant_themes_130 == 1 ~ 'society'), (relevant_themes_162 == 1 ~ 'security'), (relevant_themes_268 == 1 ~ 'eu_integration'), (relevant_themes_269 == 1 ~ 'immigration'), (relevant_themes_271 == 1 ~ 'work_incomes'), (relevant_themes_273 == 1 ~ 'economy'), (relevant_themes_274 == 1 ~ 'environment')), sure = recode (popup_answer_1037, '0' = 1, '1' = 0, '2' = 0, '3' = 0, '4' = 0), type = 'voter', index = paste(sex, age_cat3, edu_cat3, sep='.'), pindex = paste(vote2017, sex, age_cat3, edu_cat3, sep='.'), eu.bad = statement_answer_4907 - 3, privacy.limit = statement_answer_4928 - 3, eu.money.climate = statement_answer_4923 - 3, many.immigrants = statement_answer_4938 - 3, euro.out = statement_answer_4908 - 3, eu.army = statement_answer_4926 -3, eu.money.culture = 6 - (statement_answer_4925) -3, weed.legal = statement_answer_4913 -3, eu.enforce = statement_answer_4932 -3, eu.money.aid = statement_answer_4941 -3, greener.energy = statement_answer_4929 -3, eu.nat.borders = statement_answer_4921 -3, eu.money.sustainable = statement_answer_4927 -3, soc.sec.stable = statement_answer_4934 -3, islam.bad = statement_answer_4920 -3, market.health = statement_answer_4918 -3, homo.adopt = statement_answer_4912 -3, eu.veto = statement_answer_4911 -3, fire.easy = statement_answer_4917 -3, eu.fp = statement_answer_4910 -3, state.econ.out = statement_answer_4916 -3, eu.lost.id = statement_answer_4922 -3, eu.fin.solidarity = statement_answer_4924 -3, limit.asylum = statement_answer_4935 -3, eu.unempl = statement_answer_4937 -3, redistribute = statement_answer_4915 -3, fighters.back = statement_answer_4939 -3, eu.out = statement_answer_4940 -3, eu.tax = statement_answer_4931 -3, eu.limit.work = statement_answer_4909 -3 ) %>% filter (country == 'Netherlands', amount_visits == 1, seconds_on_statements > 30*3, seconds_on_statements < 30*30) %>% select (lr:ncol(.)) # remove those who answered everything with Neutral nl.all$zeroes <- nl.all %>% select (eu.bad:eu.limit.work) %>% abs () %>% rowSums(na.rm=T) nl.all <- nl.all %>% filter (zeroes!=0) %>% select (-zeroes) # create deductively derived positions on scales nl.all <- nl.all %>% mutate (obj.lr = (- redistribute - soc.sec.stable + fire.easy + state.econ.out + market.health)/5, obj.cp = (- islam.bad - privacy.limit - limit.asylum - many.immigrants + homo.adopt + weed.legal + greener.energy + fighters.back)/8, obj.eu = (- eu.bad - euro.out - eu.veto - eu.lost.id - eu.out + eu.army + eu.tax + eu.enforce + eu.fp + eu.unempl)/10, obj.eu2 = (- eu.bad - euro.out - eu.veto - eu.lost.id - eu.out + eu.army + eu.tax + eu.enforce + eu.fp - eu.limit.work - eu.nat.borders + eu.money.climate + eu.money.sustainable + eu.unempl + eu.fin.solidarity + eu.money.aid + eu.money.culture)/17 ) nl.all$id = seq.int(nrow(nl.all)) #index for merging later # save save (nl.all, file = './data out/nl2019_replication.RData') # Load and subset the data (start replication here) ------------------------------------------------------------------------- load (file = './data out/nl2019_replication.RData') # Subset policy items only nl.s <- nl.all %>% select(id, eu.bad:eu.limit.work) %>% filter(complete.cases(.)) # Number of factors ---------------------------------------------------------------- # get polychoric covariance tk.poly<-polychoric(select(nl.s, -id)) # determine number of factors ev <- eigen(cor(select(nl.s, -id))) ev.poly <- eigen(tk.poly$rho) #with the polychronic correlations # number of factors ap <- parallel(subject=nrow(nl.s),var=ncol(select(nl.s, -id)),rep=10,quantile=.95) nS <- nScree(x=ev$values, aparallel=ap$eigen$qevpea) plotnScree(nS) # PCA analysis (Table A5) ---------------------------------------------------------------- pca.1 <- principal(select(nl.s, -id), nfactors=5, rotate='none', cor='poly', scores=TRUE) print(pca.1, digits=2, cut=.31, sort=T) pca.out<-as.data.frame(unclass(pca.1$loadings)) pca.out <- pca.out %>% mutate_at(vars(starts_with('PC')), ~replace(., abs(.) < 0.4, 0)) %>% mutate_at(vars(starts_with('PC')), funs(round(., 2))) %>% arrange (-abs(PC1),-abs(PC2), -abs(PC3), -abs(PC4), -abs(PC5)) %>% mutate_at(vars(starts_with('PC')), ~replace(., abs(.) == 0, '')) %>% mutate ( issue = rownames(.)) %>% relocate (issue) pca.out <- rbind (pca.out, c('Proportion Variance', paste0(round(prop.table(pca.1$values)[1:5]*100,0), '%'))) pca.out pca.nl.table = htmlTable::htmlTable(pca.out, rnames = FALSE, col.columns = c('white','lightgrey'), tfoot = 'Note: Loadings smaller than +/-0.41 are not printed', align = c('l','r','r','r','r','r'), align.header = c('l','r','r','r','r','r'), header = c('Policy issue','PC1','PC2','PC3','PC4','PC5'), caption = "Table X. PCA results from the Netherlands, 2019") print(pca.nl.table, type = "html", file = './tables/pca_nl_table.html') # Factor analysis I (full sample) (Table 1)---------------------------------------------------------------- fa.1 <- fa (select(nl.s, -id), nfactors=6, rotate='Promax', cor='poly', scores=TRUE) print(fa.1, digits=2, cut=.33, sort=T) fa.out<-as.data.frame(unclass(fa.1$loadings)) fa.out <- fa.out %>% mutate_at(vars(starts_with('MR')), ~replace(., abs(.) < 0.31, 0)) %>% mutate_at(vars(starts_with('MR')), funs(round(., 2))) %>% arrange (-abs(MR2),-abs(MR1), -abs(MR3), -abs(MR5), -abs(MR4), -abs(MR6)) %>% mutate_at(vars(starts_with('MR')), ~replace(., abs(.) == 0, '')) %>% mutate ( issue = rownames(.)) %>% relocate (issue) fa.out <- rbind (fa.out, c('Proportion Variance', paste0(round(fa.1$Vaccounted[2,]*100,0), '%'))) fa.out fa.nl.table = htmlTable::htmlTable(fa.out, rnames = FALSE, col.columns = c('white','lightgrey'), tfoot = 'Note: Loadings smaller than +/-0.31 are not printed', align = c('l','r','r','r','r','r','r'), align.header = c('l','r','r','r','r','r','r'), header = colnames(fa.out), caption = "Table X. Factor Analysis results from the Netherlands, 2019") print(fa.nl.table, type = "html", file = './tables/fa_nl_table.html') # Factor analysis II (restricted sample) (Table A6)---------------------------------------------------------------- # Run the same factor model but only on the subset who filled in the additional module nl.s2 <- nl.all %>% select (id, lr, progr, proeu, obj.lr, obj.cp, obj.eu, obj.eu2, eu.bad:eu.limit.work) %>% na.omit # factors with fa.poly and promax rotation fa.2 <- fa (select(nl.s2, eu.bad:eu.limit.work), nfactors=6, rotate='Promax', cor='poly', scores=TRUE) print(fa.2, digits=2, cut=.33, sort=T) fa2.out<-as.data.frame(unclass(fa.2$loadings)) fa2.out <- fa2.out %>% mutate_at(vars(starts_with('MR')), ~replace(., abs(.) < 0.31, 0)) %>% mutate_at(vars(starts_with('MR')), funs(round(., 2))) %>% arrange (-abs(MR4),-abs(MR1), -abs(MR3), -abs(MR5), -abs(MR2), -abs(MR6)) %>% mutate_at(vars(starts_with('MR')), ~replace(., abs(.) == 0, '')) %>% mutate ( issue = rownames(.)) %>% relocate (issue) fa2.out <- rbind (fa2.out, c('Proportion Variance', paste0(round(fa.2$Vaccounted[2,]*100,0), '%'))) fa2.out fa2.nl.table = htmlTable::htmlTable(fa2.out, rnames = FALSE, col.columns = c('white','lightgrey'), tfoot = 'Note: Loadings smaller than +/-0.31 are not printed', align = c('l','r','r','r','r','r','r'), align.header = c('l','r','r','r','r','r','r'), header = colnames(fa2.out), caption = "Table X. Factor Analysis results from the Netherlands, 2019 [interested sample]") print(fa2.nl.table, type = "html", file = './tables/fa2_nl_table.html') # Factor analysis III (only those who are sure) (Table A7---------------------------------------------------------------- # Run the same factor model but only on the subset who are sure who they will vote for nl.s4 <- nl.all %>% filter(sure == 1) %>% select (eu.bad:eu.limit.work) %>% na.omit # factors with fa.poly and promax rotation fa.4 <- fa (select(nl.s4, eu.bad:eu.limit.work), nfactors=6, rotate='Promax', cor='poly', scores=TRUE) print(fa.4, digits=2, cut=.33, sort=T) fa4.out<-as.data.frame(unclass(fa.4$loadings)) fa4.out <- fa4.out %>% mutate_at(vars(starts_with('MR')), ~replace(., abs(.) < 0.31, 0)) %>% mutate_at(vars(starts_with('MR')), funs(round(., 2))) %>% arrange (-abs(MR1),-abs(MR5), -abs(MR3), -abs(MR6), -abs(MR2), -abs(MR4)) %>% mutate_at(vars(starts_with('MR')), ~replace(., abs(.) == 0, '')) %>% mutate ( issue = rownames(.)) %>% relocate (issue) fa4.out <- rbind (fa4.out, c('Proportion Variance', paste0(round(fa.4$Vaccounted[2,]*100,0), '%'))) fa4.out fa4.nl.table = htmlTable::htmlTable(fa4.out, rnames = FALSE, col.columns = c('white','lightgrey'), tfoot = 'Note: Loadings smaller than +/-0.31 are not printed', align = c('l','r','r','r','r','r','r'), align.header = c('l','r','r','r','r','r','r'), header = colnames(fa4.out), caption = "Table X. Factor Analysis results from the Netherlands, 2019 [certain sample]") print(fa4.nl.table, type = "html", file = './tables/fa4_nl_table.html') # Assign factor loadings to cases ---------------------------------------------------------------- nl.s <- bind_cols(nl.s, data.frame(pca.1$scores), data.frame(fa.1$scores)) nl.all <- left_join (nl.all, select(nl.s, id, PC1:MR6), by='id') # Subset ------------------------------------------------- nl.s <- nl.all %>% filter (type=='voter') %>% select ('lr','progr','proeu') %>% filter(complete.cases(.)) # Self-placement analysis: 3D plots ------------------------------------------------- #prepare data elevation.df = data.frame(x = nl.s$lr, y = nl.s$progr, z = nl.s$proeu) elevation.loess = loess(z ~ x*y, data = elevation.df, degree=2, span=0.4) xmin<- -5 xmax<- 5 ymin<- -5 ymax<- 5 elevation.fit = expand.grid(list(x = seq(xmin, xmax, 0.1), y = seq(ymin, ymax, 0.1))) z = predict(elevation.loess, newdata = elevation.fit) elevation.fit$Height = as.numeric(z) #prepare matrix with x and y and rows/cols and z as the cell values u1<-unique(elevation.fit$x) u2<-unique(elevation.fit$y) temp.data<-matrix(NA, nrow=length(u1), ncol=length(u2)) for (i in 1:length(u1)){ for (j in 1:length(u2)){ temp.data[i,j]<-elevation.fit$Height[elevation.fit$x==u1[i] & elevation.fit$y==u2[j]] } } ix <- seq(1, nrow(temp.data), length.out = 11) iy <- seq(1, ncol(temp.data), length.out = 11) # 3D in black and white (Figure 2, left panel) png('./figures/nl 2019 eu3.png', width=900*2, height=700*2, res=96) par(mfrow = c(1, 1), mar = c(0, 0, 0, 0)) persp3D(z = temp.data, theta = 45, phi = 10, facets = T, curtain = TRUE, shade=0.8, col = ramp.col(c("white", "black")), border = "grey", colkey = FALSE, xlab="Left-Right", ylab="Conservative-Progressive", zlab="Anti-Pro EU") dev.off() # 2D plot with the plot3d package (Figure A1 top) par(mfrow = c(1, 1), mar = c(6, 8, 6, 4)) #first in 2D. quite nice actually png('./figures/nl 2019 eu.png', width=900*2, height=700*2, res=96) image2D(temp.data, clab = "Anti-Pro EU",xlab="Economic Left/Right", ylab="Conservative/Progressive", cex.lab=1.5, cex.axis=1, cex.main=1.5, cex.sub=2, pt.cex = 2) dev.off() # 4 3D plots in one figure (Figure A1 bottom) png('./figures/nl 2019 eu6.png', width=900*2, height=700*2, res=96) par(mfrow = c(2, 2), mar = c(3, 3, 3, 3)) persp3D(z = temp.data, theta = 55, phi = 20, facets = T, curtain = TRUE, shade=0.8, colkey = FALSE, xlab="Left-Right", ylab="Conservative-Progressive", zlab="Anti-Pro EU", box=T, bty='b2',ticktype = "detailed") persp3D(z = temp.data, theta = 45, phi = 10, facets = T, curtain = TRUE, shade=0.8, col = ramp.col(c("white", "black")), border = "grey", colkey = FALSE, xlab="Left-Right", ylab="Conservative-Progressive", zlab="Anti-Pro EU") ribbon3D(z = temp.data[ix, ], along = "y", cex.axis = 0.8, curtain = T, space = 0.7, theta = 25, phi = 5, shade=0.2, colkey = FALSE, xlab="Left-Right", ylab="Conservative-Progressive", zlab="Anti-Pro EU", box=T, bty='b2', ticktype = "detailed") hist3D(z = temp.data[ix,iy], theta = 55, phi = 20, shade=0.2, colkey = FALSE, box=T, bty='b2',ticktype = "detailed", xlab="Left-Right", ylab="Conservative-Progressive", zlab="Anti-Pro EU") dev.off() # with the plotly library plot_ly(x = u1, y=u2, z=temp.data) %>% layout( title = "Citizen positions on political scales, the Netherlands, 2019", scene = list( xaxis = list(title = "Economic Left-Right"), yaxis = list(title = "Conservative-Progressive"), zaxis = list(title = "Anti-Pro EU") )) %>% add_surface(contours = list( z = list( show=TRUE, usecolormap=TRUE, highlightcolor="#ff0000", project=list(z=TRUE) ) )) # now with the rayshader package (Figure A2) temp.data2 = temp.data*300 raymat = ray_shade(temp.data2) ambmat = ambient_shade(temp.data2) png('./figures/raysahder1.png', width=1000, height=800, res=96) temp.data2 %>% sphere_shade(texture = "desert") %>% add_water(detect_water(temp.data2), color="desert") %>% add_shadow(ray_shade(temp.data2,zscale=3,maxsearch = 300),0.5) %>% add_shadow(ambmat,0.5) %>% plot_3d(temp.data2,zscale=35,fov=0,theta=-197,zoom=0.75,phi=10, windowsize = c(1000,800), water=TRUE, waterdepth = 0, wateralpha = 0.5,watercolor = "lightblue", waterlinecolor = "white",waterlinealpha = 0.5) render_label(temp.data2, x=1,y=1, z=1100, zscale=35,relativez=FALSE, text = "Left, Conservative",textsize = 1.7,linewidth = 2,freetype = FALSE, linecolor = 'darkred', textcolor = "darkred") render_label(temp.data2, x=101,y=1, z=1200, zscale=35,relativez=FALSE, linecolor = 'darkblue', textcolor = "darkblue", text = "Right, Conservative",textsize = 1.7,linewidth = 2,freetype = FALSE) render_label(temp.data2, x=1,y=101, z=1400, zscale=35,relativez=FALSE, linecolor = 'darkred', textcolor = "darkred", text = "Left, Progressive",textsize = 1.7,linewidth = 2,freetype = FALSE) render_label(temp.data2, x=101,y=101, z=1200, zscale=35, relativez=FALSE, linecolor = 'darkblue',textcolor = "darkblue", text = "Right, Progresive",textsize = 1.7,linewidth = 2,freetype = FALSE) render_label(temp.data2, x=26,y=101, z=1400, zscale=35, relativez=FALSE, dashed= TRUE, linecolor = 'black', textcolor = "black", text = "Peak of EU support",textsize = 1.5,linewidth = 2,freetype = FALSE, offset=0) render_label(temp.data2, x=85,y=15, z=300, zscale=35, relativez=FALSE, dashed= TRUE, linecolor = 'black', textcolor = "black", text = "Lake of Eurosceptics",textsize = 1.5, linewidth = 2,freetype = FALSE) render_snapshot() dev.off() # Additional 2D and 3D density plots -------------------------------------------------- # a combination of 2d plots (not included) p1 <- ggplot(nl.s, aes(x=lr, y=progr) ) + stat_density_2d(aes(fill = ..level..), geom = "polygon", colour="white") + theme_bw() + xlim(-5,5) + ylim(-5,5) + xlab('Left-Right') + ylab('Conservative-Progressive') + theme(legend.position="none") p2 <- ggplot(nl.s, aes(x=lr, y=proeu) ) + stat_density_2d(aes(fill = ..level..), geom = "polygon", colour="white") + theme_bw() + xlim(-5,5) + ylim(-5,5) + xlab('Left-Right') + ylab('Anti-Pro EU') + theme(legend.position="none") p3 <- ggplot(nl.s, aes(x=progr, y=proeu) ) + stat_density_2d(aes(fill = ..level..), geom = "polygon", colour="white") + theme_bw() + xlim(-5,5) + ylim(-5,5) + ylab('Anti-Pro EU') + xlab('Conservative-Progressive') + theme(legend.position="none") (p1 + plot_spacer())/(p2 +p3) # 3D histogram with the density mapped as color (Figure A3) smooth <- MASS::kde2d(nl.s$lr, nl.s$progr, lims = c(-5, 5, -5, 5), n = 11)$z par(mfrow = c(2, 2), mar = c(3, 3, 3, 3)) color = rev(rainbow(10, start = 0/6, end = 4/6)) hist3D(z = temp.data[ix,iy], theta = 35, phi = 20, shade=0.2, colvar = smooth, col=color, colkey = FALSE, box=T, bty='b2',ticktype = "detailed",xlab="Left-Right", ylab="Conservative-Progressive", zlab="Anti-Pro EU") hist3D(z = temp.data[ix,iy], theta = 35 + 180, phi = 20, shade=0.2, colvar = smooth, col=color, box=T, bty='b2',ticktype = "detailed",xlab="Left-Right", ylab="Conservative-Progressive", zlab="Anti-Pro EU") hist3D(colvar = temp.data[ix,iy], theta = 35, phi = 20, shade=0.2, z = smooth, col=color, colkey = FALSE, box=T, bty='b2',ticktype = "detailed",xlab="Left-Right", ylab="Conservative-Progressive", zlab="Density") hist3D(colvar = temp.data[ix,iy], theta = 35 + 180, phi = 20, shade=0.2, z = smooth, col=color, box=T, bty='b2',ticktype = "detailed",xlab="Left-Right", ylab="Conservative-Progressive", zlab="Density") # Subjective placement and objective positions: correlation plots (Figure 1) ---------------------------------- ggpairs( nl.all [, c('obj.lr', 'obj.cp', 'obj.eu', 'lr','progr','proeu')], upper = list(continuous = "cor"), #upper = list(continuous = "density"), lower = list(continuous = "trends"), diag = list(continuous = "barDiag"), axisLabels = 'show', columnLabels = c('Obj. Left-Right','Obj. Cons.-Progr.','Obj. Anti-Pro EU', 'Left-Right','Conservative-Progressive','Anti-Pro EU')) + theme_bw() # Objective positions: 3D plots ---------------------------------- nl.s2 <- nl.all %>% filter (type=='voter') %>% select ('obj.lr','obj.cp','obj.eu') %>% filter(complete.cases(.)) #prepare data elevation.df = data.frame(x = nl.s2$obj.lr, y = nl.s2$obj.cp, z = nl.s2$obj.eu) elevation.loess = loess(z ~ x*y, data = elevation.df) xmin<- -2 xmax<- 2 ymin<- -2 ymax<- 2 elevation.fit = expand.grid(list(x = seq(xmin, xmax, 0.1), y = seq(ymin, ymax, 0.1))) z = predict(elevation.loess, newdata = elevation.fit) elevation.fit$Height = as.numeric(z) #prepare matrix with x and y and rows/cols and z as the cell values u1<-unique(elevation.fit$x) u2<-unique(elevation.fit$y) temp.data<-matrix(NA, nrow=length(u1), ncol=length(u2)) for (i in 1:length(u1)){ for (j in 1:length(u2)){ temp.data[i,j]<-elevation.fit$Height[elevation.fit$x==u1[i] & elevation.fit$y==u2[j]] } } ix <- seq(1, nrow(temp.data), length.out = 11) iy <- seq(1, ncol(temp.data), length.out = 11) smooth <- MASS::kde2d(nl.s2$obj.lr, nl.s2$obj.cp, lims = c(-2, 2, -2, 2), n = 11)$z # 4 3D plots (Figure A4) par(mfrow = c(2, 2), mar = c(1, 1, 1, 1)) color = rev(rainbow(10, start = 0/6, end = 4/6)) hist3D(z = temp.data[ix,iy], theta = 35, phi = 20, shade=0.2, colvar = smooth, col=color, colkey = FALSE, clab='Density', box=T, bty='b2',ticktype = "detailed",xlab="Left-Right", ylab="Conservative-Progressive", zlab="Anti-Pro EU") hist3D(z = temp.data[ix,iy], theta = 35 - 90, phi = 20, shade=0.2, colvar = smooth, col=color, clab='Density', box=T, bty='b2',ticktype = "detailed",xlab="Left-Right", ylab="Conservative-Progressive", zlab="Anti-Pro EU") hist3D(colvar = temp.data[ix,iy], theta = 35, phi = 20, shade=0.2, z = smooth, col=color, colkey = FALSE, clab='Anti-Pro EU', box=T, bty='b2',ticktype = "detailed",xlab="Left-Right", ylab="Conservative-Progressive", zlab="Density") hist3D(colvar = temp.data[ix,iy], theta = 35 - 90, phi = 20, shade=0.2, z = smooth, col=color, clab='Anti-Pro EU', box=T, bty='b2',ticktype = "detailed",xlab="Left-Right", ylab="Conservative-Progressive", zlab="Density") # 3D plots in black and white (Figure 2, right panel) png('./figures/nl 2019 obj eu3.png', width=900*2, height=700*2, res=96) par(mfrow = c(1, 1), mar = c(0, 0, 0, 0)) persp3D(z = temp.data, theta = 45, phi = 10, facets = T, curtain = TRUE, shade=0.8, col = ramp.col(c("white", "black")), border = "grey", colkey = FALSE, xlab="Obj. Left-Right", ylab="Obj. Conservative-Progressive", zlab="Obj. Anti-Pro EU") dev.off() # Preparing the LISS data ---------------------------------- # These files are downloaded from the LISS data archive. The files are public but access requires registration liss_pol <- read_sav('./data in/nl ep 2019/LISS/cv20l_EN_1.0p.sav') liss_gen <- read_sav('./data in/nl ep 2019/LISS/avars_201912_EN_1.0p.sav') liss <- left_join (liss_pol, liss_gen, by='nomem_encr') liss <- liss %>% filter (leeftijd > 17) %>% mutate (age = leeftijd, age_cat = lftdcat, age_cat3 = ifelse(age<35, 'young', ifelse(age<65, 'midage', 'old')), edu = dplyr::recode (as.character(oplzon), '1'='1.lower','2'='2.midlow','3'='2.midlow','4'='3.midhigh', '5'='3.midhigh', '6'='4.high', .default=NA_character_), edu_cat3 = ifelse(edu=='4.high', 'high_edu', ifelse(edu=='3.midhigh', 'mid_edu', 'low_edu')), sex = ifelse(geslacht==2, 'woman', 'man'), vote2017 = dplyr::recode(as.character(cv20l307), '1' = 'VVD', '2' ='PVV', '3'='CDA', '4'='D66', '5'='GL', '6'='SP','7'='PvdA', '8'='CU', '9'='PvdD', '10'='50+', '11'='SGP', '12'='DENK', '13'='FvD',.default = NA_character_), index = paste(sex, age_cat3, edu_cat3, sep='.'), pindex = paste(vote2017, sex, age_cat3, edu_cat3, sep='.') ) %>% select (age, age_cat, age_cat3, edu, edu_cat3, sex, vote2017, index) %>% filter (is.na(edu_cat3)==F) # save the weigths from LISS save (liss, file = './data out/liss.RData') # Factor analysis V (weighted based on demographics and party)---------------------------------------------------------------- # load the LISS weights load (file = './data out/liss.RData') # subset the VAA data nl.nona <- nl.all %>% select (index, pindex, vote2017, sex, age_cat3, edu_cat3, eu.bad:eu.limit.work) %>% na.omit # get the relative weights w<-round(prop.table(table(liss$index, liss$vote2017), 2), 2)/ round(prop.table(table(nl.nona$index, nl.nona$vote2017), 2), 2) w <- melt(data.frame(w)) %>% mutate (pindex = paste(Var2, Var1, sep='.')) %>% mutate_if(is.numeric, list(~na_if(., Inf))) %>% select (pindex, value) # attach back to data nl.nona<- left_join (nl.nona, w, by='pindex') # Factor analysis IV (weighted based on party choice and demographics) (Table A8)---------------------------------------------------------------- nl.s5 <- nl.nona %>% select (value, eu.bad:eu.limit.work) %>% na.omit # factors with fa.poly and promax rotation fa.5 <- fa (select(nl.s5, -value), nfactors=6, rotate='Promax', cor='poly', scores=TRUE, weight=nl.s5$value) print(fa.5, digits=2, cut=.33, sort=T) fa5.out<-as.data.frame(unclass(fa.5$loadings)) fa5.out <- fa5.out %>% mutate_at(vars(starts_with('MR')), ~replace(., abs(.) < 0.31, 0)) %>% mutate_at(vars(starts_with('MR')), funs(round(., 2))) %>% arrange (-abs(MR4),-abs(MR1), -abs(MR3), -abs(MR5), -abs(MR2), -abs(MR6)) %>% mutate_at(vars(starts_with('MR')), ~replace(., abs(.) == 0, '')) %>% mutate ( issue = rownames(.)) %>% relocate (issue) fa5.out <- rbind (fa5.out, c('Proportion Variance', paste0(round(fa.5$Vaccounted[2,]*100,0), '%'))) fa5.out fa5.nl.table = htmlTable::htmlTable(fa5.out, rnames = FALSE, col.columns = c('white','lightgrey'), tfoot = 'Note: Loadings smaller than +/-0.31 are not printed', align = c('l','r','r','r','r','r','r'), align.header = c('l','r','r','r','r','r','r'), header = colnames(fa4.out), caption = "Table X. Factor Analysis results from the Netherlands, 2019 [weighted demo+party]") print(fa5.nl.table, type = "html", file = './tables/fa5_nl_table.html') # Factor analysis V (weighted based on demographics only) (Table A9)---------------------------------------------------------------- nl.nona <- nl.all %>% select (index, pindex, vote2017, sex, age_cat3, edu_cat3, eu.bad:eu.limit.work) %>% na.omit w<-round(prop.table(table(liss$index)), 2) / round(prop.table(table(nl.nona$index)), 2) w <- melt(data.frame(w)) %>% mutate (index = paste(Var1, sep='')) %>% mutate_if(is.numeric, list(~na_if(., Inf))) %>% select (index, value) # attach back to data nl.nona<- left_join (nl.nona, w, by='index') # Run the same factor model but only on the subset who are sure who they will vote for nl.s6 <- nl.nona %>% select (value, eu.bad:eu.limit.work) %>% na.omit # factors with fa.poly and promax rotation fa.6 <- fa (select(nl.s6, -value), nfactors=6, rotate='Promax', cor='poly', scores=TRUE, weight=nl.s6$value) print(fa.6, digits=2, cut=.33, sort=T) fa6.out<-as.data.frame(unclass(fa.6$loadings)) fa6.out <- fa6.out %>% mutate_at(vars(starts_with('MR')), ~replace(., abs(.) < 0.31, 0)) %>% mutate_at(vars(starts_with('MR')), funs(round(., 2))) %>% arrange (-abs(MR4),-abs(MR1), -abs(MR3), -abs(MR5), -abs(MR2), -abs(MR6)) %>% mutate_at(vars(starts_with('MR')), ~replace(., abs(.) == 0, '')) %>% mutate ( issue = rownames(.)) %>% relocate (issue) fa6.out <- rbind (fa6.out, c('Proportion Variance', paste0(round(fa.6$Vaccounted[2,]*100,0), '%'))) fa6.out fa6.nl.table = htmlTable::htmlTable(fa6.out, rnames = FALSE, col.columns = c('white','lightgrey'), tfoot = 'Note: Loadings smaller than +/-0.31 are not printed', align = c('l','r','r','r','r','r','r'), align.header = c('l','r','r','r','r','r','r'), header = colnames(fa6.out), caption = "Table X. Factor Analysis results from the Netherlands, 2019 [weighted demo only]") print(fa6.nl.table, type = "html", file = './tables/fa6_nl_table.html') ### THE END