Sys.setenv(LANG = "en")
packages <-c('dplyr','magrittr','tidyverse','corrr',
'ggplot2','PerformanceAnalytics','ggtext','sf',
'ResourceSelection')
for(p in packages) if(p %in% rownames(installed.packages()) == F) { install.packages(p) }
for(p in packages) suppressPackageStartupMessages(library(p,quietly=T,character.only=T))
theme_set(theme_bw())
| Indicator Name | Description | Higher Values Indicate |
|---|---|---|
| Environmental Performance Index (EPI) | Measures a country’s environmental health and ecosystem vitality, but does not cover comprehensive ecological impacts like greenhouse gas emissions and material footprint. | Higher levels of environmental governance and ecological health |
| Global Peace Index (GPI) | Assesses negative peace in countries (absence of violence, conflict, crime, etc.). | Higher levels of violence, conflict, and crime (Lower GPI indicates higher peace) |
| Positive Peace Index (PPI) | Measures positive peace (constructive factors including social equity, education, economic opportunity). | Higher levels of social equity, economic opportunity, and governance stability |
| Greenhouse Gas Emissions | Direct and indirect greenhouse gas emissions by country, not included in EPI. | Higher greenhouse gas emissions (greater environmental burden) |
| Ecological Footprint | Measures a country’s resource consumption impact on the environment, including domestic and foreign resource use. | Higher resource consumption and ecological burden |
| Material Footprint | Assesses a country’s material consumption and its global supply chain impact on ecology. | Higher material resource consumption |
| Sustainable Development Goals (SDGs) | UN-established development goals that may not fully consider ecological sustainability issues. | Higher goal achievement (but may ignore ecological sustainability) |
| Wellbeing-Consumption Paradox | How to achieve social development needs without exceeding Earth’s ecological carrying capacity. | Greater contradiction (more conflict between high wellbeing and high consumption) |
| Ecologically Unequal Exchange | Global South resources heavily consumed by Global North countries, with environmental impacts transferred. | Higher levels of resource exploitation and ecological inequality |
| Negative Peace | Absence of direct violence like conflict, war, crime. | Higher negative peace (less violence and conflict) |
| Positive Peace | Long-term peace maintained through social, equity, economic, and governance means. | Higher positive peace (stronger social stability and governance) |
| ND Vulnerability (nd_v) | Degree to which a country or region is affected by environmental risks (e.g., climate change, resource depletion). | Higher vulnerability to environmental risks |
The Environmental Performance Index (EPI) is a country-level metric that ranks countries based on their environmental health and ecosystem vitality. The dataset contains EPI scores from 2010 to 2022, with scores ranging from 0 to 100, where higher scores indicate better environmental performance. The EPI was developed by researchers at Yale and Columbia Universities. In our dataset, we have both the overall EPI score and the climate change performance score (epi_cc_score) for each country.
The EPI score is calculated based on 40 performance indicators across 11 issue categories:
The formula for the overall EPI score is:
EPI = 0.4 * Environmental_Health + 0.6 * Ecosystem_Vitality
Where Environmental_Health and Ecosystem_Vitality are themselves weighted combinations of their sub-indicators.
The Climate Change score (epi_cc_score) specifically focuses on greenhouse gas emission trends and mitigation efforts, making up 36% of the Ecosystem Vitality score.
# The other dataset calculated is avaiable in the orginal csv
epi <- read.csv("https://raw.githubusercontent.com/sfield2/Peace_EnvironmentalSustainability_Paradox/refs/heads/main/DATA/EPI_2010_2022.csv",header=T,fileEncoding = 'UTF-8-BOM')%>%
select(c("country","year","EPI.new","EPI.cc"))%>%
setNames(c("country_name","year","epi_score","epi_cc_score"))
epi <- epi[,1:3]%>%
spread(year,epi_score)%>%
set_colnames(c("country_name","twoten","twotwelve","twofourteen","twosixteen",
"twoeighteen","twotwenty","twotwentytwo"))%>%
mutate(twoeleven = (twoten + twotwelve)/2)%>%
mutate(twothirteen = (twotwelve + twofourteen)/2)%>%
mutate(twofifteen = (twofourteen + twosixteen)/2)%>%
mutate(twoseventeen = (twosixteen + twoeighteen)/2)%>%
mutate(twonineteen = (twoeighteen + twotwenty)/2)%>%
mutate(twotwentyone = (twotwenty + twotwentytwo)/2)%>%
setNames(c("country_name","2010","2012","2014","2016",
"2018","2020","2022","2011","2013","2015",
"2017","2019","2021"))%>%
gather(year,epi_score,-country_name)
epi_merge <- epi %>%
rename(epi.score = epi_score)
epi %>%
filter(country_name %in% c("United States of America", "United Kingdom", "China", "Afghanistan")) %>%
filter(year == "2020")
The Global Peace Index (GPI) is a country-level metric that ranks countries based on their passive peacefulness (i.e., the absence of violence and conflict). The dataset contains GPI scores from 2008 to 2024, with scores ranging from 0 to 100, where lower scores indicate better peacefulness. The GPI was developed by the Institute for Economics and Peace (IEP). In our dataset, we have the overall GPI score for each country.
The GPI score is calculated based on 23 sub-indicators across 5 pillars:
# The other dataset calculated is avaiable in the orginal csv
gpi <- read.csv("https://raw.githubusercontent.com/sfield2/Peace_EnvironmentalSustainability_Paradox/refs/heads/main/DATA/GPI_2008_2024.csv",header=T,fileEncoding = 'UTF-8-BOM')%>%
select(c("country","year","Overall.Score"))%>%
setNames(c("country_name","year","gpi_score"))
# Print column names of gpi dataset
print(colnames(read.csv("https://raw.githubusercontent.com/sfield2/Peace_EnvironmentalSustainability_Paradox/refs/heads/main/DATA/GPI_2008_2024.csv",header=T,fileEncoding = 'UTF-8-BOM')))
## [1] "country" "geocode"
## [3] "year" "Overall.Score"
## [5] "Rank" "Access.To.Small.Arms"
## [7] "Armed.Services.Personnel.Rate" "Deaths.From.External.Conflict"
## [9] "Deaths.From.Internal.Conflict" "External.Conflicts.Fought"
## [11] "Homicide.Rate" "Incarceration.Rate"
## [13] "Intensity.Of.Internal.Conflict" "Internal.Conflicts.Fought"
## [15] "Militarisation" "Military.Expenditure....Gdp."
## [17] "Neighbouring.Countries.Relations" "Nuclear.And.Heavy.Weapons"
## [19] "Ongoing.Conflict" "Perceptions.Of.Criminality"
## [21] "Police.Rate" "Political.Instability"
## [23] "Political.Terror.Scale" "Refugees.And.Idps"
## [25] "Safety.And.Security" "Terrorism.Impact"
## [27] "Un.Peacekeeping.Funding" "Violent.Crime"
## [29] "Violent.Demonstrations" "Weapons.Exports"
## [31] "Weapons.Imports"
# The other dataset calculated is avaiable in the orginal csv
gpi_merge <- read.csv("https://raw.githubusercontent.com/sfield2/Peace_EnvironmentalSustainability_Paradox/refs/heads/main/DATA/GPI_2008_2024.csv",header=T,fileEncoding = 'UTF-8-BOM')%>%
select(c("country","year","Overall.Score",
"Access.To.Small.Arms", "Armed.Services.Personnel.Rate",
"Deaths.From.External.Conflict", "Deaths.From.Internal.Conflict",
"External.Conflicts.Fought", "Homicide.Rate", "Incarceration.Rate",
"Intensity.Of.Internal.Conflict", "Internal.Conflicts.Fought",
"Militarisation", "Military.Expenditure....Gdp.",
"Neighbouring.Countries.Relations", "Nuclear.And.Heavy.Weapons",
"Ongoing.Conflict", "Perceptions.Of.Criminality",
"Police.Rate", "Political.Instability",
"Political.Terror.Scale", "Refugees.And.Idps",
"Safety.And.Security", "Terrorism.Impact",
"Un.Peacekeeping.Funding", "Violent.Crime",
"Violent.Demonstrations", "Weapons.Exports",
"Weapons.Imports"))%>%
setNames(c("country_name","year","gpi.score",
"gpi.small.arms", "gpi.armed.services",
"gpi.external.deaths", "gpi.internal.deaths",
"gpi.external.conflicts", "gpi.homicide", "gpi.incarceration",
"gpi.internal.intensity", "gpi.internal.conflicts",
"gpi.militarisation", "gpi.military.gdp",
"gpi.neighbor.relations", "gpi.nuclear.weapons",
"gpi.ongoing.conflict", "gpi.criminality",
"gpi.police", "gpi.political.instability",
"gpi.terror.scale", "gpi.refugees",
"gpi.safety.security", "gpi.terrorism",
"gpi.peacekeeping", "gpi.violent.crime",
"gpi.demonstrations", "gpi.weapons.exports",
"gpi.weapons.imports"))
gpi %>%
filter(country_name %in% c("United States of America", "United Kingdom", "China", "Afghanistan")) %>%
filter(year == "2020")
The Positive Peace Index (PPI) is a country-level metric that ranks countries based on their Positive Peacefulness (i.e., the presence of social, economic, and governance factors that contribute to peace). The dataset contains PPI scores from 2009 to 2022, with scores ranging from 0 to 100, where lower scores indicate better peacefulness. The PPI was developed by the Institute for Economics and Peace (IEP). In our dataset, we have the overall PPI score for each country.
The PPI score is calculated based on 10 sub-indicators across 5 pillars:
# The other dataset calculated is avaiable in the orginal csv
ppi <- read.csv("https://raw.githubusercontent.com/sfield2/Peace_EnvironmentalSustainability_Paradox/refs/heads/main/DATA/PPI_2009_2022.csv",header=T,fileEncoding = 'UTF-8-BOM')%>%
select(c("Country","Year","PPI.Overall.Score"))%>%
setNames(c("country_name","year","ppi_score"))
# Print column names of gpi dataset
# print(colnames(read.csv("https://raw.githubusercontent.com/sfield2/Peace_EnvironmentalSustainability_Paradox/refs/heads/main/DATA/PPI_2009_2022.csv",header=T,fileEncoding = 'UTF-8-BOM')))
ppi_merge <- read.csv("https://raw.githubusercontent.com/sfield2/Peace_EnvironmentalSustainability_Paradox/refs/heads/main/DATA/PPI_2009_2022.csv",header=T,fileEncoding = 'UTF-8-BOM')%>%
select(c("Country","Year","PPI.Overall.Score",
"Acceptance.of.the.Rights.of.Others",
"Equitable.Distribution.of.Resources",
"Free.Flow.of.Information",
"Good.Relations.with.Neighbours",
"High.Levels.of.Human.Capital",
"Low.Levels.of.Corruption",
"Sound.Business.Environment",
"Well.Functioning.Government",
"International.tourism",
"Healthy.life.expectancy..HALE.",
"Researchers.in.R.D",
"Youth.not.in.employment..education.or.training..NEET.",
"GDP.per.capita",
"Exclusion.by.socio.economic.group",
"Quality.of.information",
"Public.sector.theft",
"Control.of.corruption",
"Government.effectiveness",
"Rule.of.law",
"Regulatory.quality",
"Group.grievance",
"Factionalised.elites",
"Government.openness.and.transparency",
"Law.to.support.equal.treatment.of.population.segments",
"Equality.of.opportunity",
"Freedom.of.the.press",
"Financial.institutions.index",
"Telecom..infrastructure.index..internet.mobile.broadband.",
"Freedom.to.Trade.Internationally",
"Inequality.adjusted.life.expectancy",
"Gender.inequality",
"Education.and.income.inequality"))%>%
setNames(c("country_name","year","ppi.score",
"ppi.rights.acceptance",
"ppi.resource.distribution",
"ppi.information.flow",
"ppi.neighbour.relations",
"ppi.human.capital",
"ppi.corruption.level",
"ppi.business.environment",
"ppi.government.function",
"ppi.tourism",
"ppi.life.expectancy",
"ppi.researchers",
"ppi.youth.neet",
"ppi.gdp.capita",
"ppi.socioeconomic.exclusion",
"ppi.info.quality",
"ppi.public.theft",
"ppi.corruption.control",
"ppi.govt.effectiveness",
"ppi.rule.of.law",
"ppi.regulatory.quality",
"ppi.group.grievance",
"ppi.elite.factionalization",
"ppi.govt.transparency",
"ppi.equal.treatment.law",
"ppi.opportunity.equality",
"ppi.press.freedom",
"ppi.financial.institutions",
"ppi.telecom.infrastructure",
"ppi.trade.freedom",
"ppi.life.expectancy.inequality",
"ppi.gender.inequality",
"ppi.education.income.inequality"))
ppi %>%
filter(country_name %in% c("United States of America", "United Kingdom", "China", "Afghanistan")) %>%
filter(year == "2020")
The ND Vulnerability (nd_v) is a country-level metric that ranks countries based on their vulnerability to environmental risks. The dataset contains ND Vulnerability scores from 1995 to 2022, with scores ranging from 0 to 100, where lower scores indicate better vulnerability. The ND Vulnerability was developed by the Notre Dame Global Adaptation Index (ND-GAIN). In our dataset, we have the overall ND Vulnerability score for each country.
# No calculated metric avaiable
nd_vuln <- read.csv("https://raw.githubusercontent.com/sfield2/Peace_EnvironmentalSustainability_Paradox/refs/heads/main/DATA/ND_Vulnerability_1995_2022.csv",header=T,fileEncoding = 'UTF-8-BOM')%>%
select(c("Name","X2010","X2011","X2012","X2013","X2014","X2015","X2016",
"X2017","X2018","X2019","X2020","X2021","X2022"))%>%
setNames(c("country_name","2010","2011","2012", "2013", "2014","2015", "2016",
"2017","2018", "2019", "2020","2021","2022"))%>%
gather(.,"year","nd_vuln_score",-country_name)
nd_vuln_merge <- read.csv("https://raw.githubusercontent.com/sfield2/Peace_EnvironmentalSustainability_Paradox/refs/heads/main/DATA/ND_Vulnerability_1995_2022.csv",header=T,fileEncoding = 'UTF-8-BOM')%>%
select(c("Name","X2010","X2011","X2012","X2013","X2014","X2015","X2016",
"X2017","X2018","X2019","X2020","X2021","X2022"))%>%
setNames(c("country_name","2010","2011","2012", "2013", "2014","2015", "2016",
"2017","2018", "2019", "2020","2021","2022"))%>%
gather(.,"year","nd_vuln.score",-country_name)
nd_vuln %>%
filter(country_name %in% c("United States of America", "United Kingdom", "China", "Afghanistan")) %>%
filter(year == "2020")
The Ecological Footprint (ecofoot_score) is a country-level metric that ranks countries based on their ecological footprint (the amount of biologically productive land and water required to produce the resources consumed and absorb the waste generated by a country’s population). The dataset contains Ecological Footprint scores from 2010 to 2022. The Ecological Footprint is measured in global hectares per capita, with higher values indicating a larger ecological footprint and thus more environmental impact. The Ecological Footprint was developed by the Global Footprint Network (GFN). In our dataset, we have the overall Ecological Footprint score for each country.
The Ecological Footprint is composed of several components:
gf_ecofoot <- read.csv("https://raw.githubusercontent.com/sfield2/Peace_EnvironmentalSustainability_Paradox/refs/heads/main/DATA/GF_ECOFOOT_PERCAPITA_2010_2022.csv",header=T,fileEncoding = 'UTF-8-BOM')%>%
select(c("Country.Name","year","Total_ecofoot"))%>%
setNames(c("country_name","year","ecofoot_score"))%>%
mutate(ecofoot_score = as.numeric(ecofoot_score))
# print(colnames(read.csv("https://raw.githubusercontent.com/sfield2/Peace_EnvironmentalSustainability_Paradox/refs/heads/main/DATA/GF_ECOFOOT_PERCAPITA_2010_2022.csv",header=T,fileEncoding = 'UTF-8-BOM')))
gf_ecofoot_merge <- read.csv("https://raw.githubusercontent.com/sfield2/Peace_EnvironmentalSustainability_Paradox/refs/heads/main/DATA/GF_ECOFOOT_PERCAPITA_2010_2022.csv",header=T,fileEncoding = 'UTF-8-BOM')%>%
select(c("Country.Name","year","Total_ecofoot",
"Built.up.Land", "Carbon", "Cropland",
"Fishing.Grounds", "Forest.Products", "Grazing.Land"))%>%
setNames(c("country_name","year","ecofoot.score",
"ecofoot.built.up", "ecofoot.carbon", "ecofoot.cropland",
"ecofoot.fishing", "ecofoot.forest", "ecofoot.grazing"))%>%
mutate(ecofoot.score = as.numeric(ecofoot.score))
gf_ecofoot %>%
filter(country_name %in% c("United States of America", "United Kingdom", "China", "Afghanistan")) %>%
filter(year == "2020")
The Material Footprint (matfoot_score) is a country-level metric that ranks countries based on their material footprint (the amount of resources consumed and waste generated by a country’s population). The dataset contains Material Footprint scores from 2008 to 2022, with scores not normalized, where lower scores indicate better material footprint. The Material Footprint was developed by the Global Footprint Network (GFN). In our dataset, we have the overall Material Footprint score for each country.
mat_foot <- read.csv("https://raw.githubusercontent.com/sfield2/Peace_EnvironmentalSustainability_Paradox/refs/heads/main/DATA/MAT_FOOT_1974_2024.csv",header=T,fileEncoding = 'UTF-8-BOM')%>%
select(c("Country","Category","X2008","X2009","X2010","X2011","X2012","X2013","X2014","X2015","X2016",
"X2017","X2018","X2019","X2020","X2021","X2022"))%>%
setNames(c("country_name","material_category","2008","2009","2010","2011","2012", "2013", "2014","2015", "2016",
"2017","2018", "2019", "2020","2021","2022"))%>%
gather(.,"year","value",-c("country_name","material_category"))
mat_foot %>%
filter(country_name %in% c("United States of America", "United Kingdom", "China", "Afghanistan")) %>%
filter(year == "2020")
Then, we need to normalize the material footprint score by the population. The formula is:
mf_score = (mf_score/1102311310)/population
where 1102311310 is the total material footprint in 2020 and population is the population of the country in 2020.
# step 1.1: read in population dataset
population <- read.csv("./DATA/POPULATIONS.csv",header=T,fileEncoding = 'UTF-8-BOM')%>%
select(c("Location","Time","TPopulation1Jan"))%>%
setNames(c("country_name","year","population_thousands"))%>%
mutate(population = population_thousands * 1000)%>%
select(c("country_name","year","population"))
# step 1.2: run calculations
mat_foot <- mat_foot%>%
group_by(country_name,year)%>%
summarize(mf_score = sum(value))%>%
merge(.,population,by=c("country_name","year"))%>%
mutate(mf_score = (mf_score/1102311310)/population)%>% ## 1102311310 in 1 gigaton
select(-c("population"))
mat_foot %>%
filter(country_name %in% c("United States of America", "United Kingdom", "China", "Afghanistan")) %>%
filter(year == "2020")
# Create mat_foot_merge with overall score and individual category scores
mat_foot_merge <- mat_foot %>%
select(c("country_name", "year", "mf_score")) %>%
rename(mat_foot.score = mf_score)
# Calculate normalized scores for each material category
mat_foot_categories <- read.csv("https://raw.githubusercontent.com/sfield2/Peace_EnvironmentalSustainability_Paradox/refs/heads/main/DATA/MAT_FOOT_1974_2024.csv",header=T,fileEncoding = 'UTF-8-BOM') %>%
select(c("Country","Category","X2008","X2009","X2010","X2011","X2012","X2013","X2014","X2015","X2016",
"X2017","X2018","X2019","X2020","X2021","X2022")) %>%
setNames(c("country_name","material_category","2008","2009","2010","2011","2012", "2013", "2014","2015", "2016",
"2017","2018", "2019", "2020","2021","2022")) %>%
gather(.,"year","value",-c("country_name","material_category")) %>%
merge(.,population,by=c("country_name","year")) %>%
mutate(value = (value/1102311310)/population) %>%
select(-c("population")) %>%
spread(material_category, value) %>%
rename(
mat_foot.biomass = Biomass,
mat_foot.fossil = `Fossil fuels`,
mat_foot.metal = `Metal ores`,
mat_foot.minerals = `Non-metallic minerals`
)
# Merge overall score with category scores
mat_foot_merge <- merge(mat_foot_merge, mat_foot_categories, by=c("country_name","year"))
Greenhouse Gas Emissions (GHG) are a key indicator of a country’s contribution to climate change. The dataset contains CO2 emissions per capita from 1990 to 2021. CO2 emissions are measured in metric tons per capita, with higher values indicating greater emissions and thus more environmental impact. The data comes from the World Bank’s World Development Indicators database. In our dataset, we have the CO2 emissions per capita for each country.
CO2 emissions are produced from: - Burning of fossil fuels (coal, oil, and natural gas) - Cement manufacturing - Gas flaring
Higher per capita emissions generally indicate: - More industrialized economies - Higher energy consumption - Greater reliance on fossil fuels - Higher standards of living (though this correlation is not absolute)
Lower per capita emissions may indicate: - Less industrialized economies - Greater use of renewable energy - More efficient energy use - Lower population density
ghg <- read.csv("https://raw.githubusercontent.com/sfield2/Peace_EnvironmentalSustainability_Paradox/refs/heads/main/DATA/co2_ghg.csv",header=T,fileEncoding = 'UTF-8-BOM')%>%
select(c("country","year","co2_per_capita"))%>%
setNames(c("country_name","year","ghg"))%>%
mutate(ghg = as.numeric(ghg))%>%
filter(complete.cases(.))
ghg_merge <- ghg%>%
select(c("country_name","year","ghg"))%>%
rename(ghg.score = ghg)
ghg %>%
filter(country_name %in% c("United States of America", "United Kingdom", "China", "Afghanistan")) %>%
filter(year == "2020")
################################################################################
### Step 3: Build full dataset ###
################################################################################
full_merge <- merge(epi_merge,ghg_merge,by=c("country_name","year"))%>%
merge(.,gpi_merge,by=c("country_name","year"))%>%
merge(.,ppi_merge,by=c("country_name","year"))%>%
merge(.,nd_vuln_merge,by=c("country_name","year"))%>%
merge(.,gf_ecofoot_merge,by=c("country_name","year"))%>%
merge(.,mat_foot_merge,by=c("country_name","year"))
full_merge%>%
filter(country_name %in% c("United States of America", "United Kingdom", "China", "Afghanistan")) %>%
filter(year == "2020")
full <- merge(epi,ghg,by=c("country_name","year"))%>%
merge(.,gpi,by=c("country_name","year"))%>%
merge(.,ppi,by=c("country_name","year"))%>%
merge(.,nd_vuln,by=c("country_name","year"))%>%
merge(.,gf_ecofoot,by=c("country_name","year"))%>%
merge(.,mat_foot,by=c("country_name","year"))
full%>%
filter(country_name %in% c("United States of America", "United Kingdom", "China", "Afghanistan")) %>%
filter(year == "2020")
full_scale <- full%>%
mutate_at(c("epi_score","ghg","gpi_score","ppi_score",
"nd_vuln_score","ecofoot_score",
"mf_score"), ~(scale(.) %>% as.vector))
full_scale%>%
filter(country_name %in% c("United States of America", "United Kingdom", "China", "Afghanistan")) %>%
filter(year == "2020")
full_mapping <- full_scale%>%
mutate(environmental_sustainability = ghg+ecofoot_score+mf_score)%>%
mutate(peace = nd_vuln_score+ppi_score+gpi_score)%>%
select(c("country_name","year","environmental_sustainability","peace"))%>%
subset(.,year==c("2022"))%>%
mutate(peace=peace*-1)%>%
mutate(tot=environmental_sustainability+peace)%>%
mutate(tot_rank = rank(-tot))%>%
setNames(c("name","year","environmental_sustainability","peace","tot","tot_rank"))
country <- sf::read_sf("./DATA/world-administrative-boundaries.shp")
country <- merge(country,full_mapping,by="name")
# Plot the world map with environmental sustainability and peace scores
ggplot(data = country) +
geom_sf(aes(fill = tot)) +
scale_fill_gradient2(low = "red", mid = "yellow", high = "green",
name = "Combined Score\n(Environmental Sustainability + Peace)") +
theme_minimal() +
labs(title = "Global Environmental Sustainability and Peace Scores",
subtitle = "2022") +
theme(plot.title = element_text(hjust = 0.5),
plot.subtitle = element_text(hjust = 0.5))
Correlation for 2020
library(corrplot)
library("PerformanceAnalytics")
# Create dataset for 2020
data_2020 <- full %>%
filter(year == "2020") %>%
select(epi_score, ghg, gpi_score, ppi_score, nd_vuln_score, ecofoot_score, mf_score) %>%
setNames(c("EPI", "GHG", "GPI", "PPI", "ND Vuln", "Eco Foot", "Mat Foot"))
# Create correlation plot with histograms for 2020
chart.Correlation(data_2020,
histogram = TRUE,
pch = 19,
method = "pearson")
data_2020
Correlation for Yearly Average
# Create dataset for yearly averages by country
yearly_country_avg <- full %>%
group_by(year, country_name) %>% # Group by both year and country
summarise(across(c(epi_score, ghg, gpi_score, ppi_score, nd_vuln_score, ecofoot_score, mf_score),
mean, na.rm = TRUE), .groups = "drop") %>%
select(-year, -country_name) %>% # Remove grouping columns
setNames(c("EPI", "GHG", "GPI", "PPI", "ND Vuln", "Eco Foot", "Mat Foot"))
# Create correlation plot with histograms for yearly country averages
chart.Correlation(yearly_country_avg,
histogram = TRUE,
pch = 19,
method = "pearson")
Conclusion
Based on the previous analysis and these two correlation plots, we can confirm that there is indeed a negative correlation between ecological sustainability scores and peace scores. While our earlier research showed some slight temporal variations in this pattern, the overall relationship/pattern remains consistent across time periods.
Our spatial analysis examines how geographic features affect peace and sustainability metrics by analyzing:
The following visualizations merge our sustainability and peace indices with geographic boundary data to explore these spatial patterns.
# Create yearly averages for full_mapping
full_mapping_avg <- full_scale %>%
mutate(environmental_sustainability = ghg + ecofoot_score + mf_score) %>%
mutate(peace = nd_vuln_score + ppi_score + gpi_score) %>%
group_by(country_name) %>%
summarise(
environmental_sustainability = mean(environmental_sustainability, na.rm = TRUE),
peace = mean(peace, na.rm = TRUE)
) %>%
setNames(c("name", "environmental_sustainability", "peace"))
# Load required libraries
library(sf)
library(tmap)
library(spdep)
library(igraph)
library(ggplot2)
library(dplyr)
library(GGally)
library(gridExtra)
# Read geographic data
country <- sf::read_sf("./DATA/world-administrative-boundaries.shp")
# Calculate basic geographic features
geo_features <- country %>%
mutate(
area = as.numeric(st_area(.)),
n_neighbors = lengths(st_touches(.)),
perimeter = as.numeric(st_length(st_cast(., "MULTILINESTRING")))
)
# Merge geographic features with average metrics
geo_analysis <- geo_features %>%
left_join(full_mapping_avg, by = c("name")) %>%
mutate(
neighbor_groups = cut(n_neighbors,
breaks = c(-Inf, 0, 2, 4, 6, 8, Inf),
labels = c("0", "1-2", "3-4", "5-6", "7-8", "8+"))
)
# Create plots
ggplot(geo_analysis, aes(x = n_neighbors)) +
geom_histogram(binwidth = 1,
fill = "#73B38D",
color = "black",
alpha = 0.7) +
labs(title = "Distribution of Number of Neighboring Countries",
subtitle = "Analysis of Global Border Relationships (Period Average)",
x = "Number of Neighbors",
y = "Number of Countries") +
theme_minimal() +
theme(
plot.title = element_text(size = 14, face = "bold"),
plot.subtitle = element_text(size = 12)
)
tm_shape(geo_features) +
tm_fill("n_neighbors",
title = "Number of Neighboring Countries",
palette = "viridis") +
tm_borders() +
tm_layout(main.title = "Country Neighborhood Analysis",
legend.position = c("right", "bottom"))
# Scatter plot: Environmental Sustainability vs Peace Index
ggplot(geo_analysis,
aes(x = environmental_sustainability,
y = peace,
size = as.numeric(area))) +
# Add reference lines
geom_hline(yintercept = 0, linetype = "dashed", color = "gray80") +
geom_vline(xintercept = mean(geo_analysis$environmental_sustainability, na.rm = TRUE),
linetype = "dashed", color = "gray80") +
# Main scatter plot
geom_point(aes(color = n_neighbors),
alpha = 0.8) +
# Overall trend line
geom_smooth(method = "lm",
color = "#204645",
se = TRUE,
alpha = 0.2) +
# Adjust size scale
scale_size_continuous(name = "Country Size",
range = c(0.5, 8),
breaks = c(1e11, 5e11, 1e12),
labels = c("Small", "Medium", "Large")) +
# Modified color gradient for better distinction of small values
scale_color_gradientn(
colors = c("#73B38D", "#5B9B7C", "#44836B", "#2D6B5A", "#155349", "#204645"),
values = scales::rescale(c(0, 2, 4, 6, 8, max(geo_analysis$n_neighbors, na.rm = TRUE))),
name = "Number of\nNeighboring\nCountries"
) +
# Labels and titles
labs(title = "Environmental Sustainability and Peace Relationship",
subtitle = "Analysis of country size and border connectivity effects",
x = "Environmental Sustainability Score",
y = "Peace Index Score",
caption = "Note: Dashed lines show trends for different country sizes; solid line shows overall trend") +
# Theme customization
theme_minimal() +
theme(
plot.title = element_text(size = 14, face = "bold"),
plot.subtitle = element_text(size = 12),
legend.position = "right",
legend.box = "vertical",
legend.key.size = unit(0.8, "lines"),
panel.grid.minor = element_line(color = "gray90"),
panel.grid.major = element_line(color = "gray85")
) +
# Add correlation annotations in bottom right
annotate("text",
x = max(geo_analysis$environmental_sustainability, na.rm = TRUE) * 0.95,
y = min(geo_analysis$peace, na.rm = TRUE) * 0.5,
label = paste("Correlations:\n",
"Small countries: ",
round(cor(subset(geo_analysis, as.numeric(area) <= 1e11)$environmental_sustainability,
subset(geo_analysis, as.numeric(area) <= 1e11)$peace,
use = "complete.obs"), 2), "\n",
"Medium countries: ",
round(cor(subset(geo_analysis, as.numeric(area) > 1e11 & as.numeric(area) <= 5e11)$environmental_sustainability,
subset(geo_analysis, as.numeric(area) > 1e11 & as.numeric(area) <= 5e11)$peace,
use = "complete.obs"), 2), "\n",
"Large countries: ",
round(cor(subset(geo_analysis, as.numeric(area) > 5e11)$environmental_sustainability,
subset(geo_analysis, as.numeric(area) > 5e11)$peace,
use = "complete.obs"), 2), "\n",
"Global: ",
round(cor(geo_analysis$environmental_sustainability,
geo_analysis$peace,
use = "complete.obs"), 2)),
size = 3.5,
hjust = 1, # Right alignment
color = "#204645")
This section explores how natural resources and environmental conditions affect the relationship between peace and sustainability:
The temporal dimension provides another crucial analytical perspective. We can combine this with spatiotemporal analysis to examine: column include: year, country_name, environmental sustainability related: ghg+ecofoot_score+mf_score peace related: nd_vuln_score ppi_score gpi_score
# Fit linear model to examine relationship between peace and sustainability over time
# First calculate peace and environmental sustainability scores
full <- full %>%
mutate(environmental_sustainability = ghg + ecofoot_score + mf_score,
peace = nd_vuln_score + ppi_score + gpi_score)
# Fit linear model to examine relationship between peace and sustainability over time
temporal_model <- lm(peace ~ environmental_sustainability + year + environmental_sustainability:year, data = full)
# Create predictions for visualization
pred_data <- expand.grid(
environmental_sustainability = seq(min(full$environmental_sustainability, na.rm=TRUE),
max(full$environmental_sustainability, na.rm=TRUE),
length=100),
year = unique(full$year)
)
pred_data$peace_pred <- predict(temporal_model, newdata=pred_data)
# Plot the relationship with interaction effects
ggplot(full, aes(x=environmental_sustainability, y=peace)) +
geom_point(alpha=0.3) +
geom_line(data=pred_data, aes(x=environmental_sustainability, y=peace_pred, color=factor(year))) +
facet_wrap(~year) +
labs(title="Peace-Sustainability Relationship Evolution",
subtitle="Showing yearly trends with interaction effects",
x="Environmental Sustainability Score",
y="Peace Score",
color="Year") +
theme_minimal() +
theme(legend.position="bottom",
plot.title = element_text(hjust=0.5),
plot.subtitle = element_text(hjust=0.5)) +
scale_color_viridis_d()
# Add statistical tests
summary(temporal_model)
##
## Call:
## lm(formula = peace ~ environmental_sustainability + year + environmental_sustainability:year,
## data = full)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.47082 -0.63527 0.01846 0.57523 2.98380
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.2349176 0.1141472 54.622 < 2e-16
## environmental_sustainability -0.0757447 0.0094064 -8.053 1.41e-15
## year2011 -0.0583580 0.1618223 -0.361 0.718
## year2012 -0.0627309 0.1604420 -0.391 0.696
## year2013 -0.0266080 0.1622201 -0.164 0.870
## year2014 -0.0299581 0.1624499 -0.184 0.854
## year2015 -0.0372815 0.1618641 -0.230 0.818
## year2016 0.0016267 0.1630946 0.010 0.992
## year2017 0.0037204 0.1631925 0.023 0.982
## year2018 -0.0118403 0.1633884 -0.072 0.942
## year2019 -0.0342915 0.1635418 -0.210 0.834
## year2020 -0.0315482 0.1629288 -0.194 0.846
## year2021 -0.0037102 0.1633818 -0.023 0.982
## year2022 0.0002469 0.1627866 0.002 0.999
## environmental_sustainability:year2011 0.0033335 0.0133551 0.250 0.803
## environmental_sustainability:year2012 0.0064560 0.0131909 0.489 0.625
## environmental_sustainability:year2013 -0.0019359 0.0136483 -0.142 0.887
## environmental_sustainability:year2014 -0.0043478 0.0138330 -0.314 0.753
## environmental_sustainability:year2015 -0.0025099 0.0137624 -0.182 0.855
## environmental_sustainability:year2016 -0.0074902 0.0140751 -0.532 0.595
## environmental_sustainability:year2017 -0.0062135 0.0139822 -0.444 0.657
## environmental_sustainability:year2018 -0.0061834 0.0140659 -0.440 0.660
## environmental_sustainability:year2019 -0.0065140 0.0141834 -0.459 0.646
## environmental_sustainability:year2020 -0.0098474 0.0146693 -0.671 0.502
## environmental_sustainability:year2021 -0.0099873 0.0144659 -0.690 0.490
## environmental_sustainability:year2022 -0.0094091 0.0143906 -0.654 0.513
##
## (Intercept) ***
## environmental_sustainability ***
## year2011
## year2012
## year2013
## year2014
## year2015
## year2016
## year2017
## year2018
## year2019
## year2020
## year2021
## year2022
## environmental_sustainability:year2011
## environmental_sustainability:year2012
## environmental_sustainability:year2013
## environmental_sustainability:year2014
## environmental_sustainability:year2015
## environmental_sustainability:year2016
## environmental_sustainability:year2017
## environmental_sustainability:year2018
## environmental_sustainability:year2019
## environmental_sustainability:year2020
## environmental_sustainability:year2021
## environmental_sustainability:year2022
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.011 on 1909 degrees of freedom
## (13 observations deleted due to missingness)
## Multiple R-squared: 0.2931, Adjusted R-squared: 0.2838
## F-statistic: 31.66 on 25 and 1909 DF, p-value: < 2.2e-16
# Test for non-linear relationships
temporal_model_quad <- lm(peace ~ environmental_sustainability + I(environmental_sustainability^2) +
year + environmental_sustainability:year, data=full)
anova(temporal_model, temporal_model_quad)
# 1. Create a copy of the data for PCA analysis
pca_data <- full %>%
# Remove any rows with NA values before PCA
drop_na(ghg, ecofoot_score, mf_score, nd_vuln_score, ppi_score, gpi_score)
# Standardize variables
pca_data <- pca_data %>%
mutate(across(c(ghg, ecofoot_score, mf_score, nd_vuln_score, ppi_score, gpi_score),
scale))
# Perform PCA separately
env_pca <- prcomp(~ ghg + ecofoot_score + mf_score, data = pca_data)
peace_pca <- prcomp(~ nd_vuln_score + ppi_score + gpi_score, data = pca_data)
# Add PCA scores to dataframe
pca_data <- pca_data %>%
mutate(
environmental_sustainability = as.vector(env_pca$x[,1]),
peace = as.vector(peace_pca$x[,1])
)
# Load required libraries
library(mgcv)
library(nlme)
# Panel data analysis using lme4 instead of plm
library(lme4)
panel_model <- lmer(peace ~ environmental_sustainability + (1|country_name) + (1|year),
data = pca_data)
# Time-varying coefficient model with fixed bandwidth
library(tvReg)
tvc_model <- tvLM(peace ~ environmental_sustainability,
data = pca_data,
bw = 0.5) # Using fixed bandwidth instead of cv.mspe
# Visualization
ggplot(pca_data, aes(x = environmental_sustainability, y = peace)) +
geom_point(alpha = 0.3) +
geom_smooth(method = "gam", aes(color = factor(year)), se = TRUE) +
facet_wrap(~year) +
labs(title = "Dynamic Peace-Sustainability Relationship",
subtitle = "Non-linear patterns across time (GAM fit)",
x = "Environmental Sustainability (PC1)",
y = "Peace Index (PC1)") +
theme_minimal() +
scale_color_viridis_d()
# Granger causality test using vars package
library(vars)
ts_data <- ts(data.frame(peace = pca_data$peace,
sustainability = pca_data$environmental_sustainability),
start = min(pca_data$year))
var_model <- VAR(ts_data, p = 2)
granger_test <- causality(var_model, cause = "sustainability")
yearly_models <- pca_data %>%
group_by(year) %>%
summarise(
correlation = cor(environmental_sustainability, peace, use = "complete.obs"),
n = n()
)
ggplot(yearly_models, aes(x = year, y = correlation)) +
geom_line() +
geom_point() +
labs(title = "Yearly Correlation between Peace and Environmental Sustainability",
x = "Year",
y = "Correlation Coefficient") +
theme_minimal()