### R-code for replication
# Article: Comparing Refugee Dispersal Policies: Varieties of Responsibility-Sharing in Europe
# Authors: Walter Bartl & Philipp Lutz
# Forthcoming in 'Governance'

## load packages
library(foreign)
library(readxl)
library(car)
library(ggplot2)
library(ggpubr)
library(hrbrthemes)
library(kableExtra)
library(dplyr)
library(scales)
library(RColorBrewer)
library(eurostat)
library(sf)
library(countrycode)
library(coefplot)
library(readstata13)

## load data
dispersal <- read.dta("data_replication.dta")

## build dispersal indeces

# equal weights index
dispersal$index <- dispersal$bindingness  + dispersal$criteria + dispersal$governance + dispersal$sanction + dispersal$monitoring
dispersal <- dispersal[order(-dispersal$index),]
dispersal <- subset(dispersal, dispersal$country!="country")

# Create the refugee restrictiveness index
dispersal$refugee_restrictiveness <- with(dispersal,
                                          0.6 * bindingness +
                                            0.4 * sanction 
)

# Create the subnational authority restrictiveness index
dispersal$subnational_restrictiveness <- with(dispersal,
                                              0.5 * governance +
                                                0.3 * criteria +
                                                0.2 * monitoring
)

cor(dispersal$refugee_restrictiveness, dispersal$subnational_restrictiveness) #

# 2. Scatterplot with regression line
jpeg("figureA3.jpeg", width = 6, height = 6, units = "in", res=800)
par(mar=c(4,4,1,1), mgp=c(2,.7,0), tck=-.01, las=1, bty="l")
set.seed(123)
plot(jitter(dispersal$refugee_restrictiveness), 
     jitter(dispersal$subnational_restrictiveness),
     xlab = "Restrictiveness toward refugees",
     ylab = "Restrictiveness toward subnational authorities",
     pch = 19, 
     col = rgb(0.5, 0.5, 0.5, alpha = 0.6)) 
abline(lm(subnational_restrictiveness ~ refugee_restrictiveness, data = dispersal), 
       col = "red", lwd = 2)
dev.off()

# Optional: Create an overall restrictiveness index (average of the two)
dispersal$overall_restrictiveness <- rowMeans(dispersal[, c("refugee_restrictiveness", "subnational_restrictiveness")])

cor(dispersal$index, dispersal$overall_restrictiveness, method="pearson") # 
cor(dispersal$index, dispersal$overall_restrictiveness, method="spearman") # 

# Sort data in increasing order (reversed)
sorted_order <- order(dispersal$overall_restrictiveness, decreasing = FALSE)
sorted_values <- dispersal$overall_restrictiveness[sorted_order]
sorted_countries <- dispersal$country[sorted_order]

# create a data.frame for the barplot
df_bar <- data.frame(
  country = factor(sorted_countries, levels = rev(sorted_countries)), # preserve order
  value   = sorted_values
)

library(ggplot2)
library(hrbrthemes)
library(forcats)
p1 <- ggplot(df_bar, aes(x = value, y = fct_rev(country), fill = value)) +
  geom_col() +
  scale_x_continuous(limits = c(0, 2), expand = c(0,0)) +
  scale_fill_gradient2(
    low = "yellow", mid = "orange", high = "black",
    midpoint = 1.0, space = "Lab",
    name = "Restrictiveness",
    trans = mild_high_trans, limits = c(0,2),
    breaks = legend_breaks, labels = legend_labels
  ) +
  labs(title = "A. Country ranking", x = "Restrictiveness of refugee dispersal policies", y = "") +
  theme_ipsum(base_size = 14) +
  theme(
    legend.position = "none",
    plot.title = element_text(
      family = "Roboto Condensed", # match map
      face = "bold",
      size = 14,
      hjust = 0
    ),
    axis.text.y = element_text(size = 8),
    panel.grid.major.y = element_blank(),
    panel.grid.minor.x = element_blank(),
    plot.margin = unit(rep(5,4), "mm")
  )

