Road Density & The Roman Empire
Joona Kähkönen
2024-05-24
Introduction
This analysis explores the borders of the Roman Empire and their relationship with various socioeconomic factors and geographical features.
Libraries
# Load necessary libraries
library(cawd)
library(ggplot2)
library(rnaturalearthdata)
library(tidyverse)
library(sf)
library(dplyr)
library(terra)
library(raster)
library(elevatr)
library(rdrobust)
library(geodata)
library(lwgeom)
library(stars)
library(starsExtra)
library(modelsummary)
library(webshot2)
library(kableExtra)Download and Process Maps
# Download map of Europe
europe <- rnaturalearthdata::sovereignty50 |>
st_as_sf() |>
filter(continent == "Europe")
# Download map of Roman Empire
roman_empire <- cawd::awmc.roman.empire.117.sp |> st_as_sf()
# Extract Border of the Roman Empire
boundary_re <- st_boundary(roman_empire)
# Extract Northern Border of the Roman Empire
north_countries <- europe |>
filter(is.element(sovereignt, c("Netherlands", "Germany", "Croatia",
"Austria", "Hungary", "Slovakia",
"Ukraine", "Romania", "Bulgaria",
"Republic of Serbia")))
# Filter the LINESTRING(s) with the maximum Y-coordinate
northernmost_lines <- st_intersection(boundary_re, north_countries)## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries# Filter out south border of Croatia
coords <- matrix(c(15, 46, 20, 42, 10, 42, 10, 46, 15, 46), ncol = 2, byrow = TRUE)
polygon <- st_polygon(list(coords)) |> st_sfc(crs = 4326)
northernmost_border <- northernmost_lines |> st_difference(polygon)## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries# Build a 200km Buffer Around Northern Border and Intersection With Europe
sf_use_s2(TRUE)
buffer_northern_border <- st_buffer(northernmost_border, 200000)
sf_use_s2(FALSE)## Spherical geometry (s2) switched offborder_area <- europe |> st_intersection(buffer_northern_border)## although coordinates are longitude/latitude, st_intersection assumes that they
## are planar## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries# Visualize the area
border_area |>
ggplot() + geom_sf()
Build Grid and Save Raster
# Build a grid on this basis
grid <- st_make_grid(border_area, n = 100)
filter <- st_within(grid, st_union(border_area)) |> lengths() > 0## although coordinates are longitude/latitude, st_union assumes that they are
## planar## although coordinates are longitude/latitude, st_within assumes that they are
## planarfinal_grid <- grid[filter]
# Turn grid into an empty raster and save
raster_grid <- final_grid |> stars::st_as_stars() ## although coordinates are longitude/latitude, st_intersects assumes that they
## are planar## although coordinates are longitude/latitude, st_intersects assumes that they
## are planarstars::write_stars(raster_grid, 'raster_roman_empire.tiff', overwrite=TRUE)Load and Provess Road Density Data
# Load the data
road_density_raster <- rast("C:/Users/jookk/OneDrive/Työpöytä/R_Thesis/grip4_tp1_dens_m_km2.asc")
raster_re <- rast("C:/Users/jookk/OneDrive/Työpöytä/R_Thesis/raster_re.tiff")
# Check if rasters are loaded by examining their properties
print(raster_re)## class : SpatRaster
## dimensions : 164, 398, 1 (nrow, ncol, nlyr)
## resolution : 0.07593909, 0.07593909 (x, y)
## extent : 0.8761778, 31.09993, 41.10402, 53.55803 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (EPSG:4326)
## source : raster_re.tiff
## name : raster_reprint(road_density_raster)## class : SpatRaster
## dimensions : 2160, 4320, 1 (nrow, ncol, nlyr)
## resolution : 0.08333333, 0.08333333 (x, y)
## extent : -180, 180, -90, 90 (xmin, xmax, ymin, ymax)
## coord. ref. : lon/lat WGS 84 (CRS84) (OGC:CRS84)
## source : grip4_tp1_dens_m_km2.asc
## name : grip4_tp1_dens_m_km2if (crs(raster_re) != crs(road_density_raster)) {
road_density_raster <- project(road_density_raster, crs(raster_re))
}
# Crop the global road density raster to the extent of the raster_re to reduce computation load
road_density_cropped <- crop(road_density_raster, ext(raster_re))
# Resample the cropped road density raster to match the grid resolution of raster_re
road_density_resampled <- resample(road_density_cropped, raster_re, method = "bilinear")
# Perform the overlay operation by multiplying corresponding cells
overlay_result <- raster_re * road_density_resampled
# Plot the result to verify the overlay
plot(overlay_result, main = "Overlay of Road Density with Grid")
# Save the overlay result to a new file
writeRaster(overlay_result, "overlay_density_grid.tiff", overwrite=TRUE)Ensure CRS, Create Intersection & The RE Dummy
# Ensure both have the same CRS
roman_empire <- st_transform(roman_empire, st_crs(raster_re))
# Create an intersection
intersection <- st_intersection(roman_empire, st_as_sfc(st_bbox(raster_re)))## although coordinates are longitude/latitude, st_intersection assumes that they
## are planar## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries# Rasterize the intersection using the template
roman_empire_raster <- rasterize(intersection, raster_re, field = 1, fun = max, background = 0)
roman_empire_clipped <- mask(roman_empire_raster, raster_re)
# Plot to verify the clipped raster
plot(roman_empire_clipped, main = "Clipped Roman Empire Dummy Variable")
Load Socioeconomic Data
# Load socioeconomic data
gdp_data <- rast("C:/Users/jookk/OneDrive/Työpöytä/R_Thesis/GDP_PPP_1990_2015_5arcmin_v2.nc")
hdi_data <- rast("C:/Users/jookk/OneDrive/Työpöytä/R_Thesis/HDI_1990_2015_v2.nc")
# Resample and align GDP and HDI data with the road density raster
gdp_data_aligned <- resample(gdp_data, road_density_resampled, method = "bilinear")
hdi_data_aligned <- resample(hdi_data, road_density_resampled, method = "bilinear")
# Create dummy variables for areas within the Roman Empire and mask the GDP and HDI data
gdp_masked <- mask(gdp_data_aligned, roman_empire_clipped)
hdi_masked <- mask(hdi_data_aligned, roman_empire_clipped)
# Extract data for a specific year, I use 2015
gdp_2015 <- gdp_data_aligned[["GDP_PPP_26"]]
hdi_2015 <- hdi_data_aligned[["HDI_26"]]
# Align, mask, and prepare specific layers for analysis
gdp_2015_aligned <- resample(gdp_2015, road_density_resampled, method = "bilinear")
hdi_2015_aligned <- resample(hdi_2015, road_density_resampled, method = "bilinear")
gdp_masked_2015 <- mask(gdp_2015_aligned, roman_empire_clipped)
hdi_masked_2015 <- mask(hdi_2015_aligned, roman_empire_clipped)Load and Process Elevation Data
# Get the extent of raster_re to define the area for elevation data extraction
raster_extent <- ext(raster_re)
print(raster_extent)## SpatExtent : 0.876177769909369, 31.0999339127773, 41.1040216726612, 53.5580317415816 (xmin, xmax, ymin, ymax)# Define a grid over the extent with suitable resolution
points <- expand.grid(x = seq(0.876177769909369, 31.0999339127773, length.out = 100),
y = seq(41.1040216726612, 53.5580317415816, length.out = 100))
# Get elevation data for these points
elevation_data <- get_elev_point(points, prj = "+proj=longlat +datum=WGS84", src = "aws")## Mosaicing & Projecting## Note: Elevation units are in meters# Convert the 'sf' object to a regular dataframe with x, y, and elevation
coordinates <- st_coordinates(elevation_data)
elevation_df <- data.frame(
x = coordinates[, "X"],
y = coordinates[, "Y"],
elevation = elevation_data$elevation
)
# Convert the elevation data frame into a raster using the extent and resolution of raster_re
elevation_raster <- rasterFromXYZ(elevation_df, crs = crs(raster_re))
elevation_raster_terra <- terra::rast(elevation_raster)
elevation_aligned <- terra::resample(elevation_raster_terra, road_density_resampled, method = "bilinear")
# Crop elevation_aligned to match the extent of raster_re
elevation_cropped <- crop(elevation_aligned, ext(raster_re))
# Align it if there is any slight misalignment
elevation_cropped <- resample(elevation_cropped, raster_re, method = "bilinear")Load and Process Population Data
# Load the TIF file
pop_density_raster <- rast("C:/Users/jookk/OneDrive/Työpöytä/R_Thesis/gpw_v4_population_density_rev11_2015_30_sec.tif")
# Reproject the population density raster
if (crs(pop_density_raster) != crs(raster_re)) {
pop_density_raster <- projectRaster(pop_density_raster, crs = crs(raster_re))
}
# Crop the population density raster to the extent of the study area
pop_density_cropped <- crop(pop_density_raster, ext(raster_re))
# Resample the population density raster to match the resolution of raster_re
pop_density_resampled <- resample(pop_density_cropped, raster_re, method = "bilinear")
# Mask the resampled population density raster
pop_density_masked <- mask(pop_density_resampled, raster_re)
# Plot the masked population density raster
plot(pop_density_masked, main = "Masked Population Density 2015")
Filter Data for Analysis
# Extract the values from relevant rasters
road_density_values <- values(road_density_resampled)
gdp_values <- values(gdp_masked_2015)
hdi_values <- values(hdi_masked_2015)
roman_empire_values <- values(roman_empire_clipped)
pop_density_values <- values(pop_density_masked)
elevation_values <- values(elevation_cropped)
# Ensure all values are of the same length
if (!all(length(road_density_values) == length(gdp_values),
length(gdp_values) == length(hdi_values),
length(hdi_values) == length(roman_empire_values),
length(roman_empire_values) == length(pop_density_values),
length(pop_density_values) == length(elevation_values))) {
stop("Length of values extracted from rasters do not match.")
}
# Create a logical mask to filter out NA values across all variables
valid_mask <- !is.na(road_density_values) & !is.na(gdp_values) &
!is.na(hdi_values) & !is.na(roman_empire_values) &
!is.na(pop_density_values) & !is.na(elevation_values)
# Apply the mask to filter out NA values
filtered_road_density_values <- road_density_values[valid_mask]
filtered_gdp_values <- gdp_values[valid_mask]
filtered_hdi_values <- hdi_values[valid_mask]
filtered_roman_empire_values <- roman_empire_values[valid_mask]
filtered_pop_density_values <- pop_density_values[valid_mask]
filtered_elevation_values <- elevation_values[valid_mask]
# Create a new dataframe with filtered data
analysis_data <- data.frame(
RoadDensity = filtered_road_density_values,
GDP = filtered_gdp_values,
HDI = filtered_hdi_values,
RomanEmpire = filtered_roman_empire_values,
PopDensity2015 = filtered_pop_density_values,
Elevation = filtered_elevation_values
)
# Remove any rows with missing values
analysis_data <- na.omit(analysis_data)Create the Running Variable
# Calculate the distance to the nearest Roman Empire border
nearest_feat <- starsExtra::dist_to_nearest(raster_grid, northernmost_border)## lines or polygons## | | | 0% | | | 1% | |= | 1% | |= | 2% | |== | 2% | |== | 3% | |== | 4% | |=== | 4% | |=== | 5% | |==== | 5% | |==== | 6% | |===== | 6% | |===== | 7% | |===== | 8% | |====== | 8% | |====== | 9% | |======= | 9% | |======= | 10% | |======= | 11% | |======== | 11% | |======== | 12% | |========= | 12% | |========= | 13% | |========= | 14% | |========== | 14% | |========== | 15% | |=========== | 15% | |=========== | 16% | |============ | 16% | |============ | 17% | |============ | 18% | |============= | 18% | |============= | 19% | |============== | 19% | |============== | 20% | |============== | 21% | |=============== | 21% | |=============== | 22% | |================ | 22% | |================ | 23% | |================ | 24% | |================= | 24% | |================= | 25% | |================== | 25% | |================== | 26% | |=================== | 26% | |=================== | 27% | |=================== | 28% | |==================== | 28% | |==================== | 29% | |===================== | 29% | |===================== | 30% | |===================== | 31% | |====================== | 31% | |====================== | 32% | |======================= | 32% | |======================= | 33% | |======================= | 34% | |======================== | 34% | |======================== | 35% | |========================= | 35% | |========================= | 36% | |========================== | 36% | |========================== | 37% | |========================== | 38% | |=========================== | 38% | |=========================== | 39% | |============================ | 39% | |============================ | 40% | |============================ | 41% | |============================= | 41% | |============================= | 42% | |============================== | 42% | |============================== | 43% | |============================== | 44% | |=============================== | 44% | |=============================== | 45% | |================================ | 45% | |================================ | 46% | |================================= | 46% | |================================= | 47% | |================================= | 48% | |================================== | 48% | |================================== | 49% | |=================================== | 49% | |=================================== | 50% | |=================================== | 51% | |==================================== | 51% | |==================================== | 52% | |===================================== | 52% | |===================================== | 53% | |===================================== | 54% | |====================================== | 54% | |====================================== | 55% | |======================================= | 55% | |======================================= | 56% | |======================================== | 56% | |======================================== | 57% | |======================================== | 58% | |========================================= | 58% | |========================================= | 59% | |========================================== | 59% | |========================================== | 60% | |========================================== | 61% | |=========================================== | 61% | |=========================================== | 62% | |============================================ | 62% | |============================================ | 63% | |============================================ | 64% | |============================================= | 64% | |============================================= | 65% | |============================================== | 65% | |============================================== | 66% | |=============================================== | 66% | |=============================================== | 67% | |=============================================== | 68% | |================================================ | 68% | |================================================ | 69% | |================================================= | 69% | |================================================= | 70% | |================================================= | 71% | |================================================== | 71% | |================================================== | 72% | |=================================================== | 72% | |=================================================== | 73% | |=================================================== | 74% | |==================================================== | 74% | |==================================================== | 75% | |===================================================== | 75% | |===================================================== | 76% | |====================================================== | 76% | |====================================================== | 77% | |====================================================== | 78% | |======================================================= | 78% | |======================================================= | 79% | |======================================================== | 79% | |======================================================== | 80% | |======================================================== | 81% | |========================================================= | 81% | |========================================================= | 82% | |========================================================== | 82% | |========================================================== | 83% | |========================================================== | 84% | |=========================================================== | 84% | |=========================================================== | 85% | |============================================================ | 85% | |============================================================ | 86% | |============================================================= | 86% | |============================================================= | 87% | |============================================================= | 88% | |============================================================== | 88% | |============================================================== | 89% | |=============================================================== | 89% | |=============================================================== | 90% | |=============================================================== | 91% | |================================================================ | 91% | |================================================================ | 92% | |================================================================= | 92% | |================================================================= | 93% | |================================================================= | 94% | |================================================================== | 94% | |================================================================== | 95% | |=================================================================== | 95% | |=================================================================== | 96% | |==================================================================== | 96% | |==================================================================== | 97% | |==================================================================== | 98% | |===================================================================== | 98% | |===================================================================== | 99% | |======================================================================| 99% | |======================================================================| 100%final_raster_grid <- raster_grid[raster_grid['values'] ==1]
final_raster_grid$distance <- nearest_feat$d
# Create indicator for being inside the Roman Empire
grid_sf <- st_as_sfc(final_raster_grid, as_points = TRUE)
logical <- as.logical(st_within(grid_sf, roman_empire))## although coordinates are longitude/latitude, st_within assumes that they are
## planarfinal_raster_grid <- final_raster_grid |>
mutate(inside_re = case_when(
!is.na(logical) & values == 1 ~ 1,
is.na(logical) & values == 1 ~ 0,
TRUE ~ NA
))
# Create the running variable
final_raster_grid <- final_raster_grid |>
mutate(running = if_else(inside_re == 1, distance, -distance))
# Plot the running variable
ggplot() +
geom_stars(data = final_raster_grid['running']) +
scale_fill_viridis_c()
Get the Coordinates
# Convert final_raster_grid to a data frame with coordinates
final_raster_grid_df <- as.data.frame(final_raster_grid, xy = TRUE)
# Add the distance to boundary variable to the data frame
final_raster_grid_df$DistanceToBoundary <- final_raster_grid_df$running
# Ensure analysis_data has the coordinates from the rasters
coordinates_df <- as.data.frame(road_density_resampled, xy = TRUE)
# Add values from other rasters to the coordinates_df
road_density_df <- as.data.frame(road_density_resampled, xy = TRUE)
gdp_df <- as.data.frame(gdp_masked_2015, xy = TRUE)
hdi_df <- as.data.frame(hdi_masked_2015, xy = TRUE)
roman_empire_df <- as.data.frame(roman_empire_clipped, xy = TRUE)
pop_density_df <- as.data.frame(pop_density_masked, xy = TRUE)
elevation_df <- as.data.frame(elevation_cropped, xy = TRUE)Merge the Data Frames
# Rename columns for clarity
names(road_density_df)[3] <- "RoadDensity"
names(gdp_df)[3] <- "GDP"
names(hdi_df)[3] <- "HDI"
names(roman_empire_df)[3] <- "RomanEmpire"
names(pop_density_df)[3] <- "PopDensity"
names(elevation_df)[3] <- "Elevation"
# Merge all data frames by coordinates
analysis_data <- coordinates_df %>%
left_join(road_density_df, by = c("x", "y")) %>%
left_join(gdp_df, by = c("x", "y")) %>%
left_join(hdi_df, by = c("x", "y")) %>%
left_join(roman_empire_df, by = c("x", "y")) %>%
left_join(pop_density_df, by = c("x", "y")) %>%
left_join(elevation_df, by = c("x", "y"))
# Merge analysis_data with final_raster_grid_df using coordinates
analysis_data <- analysis_data %>%
left_join(final_raster_grid_df, by = c("x", "y"))
# Remove rows with NA values
analysis_data <- na.omit(analysis_data)Clean the dataset and Create Logs
# Ensure no NA/NaN/Inf values in covariates and outcome variable
analysis_data <- analysis_data %>%
filter(!is.na(RoadDensity) & !is.na(DistanceToBoundary) & !is.na(GDP) & !is.na(HDI) & !is.na(Elevation) & !is.na(PopDensity) &
is.finite(RoadDensity) & is.finite(DistanceToBoundary) & is.finite(GDP) & is.finite(HDI) & is.finite(Elevation) & is.finite(PopDensity))
# Ensure positive and non-zero values for log transformation
analysis_data <- analysis_data %>%
filter(GDP > 0, HDI > 0, PopDensity > 0)
# Log transform the covariates and check for any problematic values
analysis_data <- analysis_data %>%
mutate(LogGDP = log(GDP),
LogHDI = log(HDI),
LogPopDensity = log(PopDensity))Analysis
Perform the Linear Regression Analysis
# Perform the Linear Regression
linear_reg <- lm(RoadDensity ~ RomanEmpire + GDP + HDI + Elevation + PopDensity, data = analysis_data)
# Summarize the results
summary(linear_reg)##
## Call:
## lm(formula = RoadDensity ~ RomanEmpire + GDP + HDI + Elevation +
## PopDensity, data = analysis_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1209.26 -34.55 -14.98 16.62 679.48
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.717e+02 7.807e+00 -21.993 < 2e-16 ***
## RomanEmpire 1.278e+01 1.198e+00 10.667 < 2e-16 ***
## GDP 9.471e-09 1.341e-09 7.062 1.72e-12 ***
## HDI 2.331e+02 9.283e+00 25.106 < 2e-16 ***
## Elevation -3.148e-02 1.900e-03 -16.572 < 2e-16 ***
## PopDensity 1.506e-01 3.873e-03 38.885 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 71.49 on 14936 degrees of freedom
## Multiple R-squared: 0.4677, Adjusted R-squared: 0.4675
## F-statistic: 2625 on 5 and 14936 DF, p-value: < 2.2e-16Perform the RDD analysis
# Perform RDD analysis using the rdrobust package
rdd_result_mserd <- rdrobust(y = analysis_data$RoadDensity,
x = analysis_data$DistanceToBoundary,
c = 0, # threshold for treatment assignment
covs = cbind(analysis_data$LogGDP, analysis_data$LogHDI, analysis_data$Elevation, analysis_data$LogPopDensity),
covs_drop = TRUE, # Drop covariates that are constant within each side of the cut-off
p = 1, # order of polynomial
kernel = "triangular") # kernel type
# Summarize the results
summary(rdd_result_mserd)## Covariate-adjusted Sharp RD estimates using local polynomial regression.
##
## Number of Obs. 14942
## BW type mserd
## Kernel Triangular
## VCE method NN
##
## Number of Obs. 7847 7095
## Eff. Number of Obs. 2951 2961
## Order est. (p) 1 1
## Order bias (q) 2 2
## BW est. (h) 64053.158 64053.158
## BW bias (b) 105645.627 105645.627
## rho (h/b) 0.606 0.606
## Unique Obs. 7846 7095
##
## =============================================================================
## Method Coef. Std. Err. z P>|z| [ 95% C.I. ]
## =============================================================================
## Conventional 21.656 6.138 3.528 0.000 [9.626 , 33.686]
## Robust - - 3.385 0.001 [10.264 , 38.501]
## =============================================================================Perform Robustness checks with different bandwidths using bwselect
# Compute the bandwidth using the msetwo method from the rdrobust package
rdd_result_msetwo <- rdrobust(y = analysis_data$RoadDensity,
x = analysis_data$DistanceToBoundary,
c = 0,
covs = cbind(analysis_data$LogGDP, analysis_data$LogHDI, analysis_data$Elevation, analysis_data$LogPopDensity),
bwselect = "msetwo")
# Summarize the results
summary(rdd_result_msetwo)## Covariate-adjusted Sharp RD estimates using local polynomial regression.
##
## Number of Obs. 14942
## BW type msetwo
## Kernel Triangular
## VCE method NN
##
## Number of Obs. 7847 7095
## Eff. Number of Obs. 3512 2205
## Order est. (p) 1 1
## Order bias (q) 2 2
## BW est. (h) 77067.969 45517.391
## BW bias (b) 109520.868 82616.052
## rho (h/b) 0.704 0.551
## Unique Obs. 7846 7095
##
## =============================================================================
## Method Coef. Std. Err. z P>|z| [ 95% C.I. ]
## =============================================================================
## Conventional 21.419 6.438 3.327 0.001 [8.800 , 34.037]
## Robust - - 3.140 0.002 [8.968 , 38.756]
## =============================================================================# Compute the bandwidth using the msesum method from the rdrobust package
rdd_result_msesum <- rdrobust(y = analysis_data$RoadDensity,
x = analysis_data$DistanceToBoundary,
c = 0,
covs = cbind(analysis_data$LogGDP, analysis_data$LogHDI, analysis_data$Elevation, analysis_data$LogPopDensity),
bwselect = "msesum")
# Summarize the results
summary(rdd_result_msesum)## Covariate-adjusted Sharp RD estimates using local polynomial regression.
##
## Number of Obs. 14942
## BW type msesum
## Kernel Triangular
## VCE method NN
##
## Number of Obs. 7847 7095
## Eff. Number of Obs. 2261 2382
## Order est. (p) 1 1
## Order bias (q) 2 2
## BW est. (h) 49469.948 49469.948
## BW bias (b) 87457.597 87457.597
## rho (h/b) 0.566 0.566
## Unique Obs. 7846 7095
##
## =============================================================================
## Method Coef. Std. Err. z P>|z| [ 95% C.I. ]
## =============================================================================
## Conventional 20.595 7.103 2.899 0.004 [6.672 , 34.517]
## Robust - - 2.815 0.005 [6.953 , 38.821]
## =============================================================================# Compute the bandwidth using the cerrd method from the rdrobust package
rdd_result_cerrd <- rdrobust(y = analysis_data$RoadDensity,
x = analysis_data$DistanceToBoundary,
c = 0,
covs = cbind(analysis_data$LogGDP, analysis_data$LogHDI, analysis_data$Elevation, analysis_data$LogPopDensity),
bwselect = "cerrd")
# Summarize the results
summary(rdd_result_cerrd)## Covariate-adjusted Sharp RD estimates using local polynomial regression.
##
## Number of Obs. 14942
## BW type cerrd
## Kernel Triangular
## VCE method NN
##
## Number of Obs. 7847 7095
## Eff. Number of Obs. 1814 1939
## Order est. (p) 1 1
## Order bias (q) 2 2
## BW est. (h) 39611.392 39611.392
## BW bias (b) 105645.627 105645.627
## rho (h/b) 0.375 0.375
## Unique Obs. 7846 7095
##
## =============================================================================
## Method Coef. Std. Err. z P>|z| [ 95% C.I. ]
## =============================================================================
## Conventional 16.305 7.999 2.038 0.042 [0.628 , 31.982]
## Robust - - 2.058 0.040 [0.830 , 33.877]
## =============================================================================Perform Robustness checks with different bandwidths using specific bandwidths
# Perform RDD analysis for specific bandwidths without covariates first
bandwidths <- c(20000, 40000, 60000, 80000, 100000)
for (bw in bandwidths) {
rdd_result <- rdrobust(y = analysis_data$RoadDensity,
x = analysis_data$DistanceToBoundary,
c = 0, # threshold for treatment assignment
h = bw, # specific bandwidth
p = 1, # order of polynomial
kernel = "triangular") # kernel type
cat("Results for bandwidth", bw, "m (without covariates):\n")
print(summary(rdd_result))
}## Results for bandwidth 20000 m (without covariates):
## Sharp RD estimates using local polynomial regression.
##
## Number of Obs. 14942
## BW type Manual
## Kernel Triangular
## VCE method NN
##
## Number of Obs. 7847 7095
## Eff. Number of Obs. 901 991
## Order est. (p) 1 1
## Order bias (q) 2 2
## BW est. (h) 20000.000 20000.000
## BW bias (b) 20000.000 20000.000
## rho (h/b) 1.000 1.000
## Unique Obs. 7847 7095
##
## =============================================================================
## Method Coef. Std. Err. z P>|z| [ 95% C.I. ]
## =============================================================================
## Conventional 15.544 14.660 1.060 0.289 [-13.189 , 44.277]
## Robust - - 0.773 0.439 [-24.969 , 57.516]
## =============================================================================
## NULL
## Results for bandwidth 40000 m (without covariates):
## Sharp RD estimates using local polynomial regression.
##
## Number of Obs. 14942
## BW type Manual
## Kernel Triangular
## VCE method NN
##
## Number of Obs. 7847 7095
## Eff. Number of Obs. 1831 1955
## Order est. (p) 1 1
## Order bias (q) 2 2
## BW est. (h) 40000.000 40000.000
## BW bias (b) 40000.000 40000.000
## rho (h/b) 1.000 1.000
## Unique Obs. 7847 7095
##
## =============================================================================
## Method Coef. Std. Err. z P>|z| [ 95% C.I. ]
## =============================================================================
## Conventional 25.369 10.392 2.441 0.015 [5.002 , 45.736]
## Robust - - 1.006 0.314 [-14.347 , 44.637]
## =============================================================================
## NULL
## Results for bandwidth 60000 m (without covariates):
## Sharp RD estimates using local polynomial regression.
##
## Number of Obs. 14942
## BW type Manual
## Kernel Triangular
## VCE method NN
##
## Number of Obs. 7847 7095
## Eff. Number of Obs. 2760 2812
## Order est. (p) 1 1
## Order bias (q) 2 2
## BW est. (h) 60000.000 60000.000
## BW bias (b) 60000.000 60000.000
## rho (h/b) 1.000 1.000
## Unique Obs. 7847 7095
##
## =============================================================================
## Method Coef. Std. Err. z P>|z| [ 95% C.I. ]
## =============================================================================
## Conventional 30.309 8.332 3.638 0.000 [13.979 , 46.638]
## Robust - - 1.866 0.062 [-1.166 , 47.264]
## =============================================================================
## NULL
## Results for bandwidth 80000 m (without covariates):
## Sharp RD estimates using local polynomial regression.
##
## Number of Obs. 14942
## BW type Manual
## Kernel Triangular
## VCE method NN
##
## Number of Obs. 7847 7095
## Eff. Number of Obs. 3641 3634
## Order est. (p) 1 1
## Order bias (q) 2 2
## BW est. (h) 80000.000 80000.000
## BW bias (b) 80000.000 80000.000
## rho (h/b) 1.000 1.000
## Unique Obs. 7847 7095
##
## =============================================================================
## Method Coef. Std. Err. z P>|z| [ 95% C.I. ]
## =============================================================================
## Conventional 26.515 7.048 3.762 0.000 [12.700 , 40.329]
## Robust - - 2.967 0.003 [10.768 , 52.683]
## =============================================================================
## NULL
## Results for bandwidth 1e+05 m (without covariates):
## Sharp RD estimates using local polynomial regression.
##
## Number of Obs. 14942
## BW type Manual
## Kernel Triangular
## VCE method NN
##
## Number of Obs. 7847 7095
## Eff. Number of Obs. 4446 4406
## Order est. (p) 1 1
## Order bias (q) 2 2
## BW est. (h) 100000.000 100000.000
## BW bias (b) 100000.000 100000.000
## rho (h/b) 1.000 1.000
## Unique Obs. 7847 7095
##
## =============================================================================
## Method Coef. Std. Err. z P>|z| [ 95% C.I. ]
## =============================================================================
## Conventional 23.096 6.164 3.747 0.000 [11.015 , 35.177]
## Robust - - 3.379 0.001 [13.464 , 50.652]
## =============================================================================
## NULL# Perform RDD analysis for specific bandwidths with covariates
for (bw in bandwidths) {
# Add a small constant to avoid exact zeros in covariates
adjusted_data <- analysis_data %>%
mutate(LogGDP = log(GDP + 1e-6),
LogHDI = log(HDI + 1e-6),
LogPopDensity = log(PopDensity + 1e-6))
rdd_result <- rdrobust(y = adjusted_data$RoadDensity,
x = adjusted_data$DistanceToBoundary,
c = 0, # threshold for treatment assignment
h = bw, # specific bandwidth
covs = cbind(adjusted_data$LogGDP, adjusted_data$LogHDI, adjusted_data$Elevation, adjusted_data$LogPopDensity),
p = 1, # order of polynomial
kernel = "triangular") # kernel type
cat("Results for bandwidth", bw, "m (with covariates):\n")
print(summary(rdd_result))
}## Results for bandwidth 20000 m (with covariates):
## Covariate-adjusted Sharp RD estimates using local polynomial regression.
##
## Number of Obs. 14942
## BW type Manual
## Kernel Triangular
## VCE method NN
##
## Number of Obs. 7847 7095
## Eff. Number of Obs. 901 991
## Order est. (p) 1 1
## Order bias (q) 2 2
## BW est. (h) 20000.000 20000.000
## BW bias (b) 20000.000 20000.000
## rho (h/b) 1.000 1.000
## Unique Obs. 7847 7095
##
## =============================================================================
## Method Coef. Std. Err. z P>|z| [ 95% C.I. ]
## =============================================================================
## Conventional 8.914 11.339 0.786 0.432 [-13.311 , 31.139]
## Robust - - 0.327 0.744 [-26.868 , 37.634]
## =============================================================================
## NULL
## Results for bandwidth 40000 m (with covariates):
## Covariate-adjusted Sharp RD estimates using local polynomial regression.
##
## Number of Obs. 14942
## BW type Manual
## Kernel Triangular
## VCE method NN
##
## Number of Obs. 7847 7095
## Eff. Number of Obs. 1831 1955
## Order est. (p) 1 1
## Order bias (q) 2 2
## BW est. (h) 40000.000 40000.000
## BW bias (b) 40000.000 40000.000
## rho (h/b) 1.000 1.000
## Unique Obs. 7847 7095
##
## =============================================================================
## Method Coef. Std. Err. z P>|z| [ 95% C.I. ]
## =============================================================================
## Conventional 16.409 7.958 2.062 0.039 [0.812 , 32.005]
## Robust - - 0.650 0.516 [-15.266 , 30.420]
## =============================================================================
## NULL
## Results for bandwidth 60000 m (with covariates):
## Covariate-adjusted Sharp RD estimates using local polynomial regression.
##
## Number of Obs. 14942
## BW type Manual
## Kernel Triangular
## VCE method NN
##
## Number of Obs. 7847 7095
## Eff. Number of Obs. 2760 2812
## Order est. (p) 1 1
## Order bias (q) 2 2
## BW est. (h) 60000.000 60000.000
## BW bias (b) 60000.000 60000.000
## rho (h/b) 1.000 1.000
## Unique Obs. 7847 7095
##
## =============================================================================
## Method Coef. Std. Err. z P>|z| [ 95% C.I. ]
## =============================================================================
## Conventional 21.812 6.372 3.423 0.001 [9.324 , 34.301]
## Robust - - 1.524 0.128 [-4.157 , 33.190]
## =============================================================================
## NULL
## Results for bandwidth 80000 m (with covariates):
## Covariate-adjusted Sharp RD estimates using local polynomial regression.
##
## Number of Obs. 14942
## BW type Manual
## Kernel Triangular
## VCE method NN
##
## Number of Obs. 7847 7095
## Eff. Number of Obs. 3641 3634
## Order est. (p) 1 1
## Order bias (q) 2 2
## BW est. (h) 80000.000 80000.000
## BW bias (b) 80000.000 80000.000
## rho (h/b) 1.000 1.000
## Unique Obs. 7847 7095
##
## =============================================================================
## Method Coef. Std. Err. z P>|z| [ 95% C.I. ]
## =============================================================================
## Conventional 17.847 5.393 3.309 0.001 [7.276 , 28.418]
## Robust - - 2.846 0.004 [7.282 , 39.505]
## =============================================================================
## NULL
## Results for bandwidth 1e+05 m (with covariates):
## Covariate-adjusted Sharp RD estimates using local polynomial regression.
##
## Number of Obs. 14942
## BW type Manual
## Kernel Triangular
## VCE method NN
##
## Number of Obs. 7847 7095
## Eff. Number of Obs. 4446 4406
## Order est. (p) 1 1
## Order bias (q) 2 2
## BW est. (h) 100000.000 100000.000
## BW bias (b) 100000.000 100000.000
## rho (h/b) 1.000 1.000
## Unique Obs. 7847 7095
##
## =============================================================================
## Method Coef. Std. Err. z P>|z| [ 95% C.I. ]
## =============================================================================
## Conventional 15.249 4.732 3.223 0.001 [5.975 , 24.522]
## Robust - - 3.117 0.002 [8.428 , 36.972]
## =============================================================================
## NULL# Perform RDD without any controls using the mserd bandwidth selection method
rdd_result_no_controls <- rdrobust(y = analysis_data$RoadDensity,
x = analysis_data$DistanceToBoundary,
c = 0, # threshold for treatment assignment
bwselect = "mserd", # bandwidth selection method
p = 1, # order of polynomial
kernel = "triangular") # kernel type
cat("Results for RDD without any controls using mserd bandwidth selection method:\n")## Results for RDD without any controls using mserd bandwidth selection method:print(summary(rdd_result_no_controls))## Sharp RD estimates using local polynomial regression.
##
## Number of Obs. 14942
## BW type mserd
## Kernel Triangular
## VCE method NN
##
## Number of Obs. 7847 7095
## Eff. Number of Obs. 1950 2109
## Order est. (p) 1 1
## Order bias (q) 2 2
## BW est. (h) 42730.967 42730.967
## BW bias (b) 63719.525 63719.525
## rho (h/b) 0.671 0.671
## Unique Obs. 7846 7095
##
## =============================================================================
## Method Coef. Std. Err. z P>|z| [ 95% C.I. ]
## =============================================================================
## Conventional 26.840 10.040 2.673 0.008 [7.162 , 46.518]
## Robust - - 2.051 0.040 [1.100 , 48.255]
## =============================================================================
## NULLPerform Balance Tests for HDI, GDP & Population Density
# Perform a Balance test for HDI
rdd_hdi <- rdrobust(y = analysis_data$LogHDI,
x = analysis_data$DistanceToBoundary,
c = 0,
p = 1,
kernel = "triangular")
# Summarize the results for HDI
summary(rdd_hdi)## Sharp RD estimates using local polynomial regression.
##
## Number of Obs. 14942
## BW type mserd
## Kernel Triangular
## VCE method NN
##
## Number of Obs. 7847 7095
## Eff. Number of Obs. 2657 2720
## Order est. (p) 1 1
## Order bias (q) 2 2
## BW est. (h) 57753.102 57753.102
## BW bias (b) 90001.697 90001.697
## rho (h/b) 0.642 0.642
## Unique Obs. 7846 7095
##
## =============================================================================
## Method Coef. Std. Err. z P>|z| [ 95% C.I. ]
## =============================================================================
## Conventional 0.007 0.004 1.464 0.143 [-0.002 , 0.015]
## Robust - - 1.092 0.275 [-0.005 , 0.016]
## =============================================================================# Perform a Balance test for GDP
rdd_gdp <- rdrobust(y = analysis_data$LogGDP,
x = analysis_data$DistanceToBoundary,
c = 0,
p = 1,
kernel = "triangular")
# Summarize the results for log(GDP)
summary(rdd_gdp)## Sharp RD estimates using local polynomial regression.
##
## Number of Obs. 14942
## BW type mserd
## Kernel Triangular
## VCE method NN
##
## Number of Obs. 7847 7095
## Eff. Number of Obs. 1388 1510
## Order est. (p) 1 1
## Order bias (q) 2 2
## BW est. (h) 30623.426 30623.426
## BW bias (b) 54982.059 54982.059
## rho (h/b) 0.557 0.557
## Unique Obs. 7846 7095
##
## =============================================================================
## Method Coef. Std. Err. z P>|z| [ 95% C.I. ]
## =============================================================================
## Conventional 0.252 0.165 1.521 0.128 [-0.073 , 0.576]
## Robust - - 0.993 0.321 [-0.182 , 0.556]
## =============================================================================# Perform a Balance test for Population Density
rdd_popdensity <- rdrobust(y = analysis_data$LogPopDensity,
x = analysis_data$DistanceToBoundary,
c = 0,
p = 1,
kernel = "triangular")
# Summarize the results for Population Density
summary(rdd_popdensity)## Sharp RD estimates using local polynomial regression.
##
## Number of Obs. 14942
## BW type mserd
## Kernel Triangular
## VCE method NN
##
## Number of Obs. 7847 7095
## Eff. Number of Obs. 1794 1915
## Order est. (p) 1 1
## Order bias (q) 2 2
## BW est. (h) 39007.354 39007.354
## BW bias (b) 56839.240 56839.240
## rho (h/b) 0.686 0.686
## Unique Obs. 7846 7095
##
## =============================================================================
## Method Coef. Std. Err. z P>|z| [ 95% C.I. ]
## =============================================================================
## Conventional 0.126 0.099 1.272 0.203 [-0.068 , 0.321]
## Robust - - 1.124 0.261 [-0.100 , 0.369]
## =============================================================================Perform Placebo Tests for Robustness
# Perform Placebo Tests with different cutoffs
placebo_cutoffs <- c(-50000, 50000, 100000) # Placebo cutoffs
for (cutoff in placebo_cutoffs) {
placebo_rdd <- rdrobust(y = analysis_data$RoadDensity,
x = analysis_data$DistanceToBoundary,
c = cutoff,
covs = cbind(analysis_data$LogGDP, analysis_data$LogHDI, analysis_data$Elevation, analysis_data$LogPopDensity),
kernel = "triangular")
print(summary(placebo_rdd))
}## Covariate-adjusted Sharp RD estimates using local polynomial regression.
##
## Number of Obs. 14942
## BW type mserd
## Kernel Triangular
## VCE method NN
##
## Number of Obs. 5564 9378
## Eff. Number of Obs. 1967 2043
## Order est. (p) 1 1
## Order bias (q) 2 2
## BW est. (h) 44970.121 44970.121
## BW bias (b) 75398.297 75398.297
## rho (h/b) 0.596 0.596
## Unique Obs. 5563 9378
##
## =============================================================================
## Method Coef. Std. Err. z P>|z| [ 95% C.I. ]
## =============================================================================
## Conventional 1.522 5.185 0.294 0.769 [-8.641 , 11.686]
## Robust - - 0.545 0.586 [-8.491 , 15.033]
## =============================================================================
## NULL
## Covariate-adjusted Sharp RD estimates using local polynomial regression.
##
## Number of Obs. 14942
## BW type mserd
## Kernel Triangular
## VCE method NN
##
## Number of Obs. 10246 4696
## Eff. Number of Obs. 1060 962
## Order est. (p) 1 1
## Order bias (q) 2 2
## BW est. (h) 23227.371 23227.371
## BW bias (b) 45951.911 45951.911
## rho (h/b) 0.505 0.505
## Unique Obs. 10245 4696
##
## =============================================================================
## Method Coef. Std. Err. z P>|z| [ 95% C.I. ]
## =============================================================================
## Conventional 7.451 7.742 0.962 0.336 [-7.722 , 22.625]
## Robust - - 0.490 0.624 [-12.753 , 21.254]
## =============================================================================
## NULL
## Covariate-adjusted Sharp RD estimates using local polynomial regression.
##
## Number of Obs. 14942
## BW type mserd
## Kernel Triangular
## VCE method NN
##
## Number of Obs. 12253 2689
## Eff. Number of Obs. 1109 870
## Order est. (p) 1 1
## Order bias (q) 2 2
## BW est. (h) 28214.108 28214.108
## BW bias (b) 41329.816 41329.816
## rho (h/b) 0.683 0.683
## Unique Obs. 12252 2689
##
## =============================================================================
## Method Coef. Std. Err. z P>|z| [ 95% C.I. ]
## =============================================================================
## Conventional 11.060 8.644 1.280 0.201 [-5.881 , 28.001]
## Robust - - 1.188 0.235 [-7.856 , 32.029]
## =============================================================================
## NULLPerform Sensitivity analysis using a quadratic polynomial
# Perform Sensitivity analysis with a quadratic polynomial
rdd_sensitivity <- rdrobust(y = analysis_data$RoadDensity,
x = analysis_data$DistanceToBoundary,
c = 0,
p = 2,
covs = cbind(analysis_data$LogGDP, analysis_data$LogHDI, analysis_data$Elevation, analysis_data$LogPopDensity),
kernel = "triangular")
print(summary(rdd_sensitivity))## Covariate-adjusted Sharp RD estimates using local polynomial regression.
##
## Number of Obs. 14942
## BW type mserd
## Kernel Triangular
## VCE method NN
##
## Number of Obs. 7847 7095
## Eff. Number of Obs. 3099 3113
## Order est. (p) 2 2
## Order bias (q) 3 3
## BW est. (h) 67221.096 67221.096
## BW bias (b) 93304.722 93304.722
## rho (h/b) 0.720 0.720
## Unique Obs. 7846 7095
##
## =============================================================================
## Method Coef. Std. Err. z P>|z| [ 95% C.I. ]
## =============================================================================
## Conventional 18.932 8.996 2.104 0.035 [1.300 , 36.564]
## Robust - - 1.654 0.098 [-3.079 , 36.423]
## =============================================================================
## NULLCreate tables
Linear regression model Results
# Extract summary statistics
summary_stats <- summary(linear_reg)
# Create a dataframe for the results
results <- data.frame(
Metric = c("Coefficient", "Std. Error"),
Intercept = c(
paste0(round(summary_stats$coefficients[1, 1], 3),
ifelse(summary_stats$coefficients[1, 4] < 0.001, "***",
ifelse(summary_stats$coefficients[1, 4] < 0.01, "**",
ifelse(summary_stats$coefficients[1, 4] < 0.05, "*", "")))),
paste0("(", round(summary_stats$coefficients[1, 2], 3), ")")
),
RomanEmpire = c(
paste0(round(summary_stats$coefficients[2, 1], 3),
ifelse(summary_stats$coefficients[2, 4] < 0.001, "***",
ifelse(summary_stats$coefficients[2, 4] < 0.01, "**",
ifelse(summary_stats$coefficients[2, 4] < 0.05, "*", "")))),
paste0("(", round(summary_stats$coefficients[2, 2], 3), ")")
),
GDP = c(
paste0(round(summary_stats$coefficients[3, 1], 3),
ifelse(summary_stats$coefficients[3, 4] < 0.001, "***",
ifelse(summary_stats$coefficients[3, 4] < 0.01, "**",
ifelse(summary_stats$coefficients[3, 4] < 0.05, "*", "")))),
paste0("(", round(summary_stats$coefficients[3, 2], 3), ")")
),
HDI = c(
paste0(round(summary_stats$coefficients[4, 1], 3),
ifelse(summary_stats$coefficients[4, 4] < 0.001, "***",
ifelse(summary_stats$coefficients[4, 4] < 0.01, "**",
ifelse(summary_stats$coefficients[4, 4] < 0.05, "*", "")))),
paste0("(", round(summary_stats$coefficients[4, 2], 3), ")")
),
Elevation = c(
paste0(round(summary_stats$coefficients[5, 1], 3),
ifelse(summary_stats$coefficients[5, 4] < 0.001, "***",
ifelse(summary_stats$coefficients[5, 4] < 0.01, "**",
ifelse(summary_stats$coefficients[5, 4] < 0.05, "*", "")))),
paste0("(", round(summary_stats$coefficients[5, 2], 3), ")")
),
PopDensity = c(
paste0(round(summary_stats$coefficients[6, 1], 3),
ifelse(summary_stats$coefficients[6, 4] < 0.001, "***",
ifelse(summary_stats$coefficients[6, 4] < 0.01, "**",
ifelse(summary_stats$coefficients[6, 4] < 0.05, "*", "")))),
paste0("(", round(summary_stats$coefficients[6, 2], 3), ")")
)
)
# Transpose the table
results_transposed <- as.data.frame(t(results))
# Add column names
colnames(results_transposed) <- results$Metric
results_transposed <- results_transposed[-1, ] # Remove the first row with metric names
# Add coefficient names as a column
results_transposed <- cbind(Variable = rownames(results_transposed), results_transposed)
rownames(results_transposed) <- NULL
# Create a table using modelsummary
datasummary_df(results_transposed, title = "Table - Linear Regression Results", output = "kableExtra")| Variable | Coefficient | Std. Error |
|---|---|---|
| Intercept | -171.706*** | (7.807) |
| RomanEmpire | 12.784*** | (1.198) |
| GDP | 0*** |
|
| HDI | 233.056*** | (9.283) |
| Elevation | -0.031*** | (0.002) |
| PopDensity | 0.151*** | (0.004) |
RDD Results
# Perform RDD analyses and store the results
rdd_results <- list()
# Pre-specified bandwidths
bandwidths <- c(20000, 40000, 60000, 80000, 100000)
for (bw in bandwidths) {
rdd_results[[paste0(bw / 1000, "km without covariates")]] <- rdrobust(
y = analysis_data$RoadDensity,
x = analysis_data$DistanceToBoundary,
c = 0,
h = bw,
p = 1,
kernel = "triangular"
)
rdd_results[[paste0(bw / 1000, "km with covariates")]] <- rdrobust(
y = analysis_data$RoadDensity,
x = analysis_data$DistanceToBoundary,
c = 0,
h = bw,
covs = cbind(analysis_data$LogGDP, analysis_data$LogHDI, analysis_data$Elevation, analysis_data$LogPopDensity),
p = 1,
kernel = "triangular"
)
}
# Bandwidth selection methods
bwselect_methods <- c("mserd", "msesum", "msetwo", "cerrd")
for (method in bwselect_methods) {
rdd_results[[paste0(method, " without covariates")]] <- rdrobust(
y = analysis_data$RoadDensity,
x = analysis_data$DistanceToBoundary,
c = 0,
bwselect = method,
p = 1,
kernel = "triangular"
)
rdd_results[[paste0(method, " with covariates")]] <- rdrobust(
y = analysis_data$RoadDensity,
x = analysis_data$DistanceToBoundary,
c = 0,
bwselect = method,
covs = cbind(analysis_data$LogGDP, analysis_data$LogHDI, analysis_data$Elevation, analysis_data$LogPopDensity),
p = 1,
kernel = "triangular"
)
}
# Function to extract RDD results
make_table_row <- function(estimates) {
est <- round(estimates$Estimate[1], 3)
se <- round(estimates$Estimate[3], 3)
robust_p_value <- round(estimates$pv[3], 3)
stars <- case_when(
robust_p_value < 0.01 ~ "***",
robust_p_value < 0.05 ~ "**",
robust_p_value < 0.10 ~ "*",
TRUE ~ ""
)
n_control <- estimates$N[1]
n_treat <- estimates$N[2]
bw_left <- round(estimates$bws[1, 1], 3)
bw_right <- round(estimates$bws[1, 2], 3)
data.frame(
Method = "",
Coefficient = paste0(est, stars),
`SE (BC)` = paste0("(", se, ")"),
`N (Treated)` = as.character(n_treat),
`N (Control)` = as.character(n_control),
Bandwidth = paste0(bw_left, " / ", bw_right)
)
}
# Create a combined summary table
summary_table <- map_df(rdd_results, make_table_row, .id = "Method")
# Reformat the table to have methods as rows and metrics as columns
summary_table_transposed <- summary_table %>%
pivot_longer(cols = -Method, names_to = "Metric", values_to = "Value") %>%
pivot_wider(names_from = Metric, values_from = Value)
# Create a table using modelsummary
datasummary_df(summary_table_transposed, title = "Table 3 - RDD Results", output = "kableExtra")| Method | Coefficient | SE..BC. | N..Treated. | N..Control. | Bandwidth |
|---|---|---|---|---|---|
| 20km without covariates | 15.544 | (14.66) | 7095 | 7847 | 20000 / 20000 |
| 20km with covariates | 8.914 | (11.339) | 7095 | 7847 | 20000 / 20000 |
| 40km without covariates | 25.369 | (10.392) | 7095 | 7847 | 40000 / 40000 |
| 40km with covariates | 16.488 | (7.959) | 7095 | 7847 | 40000 / 40000 |
| 60km without covariates | 30.309* | (8.332) | 7095 | 7847 | 60000 / 60000 |
| 60km with covariates | 21.919 | (6.373) | 7095 | 7847 | 60000 / 60000 |
| 80km without covariates | 26.515*** | (7.048) | 7095 | 7847 | 80000 / 80000 |
| 80km with covariates | 17.948*** | (5.394) | 7095 | 7847 | 80000 / 80000 |
| 100km without covariates | 23.096*** | (6.164) | 7095 | 7847 | 1e+05 / 1e+05 |
| 100km with covariates | 15.334*** | (4.732) | 7095 | 7847 | 1e+05 / 1e+05 |
| mserd without covariates | 26.84** | (10.04) | 7095 | 7847 | 42730.967 / 42730.967 |
| mserd with covariates | 21.656*** | (6.138) | 7095 | 7847 | 64053.158 / 64053.158 |
| msesum without covariates | 21.354 | (11.931) | 7095 | 7847 | 30384.082 / 30384.082 |
| msesum with covariates | 20.595*** | (7.103) | 7095 | 7847 | 49469.948 / 49469.948 |
| msetwo without covariates | 21.539 | (11.503) | 7095 | 7847 | 30258.83 / 36229.785 |
| msetwo with covariates | 21.419*** | (6.438) | 7095 | 7847 | 77067.969 / 45517.391 |
| cerrd without covariates | 19.729 | (12.788) | 7095 | 7847 | 26425.443 / 26425.443 |
| cerrd with covariates | 16.305** | (7.999) | 7095 | 7847 | 39611.392 / 39611.392 |
Variables Summary
# Create a summary table for the variables
datasummary(
RoadDensity + RomanEmpire + GDP + HDI + Elevation + PopDensity ~ N + Mean + SD + Min + Max,
data = analysis_data,
title = "Table 1 - Summary of the Variables"
)| N | Mean | SD | Min | Max | |
|---|---|---|---|---|---|
| RoadDensity | 14942 | 49.73 | 97.97 | 0.00 | 1304.57 |
| RomanEmpire | 14942 | 0.47 | 0.50 | 0.00 | 1.00 |
| GDP | 14942 | 344856703.99 | 954220457.11 | 0.00 | 30173835264.00 |
| HDI | 14942 | 0.85 | 0.07 | 0.70 | 0.95 |
| Elevation | 14942 | 362.48 | 319.91 | -14.85 | 2348.56 |
| PopDensity | 14942 | 171.29 | 329.68 | 1.51 | 8104.50 |
Balance Test Results
# Define the balance test results
balance_tests <- list(
`Log HDI` = rdrobust(y = analysis_data$LogHDI, x = analysis_data$DistanceToBoundary, c = 0, p = 1, kernel = "triangular"),
`Log GDP` = rdrobust(y = analysis_data$LogGDP, x = analysis_data$DistanceToBoundary, c = 0, p = 1, kernel = "triangular"),
`Log Population Density` = rdrobust(y = analysis_data$LogPopDensity, x = analysis_data$DistanceToBoundary, c = 0, p = 1, kernel = "triangular")
)
# Function to make table columns from balance test results
make_balance_table_column <- function(estimates) {
est <- estimates$coef[1, 1] |> round(3)
se <- estimates$se[1, 1] |> round(3)
robust_p_value <- estimates$pv[1, 1] |> round(3)
stars <- case_when(
robust_p_value < 0.01 ~ "***",
robust_p_value < 0.05 ~ "**",
robust_p_value < 0.10 ~ "*",
TRUE ~ ""
)
n_control <- estimates$N[2]
n_treat <- estimates$N[1]
bw <- estimates$bws[1, 1] |> round(3)
data.frame(out = c(
paste0(est, stars),
paste0("(", se, ")"),
n_treat,
n_control,
bw
))
}
# Function to summarize balance test models in a neat table
modelsummary_balance <- function(list_of_balance_tests, row_names = NULL, row_names_add = NULL, ...) {
if (is.null(row_names)) {
row_names <- c("Coefficient", "SE (BC)", "N (Treated)", "N (Control)", "Bandwidth")
if (!is.null(row_names_add)) {
row_names <- c(row_names, row_names_add)
}
}
names <- names(list_of_balance_tests)
df <- list_of_balance_tests |>
map(~make_balance_table_column(.x)) |>
discard(is.null) |>
reduce(bind_cols)
if (nrow(df) == 0) {
return(NULL)
}
out <- bind_cols(" " = row_names, df)
names(out) <- c(" ", names)
datasummary_df(out, ...)
}
# Create the summary table for balance tests
modelsummary_balance(balance_tests, title = "Table 2 - Balance Tests Results")## New names:
## New names:
## New names:
## • `out` -> `out...1`
## • `out` -> `out...2`| Log HDI | Log GDP | Log Population Density | |
|---|---|---|---|
| Coefficient | 0.007 | 0.252 | 0.126 |
| SE (BC) | (0.004) | (0.165) | (0.099) |
| N (Treated) | 7847 | 7847 | 7847 |
| N (Control) | 7095 | 7095 | 7095 |
| Bandwidth | 57753.102 | 30623.426 | 39007.354 |
Placebo Test Results
# Perform Placebo RDD analyses and store the results
placebo_rdd_results <- list()
# Loop over placebo cutoffs and run rdrobust
for (cutoff in placebo_cutoffs) {
placebo_rdd_results[[paste0("Placebo cutoff at ", cutoff, " meters")]] <- rdrobust(
y = analysis_data$RoadDensity,
x = analysis_data$DistanceToBoundary,
c = cutoff,
covs = cbind(analysis_data$LogGDP, analysis_data$LogHDI, analysis_data$Elevation, analysis_data$LogPopDensity),
kernel = "triangular"
)
}
# Function to make table columns from RDD estimates
make_placebo_table_column <- function(estimates, margin_dv = 0.05, extra_rows = NULL) {
est <- estimates$Estimate[1] |> round(3)
se <- estimates$Estimate[3] |> round(3)
robust_p_value <- estimates$pv[3] |> round(3)
stars <- case_when(
robust_p_value < 0.01 ~ "***",
robust_p_value < 0.05 ~ "**",
robust_p_value < 0.10 ~ "*",
TRUE ~ ""
)
n_control <- estimates$N[1]
n_treat <- estimates$N[2]
bw <- estimates$bws[1, 1] |> round(3)
data.frame(out = c(
paste0(est, stars),
paste0("(", se, ")"),
n_treat,
n_control,
bw,
extra_rows
))
}
# Function to summarize RDD models in a neat table
modelsummary_placebo_rdrobust <- function(list_of_rdrobust_objects, row_names = NULL, row_names_add = NULL, ...) {
if (is.null(row_names)) {
row_names <- c("Coefficient", "SE (BC)", "N (Treated)", "N (Control)", "Bandwidth")
if (!is.null(row_names_add)) {
row_names <- c(row_names, row_names_add)
}
}
names <- names(list_of_rdrobust_objects)
df <- list_of_rdrobust_objects |>
map(make_placebo_table_column) |>
reduce(bind_cols)
names(df) <- names
out <- bind_cols(" " = row_names, df)
datasummary_df(out, ...)
}
# Create a table for placebo tests
modelsummary_placebo_rdrobust(placebo_rdd_results, title = "Table 5 - Placebo Test Results")## New names:
## New names:
## • `out` -> `out...1`
## • `out` -> `out...2`| Placebo cutoff at -50000 meters | Placebo cutoff at 50000 meters | Placebo cutoff at 1e+05 meters | |
|---|---|---|---|
| Coefficient | 1.522 | 7.451 | 11.06 |
| SE (BC) | (5.185) | (7.742) | (8.644) |
| N (Treated) | 9378 | 4696 | 2689 |
| N (Control) | 5564 | 10246 | 12253 |
| Bandwidth | 44970.121 | 23227.371 | 28214.108 |
# Perform sensitivity analyses with different polynomial orders and store the results
sensitivity_results <- list()
# Polynomial orders to test
polynomial_orders <- c(1, 2) # Example polynomial orders
# Loop over polynomial orders and run rdrobust
for (p_order in polynomial_orders) {
sensitivity_results[[paste0("Polynomial order ", p_order)]] <- rdrobust(
y = analysis_data$RoadDensity,
x = analysis_data$DistanceToBoundary,
c = 0,
p = p_order,
covs = cbind(analysis_data$LogGDP, analysis_data$LogHDI, analysis_data$Elevation, analysis_data$LogPopDensity),
kernel = "triangular"
)
}
# Function to make table columns from RDD estimates
make_sensitivity_table_column <- function(estimates, margin_dv = 0.05, extra_rows = NULL) {
est <- estimates$Estimate[1] |> round(3)
se <- estimates$Estimate[3] |> round(3)
robust_p_value <- estimates$pv[3] |> round(3)
stars <- case_when(
robust_p_value < 0.01 ~ "***",
robust_p_value < 0.05 ~ "**",
robust_p_value < 0.10 ~ "*",
TRUE ~ ""
)
n_control <- estimates$N[1]
n_treat <- estimates$N[2]
bw <- estimates$bws[1, 1] |> round(3)
data.frame(out = c(
paste0(est, stars),
paste0("(", se, ")"),
n_treat,
n_control,
bw,
extra_rows
))
}
# Function to summarize RDD models in a neat table
modelsummary_sensitivity_rdrobust <- function(list_of_rdrobust_objects, row_names = NULL, row_names_add = NULL, ...) {
if (is.null(row_names)) {
row_names <- c("Coefficient", "SE (BC)", "N (Treated)", "N (Control)", "Bandwidth")
if (!is.null(row_names_add)) {
row_names <- c(row_names, row_names_add)
}
}
names <- names(list_of_rdrobust_objects)
df <- list_of_rdrobust_objects |>
map(make_sensitivity_table_column) |>
reduce(bind_cols)
names(df) <- names
out <- bind_cols(" " = row_names, df)
datasummary_df(out, ...)
}
# Create the summary table for sensitivity analysis
modelsummary_sensitivity_rdrobust(sensitivity_results, title = "Table 6 - Sensitivity Analysis Results")## New names:
## • `out` -> `out...1`
## • `out` -> `out...2`| Polynomial order 1 | Polynomial order 2 | |
|---|---|---|
| Coefficient | 21.656*** | 18.932* |
| SE (BC) | (6.138) | (8.996) |
| N (Treated) | 7095 | 7095 |
| N (Control) | 7847 | 7847 |
| Bandwidth | 64053.158 | 67221.096 |
Visualization
Visualize RDD results
# Bin the data
binned_data <- analysis_data %>%
mutate(bin = cut(DistanceToBoundary, breaks = seq(-200000, 200000, by = 10000))) %>%
group_by(bin) %>%
summarise(RoadDensity = mean(RoadDensity, na.rm = TRUE),
DistanceToBoundary = mean(DistanceToBoundary, na.rm = TRUE))
# Plot the binned scatterplot with regression lines
ggplot() +
geom_point(data = binned_data, aes(x = DistanceToBoundary, y = RoadDensity), color = "blue") +
geom_smooth(data = analysis_data %>% filter(DistanceToBoundary < 0),
aes(x = DistanceToBoundary, y = RoadDensity),
method = "lm", formula = y ~ x, color = "red") +
geom_smooth(data = analysis_data %>% filter(DistanceToBoundary >= 0),
aes(x = DistanceToBoundary, y = RoadDensity),
method = "lm", formula = y ~ x, color = "green") +
geom_vline(xintercept = 0, linetype = "dashed") +
labs(title = "RDD Plot: Road Density vs. Distance to Roman Empire Boundary",
x = "Distance to Roman Empire Boundary (meters)",
y = "Road Density (m/km²)") +
theme_minimal()