################################# # United or Divided in Diversity? # Dominik Schraff & Ronja Sczepanski # European Union Politics, 2021 # R Version 4.1.0 # Windows 64bit # set working directory to folder containing the datasets setwd("...") # Install packages library(foreign) library(tidyverse) library(ggplot2) library(rgdal) library(sp) library(rmapshaper) library(lme4) library(sjPlot) library(mediation) library(stargazer) library(htmlTable) # Load survey data # The survey data can be accessed through # Centerdata (https://www.centerdata.nl/en/liss-panel), where you have to # purchase a remote access dat <- read.spss("L_location_2p_incl_zipcode.sav", to.data.frame = T) head(dat) euscept <- dat[,c("prov","v9","v10","v11","zipcode","nomem_encr", "geslacht","leeftijd","oplcat","brutohh_f", "nettohh_f","belbezig","v18a","v18b","v18c","v6", "v14")] euscept$eusup <- as.character(euscept$v9) euscept$eusup[euscept$eusup=="0 De inpassing is te ver gegaan"] <- 0 euscept$eusup[euscept$eusup=="10 De inpassing moet nog verder gaan"] <- 10 euscept$eusup <- as.integer(euscept$eusup) euscept$eu.member.good <- as.character(euscept$v10) euscept$eu.member.good[euscept$eu.member.good=="een goede zaak"] <- 1 euscept$eu.member.good[euscept$eu.member.good=="een slechte zaak" | euscept$eu.member.good=="geen goede en geen slechte zaak"] <- 0 euscept$global.good <- as.character(euscept$v11) euscept$global.good[euscept$global.good=="een goede zaak"] <- 1 euscept$global.good[euscept$global.good=="geen goede en geen slechte zaak" | euscept$global.good=="een slechte zaak"] <- 0 euscept$reg.identity <- as.character(euscept$v18a) euscept$reg.identity[euscept$reg.identity=="helemaal niet verbonden"] <- 0 euscept$reg.identity[euscept$reg.identity=="heel erg verbonden"] <- 10 euscept$reg.identity <- as.integer(euscept$reg.identity) euscept$nat.identity <- as.character(euscept$v18b) euscept$nat.identity[euscept$nat.identity=="helemaal niet verbonden"] <- 0 euscept$nat.identity[euscept$nat.identity=="heel erg verbonden"] <- 10 euscept$nat.identity <- as.integer(euscept$nat.identity) euscept$eu.identity <- as.character(euscept$v18c) euscept$eu.identity[euscept$eu.identity=="helemaal niet verbonden"] <- 0 euscept$eu.identity[euscept$eu.identity=="heel erg verbonden"] <- 10 euscept$eu.identity <- as.integer(euscept$eu.identity) euscept$soctrst <- as.character(euscept$v6) euscept$soctrst[euscept$soctrst=="0 Je kunt niet voorzichtig genoeg zijn"] <- 0 euscept$soctrst[euscept$soctrst=="10 De meeste mensen zijn te vertrouwen"] <- 10 euscept$soctrst <- as.integer(euscept$soctrst) euscept$stfeco <- as.character(euscept$v14) euscept$stfeco[euscept$stfeco=="0 Heel ontevreden"] <- 0 euscept$stfeco[euscept$stfeco=="10 Heel tevreden"] <- 10 euscept$stfeco <- as.integer(euscept$stfeco) #Neighborhood data nn.dat <- readOGR("wijkenbuurten2019/buurt_2019_v1.shp") nn.comp <- ms_simplify(nn.dat, keep = 0.01, keep_shapes = T) summary(nn.comp) nn.comp$POSTCODE <- as.integer(nn.comp$POSTCODE) nn.comp$POSTCODE[nn.comp$POSTCODE==-99999999] <- NA nn.comp$P_WEST_AL <- as.numeric(nn.comp$P_WEST_AL) nn.comp$P_WEST_AL[nn.comp$P_WEST_AL==-99999999] <- NA nn.comp$P_N_W_AL <- as.numeric(nn.comp$P_N_W_AL) nn.comp$P_N_W_AL[nn.comp$P_N_W_AL==-99999999] <- NA nn.comp$AANT_INW <- as.numeric(nn.comp$AANT_INW) nn.comp$AANT_INW[nn.comp$AANT_INW==-99999999] <- NA nn.comp$P_MAROKKO <- as.numeric(nn.comp$P_MAROKKO) nn.comp$P_MAROKKO[nn.comp$P_MAROKKO==-99999999] <- NA nn.comp$P_TURKIJE <- as.numeric(nn.comp$P_TURKIJE) nn.comp$P_TURKIJE[nn.comp$P_TURKIJE==-99999999] <- NA nn.comp$P_SURINAM <- as.numeric(nn.comp$P_SURINAM) nn.comp$P_SURINAM[nn.comp$P_SURINAM==-99999999] <- NA #Merge nn.rawdat <- nn.comp@data nn.rawdat <- nn.rawdat[,c("POSTCODE", "P_WEST_AL", "P_N_W_AL", "AANT_INW", "P_MAROKKO","P_TURKIJE","P_SURINAM")] names(nn.rawdat) <- c("zipcode","P_WEST_AL","P_N_W_AL", "AANT_INW", "P_MAROKKO","P_TURKIJE","P_SURINAM") nn.zip <- nn.rawdat %>% group_by(zipcode) %>% summarise(western = sum(P_WEST_AL*AANT_INW, na.rm = T)/sum(AANT_INW, na.rm = T), non.western = sum(P_N_W_AL*AANT_INW, na.rm = T)/sum(AANT_INW, na.rm = T), turks = sum(P_TURKIJE*AANT_INW, na.rm = T)/sum(AANT_INW, na.rm = T), marok = sum(P_MAROKKO*AANT_INW, na.rm = T)/sum(AANT_INW, na.rm = T), suriname = sum(P_SURINAM*AANT_INW, na.rm = T)/sum(AANT_INW, na.rm = T), pop = sum(AANT_INW, na.rm = T)) eu.zip <- merge(euscept, nn.zip, by = "zipcode", all.x = T) #Merge 2010 local immigration data nn.dat10 <- readOGR("shape 2010 versie3.0/buurt_2010_v3.shp") nn.comp10 <- ms_simplify(nn.dat10, keep = 0.01, keep_shapes = T) summary(nn.comp10) nn.comp10$POSTCODE <- as.integer(nn.comp10$POSTCODE) nn.comp10$POSTCODE[nn.comp10$POSTCODE<0] <- NA nn.comp10$P_WEST_AL <- as.numeric(nn.comp10$P_WEST_AL) nn.comp10$P_WEST_AL[nn.comp10$P_WEST_AL<0] <- NA nn.comp10$P_N_W_AL <- as.numeric(nn.comp10$P_N_W_AL) nn.comp10$P_N_W_AL[nn.comp10$P_N_W_AL<0] <- NA nn.comp10$AANT_INW <- as.numeric(nn.comp10$AANT_INW) nn.comp10$AANT_INW[nn.comp10$AANT_INW<0] <- NA nn.comp10$P_MAROKKO <- as.numeric(nn.comp10$P_MAROKKO) nn.comp10$P_MAROKKO[nn.comp10$P_MAROKKO<0] <- NA nn.comp10$P_TURKIJE <- as.numeric(nn.comp10$P_TURKIJE) nn.comp10$P_TURKIJE[nn.comp10$P_TURKIJE<0] <- NA nn.comp10$P_SURINAM <- as.numeric(nn.comp10$P_SURINAM) nn.comp10$P_SURINAM[nn.comp10$P_SURINAM<0] <- NA nn.rawdat10 <- nn.comp10@data nn.rawdat10 <- nn.rawdat10[,c("POSTCODE", "P_WEST_AL", "P_N_W_AL", "AANT_INW", "P_MAROKKO","P_TURKIJE","P_SURINAM")] names(nn.rawdat10) <- c("zipcode","P_WEST_AL","P_N_W_AL", "AANT_INW", "P_MAROKKO","P_TURKIJE","P_SURINAM") nn.zip10 <- nn.rawdat10 %>% group_by(zipcode) %>% summarise(western.10 = sum(P_WEST_AL*AANT_INW, na.rm = T)/sum(AANT_INW, na.rm = T), non.western.10 = sum(P_N_W_AL*AANT_INW, na.rm = T)/sum(AANT_INW, na.rm = T), turks.10 = sum(P_TURKIJE*AANT_INW, na.rm = T)/sum(AANT_INW, na.rm = T), marok.10 = sum(P_MAROKKO*AANT_INW, na.rm = T)/sum(AANT_INW, na.rm = T), suriname.10 = sum(P_SURINAM*AANT_INW, na.rm = T)/sum(AANT_INW, na.rm = T), pop.10 = sum(AANT_INW, na.rm = T)) eu.zip <- merge(eu.zip, nn.zip10, by = "zipcode", all.x = T) #Merge local HH income hhinc.zip <- read.csv("hh_inc_postcode.csv", sep = ";") hhinc.zip <- hhinc.zip[,-2] names(hhinc.zip) <- c("zipcode","hhinc.zip") hhinc.zip$hhinc.zip <- gsub(" ", "", hhinc.zip$hhinc.zip) hhinc.zip$hhinc.zip[hhinc.zip$hhinc.zip=="."] <- NA hhinc.zip$hhinc.zip <- as.integer(hhinc.zip$hhinc.zip) eu.zip <- merge(eu.zip, hhinc.zip, by = "zipcode", all.x = T) #Zipcode distance to border library(sf) nn19 <- st_as_sf(nn.comp) nl.borders <- st_as_sf(st_union(nn19)) #plot(nl.borders) zip.borders <- nn19 %>% group_by(POSTCODE) %>% summarise() %>% ungroup %>% st_as_sf() #plot(zip.borders) nl.borders <- st_cast(nl.borders, "MULTILINESTRING") #takes a few minutes dist <- st_distance(nl.borders, zip.borders, by_element = T) zip.border.dist <- data.frame(list(zip.borders$POSTCODE), dist) names(zip.border.dist) <- c("zipcode", "border.dist") eu.zip <- merge(eu.zip, zip.border.dist, by = "zipcode", all.x = T) eu.zip$bdist.th <- eu.zip$border.dist/1000 #Immigrant background backdat <- read.spss("avars_201912_EN_1.0p.sav", to.data.frame = T) backdat <- backdat[,c("nomem_encr","herkomstgroep")] eu.zip <- merge(eu.zip, backdat, by = "nomem_encr", all.x = T) eu.zip$migback[eu.zip$herkomstgroep=="Dutch background"] <- 0 eu.zip$migback[eu.zip$herkomstgroep!="Dutch background"] <- 1 #Contiguity neighborhood immigration shares library(spdep) nbs <- poly2nb(zip.borders, row.names = as.character(zip.borders$POSTCODE)) summary(nbs) nn.dat <- enframe(nbs) nn.dat$POSTCODE <- zip.borders$POSTCODE nn.zip$rnames <- row.names(nn.zip) nmean.west <- NA nmean.nonwest <- NA for (i in nn.dat$value) { temp.sub <- subset(nn.zip, subset = rnames %in% i) nn.mean.western <- mean(temp.sub$western, na.rm = T) nn.mean.nonwestern <- mean(temp.sub$non.western, na.rm = T) nmean.west <- rbind(nmean.west, nn.mean.western) nmean.nonwest <- rbind(nmean.nonwest, nn.mean.nonwestern) } nn.averages <- data.frame(nmean.west,nmean.nonwest) nn.averages <- nn.averages[-1,] nn.averages$zipcode <- nn.zip$zipcode eu.zip <- merge(eu.zip, nn.averages, by = "zipcode", all.x = T) #Prep variables eu.zip$global.good <- as.integer(eu.zip$global.good) eu.zip$eu.member.good <- as.integer(eu.zip$eu.member.good) eu.zip$pop.th <- eu.zip$pop/1000 eu.zip$twodig <- substring(as.character(eu.zip$zipcode),1,2) eu.zip$excl.nonw <- eu.zip$non.western - eu.zip$western eu.zip$std.netinc <- as.numeric(scale(eu.zip$nettohh_f)) eu.zip$std.brutinc <- scale(eu.zip$brutohh_f) eu.zip$std.age <- as.numeric(scale(eu.zip$leeftijd)) eu.zip$std.age.sq <- as.numeric(scale(eu.zip$leeftijd*eu.zip$leeftijd)) eu.zip$excl.natid <- eu.zip$nat.identity - eu.zip$eu.identity eu.zip$std.exnatid <- scale(eu.zip$excl.natid) eu.zip$std.euid <- scale(eu.zip$eu.identity) eu.zip$std.natid <- scale(eu.zip$nat.identity) eu.zip$std.hhinczip <- as.numeric(scale(eu.zip$hhinc.zip)) eu.zip$turks_marok <- eu.zip$turks + eu.zip$marok eu.zip$ln.wester <- log(eu.zip$western + 1) eu.zip$ln.nonwester <- log(eu.zip$non.western + 1) eu.zip$ln.turks <- log(eu.zip$turks + 1) eu.zip$foreign <- eu.zip$western + eu.zip$non.western eu.zip$ln.foreign <- log(eu.zip$foreign + 1) eu.zip$high.edu <- 0 eu.zip$high.edu[eu.zip$oplcat=="wo" | eu.zip$oplcat=="hbo"] <- 1 eu.zip$west.ch <- eu.zip$western - eu.zip$western.10 eu.zip$non.west.ch <- eu.zip$non.western - eu.zip$non.western.10 eu.zip$turks_marok.ch <- eu.zip$turks_marok - (eu.zip$turks.10 + eu.zip$marok.10) eu.zip$ln.nn.wester <- log(eu.zip$nmean.west + 1) eu.zip$ln.nn.nonwester <- log(eu.zip$nmean.nonwest + 1) eu.zip$bdist.th <- as.numeric(eu.zip$bdist.th) ### #Analysis setwd("D:/backup Dominik/Extract_ID_Paper") stargazer(round(cor(eu.zip[,c("ln.wester","ln.nonwester","std.hhinczip","west.ch", "non.west.ch","bdist.th")], use = "complete.obs"), 3), type = "html", out = "corrtab_postcode.html") zip.n <- eu.zip %>% group_by(zipcode) %>% summarise(n()) summary(zip.n$`n()`, na.rm = T) ### #Exclusive national identity ### # Table 3 lm.ex1 <- lm(excl.natid ~ geslacht + std.age + std.age.sq + high.edu + std.netinc + migback + ln.wester + ln.nonwester + std.hhinczip, data = eu.zip) summary(lm.ex1) lm.ex2 <- lm(excl.natid ~ geslacht + std.age + std.age.sq + high.edu + std.netinc + migback + ln.nonwester*ln.wester + std.hhinczip, data = eu.zip) summary(lm.ex2) lm.ex3 <- lm(excl.natid ~ geslacht + std.age + std.age.sq + high.edu + std.netinc + migback + ln.nonwester*ln.wester + std.hhinczip + bdist.th, data = eu.zip) summary(lm.ex3) lm.ex4 <- lm(excl.natid ~ geslacht + std.age + std.age.sq + high.edu + std.netinc + migback + ln.nonwester*ln.wester + std.hhinczip + bdist.th + non.west.ch + west.ch, data = eu.zip) summary(lm.ex4) lm.ex5 <- lm(excl.natid ~ geslacht + std.age + std.age.sq + high.edu + std.netinc + migback + ln.nonwester*ln.wester + std.hhinczip + bdist.th + non.west.ch*west.ch, data = eu.zip) summary(lm.ex5) stargazer(lm.ex1,lm.ex2,lm.ex3,lm.ex4,lm.ex5, type = "html", out = "tab3.html") ### # Figure 3 tiff("inter1.tiff", units = "in", width=6, height=5, res=800, compression = "lzw") plot_model(lm.ex2, type = "int", axis.title = c("Log Non-Western", "Exclusive National Identity"), title = "") + aes(shape = group, color = group, linetype = group) + labs(color="Log Western", shape="Log Western", linetype="Log Western") dev.off() ### # Table 4 ### #Effect on national & EU ID lm.nat <- lm(nat.identity ~ geslacht + std.age + std.age.sq + high.edu + std.netinc + migback + ln.nonwester*ln.wester + std.hhinczip + bdist.th + non.west.ch + west.ch, data = eu.zip) summary(lm.nat) lm.eu <- lm(eu.identity ~ geslacht + std.age + std.age.sq + high.edu + std.netinc + migback + ln.nonwester*ln.wester + std.hhinczip + bdist.th + non.west.ch + west.ch, data = eu.zip) summary(lm.eu) stargazer(lm.nat, lm.eu, type = "html", out = "tab4.html") ### # Figure 4 tiff("inter_euid.tiff", units = "in", width=6, height=5, res=800, compression = "lzw") plot_model(lm.eu, type = "int", axis.title = c("Log Non-Western", "EU Identity"), title = "") + aes(shape = group, color = group, linetype = group) + labs(color="Log Western", shape="Log Western", linetype="Log Western") dev.off() ### # Table 5 cor(eu.zip[,c("excl.natid","eu.member.good", "global.good","eusup")], use = "complete.obs") eu.memb1 <- glm(eu.member.good ~ geslacht + std.age + std.age.sq + high.edu + std.netinc + migback + scale(stfeco) + scale(excl.natid), data = eu.zip, family = "binomial") summary(eu.memb1) glob1 <- glm(global.good ~ geslacht + std.age + std.age.sq + high.edu + std.netinc + migback + scale(stfeco) + scale(excl.natid), data = eu.zip, family = "binomial") summary(glob1) stargazer(eu.memb1, glob1, type = "html", out = "cosm_prefs.html") ### #Appendix 1 #Descriptive Stats stargazer(eu.zip[,c("geslacht","leeftijd","high.edu","nettohh_f","migback","stfeco", "excl.natid", "nat.identity","eu.identity", "bdist.th","non.west.ch","west.ch", "western", "non.western","pop.th","nmean.west","nmean.nonwest")], type = "html", out = "descriptives_rep.html") hist(eu.zip$western) hist(eu.zip$ln.wester) hist(eu.zip$non.western) hist(eu.zip$ln.nonwester) ### # Appendix 2 lm.ap1 <- lm(excl.natid ~ geslacht + std.age + std.age.sq + high.edu + std.netinc + migback + ln.nonwester*ln.wester + std.hhinczip + bdist.th + non.west.ch + west.ch + stfeco, data = eu.zip) summary(lm.ap1) lm.ap2 <- lm(excl.natid ~ geslacht + std.age + std.age.sq + high.edu + std.netinc + migback + ln.nonwester*ln.wester + std.hhinczip + bdist.th + non.west.ch + west.ch + stfeco + prov, data = eu.zip) summary(lm.ap2) lm.ap3 <- lm(excl.natid ~ geslacht + std.age + std.age.sq + high.edu + std.netinc + migback + ln.nonwester*ln.wester + std.hhinczip + bdist.th + non.west.ch + west.ch + ln.nn.wester + ln.nn.nonwester, data = eu.zip) summary(lm.ap3) lm.ap4 <- lm(excl.natid ~ geslacht + std.age + std.age.sq + high.edu + std.netinc + migback + ln.nonwester*ln.wester + std.hhinczip + bdist.th + non.west.ch + west.ch + pop.th, data = eu.zip) summary(lm.ap4) stargazer(lm.ap1, lm.ap2, lm.ap3, lm.ap4, type = "html", out = "append_tab2.html") ##### #Maps ### #Foreigner shares by zipcode area western.sdf <- fortify(nn.comp, region = "POSTCODE") western.sdf <- merge(western.sdf, nn.zip, by.x = "id", by.y = "zipcode") obs.zips <- unique(euscept$zipcode) western.sdf$obs.zips <- 0 western.sdf$id <- as.integer(western.sdf$id) western.sdf$obs.zips <- replace(western.sdf$obs.zips, western.sdf$id%in%obs.zips, 1) ### #Bivariate map library(biscale) library(sf) data2 <- na.omit(western.sdf[,c(1,2,3,7,8,9)]) data <- bi_class(data2, x = western, y = non.western) tiff("bivar_map.tiff", units = "in", width=6, height=6, res=800, compression = "lzw") ggplot() + geom_polygon(data = data, aes(x=long, y=lat, group=group, fill = bi_class), show.legend = FALSE) + bi_scale_fill(pal = "DkBlue", dim = 3) + bi_theme() + xlab("") + ylab("") dev.off() tiff("bivar_legend.tiff", units = "in", width=4, height=4, res=800, compression = "lzw") bi_legend(pal = "DkBlue", dim = 3, xlab = "Western", ylab = "Non-Western", size = 16) dev.off() ### #Binned exclusive nat. id map dv.sdf <- fortify(nn.comp, region = "POSTCODE") d.sub <- eu.zip[,c("zipcode","excl.natid","eusup")] d.sub <- d.sub %>% group_by(zipcode) %>% summarise(eusup = mean(eusup, na.rm = T), excl.natid = mean(excl.natid, na.rm = T)) dv.sdf.m <- merge(dv.sdf, d.sub, by.x = "id", by.y = "zipcode", all.x = T) dv.sdf.m <- dv.sdf.m[order(dv.sdf.m$order),] ggplot(dv.sdf.m, aes(long,lat,group=group)) + stat_summary_2d(aes(z = excl.natid), na.rm = T, bins = 10, alpha = 0.5, geom = "tile") + geom_path() + scale_fill_gradient2(low="#d7191c", mid = "#ffffbf", high = "#1a9641", midpoint = 2.750)