# Figure 2: Distribution of RDPs by policy attributes
jpeg("figure2.jpeg", width = 8, height = 2, units = "in", res=800)
par (mar=c(4,4,1,1), mgp=c(2,.7,0), tck=-.01, las=1, bty="l")   # 
par(mfrow = c(1,5))
barplot(table(dispersal$bindingness), border=F, main="Bindingness", ylim=c(0,29))
barplot(table(dispersal$sanction), border=F, main="Sanctions", ylim=c(0,29))
barplot(table(dispersal$governance), border=F, main="Governance", ylim=c(0,29))
barplot(table(dispersal$criteria), border=F, main="Criteria", ylim=c(0,29))
barplot(table(dispersal$monitoring), border=F, main="Monitoring", ylim=c(0,29))
dev.off()

## plotting country-profiles
dispersaldata <- dispersal[,1:6]
library(reshape)
dispersaldata_df <- as.data.frame(dispersaldata)
dispersaldata1 <- melt(dispersaldata_df, id=c("country"))
names(dispersaldata1)[2] <- "criteria"

# Order countries by decreasing overall_restrictiveness
dispersal_sorted <- dispersal[order(-dispersal$overall_restrictiveness), ]
countries <- dispersal_sorted$country

# Vector of countries you want to plot
countries <- c("Germany", "Austria", "Switzerland", "Sweden", "Finland", "Norway",
               "Denmark", "Netherlands", "Belgium", "Luxembourg", "France", "Spain",
               "Portugal", "United Kingdom", "Ireland", "Italy", "Malta", "Greece",
               "Cyprus", "Turkey", "Bulgaria", "Romania", "Croatia", "Czech Republic",
               "Estonia", "Hungary", "Latvia", "Lithuania", "Poland", "Serbia",
               "Slovakia", "Slovenia")

make_country_plot <- function(cntry) {
  data_subset <- subset(dispersaldata1, country == cntry)
  ggplot(data_subset, aes(country, value, fill = criteria)) +
    geom_bar(stat = "identity", position = "dodge") +
    theme_ipsum() +
    labs(title = cntry) +
    ylab("Dispersal restrictiveness") +
    theme(axis.title.x = element_blank(),
          axis.text.x = element_blank(),
          axis.ticks.x = element_blank()) +
    guides(fill = guide_legend(title = "Indicator"))
}

plots <- lapply(countries, make_country_plot)
names(plots) <- countries

jpeg("figureA2.jpeg", width = 14, height = 21, units = "in", res = 800)
ggarrange(plotlist = plots, ncol = 4, nrow = 8)
dev.off()

# correlation matrix
dispersal_attributes <- dispersal[,c(2,3,4,5,6)]
library(corrplot)
M = cor(dispersal_attributes, method = "spearman")
par (mar=c(4,4,1,1), mgp=c(2,.7,0), tck=-.01, las=1, bty="l")   # 
par(mfrow = c(1,1))
testRes = cor.mtest(dispersal_attributes, conf.level = 0.95)

jpeg("figure3.jpeg", width = 8, height = 6.5, units = "in", res=800)
my_col <- colorRampPalette(c("blue", "white", "darkorange"))(200)
corrplot(
  M,
  p.mat = testRes$p,
  sig.level = 0.1,
  addrect = 2,
  method = "shade",
  tl.col = "black",
  col = my_col,
  col.lim = c(-1, 1))
dev.off()

# Determine Number of Factors to Extract
library(nFactors)
ev <- eigen(cor(dispersal_attributes)) # get eigenvalues
ap <- parallel(subject=nrow(dispersal_attributes),var=ncol(dispersal_attributes),
               rep=100,cent=.05)
nS <- nScree(x=ev$values, aparallel=ap$eigen$qevpea)

# Item response theory (IRT)
library(mirt)
library(ggmirt)
dispersal_irt <- as.matrix(dispersal_attributes)
mirt1 <- mirt(dispersal_irt, model=1, itemtype="gpcm", technical=list(NCYCLES=1000))
params <- coef(mirt1, IRTpars=T, simplify=T)
dispersal$mirt_index <- fscores(mirt1)

cor(dispersal$mirt_index, dispersal$overall_restrictiveness, method="pearson") # 
cor(dispersal$mirt_index, dispersal$overall_restrictiveness, method="spearman") # 

jpeg("figureA4.jpeg", width = 12, height = 6, units = "in", res=800)
par(mar=c(4,4,1,1), mgp=c(2,.7,0), tck=-.01, las=1, bty="l")
par(mfrow=c(1, 2))
set.seed(123) 
plot(jitter(dispersal$index), jitter(dispersal$overall_restrictiveness),
     xlab = "Equal weights index",
     ylab = "Theoretical weights index",
     pch = 19,
     col = rgb(0.5, 0.5, 0.5, alpha = 0.6))  
