.Rhistory

main.title.size = 0.8,
legend.bg.color = TRUE,
legend.outside = TRUE,
legend.title.size = 1,
legend.text.size = 1,
legend.outside.size = 0.2,
legend.outside.position = "bottom",
frame = FALSE
) +
tm_shape(adm_10m_fao_id_simple_rob) +
tm_borders(col = NA, lwd = 0.5)
return (index_map)
}
# ----------------------------------------------------------- Get Daniel's data
daniel_original <- read.csv("data_johannes.csv")
# Combine geometries and food loss data
country_data <- adm_10m_fao_id_simple %>%
left_join(daniel_original, by = c("ISO_3166_1" = "cntry_code"))
pal <- scico(n = 10, palette = "nuuk", direction = -1)
# Chunk 3
(plt_wheat <-
rasterize(country_data %>% filter(basketItem2 == "wheat") %>%
group_by(FAO_ID, country) %>%
summarise(foodloss_mean = mean(loss_percentage1)) ,
template_rast_5arcmin, "foodloss_mean")  %>%
create_index_map(., tocrs = "ESRI:54030",
index_main_title = "(a) Wheat",
index_label = "Average food loss (%)",
breakvals = c(0,2,4,6,8,10,15,20,Inf),
breaknames = c("<2", "2-4", "4-6","6-8", "8-10", "10-15", "15-20", ">20"),
colorpal = pal))
(plt_maize <-
rasterize(country_data %>% filter(basketItem2 == "maize") %>%
group_by(FAO_ID, country) %>%
summarise(foodloss_mean = mean(loss_percentage1)) ,
template_rast_5arcmin, "foodloss_mean")  %>%
create_index_map(., tocrs = "ESRI:54030",
index_main_title = "(b) Maize",
index_label = "Average food loss (%)",
breakvals = c(0,2,4,6,8,10,15,20,Inf),
breaknames = c("<2", "2-4", "4-6","6-8", "8-10", "10-15", "15-20", ">20"),
colorpal = pal))
(plt_rice <-
rasterize(country_data %>% filter(basketItem2 == "rice") %>%
group_by(FAO_ID, country) %>%
summarise(foodloss_mean = mean(loss_percentage1)) ,
template_rast_5arcmin, "foodloss_mean")  %>%
create_index_map(., tocrs = "ESRI:54030",
index_main_title = "(c) Rice",
index_label = "Average food loss (%)",
breakvals = c(0,2,4,6,8,10,15,20,Inf),
breaknames = c("<2", "2-4", "4-6","6-8", "8-10", "10-15", "15-20", ">20"),
colorpal = pal))
(plt_cereals <-
rasterize(country_data %>% filter(basketItem2 == "cereals") %>%
group_by(FAO_ID, country) %>%
summarise(foodloss_mean = mean(loss_percentage1)) ,
template_rast_5arcmin, "foodloss_mean")  %>%
create_index_map(., tocrs = "ESRI:54030",
index_main_title = "(d) Other cereals & pulses",
index_label = "Average food loss (%)",
breakvals = c(0,2,4,6,8,10,15,20,Inf),
breaknames = c("<2", "2-4", "4-6","6-8", "8-10", "10-15", "15-20", ">20"),
colorpal = pal))
(plt_fruits <-
rasterize(country_data %>% filter(basketItem2 == "fruits") %>%
group_by(FAO_ID, country) %>%
summarise(foodloss_mean = mean(loss_percentage1)) ,
template_rast_5arcmin, "foodloss_mean")  %>%
create_index_map(., tocrs = "ESRI:54030",
index_main_title = "(e) Fruits & vegetables",
index_label = "Average food loss (%)",
breakvals = c(0,2,4,6,8,10,15,20,Inf),
breaknames = c("<2", "2-4", "4-6","6-8", "8-10", "10-15", "15-20", ">20"),
colorpal = pal))
(plt_roots <-
rasterize(country_data %>% filter(basketItem2 == "roots") %>%
group_by(FAO_ID, country) %>%
summarise(foodloss_mean = mean(loss_percentage1)) ,
template_rast_5arcmin, "foodloss_mean")  %>%
create_index_map(., tocrs = "ESRI:54030",
index_main_title = "(f) Roots, tubers, & oil bearing crops",
index_label = "Average food loss (%)",
breakvals = c(0,2,4,6,8,10,15,20,Inf),
breaknames = c("<2", "2-4", "4-6","6-8", "8-10", "10-15", "15-20", ">20"),
colorpal = pal))
(plts_food_losses_commodities <- tmap_arrange(
plt_wheat, plt_maize,
plt_rice, plt_cereals,
plt_fruits, plt_roots,
ncol = 2))
# uncomment these!
# tmap_save(plts_food_losses_commodities,
#           "average_food_losses_commodities_rev.png",
#           dpi = 300)
#
#
# tmap_save(plts_food_losses_commodities,
#           "average_food_losses_commodities_rev.pdf",
#           dpi = 300)
# Chunk 4
(plt_farm <-
rasterize(country_data %>% filter(food_supply_stage == "Farm") %>%
group_by(FAO_ID, country) %>%
summarise(foodloss_mean = mean(loss_percentage1)) ,
template_rast_5arcmin, "foodloss_mean")  %>%
create_index_map(., tocrs = "ESRI:54030",
index_main_title = "(a) Farm",
index_label = "Average food loss (%)",
breakvals = c(0,2,4,6,8,10,15,20,Inf),
breaknames = c("<2", "2-4", "4-6","6-8", "8-10", "10-15", "15-20", ">20"),
colorpal = pal))
(plt_harvest <-
rasterize(country_data %>% filter(food_supply_stage == "Harvest") %>%
group_by(FAO_ID, country) %>%
summarise(foodloss_mean = mean(loss_percentage1)) ,
template_rast_5arcmin, "foodloss_mean")  %>%
create_index_map(., tocrs = "ESRI:54030",
index_main_title = "(b) Harvest",
index_label = "Average food loss (%)",
breakvals = c(0,2,4,6,8,10,15,20,Inf),
breaknames = c("<2", "2-4", "4-6","6-8", "8-10", "10-15", "15-20", ">20"),
colorpal = pal))
(plt_storage <-
rasterize(country_data %>% filter(food_supply_stage == "Storage") %>%
group_by(FAO_ID, country) %>%
summarise(foodloss_mean = mean(loss_percentage1)) ,
template_rast_5arcmin, "foodloss_mean")  %>%
create_index_map(., tocrs = "ESRI:54030",
index_main_title = "(c) Storage",
index_label = "Average food loss (%)",
breakvals = c(0,2,4,6,8,10,15,20,Inf),
breaknames = c("<2", "2-4", "4-6","6-8", "8-10", "10-15", "15-20", ">20"),
colorpal = pal))
(plt_transport <-
rasterize(country_data %>% filter(food_supply_stage == "Transport") %>%
group_by(FAO_ID, country) %>%
summarise(foodloss_mean = mean(loss_percentage1)) ,
template_rast_5arcmin, "foodloss_mean")  %>%
create_index_map(., tocrs = "ESRI:54030",
index_main_title = "(d) Transport",
index_label = "Average food loss (%)",
breakvals = c(0,2,4,6,8,10,15,20,Inf),
breaknames = c("<2", "2-4", "4-6","6-8", "8-10", "10-15", "15-20", ">20"),
colorpal = pal))
(plts_food_losses_supply_stage <- tmap_arrange(
plt_farm,
plt_harvest,
plt_storage,
plt_transport,
ncol = 2))
# uncomment these!
# tmap_save(plts_food_losses_supply_stage,
#           "average_food_losses_supply_stages_rev.png",
#           dpi = 300)
#
#
# tmap_save(plts_food_losses_supply_stage,
#           "average_food_losses_supply_stages_rev.pdf",
#           width = 180, height = 150, units = "mm")
# Chunk 5
# Extract data for explanatory variables, average values for 9 independent variables
explanatory_variables_averages <-
country_data %>%
group_by(FAO_ID) %>%
summarise(gni_mean = mean(gni),
agriShr_mean = mean(agriShr),
agriEmp_mean = mean(agriEmp),
years_of_schooling_mean = mean(msch),
electricityRural_mean = mean(electricityRural),
phoneSub_mean = mean(phoneSub),
polStab_mean = mean(polStab),
export_mean = mean(export),
median_incLevel_num = median(as.numeric(factor(incLevel, levels = c("L", "M", "H")))))
# plots
(plt_gni_mean <- rasterize(explanatory_variables_averages, template_rast_5arcmin, "gni_mean") %>%
create_index_map(., tocrs = "ESRI:54030",
index_main_title = "(a) GNI per capita",
index_label = "International $",
breakvals = c(0,1000,2000,4000,6000,10000,Inf),
breaknames = c("<1000", "1000-2000", "2000-4000","4000-6000", "6000-10000", ">10000"),
colorpal = pal))
(plt_agriShr_mean <- rasterize(explanatory_variables_averages, template_rast_5arcmin, "agriShr_mean") %>%
create_index_map(., tocrs = "ESRI:54030",
index_main_title = "(b) Agriculture shares of GDP",
index_label = "% of GDP",
breakvals = seq(0, 100, 10),
breaknames = c("<10", "10-20", "20-30", "30-40", "40-50",
"50-60", "60-70","70-80", "80-90", "90-100"),
colorpal = pal))
(plt_agriEmp_mean <- rasterize(explanatory_variables_averages, template_rast_5arcmin, "agriEmp_mean") %>%
create_index_map(., tocrs = "ESRI:54030",
index_main_title = "(c) Employment in agriculture",
index_label = "% of total employment",
breakvals = seq(0, 100, 10),
breaknames = c("<10", "10-20", "20-30", "30-40", "40-50",
"50-60", "60-70","70-80", "80-90", "90-100"),
colorpal = pal))
(plt_years_of_schooling_mean <- rasterize(explanatory_variables_averages, template_rast_5arcmin, "years_of_schooling_mean") %>%
create_index_map(., tocrs = "ESRI:54030",
index_main_title = "(d) Education",
index_label = "Years",
breakvals = c(0,2,4,6,8,10,Inf),
breaknames = c("<2", "2-4", "4-6", "6-8", "8-10",">10"),
colorpal = pal))
(plt_electricityRural_mean <- rasterize(explanatory_variables_averages, template_rast_5arcmin, "electricityRural_mean") %>%
create_index_map(., tocrs = "ESRI:54030",
index_main_title = "(e) Access to electricity in rural area",
index_label = "% of rural population",
breakvals = seq(0, 100, 10),
breaknames = c("<10", "10-20", "20-30", "30-40", "40-50",
"50-60", "60-70","70-80", "80-90", "90-100"),
colorpal = pal))
(plt_phoneSub_mean <- rasterize(explanatory_variables_averages, template_rast_5arcmin, "phoneSub_mean") %>%
create_index_map(., tocrs = "ESRI:54030",
index_main_title = "(f) Mobile cellular subscriptions",
index_label = "% of population",
breakvals = seq(0, 100, 10),
breaknames = c("<10", "10-20", "20-30", "30-40", "40-50",
"50-60", "60-70","70-80", "80-90", "90-100"),
colorpal = pal))
(plt_polStab_mean <- rasterize(explanatory_variables_averages, template_rast_5arcmin, "polStab_mean") %>%
"*" (100) %>%
create_index_map(., tocrs = "ESRI:54030",
index_main_title = "(g) Political stability and absence of violence",
index_label = "Index (low value=less stable; high value=more stable)",
breakvals = seq(0, 100, 10),
breaknames = c("<10", "10-20", "20-30", "30-40", "40-50",
"50-60", "60-70","70-80", "80-90", "90-100"),
colorpal = pal))
(plt_export_mean <- rasterize(explanatory_variables_averages, template_rast_5arcmin, "export_mean") %>%
create_index_map(., tocrs = "ESRI:54030",
index_main_title = "(h) Export volume index",
index_label = "Index (low value=less export; high value=more export)",
breakvals = c(0, 50 , 100, 150, 200, 250, 500, 1000, Inf),
breaknames = c("<50", "50-100", "100-150", "150-200", "200-250",
"250-500", "500-1000", ">1000"),
colorpal = pal))
# Income level
# plt_incLevel <- rasterize(explanatory_variables_averages, template_rast_5arcmin, "median_incLevel_num") %>%
#   create_index_map(., tocrs = "ESRI:54030",
#                    index_main_title = "(i) Income level",
#                    index_label = "Income level",
#                    breakvals = c(0.5, 1.5, 2.5, 3.5),
#                    breaknames = c("Low", "Middle", "High"),
#                    colorpal = pal)
#
# plt_incLevel
plts_figureAx <- tmap_arrange(plt_gni_mean, plt_agriShr_mean,
plt_agriEmp_mean, plt_years_of_schooling_mean,
plt_electricityRural_mean, plt_phoneSub_mean,
plt_polStab_mean, plt_export_mean,
ncol = 2) #  plt_incLevel removed
plts_figureAx
# # uncomment these!
# tmap_save(plts_figureAx,
#           filename = "average_explanatory_variables_rev.png",
#           dpi = 300)
#
#
# tmap_save(plts_figureAx,
#           filename = "average_explanatory_variables_rev.pdf",
#           width = 180, height = 150, units = "mm")
# Chunk 6
# Extract data for explanatory variables, average values for 9 independent variables
library(sf)
library(plotrix)
df_descriptive_using_countrydata <-
country_data %>%
st_drop_geometry() %>%
group_by(incLevel) %>%
summarise(foodloss_mean = mean(loss_percentage1),
gni_mean = mean(gni),
gni_se = std.error(gni),
agriShr_mean = mean(agriShr),
agriShr_se = std.error(agriShr),
agriEmp_mean = mean(agriEmp),
agriEmp_sd =  std.error(agriEmp),
years_of_schooling_mean = mean(msch),
years_of_schooling_se = std.error(msch),
electricityRural_mean = mean(electricityRural),
electricityRural_se =  std.error(electricityRural),
phoneSub_mean = mean(phoneSub),
phoneSub_se = std.error(phoneSub),
polStab_mean = mean(polStab),
polStab_se =  std.error(polStab),
export_mean = mean(export),
export_se =  std.error(export))
library(plotrix)
library(tidyverse)
# Get Daniel's data
daniel_original <- read.csv("data_johannes.csv")
# Explore descriptive statistics using the original data
df_se_descriptive <-
daniel_original %>%
group_by(incLevel, basketItem2, food_supply_stage) %>%
summarise(foodloss_mean = mean(loss_percentage1),
StandardError = std.error(loss_percentage1),
N = n()) %>%
mutate_if(is.numeric, round, digits = 3)
# alternatively, use data with polygons
# df_se_descriptive <-
#   country_data %>%
#   st_drop_geometry() %>%
#   group_by(incLevel, basketItem2, food_supply_stage) %>%
#   summarise(foodloss_mean = mean(loss_percentage1),
#             StandardError = std.error(loss_percentage1),
#             N = n()) %>%
#   mutate_if(is.numeric, round, digits = 3)
# Modify the table to have L and M countries side by side
df_se_descriptive_L <- df_se_descriptive %>%
filter(incLevel == "L") %>%
dplyr::select(-incLevel)
df_se_descriptive_M <- df_se_descriptive %>%
filter(incLevel == "M") %>%
dplyr::select(-incLevel)
# Join the modified tables
df_se_table <-
left_join(df_se_descriptive_L, df_se_descriptive_M,
by = c("basketItem2", "food_supply_stage"))
df_se_table
#View(df_se_table)
# Load necessary libraries
library(tidyverse)
library(fixest)
# Read the data
df_raw <- read_csv("data_johannes.csv")
# Add a new column ln_gni as the natural logarithm of gni
df_raw <- df_raw %>%
mutate(ln_gni = log(gni))
# Explore data: number of countries, crops and years
length(unique(df_raw$country)) # 81 countries
unique(df_raw$year) %>% sort() # years 2000-2020
# Group data by country and count the number of cases
df_raw %>%
group_by(country) %>%
summarise(Number_of_cases = n()) %>%
arrange(Number_of_cases)
# Create a list of countries with more than 80 observations
#my_cntrylist_over80 <-
my_countries_obs <-
df_raw %>%
group_by(country) %>%
summarise(Number_of_cases = n()) %>%
# filter(Number_of_cases > 80) %>%
arrange(Number_of_cases)# %>%   pull(country)
# Create a function to run feols excluding one country at a time
f_exclude_countries <- function(country_tobe_excluded) {
feols(loss_percentage1 ~ ln_gni + agriShr + agriEmp + msch +
electricityRural + phoneSub + polStab + export + regionalName + as.factor(year),
cluster = ~ regionalName,
data = df_raw %>% filter(country != country_tobe_excluded)) %>%
broom::tidy() %>%
dplyr::select(term, estimate, p.value) %>%
mutate(excluded_country = country_tobe_excluded)
}
# Example function calls
f_exclude_countries("Angola")
f_exclude_countries("Burundi")
# Use the function for all countries
my_countries <- unique(df_raw$country)
# Apply the function to all countries and store the results
results <- map_df(my_countries, f_exclude_countries)
# Widen the results and rename p-value columns ---> not sure if this is needed
results_wider <- results %>%
pivot_wider(names_from = term, values_from = c(estimate, p.value)) %>%
rename(
pval_intercept = `p.value_(Intercept)`,
pval_ln_gni = `p.value_ln_gni`,
pval_agriShr = `p.value_agriShr`,
pval_agriEmp = `p.value_agriEmp`,
pval_msch = `p.value_msch`,
pval_electricityRural = `p.value_electricityRural`,
pval_phoneSub = `p.value_phoneSub`,
pval_polStab = `p.value_polStab`,
pval_export = `p.value_export`
)
library(flextable)
library(modelsummary)
model <-   feols(loss_percentage1 ~ ln_gni + agriShr + agriEmp + msch +
electricityRural + phoneSub + polStab + export + regionalName + as.factor(year),
cluster = ~ regionalName,
data = df_raw)
# Create the table for the feols model
model_table <- modelsummary(model, stars = T, output = "flextable")
# uncomment this
#save_as_docx(model_table, path = ("selected_results.docx"))
results_arranged <- map_df(unique(my_countries_obs$country), f_exclude_countries)
results_wider_arranged <- results_arranged %>%
pivot_wider(names_from = term, values_from = c(estimate, p.value)) %>%
rename(
pval_intercept = `p.value_(Intercept)`,
pval_ln_gni = `p.value_ln_gni`,
pval_agriShr = `p.value_agriShr`,
pval_agriEmp = `p.value_agriEmp`,
pval_msch = `p.value_msch`,
pval_electricityRural = `p.value_electricityRural`,
pval_phoneSub = `p.value_phoneSub`,
pval_polStab = `p.value_polStab`,
pval_export = `p.value_export`
)
results_long_arranged <- results_wider_arranged %>%
pivot_longer(
cols = starts_with("estimate_"),
names_to = "variable",
values_to = "value"
) %>%
dplyr::select(excluded_country, variable, value) %>%
drop_na()
# Select only some variables
results_long_arranged <- results_long_arranged %>%
filter(variable %in% c("estimate_(Intercept)", "estimate_ln_gni",
"estimate_agriShr",
"estimate_agriEmp",
"estimate_msch",
"estimate_electricityRural",
"estimate_phoneSub",
"estimate_polStab",
"estimate_export"   ))
# Draw a combined plot
library(scico)
plt_combined_arranged <- ggplot(results_long_arranged, aes(x = factor(excluded_country, levels = unique(excluded_country)), y = value, color = variable, group = variable)) +
geom_line(size = 1.2) +
facet_grid(variable ~ .,
scales = "free_y",
labeller = as_labeller(c("estimate_(Intercept)" = "Intercept",
estimate_ln_gni = "Log GNI",
estimate_agriShr = "Agri Share",
estimate_agriEmp = "Agri Emp",
estimate_msch = "Mean Sch",
estimate_electricityRural = "Elec Acc",
estimate_phoneSub = "Phone Sub",
estimate_polStab = "Pol Stab",
estimate_export = "Exports"))) +
theme_gray() +
theme(axis.text.x = element_text(angle = 90, hjust = 1, size = 6),
axis.text.y = element_text(size = 8),
strip.text = element_text(size = 8),
legend.position = "none") +
scale_color_scico_d(palette = "roma") +
labs(title = "Impact of excluding a country on regression output",
x = "Excluded country countries with fewer observations         countries with more observations",
y = "Estimate")
plt_combined_arranged
# uncomment this
# ggsave("sensitivity_plot_x_axis_based_on_N_obs.png", plot = plt_combined_arranged, dpi = 300, width = 300, height = 225, units = "mm")
#
# ggsave("sensitivity_plot_x_axis_based_on_N_obs.pdf", plot = plt_combined_arranged)
# Prepare data for plotting
results_long_reg <- results_wider %>%
pivot_longer(
cols = starts_with("estimate_regional"),
names_to = "variable",
values_to = "value"
) %>%
dplyr::select(excluded_country, variable, value)
# Draw a combined plot
plt_combined_reg <- ggplot(results_long_reg, aes(x = excluded_country, y = value, color = variable)) +
geom_point() +
facet_grid(variable ~ .,
scales = "free_y",
labeller = as_labeller(c("estimate_regionalNameCentral & South Asia" = "Cent & S Asia",
"estimate_regionalNameEast & Southeast Asia" = "E & SE Asia",
"estimate_regionalNameEurope & North America" = "Eur & N Am",
"estimate_regionalNameLatin America & the Caribbean" = "Lat Am & Car",
"estimate_regionalNameMiddle East & North Africa" = "MEast & N Af",
"estimate_regionalNameSub-Saharan Africa" = "Sub-Sah Africa"))) +
theme_classic() +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "none") +
labs(title = "Impact of excluding a country on regression output",
x = "Excluded Country",
y = "Estimate")
plt_combined_reg
# ggsave("sensitivity_plot_reg.png", plot = plt_combined_reg, dpi = 300, width = 300, height = 225, units = "mm")
results_long_year <- results_wider %>%
pivot_longer(
cols = starts_with("estimate_as.factor"),
names_to = "variable",
values_to = "value"
) %>%
dplyr::select(excluded_country, variable, value)
plt_combined_year <- ggplot(results_long_year, aes(x = excluded_country, y = value, color = variable)) +
geom_point() +
facet_grid(variable ~ ., scales = "free_y", labeller = as_labeller(c("estimate_as.factor(year)2001" = "2001",
"estimate_as.factor(year)2002" = "2002",
"estimate_as.factor(year)2003" = "2003",
"estimate_as.factor(year)2004" = "2004",
"estimate_as.factor(year)2005" = "2005",
"estimate_as.factor(year)2006" = "2006",
"estimate_as.factor(year)2007" = "2007",
"estimate_as.factor(year)2008" = "2008",
"estimate_as.factor(year)2009" = "2009",
"estimate_as.factor(year)2010" = "2010",
"estimate_as.factor(year)2011" = "2011",
"estimate_as.factor(year)2012" = "2012",
"estimate_as.factor(year)2013" = "2013",
"estimate_as.factor(year)2014" = "2014",
"estimate_as.factor(year)2015" = "2015",
"estimate_as.factor(year)2016" = "2016",
"estimate_as.factor(year)2017" = "2017",
"estimate_as.factor(year)2018" = "2018",
"estimate_as.factor(year)2019" = "2019",
"estimate_as.factor(year)2020" = "2020"))) +
theme_classic() +
theme(axis.text.x = element_text(angle = 45, hjust = 1),
legend.position = "none") +
labs(title = "Impact of excluding a country on regression output",
x = "Excluded Country",
y = "Estimate")
plt_combined_year
# ggsave("sensitivity_plot_year.png", plot = plt_combined_year, dpi = 300, width = 300, height = 225, units = "mm")