Load Libraries
library(sf)
library(maps)
library(ggplot2)
library(ggmap)
library(dplyr)
library(mapdata)
library(tidyr)
library(osmdata)
library(tidyverse)
library(terra)
library(sp)
library(fixest)
library(raster)
library(scpi)
library(MatchIt)
library(knitr)
library(kableExtra)
library(sjPlot)
Data Preperation
Load Cities
#Load the 'world.cities' dataset
data(world.cities)
head(world.cities)
## name country.etc pop lat long capital
## 1 'Abasan al-Jadidah Palestine 5629 31.31 34.34 0
## 2 'Abasan al-Kabirah Palestine 18999 31.32 34.35 0
## 3 'Abdul Hakim Pakistan 47788 30.55 72.11 0
## 4 'Abdullah-as-Salam Kuwait 21817 29.36 47.98 0
## 5 'Abud Palestine 2456 32.03 35.07 0
## 6 'Abwein Palestine 3434 32.03 35.20 0
filtered_cities <- subset(world.cities, country.etc %in% c("France", "Germany", "South Africa"))
#Define the host cities for each country
host_cities_list <- list(
France = c("Paris", "Marseille", "Lyon", "Lens", "Nantes", "Toulouse", "Saint-Etienne", "Bordeaux", "Montpellier"),
`South Africa` = c("Johannesburg", "Cape Town", "Durban", "Pretoria", "Port Elizabeth", "Rustenburg",
"Pietersburg", "Nelspruit", "Bloemfontein"),
Germany = c("Berlin", "Dortmund", "Frankfurt", "Gelsenkirchen", "Hamburg", "Hanover",
"Kaiserslautern", "Cologne", "Leipzig", "Munich", "Nuremberg", "Stuttgart")
)
#Normalize city names and countries to handle potential mismatches
filtered_cities$name <- tolower(filtered_cities$name)
filtered_cities$country.etc <- tolower(filtered_cities$country.etc)
host_cities_list <- lapply(host_cities_list, tolower)
#Add a column 'host_city' and initialize it with 0
filtered_cities$host_city <- 0
#Update the 'host_city' column to 1 for host cities
for (country in names(host_cities_list)) {
cities <- host_cities_list[[country]]
country_index <- filtered_cities$country.etc == tolower(country)
city_index <- filtered_cities$name %in% cities
filtered_cities$host_city[country_index & city_index] <- 1
}
Filter Out Problematic cities or Rename them
filtered_cities<- subset(filtered_cities, !(name == "frankfurt" & pop == 64389))
filtered_cities <- filtered_cities %>%
mutate(name = ifelse(name == "emnambithi", "ladysmith", name))
filtered_cities <- filtered_cities %>%
mutate(name = ifelse(name == "spires", "speyer", name))
Filter out unnecassary cities
filtered_cities_fr <- subset(filtered_cities, country.etc == "france" & as.numeric(pop) > 20000)
filtered_cities_ge <- subset(filtered_cities, country.etc == "germany" & as.numeric(pop) > 20000)
filtered_cities_sa <- subset(filtered_cities, country.etc == "south africa" & as.numeric(pop) > 20000)
Match the cities to polygons
library(osmdata)
library(sf)
library(dplyr)
#Function to convert bounding box to polygon
bbox_to_polygon <- function(location, country) {
#Attempt to get bounding box from osmdata
bbox <- tryCatch(
{
full_location <- paste(location, country, sep = ", ")
getbb(full_location)
},
error = function(e) {
message(paste("Error retrieving bounding box for:", full_location, ":", e$message))
return(NULL)
}
)
if (is.null(bbox)) {
return(NULL)
}
min_x <- bbox["x", "min"]
max_x <- bbox["x", "max"]
min_y <- bbox["y", "min"]
max_y <- bbox["y", "max"]
coords <- matrix(c(
min_x, min_y,
max_x, min_y,
max_x, max_y,
min_x, max_y,
min_x, min_y
), ncol = 2, byrow = TRUE)
polygon <- st_polygon(list(coords))
sf_polygon <- st_sfc(polygon, crs = 4326)
return(sf_polygon)
}
#Function to create a polygon from a row
polygon_from_row <- function(row) {
loc <- row$name
country <- row$country.etc
bbox_to_polygon(loc, country)
}
Apply the polygon function to all the Countries
#Apply the function to each row
filtered_cities_fr <- filtered_cities_fr %>%
rowwise() %>%
mutate(polygon = list(polygon_from_row(cur_data()))) %>%
ungroup()
## Error retrieving bounding box for: bruay-la-brussiere, france : `place_name` 'bruay-la-brussiere, france' can't be found
## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `polygon = list(polygon_from_row(cur_data()))`.
## ℹ In row 1.
## Caused by warning:
## ! `cur_data()` was deprecated in dplyr 1.1.0.
## ℹ Please use `pick()` instead.
#Apply the function to each row
filtered_cities_ge <- filtered_cities_ge %>%
rowwise() %>%
mutate(polygon = list(polygon_from_row(cur_data()))) %>%
ungroup()
#Apply the function to each row
filtered_cities_sa <- filtered_cities_sa %>%
rowwise() %>%
mutate(polygon = list(polygon_from_row(cur_data()))) %>%
ungroup()
## Error retrieving bounding box for: bulfontein, south africa : `place_name` 'bulfontein, south africa' can't be found
## Error retrieving bounding box for: epumalanga, south africa : `place_name` 'epumalanga, south africa' can't be found
## Error retrieving bounding box for: garankuwa, south africa : `place_name` 'garankuwa, south africa' can't be found
## Error retrieving bounding box for: kwangema, south africa : `place_name` 'kwangema, south africa' can't be found
## Error retrieving bounding box for: nkowakowa, south africa : `place_name` 'nkowakowa, south africa' can't be found
filtered_cities_fr_sf <- filtered_cities_fr %>%
unnest_wider(polygon, names_sep = '_1') %>%
st_as_sf()
filtered_cities_ge_sf <- filtered_cities_ge %>%
unnest_wider(polygon, names_sep = '_1') %>%
st_as_sf()
filtered_cities_sa_sf <- filtered_cities_sa %>%
unnest_wider(polygon, names_sep = '_1') %>%
st_as_sf()
Remove bugs in Polygons
filtered_cities_ge_sf$polygon_area <- st_area(filtered_cities_ge_sf$polygon_11)
filtered_cities_fr_sf$polygon_area <- st_area(filtered_cities_fr_sf$polygon_11)
filtered_cities_sa_sf$polygon_area <- st_area(filtered_cities_sa_sf$polygon_11)
#Filter filtered_cities_ge_sf dataset
filtered_cities_ge_sf <- filtered_cities_ge_sf %>%
filter(!(host_city == 0 & as.numeric(polygon_area) > 10000000000))
#Filter filtered_cities_fr_sf dataset
filtered_cities_fr_sf <- filtered_cities_fr_sf %>%
filter(!(host_city == 0 & as.numeric(polygon_area) > 10000000000))
#Filter filtered_cities_sa_sf dataset
filtered_cities_sa_sf <- filtered_cities_sa_sf %>%
filter(!(host_city == 0 & as.numeric(polygon_area) > 10000000000))
saveRDS(filtered_cities_fr_sf, "filtered_cities_fr_sf.rds")
saveRDS(filtered_cities_ge_sf, "filtered_cities_ge_sf.rds")
saveRDS(filtered_cities_sa_sf, "filtered_cities_sa_sf.rds")
filtered_cities_ge_sf <- readRDS("filtered_cities_ge_sf.rds")
filtered_cities_fr_sf <- readRDS("filtered_cities_fr_sf.rds")
filtered_cities_sa_sf <- readRDS("filtered_cities_sa_sf.rds")
Bounding boxes to only contain Mainland
#Define bounding boxes for mainland regions
bounding_boxes <- list(
Germany = c(xmin = 5.87, xmax = 15.04, ymin = 47.27, ymax = 55.12),
South_Africa = c(xmin = 16.28, xmax = 32.89, ymin = -34.83, ymax = -22.09),
France = c(xmin = -5.25, xmax = 9.6, ymin = 41.25, ymax = 51.1)
)
#Define the CRS
crs <- st_crs("+proj=longlat +datum=WGS84")
#Convert bounding boxes to sf objects
bounding_boxes_sf <- lapply(bounding_boxes, function(bb) {
st_bbox(bb, crs = crs) %>%
st_as_sfc()
})
#Filter each dataset for mainland regions
filtered_cities_fr_sf <- st_intersection(filtered_cities_fr_sf, bounding_boxes_sf[["France"]])
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
filtered_cities_ge_sf <- st_intersection(filtered_cities_ge_sf, bounding_boxes_sf[["Germany"]])
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
filtered_cities_sa_sf <- st_intersection(filtered_cities_sa_sf, bounding_boxes_sf[["South_Africa"]])
## Warning: attribute variables are assumed to be spatially constant throughout
## all geometries
Plot city polygons and a comparison in nightlights for South
Africa
ggplot(data = filtered_cities_fr_sf) +
geom_sf(color = "black", alpha = 0.5) +
theme_minimal() +
labs(title = "City Polygons in France")
ggplot(data = filtered_cities_ge_sf) +
geom_sf(color = "black", alpha = 0.5) +
theme_minimal() +
labs(title = "City Polygons in Germany")
ggplot(data = filtered_cities_sa_sf) +
geom_sf(color = "black", alpha = 0.5) +
theme_minimal() +
labs(title = "City Polygons in South Africa")
shapefile_sa <- st_read("sa/gadm41_ZAF_0.shp")
## Reading layer `gadm41_ZAF_0' from data source
## `/Users/stevendenotter/Documents/Scriptie/SA/gadm41_ZAF_0.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 1 feature and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 16.45189 ymin: -34.83514 xmax: 32.89125 ymax: -22.12503
## Geodetic CRS: WGS 84
nightlights_sa_94 <- raster("nightlights/F121994.v4b_web.stable_lights.avg_vis.tif")
nightlights_sa_13 <- raster("nightlights/F182013.v4c_web.stable_lights.avg_vis.tif")
st_as_sf(shapefile_sa)
## Simple feature collection with 1 feature and 2 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 16.45189 ymin: -34.83514 xmax: 32.89125 ymax: -22.12503
## Geodetic CRS: WGS 84
## GID_0 COUNTRY geometry
## 1 ZAF South Africa MULTIPOLYGON (((19.66291 -3...
nl_sa_94 <- crop(nightlights_sa_94, shapefile_sa)
nl_sa_13 <- crop(nightlights_sa_13, shapefile_sa)
nl_sa_94 <- raster::mask(nl_sa_94 ,mask=shapefile_sa)
nl_sa_13 <- raster::mask(nl_sa_13 ,mask=shapefile_sa)
par(mfrow = c(1, 2))
#Plot 1: Nightlights in South Africa 1994
plot(nl_sa_94, col = gray.colors(100), main = "Nightlights In South Africa in 1994")
#Plot 2: Nightlights in South Africa 2013
plot(nl_sa_13, col = gray.colors(100), main = "Nightlights In South Africa in 2013")
par(mfrow = c(1, 1))
ggplot() +
geom_sf(data = shapefile_sa, fill = "lightgray", color = "black") +
geom_sf(data = filtered_cities_sa_sf, color = "black", alpha = 0.5) +
theme_minimal() +
labs(title = "City Polygons in South Africa", x = "Longitude", y = "Latitude")
Rename and make the dataframe tidy for analysis
# For France
tidy_fr <- fr_nili_data %>%
pivot_longer(cols = starts_with("nili_per_polygon"),
names_to = "Year",
values_to = "Value") %>%
mutate(Year = str_sub(Year, -4, nchar(Year)))
# For Germany
tidy_ge <- ge_nili_data %>%
pivot_longer(cols = starts_with("nili_per_polygon"),
names_to = "Year",
values_to = "Value") %>%
mutate(Year = str_sub(Year, -4, nchar(Year)))
# For South Africa
tidy_sa <- sa_nili_data %>%
pivot_longer(cols = starts_with("nili_per_polygon"),
names_to = "Year",
values_to = "Value") %>%
mutate(Year = str_sub(Year, -4, nchar(Year)))
saveRDS(tidy_fr, "tidy_fr.rds")
saveRDS(tidy_ge, "tidy_ge.rds")
saveRDS(tidy_sa, "tidy_sa.rds")
Compute Nightlights trend before Treatment
tidy_fr <- tidy_fr %>%
mutate(Year = as.numeric(Year))
#Compute Nightlights Trend Before Treatment
nightlights_trend_before_treatment <- tapply(tidy_fr$Value[tidy_fr$Year < 1998], tidy_fr$name[tidy_fr$Year < 1998], function(x) {
if (all(x == unique(63))) {
#Constant trend, assign a specific value (e.g., 0)
return(0)
} else {
#Non-constant trend, fit the linear model
return(lm(x ~ tidy_fr$Year[match(x, tidy_fr$Value)])$coefficients[2])
}
})
nightlights_trend_before_treatment[is.na(nightlights_trend_before_treatment)] <- 0
trend_df_fr <- data.frame(row_number = 1:length(nightlights_trend_before_treatment),
nightlights_trend = nightlights_trend_before_treatment)
tidy_ge <- tidy_ge %>%
mutate(Year = as.numeric(Year))
nightlights_trend_before_treatment_ge <- tapply(tidy_ge$Value[tidy_ge$Year < 2006], tidy_ge$name[tidy_ge$Year < 2006], function(x) {
if (all(x == unique(63))) {
return(0)
} else {
return(lm(x ~ tidy_ge$Year[match(x, tidy_ge$Value)])$coefficients[2])
}
})
nightlights_trend_before_treatment_ge[is.na(nightlights_trend_before_treatment_ge)] <- 0
trend_df_ge <- data.frame(row_number = 1:length(nightlights_trend_before_treatment_ge),
nightlights_trend = nightlights_trend_before_treatment_ge)
tidy_sa <- tidy_sa %>%
mutate(Year = as.numeric(Year))
nightlights_trend_before_treatment_sa <- tapply(tidy_sa$Value[tidy_sa$Year < 2010], tidy_sa$name[tidy_sa$Year < 2010], function(x) {
if (all(x == unique(63))) {
return(0)
} else {
return(lm(x ~ tidy_sa$Year[match(x, tidy_sa$Value)])$coefficients[2])
}
})
nightlights_trend_before_treatment_sa[is.na(nightlights_trend_before_treatment_sa)] <- 0
trend_df_sa <- data.frame(row_number = 1:length(nightlights_trend_before_treatment_sa),
nightlights_trend = nightlights_trend_before_treatment_sa)
Add trend to Data
#France
filtered_cities_fr_sf <- filtered_cities_fr_sf %>%
mutate(ID = row_number())
fr_trend <- merge(filtered_cities_fr_sf, tibble(ID = 1:nrow(trend_df_fr), nightlights_trend = nightlights_trend_before_treatment), by = "ID", all.x = TRUE)
cities_with_treatment <- fr_trend[fr_trend$host_city == 1, "name"]
fr_trend <- fr_trend %>%
filter(!is.na(nightlights_trend))
#Germany
filtered_cities_ge_sf <- filtered_cities_ge_sf %>%
mutate(ID = row_number())
ge_trend <- merge(filtered_cities_ge_sf, tibble(ID = 1:nrow(trend_df_ge), nightlights_trend = nightlights_trend_before_treatment_ge), by = "ID", all.x = TRUE)
cities_with_treatment <- ge_trend[ge_trend$host_city == 1, "name"]
ge_trend <- ge_trend %>%
filter(!is.na(nightlights_trend))
#South Africe
filtered_cities_sa_sf <- filtered_cities_sa_sf %>%
mutate(ID = row_number())
sa_trend <- merge(filtered_cities_sa_sf, tibble(ID = 1:nrow(trend_df_sa), nightlights_trend = nightlights_trend_before_treatment_sa), by = "ID", all.x = TRUE)
cities_with_treatment <- sa_trend[sa_trend$host_city == 1, "name"]
sa_trend <- sa_trend %>%
filter(!is.na(nightlights_trend))
Propensity Score Matching
Matching for France
#Calculate propensity scores using logistic regression
model_fr <- glm(host_city ~ pop + nightlights_trend, data = fr_trend, family = "binomial")
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
pscore <- predict(model_fr, type = "response")
#Match treatment and control groups using nearest neighbor matching
matched_data_fr <- matchit(host_city ~ pop + nightlights_trend, data = fr_trend, method = "nearest")
#Create control group
matched_data_with_distance_fr <- match.data(matched_data_fr)
control_group_fr <- subset(matched_data_with_distance_fr, host_city == 0)
treatment_group_fr <- subset(matched_data_with_distance_fr, host_city == 1)
Matching for Germany
#Calculate propensity scores using logistic regression
model_ge <- glm(host_city ~ pop + nightlights_trend , data = ge_trend, family = "binomial")
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
pscore_ge <- predict(model_ge, type = "response")
#Match treatment and control groups using nearest neighbor matching
matched_data_ge <- matchit(host_city ~ pop + nightlights_trend, data = ge_trend, method = "nearest")
#Create control group
matched_data_with_distance_ge <- match.data(matched_data_ge)
control_group_ge <- subset(matched_data_with_distance_ge, matched_data_with_distance_ge$host_city == 0)
control_group_ge
## Simple feature collection with 12 features and 13 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 6.625631 ymin: 49.41036 xmax: 13.96606 ymax: 54.43294
## Geodetic CRS: WGS 84
## First 10 features:
## ID name country.etc pop lat long capital host_city
## 66 66 bielefeld germany 335158 52.03 8.53 0 0
## 73 73 bochum germany 384208 51.48 7.20 0 0
## 74 74 bonn germany 315138 50.73 7.10 0 0
## 79 79 bremen germany 547915 53.08 8.81 0 0
## 132 132 dresden germany 489883 51.05 13.74 0 0
## 134 134 duisburg germany 502251 51.43 6.75 0 0
## 137 137 dusseldorf germany 573521 51.24 6.79 0 0
## 167 167 essen germany 596204 51.47 7.00 0 0
## 321 321 kiel germany 232422 54.32 10.12 0 0
## 385 385 mannheim germany 307579 49.50 8.47 0 0
## polygon_area nightlights_trend geometry
## 66 434358443 [m^2] -1.3742647 POLYGON ((8.377817 51.91487...
## 73 230170261 [m^2] -1.0246305 POLYGON ((7.102082 51.4105,...
## 74 208815341 [m^2] -1.1969309 POLYGON ((7.022535 50.63269...
## 79 2209553130 [m^2] -0.1388525 POLYGON ((8.481593 53.01104...
## 132 609222498 [m^2] -1.7446087 POLYGON ((13.57932 50.97494...
## 134 357437472 [m^2] -0.7709877 POLYGON ((6.625631 51.33338...
## 137 443348258 [m^2] -1.7230952 POLYGON ((6.688815 51.12437...
## 167 350006396 [m^2] -0.5092954 POLYGON ((6.894344 51.34757...
## 321 243902958 [m^2] -1.9664706 POLYGON ((10.03293 54.25071...
## 385 254295143 [m^2] -1.2576923 POLYGON ((8.41416 49.41036,...
## distance weights subclass
## 66 0.07317060 1 10
## 73 0.11990694 1 11
## 74 0.05461506 1 6
## 79 0.47666062 1 9
## 132 0.40004051 1 1
## 134 0.37388595 1 3
## 137 0.66508079 1 12
## 167 0.65561918 1 5
## 321 0.02363519 1 4
## 385 0.05052467 1 8
#Extract matched data
matched_data_ge_plot <- match.data(matched_data_ge)
#Set up the plotting layout
par(mfrow = c(1, 2))
#Plot 1: Before matching
boxplot(pop ~ host_city, data = ge_trend,
main = "Before Matching",
xlab = "Group", ylab = "Population")
#Plot 2: After matching
boxplot(pop ~ host_city, data = matched_data_ge_plot,
main = "After Matching",
xlab = "Group", ylab = "Population")
#Reset the plotting layout to default
par(mfrow = c(1, 1))
Matching for South Africa
# Calculate propensity scores using logistic regression
model_sa <- glm(host_city ~ pop + nightlights_trend, data = sa_trend, family = "binomial")
pscore_sa <- predict(model_sa, type = "response")
#Match treatment and control groups using nearest neighbor matching
library(MatchIt)
matched_data_sa <- matchit(host_city ~ pop + nightlights_trend, data = sa_trend, method = "nearest")
#Create control group
matched_data_with_distance_sa <- match.data(matched_data_sa)
control_group_sa <- subset(matched_data_with_distance_sa, matched_data_with_distance_sa$host_city == 0)
control_names_fr <- control_group_fr$name
#Grab the rows that match the specified names
france_control <- tidy_fr[tidy_fr$name %in% control_names_fr, ]
tidy_fr_treated <- tidy_fr %>%
filter(host_city == "1")
tidy_fr_all <- rbind(tidy_fr_treated, france_control)
control_names_ge <- control_group_ge$name
ge_control <- tidy_ge[tidy_ge$name %in% control_names_ge, ]
tidy_ge_treated <- tidy_ge %>%
filter(host_city == "1")
tidy_ge_all <- rbind(tidy_ge_treated, ge_control)
control_names_sa <- control_group_sa$name
sa_control <- tidy_sa[tidy_sa$name %in% control_names_sa, ]
tidy_sa_treated <- tidy_sa %>%
filter(host_city == "1")
tidy_sa_all <- rbind(tidy_sa_treated, sa_control)
Adding a dummy for the Treatment
#France
tidy_fr_all <- tidy_fr_all %>%
mutate(
treatment = if_else(host_city == "1" & Year > 1998, 1, 0),
year_diff = Year - 1998 # Change the reference year as needed
)
#Germany
tidy_ge_all <- tidy_ge_all %>%
mutate(
treatment = if_else(host_city == "1" & Year > 2006, 1, 0),
year_diff = Year - 2006 # Germany reference year
)
#South Africa
tidy_sa_all <- tidy_sa_all %>%
mutate(
treatment = if_else(host_city == "1" & Year > 2010, 1, 0),
year_diff = Year - 2010 # Change the reference year as needed
)
#Combine them all
tidy_all <- bind_rows(tidy_fr_all, tidy_ge_all, tidy_sa_all)
tidy_control_all <- bind_rows(france_control, ge_control, sa_control)
#Exclude years -4 and 4
tidy_all<- tidy_all %>%
filter(!(year_diff %in% c(-4, 4)))
Difference-in-Difference
did_all <- tidy_all %>%
mutate(treated = ifelse(host_city == 1, 1, 0),
post_period = ifelse(year_diff >= 0, 1, 0),
treated_post = treated * post_period)
did_model <- lm(Value ~ treated + post_period + treated_post, data = did_all)
table_output <- tab_model(did_model, show.ci = FALSE, show.se = TRUE, file = "model_table.html")
Synthetic Control Method
donors_est <- tidy_all %>%
filter(host_city == "0") %>%
distinct(name) %>%
pull(name)
#Leipzig caused problems in the data due not matching with a control
tidy_all <- tidy_all %>%
filter(!(name %in% c("leipzig")))
#Perform synthetic control analysis across time
results_all <- scdataMulti(tidy_all,
id.var = "name",
outcome.var = "Value",
treatment.var = "treatment",
time.var = "year_diff",
donors.est = list(donors_est),
features = list(c("Value")),
cov.adj = list(c("constant", "trend")),
effect = "time")
#Estimate across-time treatment effects with positive weights from the donor pool
res <- scest(results_all, w.constr = list("name" = "simplex"))
scplotMulti(res)
## $plot_out
#Estimate bootstrapped confidence intervals
respi <- scpi(results_all,
w.constr = list("name" = "simplex"),
cores = 4,
sims = 50,
e.method = "gaussian"
)
## ---------------------------------------------------------------
## Estimating Weights...
## Quantifying Uncertainty
## 5/50 iterations completed (10%) 10/50 iterations completed (20%) 15/50 iterations completed (30%) 20/50 iterations completed (40%) 25/50 iterations completed (50%) 30/50 iterations completed (60%) 35/50 iterations completed (70%) 40/50 iterations completed (80%) 45/50 iterations completed (90%) 50/50 iterations completed (100%)
scplotMulti(respi, type="series")
## $plot_out
Add City Fixed Effects
#Fit the fixed effects model
resid_model_city <- feols(Value ~ 1 | name, data = tidy_all)
#Extract the residuals from the model
residuals_city <- resid(resid_model_city)
#Add the residuals as a new column to the data frame
tidy_all$nl_city <- residuals_city
results_all_city <- scdataMulti(tidy_all,
id.var = "name",
outcome.var = "nl_city",
treatment.var = "treatment",
time.var = "year_diff",
donors.est = list(donors_est),
features = list(c("nl_city")),
cov.adj = list(c("constant", "trend")),
effect = "time")
#Estimate across-time treatment effects with positive weights from the donor pool
res_city <- scest(results_all_city, w.constr = list("name" = "simplex"))
scplotMulti(res_city)
## $plot_out
#Estimate bootstrapped confidence intervals
respi_city <- scpi(results_all_city,
w.constr = list("name" = "simplex"),
cores = 4,
sims = 50,
e.method = "gaussian")
## ---------------------------------------------------------------
## Estimating Weights...
## Quantifying Uncertainty
## 5/50 iterations completed (10%) 10/50 iterations completed (20%) 15/50 iterations completed (30%) 20/50 iterations completed (40%) 25/50 iterations completed (50%) 30/50 iterations completed (60%) 35/50 iterations completed (70%) 40/50 iterations completed (80%) 45/50 iterations completed (90%) 50/50 iterations completed (100%)
scplotMulti(respi_city, type="series")
## $plot_out
Add Time Fixed Effects
resid_model_time <- feols(Value ~ 1 | year_diff, data = tidy_all)
#Extract residuals
residuals_time <- resid(resid_model_time)
#Add residuals to tidy_all
tidy_all$nl_time <- residuals_time
results_all_year <- scdataMulti(tidy_all,
id.var = "name",
outcome.var = "nl_time",
treatment.var = "treatment",
time.var = "year_diff",
donors.est = list(donors_est),
features = list(c("nl_time")),
cov.adj = list(c("constant", "trend")),
effect = "time")
#Estimate across-time treatment effects with positive weights from the donor pool
res_time <- scest(results_all_year, w.constr = list("name" = "simplex"))
scplotMulti(res_time)
## $plot_out
#Estimate bootstrapped confidence intervals
respi_time <- scpi(results_all_year,
w.constr = list("name" = "simplex"),
cores = 4,
sims = 50,
e.method = "gaussian")
## ---------------------------------------------------------------
## Estimating Weights...
## Quantifying Uncertainty
## 5/50 iterations completed (10%) 10/50 iterations completed (20%) 15/50 iterations completed (30%) 20/50 iterations completed (40%) 25/50 iterations completed (50%) 30/50 iterations completed (60%) 35/50 iterations completed (70%) 40/50 iterations completed (80%) 45/50 iterations completed (90%) 50/50 iterations completed (100%)
scplotMulti(respi_time, type="series")
## $plot_out
Add City and Time Fixed Effects
resid_model_fe <- feols(Value ~ 1 | name + year_diff, data = tidy_all)
residuals_fe <- resid(resid_model_fe)
tidy_all$nl_fe <- residuals_fe
results_all_fe <- scdataMulti(tidy_all,
id.var = "name",
outcome.var = "nl_fe",
treatment.var = "treatment",
time.var = "year_diff",
donors.est = list(donors_est),
features = list(c("nl_fe")),
cov.adj = list(c("constant", "trend")),
effect = "time")
#Estimate across-time treatment effects with positive weights from the donor pool
res_fe <- scest(results_all_fe, w.constr = list("name" = "simplex"))
scplotMulti(res_fe)
## $plot_out
#Estimate bootstrapped confidence intervals
respi_fe <- scpi(results_all_fe,
w.constr = list("name" = "simplex"),
cores = 4,
sims = 50,
e.method = "gaussian")
## ---------------------------------------------------------------
## Estimating Weights...
## Quantifying Uncertainty
## 5/50 iterations completed (10%) 10/50 iterations completed (20%) 15/50 iterations completed (30%) 20/50 iterations completed (40%) 25/50 iterations completed (50%) 30/50 iterations completed (60%) 35/50 iterations completed (70%) 40/50 iterations completed (80%) 45/50 iterations completed (90%) 50/50 iterations completed (100%)
scplotMulti(respi_fe, type="series")
## $plot_out
Grab data to use in Table
filtered_tidy <- tidy_all[tidy_all$host_city == 1 & tidy_all$year_diff %in% c(1, 2, 3),
c("year_diff", "nl_fe", "Value", "nl_city", "nl_time")]
mean_nl_fe <- aggregate(nl_fe ~ year_diff, data = filtered_tidy, FUN = mean)
mean_Value <- aggregate(Value ~ year_diff, data = filtered_tidy, FUN = mean)
mean_nl_city <- aggregate(nl_city ~ year_diff, data = filtered_tidy, FUN = mean)
mean_nl_time <- aggregate(nl_time ~ year_diff, data = filtered_tidy, FUN = mean)
results_value <- as.data.frame(res$est.results$Y.post.fit)
results_fe <- as.data.frame(res_fe$est.results$Y.post.fit)
results_city <- as.data.frame(res_city$est.results$Y.post.fit)
results_time <- as.data.frame(res_time$est.results$Y.post.fit)
table_value <- cbind(results_value, mean_Value, as.data.frame(respi$inference.results$bounds$subgaussian) %>%
rename(lower = V1, upper = V2))
table_fe <- cbind(results_fe, mean_nl_fe,
as.data.frame(respi_fe$inference.results$bounds$subgaussian) %>%
rename(lower = V1, upper = V2))
table_city <- cbind(results_city, mean_nl_city, as.data.frame(respi_city$inference.results$bounds$subgaussian) %>%
rename(lower = V1, upper = V2))
table_time <- cbind(results_time, mean_nl_time, as.data.frame(respi_time$inference.results$bounds$subgaussian) %>%
rename(lower = V1, upper = V2))
table_value$difference <- table_value$Value - table_value$V1
table_fe$difference <- table_fe$nl_fe - table_fe$V1
table_city$difference <- table_city$nl_city - table_city$V1
table_time$difference <- table_time$nl_time - table_time$V1
# For table_city
table_city <- table_city %>%
rename(Untreated_NL = V1, Treated_nl = nl_city) %>%
rename(Year = year_diff,
Lower_Bound = lower,
Upper_Bound = upper,
Difference = difference) %>%
relocate(Year, Untreated_NL, Treated_nl, Difference)
# For table_time
table_time <- table_time %>%
rename(Untreated_NL = V1, Treated_nl = nl_time) %>%
rename(Year = year_diff,
Lower_Bound = lower,
Upper_Bound = upper,
Difference = difference) %>%
relocate(Year, Untreated_NL, Treated_nl, Difference)
# For table_fe
table_fe <- table_fe %>%
rename(Untreated_NL = V1, Treated_nl = nl_fe) %>%
rename(Year = year_diff,
Lower_Bound = lower,
Upper_Bound = upper,
Difference = difference) %>%
relocate(Year, Untreated_NL, Treated_nl, Difference)
# For table_value
table_value <- table_value %>%
rename(Untreated_NL = V1, Treated_nl = Value) %>%
rename(Year = year_diff,
Lower_Bound = lower,
Upper_Bound = upper,
Difference = difference) %>%
relocate(Year, Untreated_NL, Treated_nl, Difference)
# Load necessary libraries
library(knitr)
library(kableExtra)
# Function to format tables with a title and two decimal places
format_table <- function(df, title) {
df_formatted <- df %>%
mutate(across(where(is.numeric), round, 2)) # Format all numeric values to 2 decimals
kbl(df_formatted, caption = title, format = "html") %>%
kable_classic(full_width = FALSE)
}
# Print tables with captions and headers
cat("Table 1: Comparison in City Data\n")
## Table 1: Comparison in City Data
format_table(table_city, "Results with City Fixed Effects")
## Warning: There was 1 warning in `mutate()`.
## ℹ In argument: `across(where(is.numeric), round, 2)`.
## Caused by warning:
## ! The `...` argument of `across()` is deprecated as of dplyr 1.1.0.
## Supply arguments directly to `.fns` through an anonymous function instead.
##
## # Previously
## across(a:b, mean, na.rm = TRUE)
##
## # Now
## across(a:b, \(x) mean(x, na.rm = TRUE))
Results with City Fixed Effects
|
Year
|
Untreated_NL
|
Treated_nl
|
Difference
|
Lower_Bound
|
Upper_Bound
|
|
1
|
-0.38
|
0.45
|
0.83
|
-1.94
|
1.77
|
|
2
|
2.45
|
1.67
|
-0.79
|
-2.17
|
2.14
|
|
3
|
1.67
|
1.16
|
-0.51
|
-3.04
|
2.68
|
cat("\nTable 2: Comparison over Time\n")
##
## Table 2: Comparison over Time
format_table(table_time, "Results with Time Fixed Effects")
Results with Time Fixed Effects
|
Year
|
Untreated_NL
|
Treated_nl
|
Difference
|
Lower_Bound
|
Upper_Bound
|
|
1
|
-1.78
|
-1.13
|
0.65
|
-1.89
|
1.61
|
|
2
|
-0.37
|
-1.55
|
-1.18
|
-1.99
|
1.89
|
|
3
|
-0.73
|
-1.59
|
-0.86
|
-2.85
|
2.31
|
cat("\nTable 3: Comparison in Fixed Effects\n")
##
## Table 3: Comparison in Fixed Effects
format_table(table_fe, "Results with Time and City Fixed Effects")
Results with Time and City Fixed Effects
|
Year
|
Untreated_NL
|
Treated_nl
|
Difference
|
Lower_Bound
|
Upper_Bound
|
|
1
|
-0.66
|
0.14
|
0.80
|
-2.08
|
1.83
|
|
2
|
0.55
|
-0.28
|
-0.83
|
-2.28
|
2.14
|
|
3
|
0.21
|
-0.33
|
-0.54
|
-3.23
|
2.69
|
cat("\nTable 4: Comparison in Value\n")
##
## Table 4: Comparison in Value
format_table(table_value, "Results")
Results
|
Year
|
Untreated_NL
|
Treated_nl
|
Difference
|
Lower_Bound
|
Upper_Bound
|
|
1
|
42.47
|
43.23
|
0.77
|
-2.39
|
1.82
|
|
2
|
45.33
|
44.45
|
-0.88
|
-2.47
|
2.33
|
|
3
|
44.50
|
43.94
|
-0.55
|
-3.53
|
2.97
|