abline(lm(overall_restrictiveness ~ index, data = dispersal), col = "red", lwd = 2)
correlation <- cor(dispersal$index, dispersal$overall_restrictiveness, use = "complete.obs")
text(x = max(dispersal$index, na.rm = TRUE) * 0.05, 
     y = max(dispersal$overall_restrictiveness, na.rm = TRUE) * 0.95,  
     labels = paste("r =", round(correlation, 3)),
     pos = 4, 
     cex = 1.1, 
     font = 2)   

set.seed(123)  #
plot(jitter(dispersal$mirt_index), jitter(dispersal$overall_restrictiveness),
     xlab = "IRT weights index",
     ylab = "Theoretical weights index",
     pch = 19,
     col = rgb(0.5, 0.5, 0.5, alpha = 0.6)) 
abline(lm(overall_restrictiveness ~ mirt_index, data = dispersal), col = "red", lwd = 2)
correlation <- cor(dispersal$mirt_index, dispersal$overall_restrictiveness, use = "complete.obs")
text(x = max(dispersal$mirt_index, na.rm = TRUE) * 0.05, 
     y = max(dispersal$overall_restrictiveness, na.rm = TRUE) * 0.95, 
     labels = paste("r =", round(correlation, 3)),
     pos = 4,  
     cex = 1.1,  
     font = 2)
dev.off()

## Principal Component Analysis (for categorical data CATPCA)

library(Gifi)
## ordinal PCA
fitord <- princals(as.matrix(dispersal_attributes), levels = "ordinal")  
summary(fitord)

jpeg("figureA1.jpeg", width = 7, height = 7, units = "in", res=800)
par (mar=c(4,4,1,1), mgp=c(2,.7,0), tck=-.01, las=1, bty="l")   # 
plot(fitord, "loadplot", main = "")  ## aspect ratio = 1
dev.off()

### draw map
dispersal$geo <- countrycode(dispersal$country, origin = 'country.name', destination = 'iso2c')
dispersal$geo[dispersal$geo=="GB"] <- "UK"
dispersal$geo[dispersal$geo=="GR"] <- "EL"

brbg3 <- brewer.pal(11,"BrBG")[c(8,2,11)]

get_eurostat_geospatial(resolution = 10, 
                        nuts_level = 0, 
                        year = 2013)
SHP_0 <- get_eurostat_geospatial(resolution = 10, 
                                 nuts_level = 0, 
                                 year = 2013)

# Define breaks and custom labels
legend_breaks <- seq(0, 2, 0.2)
legend_labels <- ifelse(legend_breaks == 0.2, "", as.character(legend_breaks))

library(scales)

# Mild custom transformation to emphasize high-end variation
mild_high_trans <- trans_new(
  name = "mild_high_emphasis",
  transform = function(x) x^1.8,
  inverse = function(x) x^(1/1.8)
)

eb_shp <- dispersal %>% 
  dplyr::select(geo, overall_restrictiveness) %>% 
  inner_join(SHP_0, by = "geo") %>% 
  st_as_sf()

# Figure 1b
p2 <- eb_shp %>% 
  ggplot(aes(fill = overall_restrictiveness)) +
  geom_sf(size = 0.1, color = "#F3F3F3") +
  scale_fill_gradient2(
    low = "yellow", 
    mid = "orange", 
    high = "black",
    midpoint = 1.0,
    space = "Lab",
    name = "Restrictiveness",
    trans = mild_high_trans,
    limits = c(0, 2),
    breaks = legend_breaks,
    labels = legend_labels,  
    guide = guide_colorbar(
      direction = "vertical", 
      title.position = "top", 
      label.position = "right",  
      barwidth = unit(0.4, "cm"), 
      barheight = unit(7, "cm"),  
      ticks = TRUE
    )
  ) + 
  scale_x_continuous(limits = c(-10, 35)) +
  scale_y_continuous(limits = c(35, 65)) +
  labs(title="B. Choropleth of dispersal policies",subtitle = "",
       caption = ""
  ) +
  theme_void() +
  theme(legend.position = "left", plot.title = element_text(face = "bold", size = 14), plot.margin = unit(rep(5, 4), "mm"))

