---
title: "AreaMatch"
subtitle: "Assessing Geospatial Linking of Incongruent Units"
author: 
  - "Anne Stroppe"
  - "Stefan Jünger"
image: img/logo_areamatch.jpg
image-alt: logo
bibliography: references.bib
format: html
---


```{r, echo = F}

knitr::opts_chunk$set(
  echo = TRUE, 
  fig.align = "center", 
  message = FALSE,
  warning = FALSE
)

```


# At a glance

What do we do when we want to link overlapping but incongruent spatial units? This tool provides replicable code to guide you through the linking process and evaluate your linking strategy. The tool application highlights the following steps

- Apply three different matching methods (centroid matching, areal matching, and areal interpolation)
- Evaluate the differences, effectiveness, and accuracy of the different matching techniques
- Evaluate the impact of the matching decision on the research question 

### Table of Content


[Introduction](#introduction)

[Setup](#setup)

[Tool application](#tool-application)

[Conclusion and recommendations](#conclusion-and-recommendations)


# Introduction

Data integration of social indicators from surveys with geospatial context variables has rapidly progressed. 
Applications in the social sciences are manifold, covering issues such as conflict and migration, political participation, environmental attitudes, and inequality. 
Geospatial approaches allow researchers to introduce new perspectives in explaining societal processes and emphasize the local aspects of globally relevant questions.

As an example, this tool application addresses the discourse on the rural-urban divide, which has received considerable attention in recent years. 
Rural areas can serve as breeding grounds for political discontent because citizens’ perceptions of their living environment may foster a sense of neglect, resource deprivation, and lack of societal respect. 
Consequently, citizens may become more susceptible to far-right narratives, which could explain the rise of populist and far-right sentiments in rural areas. 
Therefore, one central research question is: Is living in a less densely populated area, such as a rural place, associated with far-right party preferences?
To answer this question, it is necessary to link survey data with information about the respondents’ living environment.

The linking process is crucial to answering this research question but is prone to certain pitfalls, such as incongruent geospatial units, e.g. overlapping geospatial units that do not share a common border.
In many (online) surveys, respondents can indicate where they live by, for example, entering the ZIP code of their home address. 
This information allows for the linking of contextual information about the respondents’ living environment, like population density, to the survey data. 
However, contextual information is often not available at the ZIP code level.
Instead, it refers to administrative units, such as municipalities. 
These two different spatial units are often not congruent, increasing uncertainty in the linking process.

```{r, echo = F}

# Load necessary libraries
library(sf)
library(ggplot2)
library(dplyr)
library(viridis)

# Function to create a realistic polygon for Municipality 1
create_municipality_polygon <- function() {
  coords <- matrix(c(
    0, 0,     # Bottom-left
    3, 0.5,   # Bottom-right
    4, 2.5,   # Middle-right
    3, 3,     # Top-right
    0, 3,     # Top-left
    -1, 1.5,  # Middle-left
    0, 0      # Closing the polygon
  ), ncol = 2, byrow = TRUE)
  
  sf::st_polygon(list(coords)) |>
    sf::st_sfc(crs = 4326) |>
    sf::st_sf()
}

# Create Municipality 1
municipality1 <- create_municipality_polygon()
municipality1$id <- "Municipality"

# Function to create a realistic polygon for Postal Code Area 1
create_zip_code_polygon1 <- function() {
  coords <- matrix(c(
    1, 0,     # Bottom-left
    3, 0,     # Bottom-right
    3, 2,     # Top-right
    1, 2,     # Top-left
    1, 0      # Closing the polygon
  ), ncol = 2, byrow = TRUE)
  
  sf::st_polygon(list(coords)) |>
    sf::st_sfc(crs = 4326) |>
    sf::st_sf()
}

# Create Postal Code Area 1
zip_code1 <- create_zip_code_polygon1()
zip_code1$id <- "ZIP Code"

# Combine the shapes into a single data frame for plotting
fake_zip_codes <- zip_code1

# Calculate centroid for the ZIP code areas
fake_zip_codes_centroid <- sf::st_centroid(fake_zip_codes)

# Create Municipality 1
municipality2 <- create_municipality_polygon()
municipality2$id <- "Municipality"
zip_code1$id <- "ZIP Code"

# Generate colors from the inferno palette
viridis_colors <- viridis::viridis(5)

# Create a named vector for fill colors
color_mapping <- c(
  "Municipality 1" = viridis_colors[1],       # Color for Municipality
  "ZIP Code" = viridis_colors[3]      # Color for Fake ZIP Codes
)

# Set desired legend order
desired_order <- c("Municipality 1", "ZIP Code")

# Plot the selected municipality and postal code areas with overlap highlighted
ggplot2::ggplot() +
  ggplot2::geom_sf(data = municipality1, aes(fill = factor("Municipality 1", levels = desired_order)), alpha = 0.4, color = "black") +
  ggplot2::geom_sf(data = fake_zip_codes, aes(fill = factor("ZIP Code", levels = desired_order)), alpha = 0.8, color = "black") +
  ggplot2::scale_fill_manual(values = color_mapping, name = "Layers") +
  ggplot2::labs(title = "Overlapping But Incongruent Boundaries",
       fill = "Layers") +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "right")



```

This tool addresses this challenge by highlighting linking techniques that can be used to successfully link the datasets and assessing the uncertainty of the linking process.

## Data Prerequisites and Description

To replicate this tool, researchers need the following:

- Survey data that includes an identifier for the target areal units where respondents reside and for which the attribute needs to be linked or estimated
- The geometries of the target areal units 
- The geometries of the source unit that include the attributes information they are interested in

For this tool application, we rely on three data sources: (synthetic) survey data including self-reported ZIP codes, the geographies of German ZIP code areas, and the geographies of German municipalities, including attributes of interest, in this case population density.

We base our analysis example on survey data from the German Longitudinal Election Study (GLES). 
The GLES Tracking consists of short cross-sectional online surveys (CAWI) conducted three times a year.
Each cross-sectional sample includes approximately 1,000 respondents [@gles_gles_2024].
For this analysis, we create a synthetic data set based on the [GLES Tracking November 2023, T56](https://search.gesis.org/research_data/ZA7714){target="_blank"}. 
Respondents were asked to indicate their approval of the German radical right party AfD on a scale from 1 ("I don't think highly of this party at all") to 11 ("I think very highly of this party") (variable t14h). 
At the end of the questionnaire, they were also asked to enter the 5-digit ZIP code of their primary residence (variable t71).

::: {.callout-important}
Due to data protection regulations, ZIP codes cannot be published in the Scientific Use Files of the survey data, but only accessed through the [Secure Data Center at GESIS](https://www.gesis.org/en/services/processing-and-analyzing-data/analysis-of-sensitive-data/secure-data-center-sdc){target="_blank"}. For this tool application, we provide simulated data that reproduces the real correlation of the GLES data without allowing conclusions to be drawn about the place of residence or the interviewee.
:::

The two additional data sources provide the necessary geographies of ZIP code areas and municipalities and their corresponding attributes. 
We rely on the Open Data portal ArcGIS Hub to access the geometries of the [German ZIP code areas](https://hub.arcgis.com/maps/esri-de-content::postleitzahlengebiete-in-deutschland){target="_blank"} and [municipalities](https://hub.arcgis.com/maps/60eb682c95f44ba7b10fee66d871859d){target="_blank"}.
Data can be downloaded directly as a shapefile or GeoPackage file. 
We work with a geospatial data type called vector data, which is organized similarly to any other data table: each row represents a geometric object (e.g., a ZIP code area or a municipality), and each column holds an attribute (e.g., population density).
For detailed information on handling geospatial data and using the package `sf`, you can refer to our course [“Introduction to Geospatial Techniques for Social Scientists in R”](https://stefanjuenger.github.io/gesis-workshop-geospatial-techniques-R-2024/){target="_blank"}.

## Tool Functions To Assess the Spatial Linking Process

This tool application provides replicable code to guide the user through the linking process and assess their linking strategy. The tutorial highlights the following steps:

- Centroid linkage: This technique involves linking data based on the central point (centroid) of a geographic area. For example, the centroid of a ZIP code area is linked to the municipality in which the centroid is located.

- Areal matching: This method matches entire areas to one another, such as linking ZIP code areas to municipalities based on the largest overlapping area.

- Areal interpolation: This technique redistributes data from one set of geographic units to another. The basic approach, area-weighted spatial interpolation, takes into account the size of the areal overlap and assigns the value of the overlapping area in proportion to the overlap. 

- Overall Assessment of Linking Techniques: Evaluates the differences, effectiveness, and accuracy of the different linking techniques to determine the best method for integrating the survey and geospatial data at hand.

- Influence of Linking Technique on Research Question: Analyzes how the choice of linking technique affects the outcomes and interpretations of the research question, ensuring that the conclusions drawn are robust and reliable.

# Setup

## Getting started

### Packages

There are several packages out there that allow geospatial data handling in R.
Since we are working mainly with vector data, we heavily rely on the package [sf](https://doi.org/10.32614/CRAN.package.sf).
The packages `dplyr`, `ggplot2`, and `tibble` are used for data manipulation and visualization.


```{r }

# Load necessary packages ----

library(dplyr)       # For data manipulation
library(ggplot2)     # For data visualization
library(sf)          # For handling spatial (geometric) data
library(tibble)      # For creating and managing tibbles (data frames)

```

### Geospatial Data

In Germany, the 5-character ZIP code areas often do not align with municipal boundaries, leading to both overlap and incongruence. This mismatch can be more pronounced in rural areas, where larger ZIP code regions may encompass multiple small municipalities. In contrast, urban areas tend to have several distinct ZIP codes that fall within a single municipality. For example, the capital Berlin, as a German 'city-state', is one municipality with 190 ZIP codes.

When selecting and loading the geospatial data, we always check for three sources of errors:

- Consistent (projected) coordinate reference systems. In our case, both shapefiles are projected in WGS 84 (EPSG: 3857).

- Timeliness of administrative boundaries and territorial reforms. For example, we use the 2022 municipality boundaries and the 2023 ZIP code areas.

- Existence, column names, and format of identifiers and variables of interest. For instance, we require the ZIP code and municipality code to be in character format in our data frame.

```{r}

# Load ZIP code data ----

zip_codes <-
  # Load spatial data from GeoPackage file
  sf::st_read("./data-raw/PLZ_Gebiete_7155512659375659703.gpkg") |>
  # Select only relevant columns: ZIP code and population
  dplyr::select(zip_code = plz, inhabitants_zip_code = einwohner)

# Load municipality data ----

municipalities <-
  # Load spatial data from GeoPackage file
  sf::st_read("./data-raw/Gemeindegrenzen_2022__mit_Einwohnerzahl_4398740898366155627.gpkg") |>
  # Select AGS (municipality code), population and area size in square kilometers columns
  dplyr::select(ags = AGS, inhabitants_municipality = EWZ, area_municipality = KFL)

```

We also advise to always plot the data for a visual inspection of the projection and completeness of the geospatial information.

```{r, echo = F}

library(patchwork)

# Create the municipalities plot
g_mun <- ggplot(data = municipalities) +
  ggplot2::geom_sf(fill = "lightgrey", color = "black", size = 0.05) +  
  ggplot2::ggtitle("Municipalities") +
  ggplot2::theme_minimal() 

# Create the ZIP codes plot
g_zip <- ggplot() +
  ggplot2::geom_sf(data = zip_codes, fill = "lightgrey", color = "black", size = 0.05) + 
# ggplot2::geom_sf(data = zip_codes_reduced, fill = "red", color = "black") +
  ggplot2::ggtitle("ZIP Codes") +
  ggplot2::theme_minimal() 

# Combine the two plots using patchwork
g_mun + g_zip



```


### Survey Data

Due to data protection regulations, we cannot use real-world survey data, but we base our tool application on simulated data that were created on  the basis of the [GLES Tracking November 2023, T56](https://search.gesis.org/research_data/ZA7714). 
To prepare the simulated data, we clean the original survey data to retain only valid ZIP codes and then link the data with geospatial data using centroid and areal matching techniques, followed by area-weighted interpolation for population and area estimates. 
Finally, we calculate correlations between the variable of interest, "afd_rating" and population densities derived from the geospatial data.
In a next step, we generate random ZIP codes, introducing some invalid entries to mimic the real-world data.
We repeat the linking process and then simulate the dependent variable based on the original correlations to derive the simulated data.
To replicate the creation of the simulated survey data, we provide the R-Script [here](./test/create_simulated_survey.R).

```{r}

survey <- readRDS("./data-raw/simulated_survey_data.rds")

```


##  Inspecting data

Before linking survey data, it is crucial to inspect the ZIP codes provided by respondents, especially when  participants entered their ZIP codes independently.
The self-reported nature of data entry can lead to several issues, like refusal to respond or typographical errors, which can result in invalid or missing ZIP codes. 
The provided code categorizes the self-reported ZIP codes into three statuses: "Non Response", "Invalid ZIP Code" (For instance, ZIP codes that are not five characters long or that do not exist in the predefined list of valid ZIP codes), and "Valid". 
By utilizing this categorization, researchers can identify and address potential data quality issues before data integration.
<!--- footnote: could also try to resolve issues with ZIP codes -->

```{r}

# Create a new column 'status' to categorize the ZIP codes
survey <- survey |>
  dplyr::mutate(
    # Use case_when to classify each ZIP_code into different statuses
    status = dplyr::case_when(
      # Condition 1: If zip_code is NA, classify as "1 Non Response"
      is.na(zip_code) ~ "1 Non Response",
      # Condition 2: If zip_code does not have exactly 5 characters, 
      # classify as "2 Invalid ZIP Code"
      nchar(zip_code) != 5 ~ "2 Invalid ZIP Code",
      # Condition 3: If zip_code is not found in the list of valid ZIP codes, 
      # classify as "2 Invalid ZIP Code"
      !zip_code %in% zip_codes$zip_code ~ "2 Invalid ZIP Code",
      # Condition 4: If zip_code is found in the list of valid ZIP codes, 
      # classify as "3 Valid"
      zip_code %in% zip_codes$zip_code ~ "3 Valid"
    )
  )

# Create a summary table to count the occurrences of each status
summary_table <- survey |>
  dplyr::group_by(status) |>  # Group the data by the 'status' column
  dplyr::summarise(
    # Count the number of occurrences for each status
    count = dplyr::n(),  
    # Calculate the frequency 
    freq = (dplyr::n() / nrow(survey)) * 100  
  ) |>
  # Convert the result into a tibble for easier viewing and manipulation
  tibble::as_tibble()  

print(summary_table)

```

In this case, `r summary_table$freq[summary_table$status == "1 Non Response"]`% of respondents did not answer the ZIP code question, and `r summary_table$freq[summary_table$status == "2 Invalid ZIP Code"]`% provided an invalid ZIP code. What we do not check here is whether respondents entered their correct ZIP code. Depending on the data at hand, it might be possible to verify whether other self-reported data, such as the federal state, or data from other sources align with this information.

For the ongoing linking process, we exclude all observations that did not enter a ZIP code or entered an invalid one, working instead with a subsample of all valid ZIP codes in the survey. Since we will use more spatial techniques in the following analysis, we rely on this reduced spatial data frame.

```{r}

# Create a reduced dataset based on the valid ZIP codes from the survey data
zip_codes_valid <-
  survey |>
  # Filter the survey data to include only rows with a valid ZIP code status
  dplyr::filter(status == "3 Valid") |>
  # Perform a left join with the zip_codes dataset on the 'zip_code' column
  dplyr::left_join(zip_codes, ., by = "zip_code") |>
  # Convert the resulting data frame to a simple features (sf) 
  sf::st_as_sf()

```

::: {.callout-important} 
If you want to create a complete matching list between all ZIP codes and all municipalities in Germany, you can also use the complete ZIP code data set. However, it might take a while to create this data set and run the spatial joins. <!-- The script and data are provided under additional material -->
:::

# Tool application

## Matching methods

### Data simulation 

We're sampling 1,000 points within each valid ZIP code area to later assess the accuracy of our linking techniques. 
By simulating 1,000 hypothetical living locations for respondents who reported living in a specific ZIP code area, we can evaluate whether our method accurately aligns each sampled point with the correct municipality.
This approach allows us to identify any discrepancies or inaccuracies in the linking process and compare the three linking methods.

```{r}

# Sample points within ZIP code areas ----
# Randomly sample points within each ZIP code area. Note: This can take a long
# time due to the number of points sampled.
points_in_zip_codes <-
  zip_codes_valid |>
  # Sample 1000 points per area
  sf::st_sample(size = c(1000, 1000), progress = TRUE, exact = FALSE) |>
  # Convert sampled points to an sf object
  sf::st_as_sf() |>
  # Spatially join sampled points with ZIP code data
  sf::st_join(zip_codes) |>
  # Arrange points by ZIP code for easier viewing
  dplyr::arrange(zip_code) |>
  # Assign unique IDs to each sampled point
  dplyr::mutate(id = 1:dplyr::n()) |>
  # Select only the ID and ZIP code columns
  dplyr::select(id, zip_code)


# Join sampled points with municipality inhabitants data ----
# Spatially join the sampled points data with municipality data to match each
# point with actual inhabitants data from municipalities
points_with_real_inhabitants_municipality <-
  points_in_zip_codes |>
  # Join points with municipality data
  sf::st_join(municipalities |> dplyr::select(-area_municipality)) |>
  # Rename population column for clarity
  dplyr::rename(real_inhabitants_municipality = inhabitants_municipality)



```

### Centroid matching

This method uses the centroid (center point) of each ZIP code area to find the corresponding municipality.
It assumes that the center point accurately represents the ZIP code's location in terms of municipality boundaries.

```{r, echo = F}

# Load necessary library
library(patchwork)

# Function to create a realistic polygon for Municipality 1
create_municipality_polygon <- function() {
  coords <- matrix(c(
    0, 0,     # Bottom-left
    3, 0.5,   # Bottom-right
    4, 2.5,   # Middle-right
    3, 3,     # Top-right
    0, 3,     # Top-left
    -1, 1.5,  # Middle-left
    0, 0      # Closing the polygon
  ), ncol = 2, byrow = TRUE)
  
  sf::st_polygon(list(coords)) |>
    sf::st_sfc(crs = 4326) |>
    sf::st_sf()
}

# Create Municipality 1
municipality1 <- create_municipality_polygon()
municipality1$id <- "Municipality"

# Function to create a realistic polygon for Postal Code Area 1
create_zip_code_polygon1 <- function() {
  coords <- matrix(c(
    1, 0,     # Bottom-left
    3, 0,     # Bottom-right
    3, 2,     # Top-right
    1, 2,     # Top-left
    1, 0      # Closing the polygon
  ), ncol = 2, byrow = TRUE)
  
  sf::st_polygon(list(coords)) |>
    sf::st_sfc(crs = 4326) |>
    sf::st_sf()
}

# Create Postal Code Area 1
zip_code1 <- create_zip_code_polygon1()
zip_code1$id <- "ZIP Code"

# Combine the shapes into a single data frame for plotting
fake_zip_codes <- zip_code1

fake_zip_codes_centroid <- sf::st_centroid(fake_zip_codes)

# Add centroids to ZIP code data

# Function to create a realistic polygon for Municipality 1
create_municipality_polygon <- function() {
  coords <- matrix(c(
    0, 0.15,     # Bottom-left
    1.5, 0.15,   # Bottom-right
    2, 1.5,
    2.2, 1,
    2.4, 0.15,   # Middle-right
    3.8, 0.4,
    3.5, 3,     # Top-right
    0, 3,     # Top-left
    1.5, 2.5,
    -1, 1.5,  # Middle-left
    0, 0.15      # Closing the polygon
  ), ncol = 2, byrow = TRUE)
  
  sf::st_polygon(list(coords)) |>
    sf::st_sfc(crs = 4326) |>
    sf::st_sf()
}

# Create Municipality 1
municipality2 <- create_municipality_polygon()
municipality2$id <- "Municipality"
zip_code1$id <- "ZIP Code"

library(viridis)
# Generate colors from the inferno palette
viridis_colors <- viridis::viridis(5)

# Create a named vector for fill colors
color_mapping1 <- c(
  "Municipality 1" = viridis_colors[1],       # Color for Municipality
  "ZIP Code" = viridis_colors[3],      # Color for Fake ZIP Codes
  "Centroid" = "black"                         # Color for Centroid
)

# Set desired legend order
desired_order1 <- c("Municipality 1", "ZIP Code", "Centroid")


# Plot for municipality1
map1 <- ggplot2::ggplot() +
  ggplot2::geom_sf(data = municipality1, aes(fill = factor("Municipality 1", levels = desired_order1)), alpha = 0.4, color = "black") +
  ggplot2::geom_sf(data = fake_zip_codes, aes(fill = factor("ZIP Code", levels = desired_order1)), alpha = 0.8, color = "black") +
  ggplot2::geom_sf(data = fake_zip_codes_centroid, aes(fill = factor("Centroid", levels = desired_order1)), shape = 20, color = "black", size = 4) +
  ggplot2::scale_fill_manual(values = color_mapping1, name = "Layers") +
  ggplot2::labs(title = "High Accuracy") +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "right", 
                 plot.margin = unit(c(0.1, 0.1, 0.5, 0.1), "cm"),
                 plot.title = element_text(size = 10))

# Plot for municipality2
map2 <- ggplot2::ggplot() +
  ggplot2::geom_sf(data = municipality2, aes(fill = factor("Municipality 1", levels = desired_order1)), alpha = 0.4, color = "black") +
  ggplot2::geom_sf(data = fake_zip_codes, aes(fill = factor("ZIP Code", levels = desired_order1)), alpha = 0.8, color = "black") +
  ggplot2::geom_sf(data = fake_zip_codes_centroid, aes(fill = factor("Centroid", levels = desired_order1)), shape = 20, color = "black", size = 4) +
  ggplot2::scale_fill_manual(values = color_mapping1, name = "Layers") +
  ggplot2::labs(title = "Low Accuracy") +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "right", 
                 plot.margin = unit(c(0.5, 0.1, 0.1, 0.1), "cm"),
                 plot.title = element_text(size = 10))

# Combine the two plots using patchwork
(map1 / map2) +  
  patchwork::plot_layout(guides = "collect", nrow = 2) +  
  patchwork::plot_annotation(
    title = "Accuracy of Centroid Linking",
    theme = theme(plot.margin = unit(c(0.1, 0, 0.1, 0), "cm"))
  )


```


One benefit of this approach is its simplicity and efficiency, allowing for a quick linking process. 
However, a drawback is that it may overlook important local variations, as the centroid may not accurately reflect the distribution of residents or land use within the ZIP code area, nor does it account for the actual overlap or shape of the units.
To do the linking process, we need to calculate centroids for each ZIP code and use the `sf::st_join` function to identify in which municipality each centroid is located.

```{r}

# Centroid linking ----
centroid_matched <-
  zip_codes_valid |>
  # Calculate the centroid for each ZIP code area
  sf::st_point_on_surface() |>
  # Spatially join centroids with municipality data
  sf::st_join(municipalities) |>
  # Select relevant columns
  dplyr::select(zip_code, inhabitants_zip_code, inhabitants_municipality, area_municipality ) |>
  # Remove any duplicate rows
  dplyr::distinct() |>
  # Arrange by ZIP code for easy viewing
  dplyr::arrange(zip_code)

```

### Areal matching

This method assigns each ZIP code area to the municipality in which the  majority of its area lies.

```{r, echo = FALSE}


# Load necessary libraries
library(sf)        # For spatial data handling
library(viridis)   # For color palettes
library(ggplot2)   # For plotting
library(ggpattern) # For pattern filling in geom_sf

# Function to create a realistic polygon for Municipality 1
create_municipality_polygon <- function() {
  coords <- matrix(c(
    0, 0,     # Bottom-left
    1, 0,
    3.05, 0.52,   # Bottom-right
    3.4, 1.5,   # Middle-right
    3, 2,     # Top-right
    0.5, 2,     # Top-left
    0.8, 1.5,  # Middle-left
    0, 0     # Closing the polygon
  ), ncol = 2, byrow = TRUE)
  
  sf::st_polygon(list(coords)) |>
    sf::st_sfc(crs = 4326) |>
    sf::st_sf()
}

# Create Municipality 1
municipality1 <- create_municipality_polygon()
municipality1$id <- "Municipality 1"


create_municipality_polygon <- function() {
  coords <- matrix(c(
    0, 0,     # Bottom-left
    1,0,
    3.05, 0.52,   # Bottom-right
    3.4, 1.5,   # Middle-right
    3.4,-0.5,
    0,0
  ), ncol = 2, byrow = TRUE)
  
  sf::st_polygon(list(coords)) |>
    sf::st_sfc(crs = 4326) |>
    sf::st_sf()
}

municipality2 <- create_municipality_polygon()
municipality2$id <- "Municipality 2"

# Function to create a realistic polygon for Postal Code Area 1
create_zip_code_polygon1 <- function() {
  coords <- matrix(c(
    1, -0.5,     # Bottom-left
    3, -0.5,    # Bottom-right
    3, 1.5,     # Top-right
    1, 1.5,     # Top-left
    1, -0.5      # Closing the polygon
  ), ncol = 2, byrow = TRUE)
  
  sf::st_polygon(list(coords)) |>
    sf::st_sfc(crs = 4326) |>
    sf::st_sf()
}

# Create Postal Code Area 1
zip_code1 <- create_zip_code_polygon1()
zip_code1$id <- "ZIP Code"


# Overlapping Area
# Function to create a realistic polygon for Postal Code Area 1
create_overlap_polygon <- function() {
  coords <- matrix(c(
    1, 0,     # Bottom-left
    3, 0.5,    # Bottom-right
    3, 1.5,     # Top-right
    1, 1.5,     # Top-left
    1, 0       # Closing the polygon
  ), ncol = 2, byrow = TRUE)
  
  sf::st_polygon(list(coords)) |>
    sf::st_sfc(crs = 4326) |>
    sf::st_sf()
}

overlap <- create_overlap_polygon()
overlap$id <- "Largest Overlap"




# Plot the selected municipality and postal code areas with overlap highlighted


# Create a named vector for fill colors
color_mapping2 <- c(
  "Municipality 1" = viridis::viridis(5)[1],  # Color for Municipality 1
  "Municipality 2" = viridis::viridis(5)[2],  # Color for Municipality 2
  "ZIP Code" = viridis::viridis(5)[3],        # Color for ZIP Code
  "Largest Overlap" = viridis::viridis(5)[4]          # Color for Overlap
)

# Set desired legend order
desired_order2 <- c("Municipality 1", "Municipality 2", "ZIP Code", "Largest Overlap")

# Create a plot
ggplot2::ggplot() +
  ggplot2::geom_sf(data = municipality1, aes(fill = factor("Municipality 1", levels = desired_order2)), alpha = 0.3, color = "black") +  # Color for Municipality 1
  ggplot2::geom_sf(data = municipality2, aes(fill = factor("Municipality 2", levels = desired_order2)), alpha = 0.3, color = "black") +  # Color for Municipality 2
  ggplot2::geom_sf(data = zip_code1, aes(fill = factor("ZIP Code", levels = desired_order2)), alpha = 0.8, color = "black") +      # Color for ZIP Code
  ggpattern::geom_sf_pattern(data = overlap, aes(fill = factor("Largest Overlap", levels = desired_order2)), 
                  pattern = "stripe", pattern_fill = "black", pattern_spacing = 0.02, pattern_density = 0.02, alpha = 0.6, color = "grey") +  # Closer stripe pattern for Overlap
  ggplot2::scale_fill_manual(values = color_mapping2, 
                    name = "Layers") +  # Use the named vector for specific colors
  ggplot2::labs(title = "Areal Matching: Largest Overlap") +
  ggplot2::theme_minimal() +
  ggplot2::theme(legend.position = "right")


```

The `sf::st_join` function identifies which municipality overlaps with each ZIP code area. 
The `largest = TRUE` argument indicates that if a ZIP code overlaps with multiple municipalities, it will only keep the attributes, e.g. the number of inhabitants, of the municipality with the largest overlapping area.

```{r}

# Areal matching method ----

areal_matched <-
  zip_codes_valid |>
  # Spatial join using the largest overlap municipality for each ZIP code area
  sf::st_join(municipalities, left = TRUE, largest = TRUE) |>
  # Select relevant columns
  dplyr::select(zip_code, inhabitants_zip_code, inhabitants_municipality, area_municipality) |>
  # Remove duplicates to ensure each ZIP code matches one municipality
  dplyr::distinct() |>
  # Sort by ZIP code
  dplyr::arrange(zip_code)

```


Compared to centroid matching, areal matching considers the full geographic shape of ZIP code areas, providing a more accurate association with municipalities based on actual boundaries. One drawback of areal matching is that it can be more computationally intensive, especially when dealing with a large number of geographic units.

Another challenge arises when multiple municipalities overlap with a single ZIP code, requiring additional decisions about which municipality to assign. In this application, we take a straightforward approach and assign the municipality with the largest overlap.

However, we can also explore the number of municipalities each ZIP code overlaps with. 
To do this, we can modify the code to calculate the number of overlapping municipalities as well as the share of area overlap.


```{r}

areal_matched <- zip_codes_valid |>
  # Perform the spatial join using the largest overlap municipality for each ZIP code
  sf::st_join(municipalities, left = TRUE, largest = TRUE) |>
  # Count the number of municipalities that overlap each ZIP code
  dplyr::group_by(zip_code) |>
  #Identify ambiguous cases (those with more than 1 municipality)
  dplyr::mutate(num_municipalities_overlap = n()) |>
  # Select relevant columns
  dplyr::select(zip_code, inhabitants_zip_code, inhabitants_municipality, area_municipality, num_municipalities_overlap) |>
  # Remove duplicates to ensure each ZIP code is listed once
  dplyr::distinct()

# Filter ambiguous cases: ZIP codes with more than one municipality overlap
ambiguous_cases <- areal_matched |>
  dplyr::filter(num_municipalities_overlap > 1)

# Spatial intersection between ZIP codes and municipalities
overlap_areas <- sf::st_intersection(ambiguous_cases, municipalities)

# Calculate overlap area and share of overlap
overlap_areas <- overlap_areas %>%
  dplyr::mutate(
    overlap_area = sf::st_area(.),  # Area of overlap
    total_zip_area = sf::st_area(
      ambiguous_cases[match(zip_code, ambiguous_cases$zip_code),]), # ZIP code area
    share_of_overlap = overlap_area / total_zip_area  # Calculate share of overlap
  )

```

We use this information to generate summary plots that highlight the uncertainty introduced by the areal matching method in the linking process.


```{r, echo = F}  

# Create summary table: count ZIP codes by number of overlapping municipalities
overlap_summary <- areal_matched %>%
  dplyr::group_by(num_municipalities_overlap) %>%
  dplyr::summarise(num_zip_codes = n()) %>%
  dplyr::arrange(num_municipalities_overlap)

# Create bar chart of ZIP code overlap counts
summary_table_plot <- 
  ggplot2::ggplot(overlap_summary,
                  aes(x = as.factor(num_municipalities_overlap), 
                      y = num_zip_codes)) +
  ggplot2::geom_col(fill = viridis::viridis(5)[1]) +  # Bar chart for overlap
  ggplot2::labs(
    title = "Summary of ZIP Code Overlaps",  # Title for the plot
    x = "Number of Overlapping Municipalities",  # X-axis label
    y = "Number of ZIP Codes"  # Y-axis label
  ) +
  ggplot2::theme_minimal() +  # Minimal plot theme
  ggplot2::theme(
    plot.title = ggplot2::element_text(hjust = 0.5),  # Center title
    axis.text.x = ggplot2::element_text(angle = 0)  # Horizontal X labels
  )


# Create density plot for share of overlap in ambiguous cases
density_plot <- 
  ggplot2::ggplot(overlap_areas, aes(x = as.numeric(share_of_overlap))) +
  ggplot2::geom_density(fill = viridis::viridis(5)[1], color = "black", 
                        alpha = 0.6) +  # Density plot
  ggplot2::labs(
    title = "Aggregated Distribution of Overlap Shares",  # Title for density plot
    x = "Share of Overlap",  # X-axis label
    y = "Density"  # Y-axis label
  ) +
  ggplot2::theme_minimal() +  # Minimal plot theme
  ggplot2::theme(
    plot.title = ggplot2::element_text(hjust = 0.5),  # Center title
    axis.text.x = ggplot2::element_text(angle = 0, hjust = 0.5), # Readable X labels
    axis.text.y = ggplot2::element_text(angle = 0)  # Clear Y labels
  )

# Combine the bar chart and density plot side by side using patchwork
combined_plot <- summary_table_plot | 
  density_plot + 
  patchwork::plot_layout(ncol = 2, widths = c(3, 1))  # Layout: 2 columns

# Display the combined plot
combined_plot


```
About half of the sampled ZIP codes overlap with more than one municipality, which we refer to as "ambiguous" cases. 
However, when examining the share of area overlaps, most of these cases are less ambiguous than initially expected.
ZIP code areas that overlap with multiple municipalities often have extreme overlap shares: either more than 90% of the area is shared with one municipality, or less than 10% overlaps with others. 
Very few cases show nearly equal overlap between two or more municipalities, making the decision to link to the largest municipality in most of these cases indisputable.

We can also map ZIP codes that overlap with one (black outline) or more municipalities (pink outline) in our sample, with the underlying map displaying the (logged) population size of each municipality. 

```{r, echo = F}


# Create a map of ambiguous cases
ambiguous_cases_map <- ggplot2::ggplot() +
  
  # Municipality population as background, with transparency
  ggplot2::geom_sf(
    data = municipalities, aes(fill = as.numeric(inhabitants_municipality)), 
    color = NA, alpha = 0.7  # Transparent fill for better visibility
  ) +
  
  # Add ZIP code borders with no fill (no gray for NA)
  ggplot2::geom_sf(
    data = zip_codes_valid, aes(color = "ZIP Code Border", fill = NA), 
    lwd = 0.5, alpha = 0  # Only borders, no fill
  ) +
  
  # Add ambiguous ZIP codes with pink borders
  ggplot2::geom_sf(
    data = ambiguous_cases, aes(color = "Ambiguous ZIP Code", fill = NA), 
    lwd = 0.5, alpha = 0  # Pink border, no fill
  ) +
  
  # Apply minimal theme for clean visualization
  ggplot2::theme_minimal() +
  
  # Color scale for municipality population (log scale)
  ggplot2::scale_fill_viridis_c(
    name = "Number of Inhabitants",  # Legend title
    trans = "log"  # Log scale for better population range display
  ) +
  
  # Manual color scale for the two layers (no fill, only borders)
  ggplot2::scale_color_manual(
    name = "Legend",  # Title for the border color legend
    values = c("Unambiguous ZIP Code" = "black", "Ambiguous ZIP Code" = "#EA4F88"),  # Set colors for borders
    guide = ggplot2::guide_legend(override.aes = list(lwd = 1))  # Ensure the lines appear clearly in the legend
  ) +
  
  # Map title
  ggplot2::labs(
    title = "Map of Ambiguous ZIP Codes"
  )

# Display the map
ambiguous_cases_map




```

ZIP codes that overlap with multiple municipalities are not concentrated in any specific region; they are spread across all federal states. 
However, there appears to be a partial correlation between a municipality's population size and the likelihood of ZIP code areas overlapping with multiple municipalities. 
In other words, in less densely populated areas, ZIP codes are more likely to span multiple municipal borders.

Overall, the clusters of overlapping ZIP codes are not very pronounced, and the distribution of overlap shares does not raise significant concerns about the performance of the areal matching approach. 
The steps we have outlined here can serve as a guide for others using areal matching.
These steps help assess whether linking to the largest municipality is justifiable based on the data or if alternative thresholds or robustness tests should be considered.

### Areal interpolation

This method uses areal interpolation to distribute municipality inhabitants data proportionally across overlapping ZIP code areas.


```{r, echo = FALSE}

# Load necessary libraries
library(sf)        # For spatial data handling
library(ggplot2)   # For plotting
library(dplyr)     # For data manipulation
library(viridis)   # For color palettes
library(ggpattern) # For pattern filling in geom_sf

# Function to create a realistic polygon for Municipality 1
create_municipality_polygon1 <- function() {
  coords <- matrix(c(
    0, 0,     # Bottom-left
    1, 0,
    3, 0.5,   # Bottom-right
    3.4, 1.5, # Middle-right
    3, 2.5,   # Top-right
    0.5, 2.5, # Top-left
    0.8, 1.5, # Middle-left
    0, 0      # Closing the polygon
  ), ncol = 2, byrow = TRUE)
  
  sf::st_polygon(list(coords)) |>
    sf::st_sfc(crs = 4326) |>
    sf::st_sf()
}

# Create Municipality 1
municipality1 <- create_municipality_polygon1()
municipality1$id <- "Municipality 1"

# Function to create a realistic polygon for Municipality 2
create_municipality_polygon2 <- function() {
  coords <- matrix(c(
    0, 0,     # Bottom-left
    1, 0,
    3, 0.5,   # Bottom-right
    3.4, 1.5, # Middle-right
    3.4, -0.8,  # Bottom-right extension
    0, 0      # Closing the polygon
  ), ncol = 2, byrow = TRUE)
  
  sf::st_polygon(list(coords)) |>
    sf::st_sfc(crs = 4326) |>
    sf::st_sf()
}

# Create Municipality 2
municipality2 <- create_municipality_polygon2()
municipality2$id <- "Municipality 2"

# Function to create a realistic polygon for Postal Code Area
create_zip_code_polygon <- function() {
  coords <- matrix(c(
    1, -0.5,  # Bottom-left
    3, -0.5,  # Bottom-right
    3, 1.5,   # Top-right
    1, 1.5,   # Top-left
    1, -0.5   # Closing the polygon
  ), ncol = 2, byrow = TRUE)
  
  sf::st_polygon(list(coords)) |>
    sf::st_sfc(crs = 4326) |>
    sf::st_sf()
}

# Create Postal Code Area
zip_code <- create_zip_code_polygon()
zip_code$id <- "ZIP Code"


# Create geometries for each municipality's share in the interpolation area
interpolated_municipality1 <- sf::st_intersection(zip_code, municipality1) # Share from Municipality 1
interpolated_municipality1$id <- "Share from Municipality 1"

interpolated_municipality2 <- sf::st_intersection(zip_code, municipality2) # Share from Municipality 2
interpolated_municipality2$id <- "Share from Municipality 2"

# Combine the shares into a single data frame
shares_combined <- dplyr::bind_rows(interpolated_municipality1, interpolated_municipality2)


# Create a named vector for fill colors
color_mapping3 <- c(
  "Municipality 1" = viridis::viridis(5)[1],     # Color for Municipality 1
  "Municipality 2" = viridis::viridis(5)[2],     # Color for Municipality 2
  "ZIP Code" = viridis::viridis(5)[3],           # Color for ZIP Code
  "Shares Municipality 1" = viridis::viridis(5)[4],  # Color for Shares of Municipality 1
  "Shares Municipality 2" = viridis::viridis(5)[5]   # Color for Shares of Municipality 2
)

# Set desired legend order
desired_order3 <- c("Municipality 1", "Municipality 2", "ZIP Code", "Shares Municipality 1", "Shares Municipality 2")


# Create the plot

ggplot2::ggplot() +
  ggplot2::geom_sf(data = municipality1, aes(fill = factor("Municipality 1", levels = desired_order3)), 
          alpha = 0.3, color = "black") +  # Color for Municipality 1
  ggplot2::geom_sf(data = municipality2, aes(fill = factor("Municipality 2", levels = desired_order3)), 
          alpha = 0.3, color = "black") +  # Color for Municipality 2
  ggplot2::geom_sf(data = zip_code, aes(fill = factor("ZIP Code", levels = desired_order3)), 
          alpha = 0.7, color = "black") +  # Color for ZIP Code
  ggpattern::geom_sf_pattern(data = interpolated_municipality1, aes(fill = factor("Shares Municipality 1", levels = desired_order3)), 
                  pattern = "stripe", pattern_fill = "black", pattern_spacing = 0.02, pattern_density = 0.02, alpha = 0.6, color = "grey") +  # Closer stripe pattern for Shares Municipality 1
  ggpattern::geom_sf_pattern(data = interpolated_municipality2, aes(fill = factor("Shares Municipality 2", levels = desired_order3)), 
                  pattern = "circle", pattern_fill = "grey", pattern_spacing = 0.05, pattern_density = 0.02, alpha = 0.6, color = "grey") +  # Adjusted spacing for Shares Municipality 2
  ggplot2::scale_fill_manual(values = color_mapping3, 
                    name = "Layers") +
  ggplot2::theme_minimal()



```

It estimates the inhabitants for each ZIP code based on the proportion of its area that overlaps with each municipality by:

1. Overlap Calculation: For each ZIP code area, the function determines how much of its area overlaps with each municipality and, as such, calculates the intersection of each unit of both areal unit layers.

2. Proportional Distribution: Once the overlaps are identified, it calculates the proportion of each municipality's area that overlaps with the ZIP code area. 

3. Population Estimation: The estimated population for each ZIP code area is then computed by multiplying the municipality's population by the overlap proportion. For example, if a municipality has 10,000 inhabitants and 40% of its area overlaps with a particular ZIP code, then 4,000 inhabitants would be allocated to that ZIP code area. If multiple municipalities overlap with a ZIP code area, it sums the estimated populations from each municipality to get the total population for that ZIP code.

As such, a key advantage of areal interpolation is its flexibility, as it can effectively manage irregular boundaries and varying sizes of geographic units. 
Nonetheless, it operates under the assumption that populations are uniformly distributed within areas, which may not reflect reality and can lead to inaccuracies if the input data is flawed or unrepresentative.
Similar to areal matching, the process can be computationally demanding, especially when dealing with multiple overlapping areas and variable populations, as well as a large number of indicators that need to be interpolated.

```{r}

# Areal interpolation matching method ----

areal_interpolation_matched <-
  sf::st_interpolate_aw(
    # Use municipality inhabitants data for interpolation
    municipalities["inhabitants_municipality"],
    # Target ZIP code areas for the interpolation
    zip_codes_valid,
    # Set to FALSE as population data is not "extensive" (not purely additive)
    extensive = FALSE
  ) |>
  # Combine interpolated results with original ZIP code data
  dplyr::bind_cols(
    zip_codes_valid |>
      # Drop geometry to avoid duplication issues in final output
      sf::st_drop_geometry() |>
      # Select only ZIP code and its inhabitants count
      dplyr::select(zip_code, inhabitants_zip_code)
  ) |>
  # Choose relevant columns for output
  dplyr::select(zip_code, inhabitants_zip_code, inhabitants_municipality) |>
  # Ensure unique rows
  dplyr::distinct() |>
  # Sort by ZIP code
  dplyr::arrange(zip_code)

```

If you have several variables of interest, you need to perform the areal interpolation separately or create a function to loop over several columns. 
In our example, we also need the interpolated area size, so we repeat the step.

```{r}

# Perform area-weighted interpolation for the area
areal_interpolation_matched <- sf::st_interpolate_aw(
  # Use municipality are data for interpolation
  municipalities["area_municipality"],
  # Target zip code areas for the interpolation
  zip_codes_valid,
  # Set to FALSE as area is not "extensive" (not purely additive)
  extensive = FALSE) |>
  # Combine interpolated results with original zip code data
  dplyr::bind_cols(
    zip_codes_valid |>
      # Drop geometry to avoid duplication issues in final output
      sf::st_drop_geometry() |>
      # Select only zip code and its inhabitants count
      dplyr::select(zip_code, inhabitants_zip_code)
  ) |>
  # Choose relevant columns for output
  dplyr::select(zip_code, area_municipality) |>
  # Ensure unique rows
  dplyr::distinct() |>
  # Drop geometry
  sf::st_drop_geometry() |>
  # Join with original
  dplyr::left_join(
    areal_interpolation_matched,
    .,
    by = "zip_code")

```

For advanced users, we would also like to point out the package [`areal`](https://github.com/chris-prener/areal){target="_blank"} that builds upon the `sf`package we use here and expands the functions for area-weighted interpolations [@prener_areal_2020].

## Assessing the matching methods


### Data preparation

We sampled nearly 1,000,000 points to simulate 1,000 hypothetical living locations within each ZIP code area (see [Data simulation](#data-simulation)). 
This allows us to assess the accuracy and consistency of the linking process and compare the three matching methods. 
The step is crucial because, for each simulated point, we know the exact municipality match, which provides the correct value for the number of inhabitants in the respondents' municipality.

It is important to note that we focus here on the differences in the number of inhabitants per municipality to evaluate the matching process because.
Specifically, if a ZIP code area overlaps with spatial clusters of municipalities that have a similar number of inhabitants - areas with high positive spatial autocorrelation - the differences will naturally be smaller, even if the matching methods yield different results.
Therefore, the evaluation of the linking process and ultimately the selection of the most appropriate method should consider the main variable(s) of interest. 

Let's calculate for each point location the difference between real and estimated number of inhabitants for each matching method (centroid matching, area matching and area-weighted interpolation).

```{r}

# 1. Difference with Centroid Matching ----
diff_real_centroid <-
  # Perform a left join between the real data and centroid-matched data
  dplyr::left_join(
    # The real inhabitants data joined to points
    points_with_real_inhabitants_municipality,
    # Centroid matched data (drop the geometry for non-spatial comparison)
    centroid_matched |>
      # Drop geometry column, keeping only tabular data
      sf::st_drop_geometry()
  ) |>
  # Create new columns for the differences
  dplyr::mutate(
    # Add a column indicating the type of comparison
    `Type` = "Difference with Centroid Matching",
    # Calculate the difference in inhabitants
    Difference = real_inhabitants_municipality - inhabitants_municipality
  )

# 2. Difference with Areal Matching ----
diff_real_areal <-
  # Left join between real inhabitants data and areal-matched data
  dplyr::left_join(
    # The real inhabitants data
    points_with_real_inhabitants_municipality,
    # Areal matched data (drop geometry for comparison)
    areal_matched |>
      sf::st_drop_geometry()
  ) |>
  # Create new columns
  dplyr::mutate(
    # Indicate the type of comparison
    `Type` = "Difference with Areal Matching",
    # Calculate difference between real and estimated inhabitants
    Difference = real_inhabitants_municipality - inhabitants_municipality
  )

# 3. Difference with Areal Interpolation Matching ----
diff_real_areal_interpolation <-
  # Join real inhabitants data with areal interpolation-matched data
  dplyr::left_join(
    # The real inhabitants data
    points_with_real_inhabitants_municipality,
    # Areal interpolation matched data (without geometry)
    areal_interpolation_matched |>
      sf::st_drop_geometry()
  ) |>
  # Create new columns
  dplyr::mutate(
    # Indicate the comparison type
    `Type` = "Difference with Areal Interpolation Matching",
    # Calculate the difference between real and interpolated inhabitants
    Difference = real_inhabitants_municipality - inhabitants_municipality
  ) |>
  # Ensure unique rows to avoid duplication
  dplyr::distinct()

```



In the next step, we combine and process the differences from three matching methods to calculate summary statistics, e.g. the median, mean, minimum, maximum, standard deviation, variance, and interquartile range.
We also calculate accuracy indicators, including binary indicators for exact matches and for differences within ±500 inhabitants. Finally, we calculate the proportion of exact and near-matches for each method.

```{r}

differences <-
  # Combine rows from all difference data frames into one
  dplyr::bind_rows(
    # Difference data from Centroid Matching
    diff_real_centroid,
    # Difference data from Areal Matching
    diff_real_areal,
    # Difference data from Areal Interpolation Matching
    diff_real_areal_interpolation
  ) |>
  # Drop spatial geometry for non-spatial analysis
  sf::st_drop_geometry() |>
  # Group by matching method type (Centroid, Areal, etc.)
  dplyr::group_by(Type) |>
  # Add indicators for accuracy evaluation
  dplyr::mutate(
    # Binary indicator if Difference is exactly 0 (perfect match)
    correct = ifelse(Difference == 0, 1, 0),
    # Indicator if Difference within ±500 inhabitants
    more_or_less_correct = ifelse(Difference > -500 & Difference < 500, 1, 0)
  ) |>
  # Summarize the data for each matching method
  dplyr::summarize(
    # Median of the Difference column
    median = median(Difference, na.rm = TRUE),
    # Mean of the Difference column
    mean = mean(Difference, na.rm = TRUE),
    # Minimum difference
    min = min(Difference, na.rm = TRUE),
    # Maximum difference
    max = max(Difference, na.rm = TRUE),
    # Standard deviation of the differences
    sd = sd(Difference, na.rm = TRUE),
    # Variance of differences divided by 1000 for scaling
    var1000 = var(Difference, na.rm = TRUE) / 1000,
    # Interquartile range (IQR) of differences
    iqr = IQR(Difference, na.rm = TRUE),
    # Proportion of exact matches (where Difference = 0)
    prop_correct = mean(correct, na.rm = TRUE),
    # Proportion of matches within ±500 inhabitants
    prop_more_or_less_correct = mean(more_or_less_correct, na.rm = TRUE)
  )

```

Lastly, we calculate the differences in estimated inhabitants between the different matching methods to better evaluate their performance. 
These calculations are stored as new variables, each labeled with the corresponding matching method.

```{r}
# Calculate the differences between inhabitants estimated by 
# Centroid Matching and Areal Matching 
diff_centroid_areal <-
  # Create a new tibble (data frame) to store the differences
  tibble::tibble(
    # Assign a label for the type of difference calculated
    `Type` = "Difference Centroid and Areal Matching",
    # Calculate the difference in inhabitants between the two matching methods
    Difference = centroid_matched$inhabitants_municipality -
      # Subtract Areal Matching inhabitants from Centroid Matching
      areal_matched$inhabitants_municipality
  )

# Calculate the differences between inhabitants estimated by Centroid Matching 
# and Areal Interpolation Matching
diff_centroid_interpolated <-
  # Create a new tibble to store the differences
  tibble::tibble(
    # Label for the difference type
    `Type` = "Difference Centroid and Areal Interpolation Matching",
    # Calculate the difference for this matching comparison
    Difference = centroid_matched$inhabitants_municipality -
      # Subtract Areal Interpolation inhabitants from Centroid Matching
      areal_interpolation_matched$inhabitants_municipality
  )

# Calculate the differences between inhabitants estimated by Areal Matching and 
# Areal Interpolation Matching
diff_areal_interpolated <-
  # Create a new tibble to store the differences
  tibble::tibble(
    # Label for the difference type
    `Type` = "Difference Areal and Areal Interpolation Matching",
    # Calculate the difference in inhabitants between the two methods
    Difference = areal_matched$inhabitants_municipality -
      # Subtract Areal Interpolation inhabitants from Areal Matching
      areal_interpolation_matched$inhabitants_municipality
  )

```


Since we are ultimately interested in how the matching compares at the ZIP code level, rather than at the (simulated) point level, we aggregate the population estimate differences at the ZIP code level. 
This process allows us to compare the results of the three different matching methods (centroid matching, areal matching, and areal interpolation matching) and identify which method produces the smallest population difference for each ZIP code.

```{r}

# Aggregate differences at the ZIP code level ----

diff_real_zip_code_aggregated <-
  diff_real_centroid |>
  # Group data by ZIP code
  dplyr::group_by(zip_code) |>
  # Drop geometry for non-spatial operations
  sf::st_drop_geometry() |>
  # Summarize by calculating mean difference
  dplyr::summarize(
    # Mean difference for centroid matching method
    mean_difference_centroid = mean(Difference, na.rm = TRUE)
  ) |>
  # Join with the areal matching differences
  dplyr::left_join(
    diff_real_areal |>
      dplyr::group_by(zip_code) |>
      sf::st_drop_geometry() |>
      dplyr::summarize(
        # Mean difference for areal matching method
        mean_difference_areal = mean(Difference, na.rm = TRUE)
      )
  ) |>
  # Join with the areal interpolation matching differences
  dplyr::left_join(
    diff_real_areal_interpolation |>
      dplyr::group_by(zip_code) |>
      sf::st_drop_geometry() |>
      dplyr::summarize(
        # Mean difference for areal interpolation
        mean_difference_areal_interpolation = mean(Difference, na.rm = TRUE)
      )
  )

# Calculate absolute mean differences ----
diff_real_zip_code_aggregated <-
  diff_real_zip_code_aggregated |>
  # Create columns for absolute mean differences
  dplyr::mutate(
    # Absolute value of mean differences (centroid method)
    mean_difference_centroid_abs = abs(mean_difference_centroid),
    # Absolute mean difference for areal method
    mean_difference_areal_abs = abs(mean_difference_areal),
    # Absolute mean difference for areal interpolation
    mean_difference_areal_interpolation_abs =
      abs(mean_difference_areal_interpolation)
  )

# Identify the method with the smallest absolute difference per ZIP code ----
diff_real_zip_code_aggregated <-
  diff_real_zip_code_aggregated |>
  # Use bind_cols to add the method with minimum difference
  dplyr::bind_cols(
    # Find column name with minimum absolute difference
    min_method = names(diff_real_zip_code_aggregated[-c(1:4)])[
      # Inverse to get min absolute values, excluding grouping cols
      max.col(-diff_real_zip_code_aggregated[-c(1:4)])
    ]
  ) |>
  # Adjust cases with ties for consistency
  dplyr::mutate(
    min_method = ifelse(
      # Tie condition check
      mean_difference_centroid_abs == mean_difference_areal_interpolation_abs &
        mean_difference_centroid_abs == mean_difference_areal_abs,
      # Set to "centroid" if all methods tie
      "mean_difference_centroid_abs",
      min_method
    )
  ) |>
  # Remove intermediate absolute difference columns
  dplyr::select(-contains("abs"))

# Join with ZIP code geometry and assign codes for each method ----
diff_real_zip_code_aggregated <-
  diff_real_zip_code_aggregated |>
  # Join back with ZIP code spatial data
  dplyr::left_join(zip_codes) |>
  # Convert back to an sf object for spatial analysis
  sf::st_as_sf() |>
  # Create a code for each method for easier visualization
  dplyr::mutate(
    min_method_code = dplyr::case_when(
      # Assign code 1 for centroid method
      min_method == "mean_difference_centroid_abs" ~ 1,
      # Code 2 for areal method
      min_method == "mean_difference_areal_abs" ~ 2,
      # Code 3 for areal interpolation method
      min_method == "mean_difference_areal_interpolation_abs" ~ 3
    )
  ) |>
  dplyr::mutate(
    difference = dplyr::case_when(
      min_method == "mean_difference_centroid_abs" ~ mean_difference_centroid,
      min_method == "mean_difference_areal_abs" ~ mean_difference_areal,
      min_method == "mean_difference_areal_interpolation_abs" ~
        mean_difference_areal_interpolation,
      TRUE ~ NA
    )
  )

```

### Accuracy of matching methods

First, we plot the differences between our variable of interest, number of inhabitants, for each matching method (centroid matching, areal matching, and area-weighted interpolation) compared with the real estimate based on the point locations. 
The histogram shows that all three matching methods differ from the real estimate, but on average, they exhibit similar levels of accuracy.

```{r}


# Create a histogram to visualize the distribution of differences for each matching method ----
difference_histogram <-
  # Combine differences from all matching methods into a single data frame
  dplyr::bind_rows(
    # Centroid method differences
    diff_real_centroid,
    # Areal method differences
    diff_real_areal,
    # Areal interpolation method differences
    diff_real_areal_interpolation
  ) |>
  # Remove spatial geometry for visualization purposes
  sf::st_drop_geometry() |>
  # Initialize ggplot with Difference as x-axis variable
  ggplot2::ggplot(aes(x = Difference)) +
  # Add histogram with 10 bins
  ggplot2::geom_histogram(bins = 10) +
  # Limit x-axis for a clear view of the distribution
  ggplot2::xlim(-2500000, 2500000) +
  # Create a separate histogram for each matching method type
  ggplot2::facet_wrap(~`Type`, labeller = ggplot2::label_wrap_gen(width = 30)) + 
   ggplot2::theme_minimal() +  # Apply minimal theme
  # Wrap the facet labels
  ggplot2::theme(
    strip.text = ggplot2::element_text(size = 8) # Adjust font size
  )

difference_histogram


```

### Consistency of matching methods

In the next step, we aim to assess whether the matching methods are consistent. 
To achieve this, we plot the densities of the scaled differences across the different matching methods.
By examining the pairwise comparisons of the matching methods, we can determine which methods produce more consistent results (indicated by smaller variance in the differences) and which methods tend to diverge more (indicated by larger variance).

Overall, we observe that the differences between the matching methods are relatively small, as all three density plots are centered around zero. This suggests that the compared matching methods produce similar results when linking the data.

While the location of the peak tells us about the magnitude of the differences, the width of the distribution reflects the variability of those differences.
The first plot indicates small and consistent differences, implying that the centroid and areal matching perform similarly across the data.
In contrast, when comparing centroid matching with the two other matching methods, we observe two slightly wider density plots.
The plots suggest that the number of inhabitants estimated by centroid matching  differs - at least in some areas - from the estimation of the areal matching and areal interpolation.

```{r}

# Create density plots to compare scaled differences across matching methods ----
difference_densities <-
  # Combine pairwise differences for each matching method
  dplyr::bind_rows(
    # Centroid and Areal method differences
    diff_centroid_areal,
    # Centroid and Interpolated method differences
    diff_centroid_interpolated,
    # Areal and Interpolated method differences
    diff_areal_interpolated
  ) |>
  # Group data by matching method type
  dplyr::group_by(`Type`) |>
  # Scale differences for comparability across methods
  dplyr::mutate(Difference = scale(Difference)) |>
  # Initialize ggplot with Difference as x-axis variable
  ggplot2::ggplot(aes(x = Difference)) +
  # Create density plot for distribution visualization
  ggplot2::geom_density() +
  # Set x-axis limits to focus on main data range
  ggplot2::xlim(-10, 10) +
  # Separate density plots by type of matching method
  ggplot2::facet_wrap(~`Type`, labeller = ggplot2::label_wrap_gen(width = 30)) +  
  ggplot2::theme_minimal() +  # Apply minimal theme
  # Wrap the facet labels
  ggplot2::theme(
    strip.text = ggplot2::element_text(size = 8) # Adjust font size
  )

difference_densities


```

To further assess the performance of the matching methods in different areas, we create a map that allows us to identify which method is most consistently successful at minimizing population differences across ZIP codes.
The map shows which method results in the smallest difference when compared to the actual population values, as estimated via the simulated points.

This map allows us to identify which method is most consistently successful at minimizing population differences across ZIP codes.
We observe that different matching methods perform better in certain regions. 
For example, in densely populated areas like cities, ZIP codes are often contained within a single municipality, so centroid matching performs particularly well in these areas. On average, areal interpolation, which takes multiple factors into account, performs the best across the entire study area.
However, this method requires (all) spatial variables of interest to be available for the matching process, which may not always be the case and can be very resource-intensive.

```{r}


# Create a map to visualize which method minimizes the difference for each ZIP code ----
diff_real_zip_code_aggregated_map <-
  # Initialize ggplot for map
  ggplot2::ggplot() +
  # Add map layers with fill representing the best matching method
  ggplot2::geom_sf(
    data = diff_real_zip_code_aggregated, aes(fill = min_method), lwd = 0
  ) +
  # Use viridis color scale to differentiate methods
  scale_fill_viridis_d() +
  ggplot2::theme_minimal()   # Apply minimal theme

# Save the map plot to a file ----
ggplot2::ggsave(
  # Specify file name and output path
  "./test/diff_real_zip_code_aggregated_map.png",
  # Reference the map plot
  diff_real_zip_code_aggregated_map,
  # High resolution for clear saved image
  dpi = 600
)

diff_real_zip_code_aggregated_map_differences <-
  # Initialize ggplot for map
  ggplot2::ggplot() +
  ggplot2::geom_sf(data = zip_codes, lwd = .01) +
  # Add map layers with fill representing the difference
  ggplot2::geom_sf(
    data = diff_real_zip_code_aggregated, aes(fill = difference), lwd = 0
  ) +
  # Use viridis color scale to differentiate differences
  scale_fill_viridis_c() +
  ggplot2::facet_wrap(~min_method) +
  ggplot2::theme_minimal()   # Apply minimal theme

diff_real_zip_code_aggregated_map


```

### Effect of matching method on bivariate correlation  

Finally, we aim to assess whether and how the different matching methods influence the answers to our research question. 
Rather than developing an elaborate model, we take a more straightforward approach by examining the correlation between our variable of interest—the AfD rating — and the matched context variable.
Since we are not interested in the absolute value of inhabitants but the population density, we calculate the population density accordingly.

```{r}

# Linking the datasets
linked_survey <-
  survey |>
  # Left join with centroid_matched, dropping geometry information
  dplyr::left_join(centroid_matched |>
                     sf::st_drop_geometry(), by = "zip_code") |>
  # Rename columns for clarity
  dplyr::rename(cent_inhabitants_mun = inhabitants_municipality,
                cent_area_mun = area_municipality) |>
  # Left join with areal_matched, dropping specific columns and geometry
  dplyr::left_join(areal_matched |>
                     sf::st_drop_geometry() |>
                     dplyr::select(-inhabitants_zip_code), by = "zip_code") |>
  # Rename columns for clarity
  dplyr::rename(areal_inhabitants_mun = inhabitants_municipality,
                areal_area_mun = area_municipality) |>
  # Left join with areal_interpolation_matched, dropping specific columns and geometry
  dplyr::left_join(areal_interpolation_matched |>
                     sf::st_drop_geometry() |>
                     dplyr::select(-inhabitants_zip_code), by = "zip_code") |>
  # Rename columns for clarity
  dplyr::rename(interpolation_inhabitants_mun = inhabitants_municipality,
                interpolation_area_mun = area_municipality)  |>
  # Calculate population densities based on the joined data
  dplyr::mutate(
    cent_pop_dens = cent_inhabitants_mun / cent_area_mun,
    areal_pop_dens = areal_inhabitants_mun / areal_area_mun,
    interpolation_pop_dens = interpolation_inhabitants_mun / interpolation_area_mun
  )


# Calculate pairwise correlations with 'afd_rating' and columns containing 'population density'
correlations <- linked_survey |>
  dplyr::select(afd_rating, contains("pop_dens")) |>
  corrr::correlate(use = "pairwise.complete.obs")

# Melt the correlation matrix into a long format
correlations_long <- correlations %>%
  # Pivot the data into long format
  tidyr::pivot_longer(
    cols = -term,        # Keep the 'term' column as is
    names_to = "y",      # New column 'y' for variable names
    values_to = "value"  # New column 'value' for correlation values
  ) %>%
  # Rename 'term' column to 'x' for consistency
  dplyr::rename(x = term) %>%
  # Filter for correlations with 'afd_rating' and not with itself
  dplyr::filter(x == "afd_rating" & x != y)

# Create a custom label mapping (this is just an example)
custom_labels <- c(
  "zip_pop_dens" = "ZIP code",
  "cent_pop_dens" = "Centroid matching",
  "areal_pop_dens" = "Areal matching",
  "interpolation_pop_dens" = "Areal interpolation"
)

# Plot the correlation coefficients as a bar chart
ggplot2::ggplot(correlations_long, ggplot2::aes(x = y, y = value, fill = value)) +
  ggplot2::geom_bar(stat = "identity", show.legend = FALSE) +
  viridis::scale_fill_viridis(option = "D") +  # Apply viridis color scale
  ggplot2::labs(title = "Correlation of Population Density Variables with AFD Rating",
       x = "Population Density Variables",
       y = "Correlation with AFD Rating") +
  # Customize the x-axis labels with custom names
  ggplot2::scale_x_discrete(labels = custom_labels) +
  ggplot2::theme_minimal() +
  ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 45, hjust = 1))


```

As shown in the plot, there are only marginal differences in the correlation coefficients between the various matching methods. 
This suggests that, as previously indicated, the linking processes perform similarly across these matching methods.

From a substantive perspective, the correlations are consistently small and close to zero. 
Specifically, we observe that living in areas with higher population density is weakly correlated with lower AfD ratings.
In other words, individuals living in more densely populated areas tend to support the AfD less — although the strength of this bivariate relationship is marginal.

In a more complex model, especially when controlling for individual-level characteristics, the relationship may change. 
Not only could the strength of the correlation differ, but the differences between matching methods could also become more pronounced. 
Similarly, the results may vary when matching other spatial units. For 
example, the analysis could yield a different picture when matching electoral districts, which by design have a similar number of inhabitants but vary significantly in area size.
In conclusion, we recommend, whenever computationally and data-wise feasible, checking the robustness of results across different matching methods.


# Conclusion and recommendations

In this tool application, we evaluated three matching methods—centroid matching, areal matching, and area-weighted interpolation—to link survey data with geospatial context variables for overlapping but incongruent spatial units. Our goal was to assess the accuracy of each method, compare their consistency, and ultimately evaluate their impact on our research question: Is there an association between population density and far-right political preferences, and is this association robust?

Our findings indicate that while all three methods show slight deviations from the actual population estimates (derived from simulated point locations), they generally produce similar results in terms of accuracy. In terms of consistency, subtle differences were observed: centroid matching performed slightly better in densely populated areas, while area-weighted interpolation provided the most reliable results on average across the entire study area. Despite these variations, the choice of matching method had minimal impact on the research outcomes, such as the correlation between AfD ratings and population density.

Key Takeaways:

- **Centroid matching**: Although it does not account for the area of overlap or local population variations, centroid matching performs surprisingly well in certain cases when linking ZIP code areas and municipalities in Germany. It offers a good trade-off between efficiency and quality, particularly in areas where the target spatial unit is known to lie within the source unit.
- **Areal matching**: We recommend investigating the ambiguity in the matching process, such as the number of source units with which the target unit overlaps, the area shares of the overlap, and the spatial clusters of ambiguous cases. Additional robustness checks and flagging of ambiguous cases may be necessary depending on the spatial units involved, particularly when the target unit overlaps significantly with (many) source units.
- **Area-weighted interpolation**: This method addresses the limitation of areal matching by interpolating data from all overlapping units. However, it has the highest computational requirements and data needs. In our example, it did not clearly outperform the other methods but provided the most consistent results across all spatial units.
- **Robustness Checks**: To ensure the reliability of findings, researchers should assess the robustness of their matching method by comparing it with at least one alternative approach. This is especially important when exploring how contextual information, such as population density or other geographical factors, is associated with individual-level attitudes or behavior. Robustness checks help confirm that the observed relationships are not influenced by the matching process but show robust associations.

Still, choosing the most suitable matching method is just one part of the equation.
While the technical aspects of matching can improve the accuracy and consistency of linking processes, they do not replace the need for careful consideration of the conceptual underpinnings: what is the geographic unit that best captures the relevant social or environmental context for the research question under investigation?
For example, the uncertain geographic context problem (UGCoP) [@kwan_uncertain_2012] highlights the importance of carefully considering spatial boundaries in social research. 
The choice of neighborhood unit can significantly affect the results, particularly when studying phenomena such as support for far-right parties.
The methods evaluated in this application provide different ways of addressing the challenges of mismatched or incongruent boundaries when defining a "neighborhood." 
Still, researchers must remain aware that the linking process should be guided by a clear understanding of how these spatial definitions shape the conclusions drawn from integrated survey and geospatial data.


# References