library(patchwork)
combined <- p1 + p2 + plot_layout(ncol = 2)
ggsave("figure1.jpeg", combined, width = 12, height = 6.5, dpi = 800)


## SUBSET datset to the variables actually used and maybe separate data preparation and data analysis
# dispersal <- subset()

## external validity tests

# IMPIC dataset
cor(dispersal$AvgS_c11, dispersal$overall_restrictiveness, use="complete.obs") # 
fit <- lm(overall_restrictiveness ~ AvgS_c11, data=dispersal)

# DWRAP dataset: move: value 0 (no free movement), 1 (permission required), 2 (free movement)
dispersal$move_reverse <- max(dispersal$move, na.rm = TRUE) - dispersal$move
cor(dispersal$move_reverse, dispersal$overall_restrictiveness, use="complete.obs") # -

## Plot association of RDPs with contextual variables
plot1 <- ggplot(data = dispersal, aes(x = total1, y = overall_restrictiveness)) + 
  geom_point() + 
  geom_smooth(color = "lightblue", fill = "lightblue") + 
  geom_smooth(method = 'lm', lwd = 2, fill = "blue") +
  geom_rug() +
  xlab("Asylum requests (log)") + 
  ylab("Dispersal policy restrictiveness") + 
  coord_cartesian(ylim = c(0, 2)) + 
  theme_ipsum()

plot2 <- ggplot(data = dispersal, aes(x = localautonomy_2020, y = overall_restrictiveness)) + 
  geom_point() + 
  geom_smooth(color="lightblue", fill="lightblue") + 
  xlab("Local autonomy index") +  
  coord_cartesian(ylim = c(0, 2)) + 
  ylab("Dispersal policy restrictiveness") + 
  geom_rug() +
  geom_smooth(method = 'lm',lwd=2, fill="blue") + 
  theme_ipsum()

plot3 <- ggplot(data = dispersal, aes(x = rrpp, y = overall_restrictiveness)) + 
  geom_point() + geom_smooth(color="lightblue", fill="lightblue") + 
  xlab("Vote share of radical-right populist parties (in%)") + 
  coord_cartesian(ylim = c(0, 2)) + 
  ylab("Dispersal policy restrictiveness") + 
  geom_rug() + 
  geom_smooth(method = 'lm',lwd=2, fill="blue") + 
  theme_ipsum()

jpeg("figure4.jpeg", width = 11, height = 4, units = "in", res=800)
ggarrange(plot1,plot2,plot3, nrow=1, ncol=3)
dev.off()

### Estimate regression models

# first scale all independent variable to a range from 0 to 1
dispersal$localautonomy_2020 <- rescale(dispersal$localautonomy_2020, to = c(0, 1))
dispersal$rrpp <- rescale(dispersal$rrpp, to = c(0, 1))
dispersal$total1 <- rescale(dispersal$total1, to = c(0, 1))

Index <- lm(overall_restrictiveness ~ total1 + localautonomy_2020 + rrpp, data=dispersal)
summary(Index) #
IndexRefugees <- lm(refugee_restrictiveness ~ total1 + localautonomy_2020 + rrpp, data=dispersal)
summary(IndexRefugees) #
IndexGovernments <- lm(subnational_restrictiveness ~ total1 + localautonomy_2020 + rrpp, data=dispersal)
summary(IndexGovernments) #

Bindingness <- lm(bindingness ~ total1 + localautonomy_2020 + rrpp, data=dispersal)
summary(Bindingness) # 
Criteria <- lm(criteria ~ total1 + localautonomy_2020 + rrpp, data=dispersal)
summary(Criteria) # 
Governance <- lm(governance ~ total1 + localautonomy_2020 + rrpp, data=dispersal)
summary(Governance) # 
Sanction <- lm(sanction ~ total1 + localautonomy_2020 + rrpp, data=dispersal)
summary(Sanction) # 
Monitoring <- lm(monitoring ~ total1 + localautonomy_2020 + rrpp, data=dispersal)
summary(Monitoring) # 

library(stargazer)
stargazer(Index, Bindingness,  Sanction, Governance, Criteria, Monitoring, 
          title="Model output", 
          align=F, 
          covariate.labels=c("Refugee arrivals", "Local autonomy",  "Radical-right strength"),  
          no.space=TRUE, 
          dep.var.caption="", 
          model.numbers=T,
          star.cutoffs = c(0.05, 0.01, 0.001))