Data Analysis - Conflict Project
David Leydet
2022-10-22
This project is an exploration of climate change and conflict in Pakistan. This project uses the following datasets:
Pakistan Conflict Data - Acquired from the Armed Conflict Location and Event Data Project (ACLED) on 21SEP2022. See https://acleddata.com/terms-of-use/
Regional Climate Data - ERA5 Dataset with monthly averaged temperature and precipitation data for January and July beginning in 2010 and ending in 2022. The spatial resolution is approximately 9 kilometers. Temperature data is calculated at 2m above the land surface. Precipitation data is total precipitation between forecast steps (monthly). Reanalysis model. See https://cds.climate.copernicus.eu/cdsapp#!/dataset/reanalysis-era5-land-monthly-means?tab=overview
Log Gross National Income Per Capital Variable - From: https://globaldatalab.org/shdi/metadata/lgnic/
Population Data acquired on 07NOV2022 from: https://datacommons.org/place/country/PAK?utm_medium=explore&mprop=count&popt=Person&hl=en
Pakistan Gross Domestic Product acquired on 07NOV2022 from: https://data.worldbank.org/indicator/NY.GDP.MKTP.CD?end=2021&locations=PK&start=2010
Initial Data Setup
## Set working directory
setwd("~/Desktop/University of Utah PhD /Research/r_code")
## Load packages we are using
library(mapview) #mapping package
library(raster) #raster data manipulation (Climate Data)
library(RColorBrewer) #color palettes for visualization
library(sf) #simple features for spatial data
library(tmap) #mapping package
library(viridis) #color palette for visualization
library(ncdf4) #working with netCDF files (Climate Data)
library(leaflet) #basemaps for mapview
library(ggplot2) #better figures
library(ggcorrplot) #Load the correlation plot package
library(plotly) #interactive figures
library(maps) #mapping
library(kableExtra) #creating better tables and outputs
library(dplyr) #count and data functions
library(reshape2) ## Package used to reformat data - wide to long
library(tidyverse) ##Formatting dataframes, merge, and join
library(stargazer) ##Formatting model outputs to tables
library(pscl) ##Used to calculate pseudo r^2 values for log regression models (poisson)
library(janitor) ##Used to count/provide summaries for dataframes
library(jtools) ##Used to produce aesthetically pleasing model output tables
library(huxtable) ##Used in conjunction with jtools to export model outputs
library(flextable) ##Needed to knit. linked to the janitor libraryPakistan Conflict Data
## Initial
## Pakistan Conflict Data
pak.con = read.csv("../data/conflict/pak_conflict.csv")
## Change data to a factor
#year don't set as a factor
## pak.con$year = factor(pak.con$year)
#interaction - this is the interaction code between actors (see ACLED codebook)
pak.con$interaction = factor(pak.con$interaction)
#inter1 - this is the type of actor for actor 1 (see ACLED codebook)
pak.con$inter1 = factor(pak.con$inter1)
#inter2 - this is the type of actor for actor 2(see ACLED codebook)
pak.con$inter2 = factor(pak.con$inter2)
##Summary of the data
table1 = tabyl(pak.con, event_type)
table1| event_type | n | percent |
|---|---|---|
| Battles | 5397 | 0.0701 |
| Explosions/Remote violence | 5585 | 0.0726 |
| Protests | 57819 | 0.751 |
| Riots | 3517 | 0.0457 |
| Strategic developments | 805 | 0.0105 |
| Violence against civilians | 3825 | 0.0497 |
##Convert Table 1 into a more user friendly figure
table1.convert = table1 %>% kbl() %>% kable_classic()
#Visualize
table1.convert| event_type | n | percent |
|---|---|---|
| Battles | 5397 | 0.0701383 |
| Explosions/Remote violence | 5585 | 0.0725815 |
| Protests | 57819 | 0.7514035 |
| Riots | 3517 | 0.0457062 |
| Strategic developments | 805 | 0.0104616 |
| Violence against civilians | 3825 | 0.0497089 |
##Save File
##To work properly load the following libraries for images**
library(magick)
library(webshot)
save_kable(table1.convert, "./figures/event_table_summary1.jpg")## Read in country shape file
countries = st_read("../data/conflict/ne_50m_admin_0_countries/ne_50m_admin_0_countries.shp")## Reading layer `ne_50m_admin_0_countries' from data source
## `/Users/davidleydet/Desktop/University of Utah PhD /Research/data/conflict/ne_50m_admin_0_countries/ne_50m_admin_0_countries.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 241 features and 94 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -180 ymin: -89.99893 xmax: 180 ymax: 83.59961
## Geodetic CRS: WGS 84## Read in administrative boundary shape file
pak.bound = st_read("../data/conflict/pak_adm_wfp_20220909_shp/pak_admbnda_adm1_wfp_20220909.shp")## Reading layer `pak_admbnda_adm1_wfp_20220909' from data source
## `/Users/davidleydet/Desktop/University of Utah PhD /Research/data/conflict/pak_adm_wfp_20220909_shp/pak_admbnda_adm1_wfp_20220909.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 7 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 60.8786 ymin: 23.69468 xmax: 77.83397 ymax: 37.08942
## Geodetic CRS: WGS 84pak.bound.adm2 = st_read("../data/conflict/pak_adm_wfp_20220909_shp/pak_admbnda_adm2_wfp_20220909.shp")## Reading layer `pak_admbnda_adm2_wfp_20220909' from data source
## `/Users/davidleydet/Desktop/University of Utah PhD /Research/data/conflict/pak_adm_wfp_20220909_shp/pak_admbnda_adm2_wfp_20220909.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 160 features and 14 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 60.8786 ymin: 23.69468 xmax: 77.83397 ymax: 37.08942
## Geodetic CRS: WGS 84## Set as a spatial feature for visualization
pak.con.sf = st_as_sf(pak.con,
coords = c("longitude", "latitude"),
crs = 4326) #WGS84## Using grepl to subset based on keywords
## keywords: water, rain, snow, river, flood
## ArcPro's output for the definition query is ~4400. ***What is the delta??***
pak.con.water.filter = dplyr::filter(pak.con.sf, grepl("water|flood|rain|river|snow", pak.con$notes, ignore.case = TRUE))
str(pak.con.water.filter)## Classes 'sf' and 'data.frame': 4761 obs. of 22 variables:
## $ data_id : int 9491475 9491460 9491461 9491470 9491741 9491742 9491743 9491775 9491394 9491739 ...
## $ event_id_cnty : chr "PAK77601" "PAK77593" "PAK77558" "PAK77585" ...
## $ event_date : chr "15-Sep-22" "13-Sep-22" "13-Sep-22" "13-Sep-22" ...
## $ year : int 2022 2022 2022 2022 2022 2022 2022 2022 2022 2022 ...
## $ event_type : chr "Protests" "Protests" "Protests" "Protests" ...
## $ sub_event_type: chr "Peaceful protest" "Peaceful protest" "Peaceful protest" "Peaceful protest" ...
## $ actor1 : chr "Protesters (Pakistan)" "Protesters (Pakistan)" "Protesters (Pakistan)" "Protesters (Pakistan)" ...
## $ assoc_actor_1 : chr "" "" "" "" ...
## $ inter1 : Factor w/ 8 levels "1","2","3","4",..: 6 6 6 6 6 6 6 6 6 6 ...
## $ actor2 : chr "" "" "" "" ...
## $ assoc_actor_2 : chr "" "" "" "" ...
## $ inter2 : Factor w/ 9 levels "0","1","2","3",..: 1 1 1 1 1 1 1 1 1 1 ...
## $ interaction : Factor w/ 43 levels "10","11","12",..: 36 36 36 36 36 36 36 36 36 36 ...
## $ admin1 : chr "Sindh" "Balochistan" "Sindh" "Sindh" ...
## $ admin2 : chr "Badin" "Sohbatpur" "Sujawal" "Umerkot" ...
## $ admin3 : chr "Badin" "Sohbatpur" "Mirpur Bathoro" "Pithoro" ...
## $ location : chr "Talhar" "Sohbatpur" "Mirpur Bathoro" "Shadi Palli" ...
## $ source : chr "Daily Regional Times (Pakistan)" "Daily Regional Times (Pakistan)" "Daily Regional Times (Pakistan)" "Daily Regional Times (Pakistan)" ...
## $ source_scale : chr "Subnational" "National" "Subnational" "Subnational" ...
## $ notes : chr "On 15 September 2022, residents held a protest demonstration in Talhar town (Badin, Sindh), demanding drainage "| __truncated__ "On 13 September 2022, flood affectees held a protest demonstration in Sohbatpur town (Sohbatpur, Sindh), demand"| __truncated__ "On 13 September 2022, flood affectees held a protest demonstration in Mirpur Bathoro town (Sujawal, Sindh), dem"| __truncated__ "On 13 September 2022, flood affectees held a protest demonstration in Shadi Palli town (Umerkot, Sindh), demand"| __truncated__ ...
## $ fatalities : int 0 0 0 0 0 0 0 0 0 0 ...
## $ geometry :sfc_POINT of length 4761; first list element: 'XY' num 68.8 24.9
## - attr(*, "sf_column")= chr "geometry"
## - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA NA NA NA NA NA NA ...
## ..- attr(*, "names")= chr [1:21] "data_id" "event_id_cnty" "event_date" "year" ...pak.con.water.filter.sf = st_as_sf(pak.con.water.filter)
st_write(pak.con.water.filter.sf,
dsn = "./figures/pak_water_con_sf.shp",
delete_layer = TRUE) ##Overrite the existing fileInitial Visualization Exploration
## Total events
ggplot(pak.con.sf) +
geom_sf(color = "firebrick3", alpha = 0.4)
con.total.map =
tm_shape(countries , bbox = st_bbox(pak.con.sf)) +
tm_borders(col = "gray") +
tm_shape(pak.con.sf) +
tm_symbols(col = "event_type",
alpha = 0.4,
size = 0.2,
title.col = "Conflicts By Event Type") +
tm_layout(main.title = "Pakistan Conflicts: 2010 - 2022",
legend.outside = TRUE)
con.total.map
pdf("./figures/pak_conflict_total")
con.total.map
dev.off()con.water.map =
tm_shape(countries , bbox = st_bbox(pak.con.sf)) +
tm_borders(col = "gray") +
tm_shape(pak.con.water.filter) +
tm_symbols(col = "event_type",
alpha = 0.4,
size = 0.2,
title.col = "Water Conflicts By Event Type") +
tm_layout(main.title = "Pakistan Water Conflicts: 2010 - 2022",
legend.outside = TRUE)
con.water.map
water.con.by.year.map =
tm_shape(countries , bbox = st_bbox(pak.con.sf)) +
tm_borders(col = "gray") +
tm_shape(pak.con.water.filter) +
tm_facets(by = "year") +
tm_symbols(col = "dodgerblue1",
alpha = 0.5,
size = 0.2) +
tm_layout(main.title = "Pakistan Water Conflicts: 2010 - 2022",
main.title.size = 0.75,
legend.outside = FALSE,
legend.position = c("right", "bottom"))
water.con.by.year.map
pdf("./figures/pak_water_con_by_year")
water.con.by.year.map
dev.off()con.by.year.map =
tm_shape(countries , bbox = st_bbox(pak.con.sf)) +
tm_borders(col = "gray") +
tm_shape(pak.con.sf) +
tm_facets(by = "year") +
tm_symbols(col = "firebrick3",
alpha = 0.5,
size = 0.2) +
tm_layout(main.title = "Pakistan Conflicts: 2010 - 2022",
legend.outside = TRUE)
con.by.year.map
pdf("./figures/pak_con_by_year")
con.by.year.map
dev.off()# Subset to exclude 0 fatality events
pak.con.sf.fat = subset(pak.con.sf, fatalities > 10)
## NOT WORKING PROPERLY*****
tm_shape(countries , bbox = st_bbox(pak.con.sf.fat)) +
tm_borders(col = "gray") +
tm_shape(pak.con.sf.fat) +
tm_symbols(col = "fatalities",
palette = "Reds",
alpha = 0.4,
title.col = "Number of Fatalities",
size = 0.2) +
tm_layout(main.title = "Pakistan Conflicts: 2010 - 2022",
legend.outside = TRUE) 
tm_shape(pak.bound) +
tm_borders() +
tm_fill(col = "ADM1_EN",
title = "Pakistan Provinces")
## Generate a new color palette
mycol = turbo(n = 8, begin = 0.5, end = 1, direction = 1)
con.by.district =
tm_shape(pak.bound) +
tm_borders() +
tm_fill(col = "ADM1_EN",
title = "Pakistan Provinces") +
tm_shape(pak.con.sf.fat) +
tm_symbols(col = "fatalities",
palette = mycol,
alpha = 0.4,
title.col = "Number of Fatalities",
size = 0.2) +
tm_layout(main.title = "Pakistan Conflicts by Province: 2010 - 2022",
legend.outside = TRUE) +
tm_compass(position = c("left", "top"))
con.by.district
pdf("./figures/con_by_district")
con.by.district
dev.off()tm_shape(pak.bound.adm2) +
tm_borders() +
tm_fill(col = "ADM2_EN",
title = "Pakistan Divisions") +
tm_layout(legend.show = FALSE)
Initial Data Analysis
## Create the counts by year
year.2022 = as.numeric(nrow(filter(pak.con, year == 2022)))
year.2021 = as.numeric(nrow(filter(pak.con, year == 2021)))
year.2020 = as.numeric(nrow(filter(pak.con, year == 2020)))
year.2019 = as.numeric(nrow(filter(pak.con, year == 2019)))
year.2018 = as.numeric(nrow(filter(pak.con, year == 2018)))
year.2017 = as.numeric(nrow(filter(pak.con, year == 2017)))
year.2016 = as.numeric(nrow(filter(pak.con, year == 2016)))
year.2015 = as.numeric(nrow(filter(pak.con, year == 2015)))
year.2014 = as.numeric(nrow(filter(pak.con, year == 2014)))
year.2013 = as.numeric(nrow(filter(pak.con, year == 2013)))
year.2012 = as.numeric(nrow(filter(pak.con, year == 2012)))
year.2011 = as.numeric(nrow(filter(pak.con, year == 2011)))
year.2010 = as.numeric(nrow(filter(pak.con, year == 2010)))
## Create a new data frame for the counts
# vector of years
years = c("2010", "2011", "2012", "2013", "2014", "2015", "2016", "2017", "2018", "2019", "2020", "2021", "2022")
# vector of years as a number for plotting purposes
years.num = seq(from = 2010, to = 2022, by = 1)
# vector of conflict counts
num.conflicts = c(year.2010, year.2011, year.2012, year.2013, year.2014, year.2015, year.2016, year.2017, year.2018, year.2019, year.2020, year.2021, year.2022)
# Combine the vectors into a data frame
con.by.year = data.frame (years, num.conflicts)
# Combine the vectors using the years as a number
con.by.year.num = data.frame(years.num, num.conflicts)
# Visualize the Count
# How does one add a trendline?
num.conflicts.sctrplot =
ggplot(data = con.by.year.num, aes(x = years.num, y = num.conflicts)) +
geom_point(shape = 18, size = 3, color = "red") +
geom_smooth(method = "lm",
se = FALSE) +
xlab("Year") +
ylab("Number of Conflicts") +
labs(title = "Pakistan Conflicts by Year") +
scale_x_discrete(limits = c(years.num)) +
theme_bw()
num.conflicts.sctrplot## `geom_smooth()` using formula = 'y ~ x'
pdf("./figures/con_by_year_scatterplot")
num.conflicts.sctrplot
dev.off()##Count the number of water conflicts by year
water.year.2022 = as.numeric(nrow(filter(pak.con.water.filter, year == 2022)))
water.year.2021 = as.numeric(nrow(filter(pak.con.water.filter, year == 2021)))
water.year.2020 = as.numeric(nrow(filter(pak.con.water.filter, year == 2020)))
water.year.2019 = as.numeric(nrow(filter(pak.con.water.filter, year == 2019)))
water.year.2018 = as.numeric(nrow(filter(pak.con.water.filter, year == 2018)))
water.year.2017 = as.numeric(nrow(filter(pak.con.water.filter, year == 2017)))
water.year.2016 = as.numeric(nrow(filter(pak.con.water.filter, year == 2016)))
water.year.2015 = as.numeric(nrow(filter(pak.con.water.filter, year == 2015)))
water.year.2014 = as.numeric(nrow(filter(pak.con.water.filter, year == 2014)))
water.year.2013 = as.numeric(nrow(filter(pak.con.water.filter, year == 2013)))
water.year.2012 = as.numeric(nrow(filter(pak.con.water.filter, year == 2012)))
water.year.2011 = as.numeric(nrow(filter(pak.con.water.filter, year == 2011)))
water.year.2010 = as.numeric(nrow(filter(pak.con.water.filter, year == 2010)))
##Automate this: function script?
##Function Works!
##Automate loop??
##mycountbyyear = function(x) {
##as.numeric(nrow(filter(x, year == 2010)))
#}
##testyear = mycountbyyear(pak.con.water.filter)
##Create a vector of water conflicts by year
water.num.conflicts = c(water.year.2010, water.year.2011, water.year.2012, water.year.2013, water.year.2014, water.year.2015, water.year.2016, water.year.2017, water.year.2018, water.year.2019, water.year.2020, water.year.2021, water.year.2022)
# Combine the vectors into a data frame. Need to use years.num in order to plot a trend line
water.con.by.year = data.frame (years.num, water.num.conflicts)
##Plot it!
num.water.conflicts.sctrplot =
ggplot(data = water.con.by.year, aes(x = years.num, y = water.num.conflicts)) +
geom_point(shape = 18, size = 3, color = "blue") +
geom_smooth(method = NULL,
se = TRUE,
color = "black",
linetype = "dashed",
size = 0.5) +
xlab("Year") +
ylab("Number of Water Conflicts") +
labs(title = "Pakistan Water Conflicts by Year") +
scale_x_discrete(limits = c(years.num)) +
theme_bw() ## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.## Warning: Continuous limits supplied to discrete scale.
## ℹ Did you mean `limits = factor(...)` or `scale_*_continuous()`?num.water.conflicts.sctrplot## `geom_smooth()` using method = 'loess' and formula = 'y ~ x'
pdf("./figures/water_con_by_year_scatterplot")
num.water.conflicts.sctrplot
dev.off()Pakistan GNI Data
## Initial Read
## ****Note - the years are from 1990 - 2019****
## Log Gross National Income per capita in thousands of US Dollars (2011 PPP)
## Ask Simon about time series data
pak.gni = read.csv("../data/conflict/pak_subnat_gdp.csv")
## Convert it from wide format to long format
pak.gni.long = melt(pak.gni, id.vars = "Region")
## Change the column name to year
colnames(pak.gni.long)[2] = "year"
## Remove the "x" from the years
pak.gni.long$year = gsub("X", "", as.factor(pak.gni.long$year))
## Subset the temporal timeframe: 2010 - 2022 (2019 in this case)
###pak.gni.time.subset = subset(pak.gni.long, Factor == "2010")pak.gni.plot =
ggplot(data = pak.gni.long, aes(x = year, y = value, color = Region)) +
geom_point() +
xlab("Year") +
ylab("GNI per Capita (Thousands USD)") +
# xlim(2010, 2019) +
theme(axis.text.x = element_text(angle = 90))
theme_bw()## List of 94
## $ line :List of 6
## ..$ colour : chr "black"
## ..$ linewidth : num 0.5
## ..$ linetype : num 1
## ..$ lineend : chr "butt"
## ..$ arrow : logi FALSE
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_line" "element"
## $ rect :List of 5
## ..$ fill : chr "white"
## ..$ colour : chr "black"
## ..$ linewidth : num 0.5
## ..$ linetype : num 1
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ text :List of 11
## ..$ family : chr ""
## ..$ face : chr "plain"
## ..$ colour : chr "black"
## ..$ size : num 11
## ..$ hjust : num 0.5
## ..$ vjust : num 0.5
## ..$ angle : num 0
## ..$ lineheight : num 0.9
## ..$ margin : 'margin' num [1:4] 0points 0points 0points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : logi FALSE
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ title : NULL
## $ aspect.ratio : NULL
## $ axis.title : NULL
## $ axis.title.x :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 1
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 2.75points 0points 0points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.title.x.top :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 0
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0points 0points 2.75points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.title.x.bottom : NULL
## $ axis.title.y :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 1
## ..$ angle : num 90
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0points 2.75points 0points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.title.y.left : NULL
## $ axis.title.y.right :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 0
## ..$ angle : num -90
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0points 0points 0points 2.75points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.text :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : chr "grey30"
## ..$ size : 'rel' num 0.8
## ..$ hjust : NULL
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.text.x :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 1
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 2.2points 0points 0points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.text.x.top :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : num 0
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0points 0points 2.2points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.text.x.bottom : NULL
## $ axis.text.y :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : num 1
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0points 2.2points 0points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.text.y.left : NULL
## $ axis.text.y.right :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : num 0
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0points 0points 0points 2.2points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ axis.ticks :List of 6
## ..$ colour : chr "grey20"
## ..$ linewidth : NULL
## ..$ linetype : NULL
## ..$ lineend : NULL
## ..$ arrow : logi FALSE
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_line" "element"
## $ axis.ticks.x : NULL
## $ axis.ticks.x.top : NULL
## $ axis.ticks.x.bottom : NULL
## $ axis.ticks.y : NULL
## $ axis.ticks.y.left : NULL
## $ axis.ticks.y.right : NULL
## $ axis.ticks.length : 'simpleUnit' num 2.75points
## ..- attr(*, "unit")= int 8
## $ axis.ticks.length.x : NULL
## $ axis.ticks.length.x.top : NULL
## $ axis.ticks.length.x.bottom: NULL
## $ axis.ticks.length.y : NULL
## $ axis.ticks.length.y.left : NULL
## $ axis.ticks.length.y.right : NULL
## $ axis.line : list()
## ..- attr(*, "class")= chr [1:2] "element_blank" "element"
## $ axis.line.x : NULL
## $ axis.line.x.top : NULL
## $ axis.line.x.bottom : NULL
## $ axis.line.y : NULL
## $ axis.line.y.left : NULL
## $ axis.line.y.right : NULL
## $ legend.background :List of 5
## ..$ fill : NULL
## ..$ colour : logi NA
## ..$ linewidth : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ legend.margin : 'margin' num [1:4] 5.5points 5.5points 5.5points 5.5points
## ..- attr(*, "unit")= int 8
## $ legend.spacing : 'simpleUnit' num 11points
## ..- attr(*, "unit")= int 8
## $ legend.spacing.x : NULL
## $ legend.spacing.y : NULL
## $ legend.key :List of 5
## ..$ fill : chr "white"
## ..$ colour : logi NA
## ..$ linewidth : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ legend.key.size : 'simpleUnit' num 1.2lines
## ..- attr(*, "unit")= int 3
## $ legend.key.height : NULL
## $ legend.key.width : NULL
## $ legend.text :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : 'rel' num 0.8
## ..$ hjust : NULL
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ legend.text.align : NULL
## $ legend.title :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : num 0
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ legend.title.align : NULL
## $ legend.position : chr "right"
## $ legend.direction : NULL
## $ legend.justification : chr "center"
## $ legend.box : NULL
## $ legend.box.just : NULL
## $ legend.box.margin : 'margin' num [1:4] 0cm 0cm 0cm 0cm
## ..- attr(*, "unit")= int 1
## $ legend.box.background : list()
## ..- attr(*, "class")= chr [1:2] "element_blank" "element"
## $ legend.box.spacing : 'simpleUnit' num 11points
## ..- attr(*, "unit")= int 8
## $ panel.background :List of 5
## ..$ fill : chr "white"
## ..$ colour : logi NA
## ..$ linewidth : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ panel.border :List of 5
## ..$ fill : logi NA
## ..$ colour : chr "grey20"
## ..$ linewidth : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ panel.spacing : 'simpleUnit' num 5.5points
## ..- attr(*, "unit")= int 8
## $ panel.spacing.x : NULL
## $ panel.spacing.y : NULL
## $ panel.grid :List of 6
## ..$ colour : chr "grey92"
## ..$ linewidth : NULL
## ..$ linetype : NULL
## ..$ lineend : NULL
## ..$ arrow : logi FALSE
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_line" "element"
## $ panel.grid.major : NULL
## $ panel.grid.minor :List of 6
## ..$ colour : NULL
## ..$ linewidth : 'rel' num 0.5
## ..$ linetype : NULL
## ..$ lineend : NULL
## ..$ arrow : logi FALSE
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_line" "element"
## $ panel.grid.major.x : NULL
## $ panel.grid.major.y : NULL
## $ panel.grid.minor.x : NULL
## $ panel.grid.minor.y : NULL
## $ panel.ontop : logi FALSE
## $ plot.background :List of 5
## ..$ fill : NULL
## ..$ colour : chr "white"
## ..$ linewidth : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ plot.title :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : 'rel' num 1.2
## ..$ hjust : num 0
## ..$ vjust : num 1
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0points 0points 5.5points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ plot.title.position : chr "panel"
## $ plot.subtitle :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : num 0
## ..$ vjust : num 1
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 0points 0points 5.5points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ plot.caption :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : 'rel' num 0.8
## ..$ hjust : num 1
## ..$ vjust : num 1
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 5.5points 0points 0points 0points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ plot.caption.position : chr "panel"
## $ plot.tag :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : 'rel' num 1.2
## ..$ hjust : num 0.5
## ..$ vjust : num 0.5
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ plot.tag.position : chr "topleft"
## $ plot.margin : 'margin' num [1:4] 5.5points 5.5points 5.5points 5.5points
## ..- attr(*, "unit")= int 8
## $ strip.background :List of 5
## ..$ fill : chr "grey85"
## ..$ colour : chr "grey20"
## ..$ linewidth : NULL
## ..$ linetype : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_rect" "element"
## $ strip.background.x : NULL
## $ strip.background.y : NULL
## $ strip.clip : chr "inherit"
## $ strip.placement : chr "inside"
## $ strip.text :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : chr "grey10"
## ..$ size : 'rel' num 0.8
## ..$ hjust : NULL
## ..$ vjust : NULL
## ..$ angle : NULL
## ..$ lineheight : NULL
## ..$ margin : 'margin' num [1:4] 4.4points 4.4points 4.4points 4.4points
## .. ..- attr(*, "unit")= int 8
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ strip.text.x : NULL
## $ strip.text.y :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : NULL
## ..$ angle : num -90
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## $ strip.switch.pad.grid : 'simpleUnit' num 2.75points
## ..- attr(*, "unit")= int 8
## $ strip.switch.pad.wrap : 'simpleUnit' num 2.75points
## ..- attr(*, "unit")= int 8
## $ strip.text.y.left :List of 11
## ..$ family : NULL
## ..$ face : NULL
## ..$ colour : NULL
## ..$ size : NULL
## ..$ hjust : NULL
## ..$ vjust : NULL
## ..$ angle : num 90
## ..$ lineheight : NULL
## ..$ margin : NULL
## ..$ debug : NULL
## ..$ inherit.blank: logi TRUE
## ..- attr(*, "class")= chr [1:2] "element_text" "element"
## - attr(*, "class")= chr [1:2] "theme" "gg"
## - attr(*, "complete")= logi TRUE
## - attr(*, "validate")= logi TRUEpak.gni.plot
## Save as a .pdf
ggsave("./figures/gni_by_year_district.pdf", pak.gni.plot)##2022 GDP has not been published yet, so we estimated based on an average gdp growth of 3.7%
pak.gdp = read.csv("../data/conflict/pak_gdp_2010_2022.csv")
pak.gdp.plot =
ggplot(data = pak.gdp, aes(x = year, y = gdp)) +
geom_point(col = "forestgreen") +
geom_line(linetype = "dashed",
col = "forestgreen") +
labs(title = "Pakistan GDP (2010-2022)") +
xlab("Year") +
ylab("GDP USD($)")
pak.gdp.plot
#May want to scale the GDPpdf("./figures/pak_gdp_by_year")
pak.gdp.plot
dev.off()Pakistan Climate Data
Temperature Data
## Read in the ERA5 climate data
pak.temp = raster("../data/conflict/pak_clim.nc", varname = "t2m")
## Look at the data
##26 bands (layers)
##13 years worth of data * 2 months(Jan, Jul) = 26 observations (bands?)
## Syntax variable = raster("filepath", band = #) to extract the specific band you need
## library(ncdf4) this package allows you to read into the metadata
## Could use stack() to look at all of the bands
pak.temp## class : RasterLayer
## band : 1 (of 26 bands)
## dimensions : 201, 251, 50451 (nrow, ncol, ncell)
## resolution : 0.1, 0.1 (x, y)
## extent : 54.95, 80.05, 19.95, 40.05 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : pak_clim.nc
## names : X2.metre.temperature
## z-value : 2010-01-01
## zvar : t2m## CRS Check
## Coordinate Reference System:
## Deprecated Proj.4 representation: +proj=longlat +datum=WGS84 +no_defs ## Temperature is in Kelvin (how do I convert the scale to Celsius? Kelvin Temp - 273.15)
## January Temperature. How to plot July?
mytemp.pal = brewer.pal(n = 9, name = "OrRd")
plot(pak.temp,
main = "January Temperature - 2010",
col = mytemp.pal)
plot(st_geometry(pak.bound),
border = "black",
add = TRUE)
##Not working correctly now? 20221107
pak.temp.masked = mask(pak.temp, mask = pak.bound)
plot(pak.temp.masked,
main = "Pakistan Jan 2010 Temperature - Masked",
col = mytemp.pal)
plot(st_geometry(pak.bound), add = TRUE)
##Use stack() to load all of the temperature bands
temp.stack = stack("../data/conflict/pak_clim.nc", varname = "t2m")
##View the stack data
temp.stack## class : RasterStack
## dimensions : 201, 251, 50451, 26 (nrow, ncol, ncell, nlayers)
## resolution : 0.1, 0.1 (x, y)
## extent : 54.95, 80.05, 19.95, 40.05 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## names : X2010.01.01, X2010.07.01, X2011.01.01, X2011.07.01, X2012.01.01, X2012.07.01, X2013.01.01, X2013.07.01, X2014.01.01, X2014.07.01, X2015.01.01, X2015.07.01, X2016.01.01, X2016.07.01, X2017.01.01, ...
## Date/time : 2010-01-01 - 2022-07-01 (range)##January 2010 stack
temp.stack.jan2010 = temp.stack[[1]]
temp.stack.jan2010## class : RasterLayer
## band : 1 (of 26 bands)
## dimensions : 201, 251, 50451 (nrow, ncol, ncell)
## resolution : 0.1, 0.1 (x, y)
## extent : 54.95, 80.05, 19.95, 40.05 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : pak_clim.nc
## names : X2010.01.01
## Date/time : 2010-01-01
## zvar : t2m##Crop the temperature stack by the Pakistan Boundary
temp.stack.crop = crop(temp.stack, pak.bound)
temp.stack.crop## class : RasterBrick
## dimensions : 134, 170, 22780, 26 (nrow, ncol, ncell, nlayers)
## resolution : 0.1, 0.1 (x, y)
## extent : 60.85, 77.85, 23.65, 37.05 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : memory
## names : X2010.01.01, X2010.07.01, X2011.01.01, X2011.07.01, X2012.01.01, X2012.07.01, X2013.01.01, X2013.07.01, X2014.01.01, X2014.07.01, X2015.01.01, X2015.07.01, X2016.01.01, X2016.07.01, X2017.01.01, ...
## min values : 242.9010, 262.5190, 239.7853, 264.1235, 238.9400, 261.9992, 239.1240, 262.7146, 240.7158, 262.4167, 240.7767, 265.0068, 240.8986, 265.3265, 242.3812, ...
## max values : 294.2676, 310.3270, 293.8248, 310.1821, 293.5142, 311.2782, 293.4889, 311.3368, 292.2054, 311.5416, 293.4717, 311.2161, 295.2015, 311.4714, 293.9168, ...
## time : 2010-01-01, 2022-07-01 (min, max)##Color Template for Temperature
##Create a function to add the Pakistan Level 2 boundaries to the figure
addborder = function(){
plot(as_Spatial(pak.bound),
add = TRUE)
}
plot(temp.stack.crop[[1]],
col = mytemp.pal,
zlim = c(250, 320), ##zlim makes sure the scale is the same
addfun = addborder)
pak.temp.stack.mask = mask(temp.stack, mask = pak.bound)
plot(pak.temp.stack.mask,
col = mytemp.pal,
zlim = c(250, 320), ##zlim makes sure the scale is the same
addfun = addborder)
##Use the temperature stack (masked) to calculate the mean average temperature for the country
temp.avg.all = cellStats(pak.temp.stack.mask, mean)
##Visualize the data - this differs slighlt from running the entire raster data (without the mask)
#plot(temp.avg.all,
#type = "p")
temp.avg.all## X2010.01.01 X2010.07.01 X2011.01.01 X2011.07.01 X2012.01.01 X2012.07.01
## 283.0038 300.5736 281.0138 300.7901 279.7675 301.7816
## X2013.01.01 X2013.07.01 X2014.01.01 X2014.07.01 X2015.01.01 X2015.07.01
## 280.9438 301.1685 280.7888 301.3565 281.6683 300.4983
## X2016.01.01 X2016.07.01 X2017.01.01 X2017.07.01 X2018.01.01 X2018.07.01
## 282.7455 301.3774 281.1049 300.9322 282.4930 301.4175
## X2019.01.01 X2019.07.01 X2020.01.01 X2020.07.01 X2021.01.01 X2021.07.01
## 281.5669 301.5670 279.5454 301.7934 280.6591 301.7339
## X2022.01.01 X2022.07.01
## 280.8029 299.3014##Shape the average temperature values into a data frame
temp.avg.df = as.data.frame(temp.avg.all)
temp.avg.df| temp.avg.all |
|---|
| 283 |
| 301 |
| 281 |
| 301 |
| 280 |
| 302 |
| 281 |
| 301 |
| 281 |
| 301 |
| 282 |
| 300 |
| 283 |
| 301 |
| 281 |
| 301 |
| 282 |
| 301 |
| 282 |
| 302 |
| 280 |
| 302 |
| 281 |
| 302 |
| 281 |
| 299 |
##Write to csv to transform....need to find a better way to do this!!!!***
write.csv(temp.avg.df,
"../data/conflict/pak_avg_temp")##Read in the converted csv
pak.temp.avg = read.csv("../data/conflict/pak_avg_temp_converted.csv")
pak.temp.avg| year | avg_jan_temp | avg_jul_temp | avg_jan_celsius | avg_jul_celsius |
|---|---|---|---|---|
| 2010 | 283 | 301 | 9.85 | 27.4 |
| 2011 | 281 | 301 | 7.86 | 27.6 |
| 2012 | 280 | 302 | 6.62 | 28.6 |
| 2013 | 281 | 301 | 7.79 | 28 |
| 2014 | 281 | 301 | 7.64 | 28.2 |
| 2015 | 282 | 300 | 8.52 | 27.3 |
| 2016 | 283 | 301 | 9.6 | 28.2 |
| 2017 | 281 | 301 | 7.95 | 27.8 |
| 2018 | 282 | 301 | 9.34 | 28.3 |
| 2019 | 282 | 302 | 8.42 | 28.4 |
| 2020 | 280 | 302 | 6.4 | 28.6 |
| 2021 | 281 | 302 | 7.51 | 28.6 |
| 2022 | 281 | 299 | 7.65 | 26.2 |
avg.temp.plot =
ggplot(data = pak.temp.avg) +
geom_line(aes(x=year, y = avg_jan_celsius), col = "blue") +
geom_line(aes(x=year, y = avg_jul_celsius), col = "red") +
xlab("Year") +
ylab("Temperature (Celsius)")
avg.temp.plot
Precipitation Data
## Read in the ERA5 climate data
## Value is in meters
##26 bands (layers)
##13 years worth of data * 2 months(Jan and Jul) = 26 observations (bands?)
## Syntax variable = raster("filepath", band = #) to extract the specific band you need
## library(ncdf4) this package allows you to read into the metadata
## Could use stack() to look at all of the bands
pak.precip = raster("../data/conflict/pak_clim.nc", varname = "tp")
## Look at the data
pak.precip## class : RasterLayer
## band : 1 (of 26 bands)
## dimensions : 201, 251, 50451 (nrow, ncol, ncell)
## resolution : 0.1, 0.1 (x, y)
## extent : 54.95, 80.05, 19.95, 40.05 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## source : pak_clim.nc
## names : Total.precipitation
## z-value : 2010-01-01
## zvar : tp##January precipitation
myprecip.pal = brewer.pal(n = 9, name = "Blues")
plot(pak.precip,
main = "January Precipitation - 2010",
col = myprecip.pal)
plot(st_geometry(pak.bound.adm2),
border = "black",
add = TRUE)
##Use stack() to load all of the precipitation bands
##REMINDER This is total precipitation in meters
precip.stack = stack("../data/conflict/pak_clim.nc", varname = "tp")
##View the stack data
precip.stack## class : RasterStack
## dimensions : 201, 251, 50451, 26 (nrow, ncol, ncell, nlayers)
## resolution : 0.1, 0.1 (x, y)
## extent : 54.95, 80.05, 19.95, 40.05 (xmin, xmax, ymin, ymax)
## crs : +proj=longlat +datum=WGS84 +no_defs
## names : X2010.01.01, X2010.07.01, X2011.01.01, X2011.07.01, X2012.01.01, X2012.07.01, X2013.01.01, X2013.07.01, X2014.01.01, X2014.07.01, X2015.01.01, X2015.07.01, X2016.01.01, X2016.07.01, X2017.01.01, ...
## Date/time : 2010-01-01 - 2022-07-01 (range)pak.precip.stack.mask = mask(precip.stack, mask = pak.bound)
plot(pak.precip.stack.mask,
col = myprecip.pal,
zlim = c(0, 0.0030), ##zlim makes sure the scale is the same
addfun = addborder)
##Use the precipitation stack (masked) to calculate the mean average temperature for the country
precip.avg.all = cellStats(pak.precip.stack.mask, mean)
##Visualize the data - this differs slighlt from running the entire raster data (without the mask)
#plot(precip.avg.all,
#type = "p")
precip.avg.all## X2010.01.01 X2010.07.01 X2011.01.01 X2011.07.01 X2012.01.01 X2012.07.01
## 0.0005839226 0.0043177244 0.0003963835 0.0020611769 0.0009019366 0.0010772682
## X2013.01.01 X2013.07.01 X2014.01.01 X2014.07.01 X2015.01.01 X2015.07.01
## 0.0003545558 0.0018138459 0.0003568661 0.0015595845 0.0007123514 0.0038597275
## X2016.01.01 X2016.07.01 X2017.01.01 X2017.07.01 X2018.01.01 X2018.07.01
## 0.0006265805 0.0022160874 0.0019144569 0.0026588202 0.0002786720 0.0019427409
## X2019.01.01 X2019.07.01 X2020.01.01 X2020.07.01 X2021.01.01 X2021.07.01
## 0.0012742865 0.0020474729 0.0018617504 0.0016054658 0.0003895901 0.0021279386
## X2022.01.01 X2022.07.01
## 0.0017411640 0.0063039410##Shape the average temperature values into a data frame
precip.avg.df = as.data.frame(precip.avg.all)
precip.avg.df| precip.avg.all |
|---|
| 0.000584 |
| 0.00432 |
| 0.000396 |
| 0.00206 |
| 0.000902 |
| 0.00108 |
| 0.000355 |
| 0.00181 |
| 0.000357 |
| 0.00156 |
| 0.000712 |
| 0.00386 |
| 0.000627 |
| 0.00222 |
| 0.00191 |
| 0.00266 |
| 0.000279 |
| 0.00194 |
| 0.00127 |
| 0.00205 |
| 0.00186 |
| 0.00161 |
| 0.00039 |
| 0.00213 |
| 0.00174 |
| 0.0063 |
##Write to csv to transform....need to find a better way to do this!!!!***
write.csv(precip.avg.df,
"../data/conflict/pak_avg_precip")##Read in the converted csv
pak.precip.avg = read.csv("../data/conflict/pak_avg_precip_converted.csv")
pak.precip.avg| year | avg_jan_precip | avg_jul_precip | avg_jan_precip_mm | avg_jul_precip_mm |
|---|---|---|---|---|
| 2010 | 0.000584 | 0.00432 | 0.584 | 4.32 |
| 2011 | 0.000396 | 0.00206 | 0.396 | 2.06 |
| 2012 | 0.000902 | 0.00108 | 0.902 | 1.08 |
| 2013 | 0.000355 | 0.00181 | 0.355 | 1.81 |
| 2014 | 0.000357 | 0.00156 | 0.357 | 1.56 |
| 2015 | 0.000712 | 0.00386 | 0.712 | 3.86 |
| 2016 | 0.000627 | 0.00222 | 0.627 | 2.22 |
| 2017 | 0.00191 | 0.00266 | 1.91 | 2.66 |
| 2018 | 0.000279 | 0.00194 | 0.279 | 1.94 |
| 2019 | 0.00127 | 0.00205 | 1.27 | 2.05 |
| 2020 | 0.00186 | 0.00161 | 1.86 | 1.61 |
| 2021 | 0.00039 | 0.00213 | 0.39 | 2.13 |
| 2022 | 0.00174 | 0.0063 | 1.74 | 6.3 |
avg.precip.plot =
ggplot(data = pak.precip.avg) +
geom_line(aes(x=year, y = avg_jan_precip_mm, color = "avg_jan_precip")) +
geom_line(aes(x=year, y = avg_jul_precip_mm, color = "avg_jul_precip")) +
xlab("Year") +
ylab("Total Precipitation (mm)") +
scale_color_manual(name = "Seasonal Precipitation",
values = c("avg_jan_precip" = "blue",
"avg_jul_precip" = "red"))
avg.precip.plot
Population Data
##Data acquired from datacommons.org (world bank)
##2022 is estimated based on the 10-year average growth rate
pak.pop = read.csv("../data/conflict/pak_pop.csv")
ggplot(data = pak.pop, aes(x = year, y = pak_pop)) +
geom_point(col = "darkorange2") +
geom_line(linetype = "dashed",
col = "darkorange2") +
labs(title = "Pakistan Population (2010-2022)") +
xlab("Year") +
ylab("Population")
Models Build
##Test: Merge the GDP and Number of Conflicts together
##Change the column header in order to merge
colnames(con.by.year)[1] = "year"
##Merge
pak.model.df = merge(con.by.year, pak.gdp, by = "year")
##It worked!
##Add water conflicts
##use by.x and by.y to merge using different variable names
pak.model.df = merge(pak.model.df, water.con.by.year, by.x ="year", by.y = as.character("years.num"))
##Now to add population
pak.model.df = merge(pak.model.df, pak.pop, by = "year")
##Convert year to number
pak.model.df$year = as.numeric(pak.model.df$year)
##Merge temperature data
pak.model.df = merge(pak.model.df, pak.temp.avg, by = "year")
##Merge precipitation data
pak.model.df = merge(pak.model.df, pak.precip.avg, by = "year")##Write to csv to save***
write.csv(pak.model.df,
"../data/conflict/pak_model_df")#pairs(pak.model.df)
corr.plot = ggcorrplot(cor(pak.model.df),
method = "square",
type = "lower",
lab = TRUE,
hc.order = TRUE,
colors = c("blue", "darksalmon", "firebrick"))
#p.mat = cor_pmat(pak.model.df)) #x out non-significant p-values
corr.plot
ggsave("./figures/corr_plot.jpeg", corr.plot)OLS LM
pak.nat.lm = lm(num.conflicts ~ year + gdp + pak_pop + avg_jan_celsius + avg_jul_celsius + avg_jan_precip_mm + avg_jul_precip_mm,
data = pak.model.df)
summary(pak.nat.lm)##
## Call:
## lm(formula = num.conflicts ~ year + gdp + pak_pop + avg_jan_celsius +
## avg_jul_celsius + avg_jan_precip_mm + avg_jul_precip_mm,
## data = pak.model.df)
##
## Residuals:
## 1 2 3 4 5 6 7 8
## -38.71 -429.17 351.84 1023.76 -1036.99 451.84 -933.68 -338.94
## 9 10 11 12 13
## 977.89 428.42 122.66 -559.13 -19.79
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.486e+06 9.936e+06 0.955 0.384
## year -4.869e+03 5.047e+03 -0.965 0.379
## gdp -7.691e-09 1.958e-08 -0.393 0.711
## pak_pop 1.214e-03 1.186e-03 1.024 0.353
## avg_jan_celsius -2.114e+02 5.503e+02 -0.384 0.717
## avg_jul_celsius 3.232e+03 1.667e+03 1.939 0.110
## avg_jan_precip_mm 1.838e+02 6.734e+02 0.273 0.796
## avg_jul_precip_mm 9.445e+02 9.808e+02 0.963 0.380
##
## Residual standard error: 1009 on 5 degrees of freedom
## Multiple R-squared: 0.827, Adjusted R-squared: 0.5847
## F-statistic: 3.413 on 7 and 5 DF, p-value: 0.09764##Scale the dataframe
#pak.model.df.scaled = as.data.frame(scale(pak.model.df, center = FALSE))
##Simon's Code to scale the data
library(dplyr)
pak.model.df.scaled = pak.model.df %>%
mutate(year = scale(year),
gdp = scale(gdp),
pak_pop = scale(pak_pop),
avg_jan_celsius = scale(avg_jan_celsius),
avg_jul_celsius = scale(avg_jul_celsius),
avg_jan_precip_mm = scale(avg_jan_precip_mm),
avg_jul_precip_mm = scale(avg_jul_precip_mm))GLM - Poisson Intercept Model
pak.nat.glm0 = glm(num.conflicts ~ 1,
family = poisson(link = 'log'),
data = pak.model.df.scaled)
summary(pak.nat.glm0)##
## Call:
## glm(formula = num.conflicts ~ 1, family = poisson(link = "log"),
## data = pak.model.df.scaled)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -30.535 -15.742 -4.689 13.344 34.658
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 8.685936 0.003605 2409 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 4909.8 on 12 degrees of freedom
## Residual deviance: 4909.8 on 12 degrees of freedom
## AIC: 5048.2
##
## Number of Fisher Scoring iterations: 4##Convert the coefficients
exp(coef(pak.nat.glm0))## (Intercept)
## 5919.077##Use jtools to make a more aesthetically pleasing table
##Use exp=TRUE syntax to report back on the transformed coefficients
pak.nat.glm0.sum = summ(pak.nat.glm0,
exp = TRUE)
pak.nat.glm0.sum| Observations | 13 |
| Dependent variable | num.conflicts |
| Type | Generalized linear model |
| Family | poisson |
| Link | log |
| 𝛘²(0) | -0.00 |
| Pseudo-R² (Cragg-Uhler) | 0.00 |
| Pseudo-R² (McFadden) | 0.00 |
| AIC | 5048.18 |
| BIC | 5048.74 |
| exp(Est.) | 2.5% | 97.5% | z val. | p | |
|---|---|---|---|---|---|
| (Intercept) | 5919.08 | 5877.40 | 5961.05 | 2409.44 | 0.00 |
| Standard errors: MLE |
##Visualize Residuals
plot(pak.nat.glm0)
## hat values (leverages) are all = 0.07692308
## and there are no factor predictors; no plot no. 5
##Pseudo R2
pR2(pak.nat.glm0)## fitting null model for pseudo-r2## llh llhNull G2 McFadden r2ML r2CU
## -2523.089 -2523.089 0.000 0.000 0.000 0.000GLM - Poisson - Full Model
##Starting with a generalized linear model - Poisson - as the number of conflicts is count data
## Try running it without year
## Running the Z-score for the variables (the interpretation would be the change in one standard deviation) - use function scale()
pak.nat.glm = glm(num.conflicts ~ year + gdp + pak_pop + avg_jan_celsius + avg_jul_celsius + avg_jan_precip_mm + avg_jul_precip_mm,
family = poisson(link = 'log'),
data = pak.model.df.scaled)
summary(pak.nat.glm)##
## Call:
## glm(formula = num.conflicts ~ year + gdp + pak_pop + avg_jan_celsius +
## avg_jul_celsius + avg_jan_precip_mm + avg_jul_precip_mm,
## family = poisson(link = "log"), data = pak.model.df.scaled)
##
## Deviance Residuals:
## 1 2 3 4 5 6 7 8
## -0.5569 -9.4196 6.0944 14.9581 -12.9770 4.7249 -13.0862 -5.8290
## 9 10 11 12 13
## 13.3055 6.3542 0.1537 -7.4487 1.2489
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 8.661796 0.003688 2348.690 < 2e-16 ***
## year -3.170094 0.280344 -11.308 < 2e-16 ***
## gdp -0.045146 0.014938 -3.022 0.00251 **
## pak_pop 3.278815 0.277629 11.810 < 2e-16 ***
## avg_jan_celsius -0.033123 0.007115 -4.655 3.23e-06 ***
## avg_jul_celsius 0.378056 0.015330 24.661 < 2e-16 ***
## avg_jan_precip_mm 0.022933 0.004929 4.653 3.27e-06 ***
## avg_jul_precip_mm 0.237225 0.018737 12.661 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for poisson family taken to be 1)
##
## Null deviance: 4909.8 on 12 degrees of freedom
## Residual deviance: 1020.4 on 5 degrees of freedom
## AIC: 1172.7
##
## Number of Fisher Scoring iterations: 4##Need to transform the coefficients ##Summ Visualization of the Model Output
summ(pak.nat.glm,
exp = TRUE)| Observations | 13 |
| Dependent variable | num.conflicts |
| Type | Generalized linear model |
| Family | poisson |
| Link | log |
| 𝛘²(7) | 3889.43 |
| Pseudo-R² (Cragg-Uhler) | 1.00 |
| Pseudo-R² (McFadden) | 0.77 |
| AIC | 1172.74 |
| BIC | 1177.26 |
| exp(Est.) | 2.5% | 97.5% | z val. | p | |
|---|---|---|---|---|---|
| (Intercept) | 5777.90 | 5736.29 | 5819.82 | 2348.69 | 0.00 |
| year | 0.04 | 0.02 | 0.07 | -11.31 | 0.00 |
| gdp | 0.96 | 0.93 | 0.98 | -3.02 | 0.00 |
| pak_pop | 26.54 | 15.40 | 45.74 | 11.81 | 0.00 |
| avg_jan_celsius | 0.97 | 0.95 | 0.98 | -4.66 | 0.00 |
| avg_jul_celsius | 1.46 | 1.42 | 1.50 | 24.66 | 0.00 |
| avg_jan_precip_mm | 1.02 | 1.01 | 1.03 | 4.65 | 0.00 |
| avg_jul_precip_mm | 1.27 | 1.22 | 1.32 | 12.66 | 0.00 |
| Standard errors: MLE |
##Coefficient transformation as the output is in the log values (link = log)
##Reminder Interpretation Note:
##Log odds scale to odds scale - it becomes multiplicative - a rate of change (bigger than 1 it is a positive change, if it is less than 1 it is a negative change)
##Scale the gdp coefficient and population
pak.nat.glm.coef = exp(coef(pak.nat.glm))
pak.nat.glm.coef## (Intercept) year gdp pak_pop
## 5.777904e+03 4.199965e-02 9.558577e-01 2.654430e+01
## avg_jan_celsius avg_jul_celsius avg_jan_precip_mm avg_jul_precip_mm
## 9.674196e-01 1.459444e+00 1.023198e+00 1.267726e+00##Calculate McFadden's Pseudo R^2
##References
##“Conditional logit analysis of qualitative choice behavior.” Pp. 105-142 in P. Zarembka (ed.), Frontiers in Econometrics. Academic Press. http://eml.berkeley.edu/~mcfadden/travel.html
##Bahvioural Travel Modelling. Edited by David Hensher and Peter Stopher. 1979. McFadden contributed Ch. 15 "Quantitative Methods for Analyzing Travel Behaviour on Individuals: Some Recent Developments"
pR2(pak.nat.glm)## fitting null model for pseudo-r2## llh llhNull G2 McFadden r2ML
## -578.3723540 -2523.0893672 3889.4340265 0.7707682 1.0000000
## r2CU
## 1.0000000##Visualization to improve the model output
export_summs(pak.nat.glm0, pak.nat.glm,
model.names = c("Intercept Model", "Full Model"),
coefs = c("Intercept" = "(Intercept)",
"Year" = "year", "GDP" = "gdp",
"Population" = "pak_pop" ,
"January Temperature" = "avg_jan_celsius" ,
"July Temperature" = "avg_jul_celsius",
"January Precipitation" = "avg_jan_precip_mm",
"July Precipitation" = "avg_jul_precip_mm"),
exp = TRUE)| Intercept Model | Full Model | |
|---|---|---|
| Intercept | 5919.08 *** | 5777.90 *** |
| (0.00) | (0.00) | |
| Year | 0.04 *** | |
| (0.28) | ||
| GDP | 0.96 ** | |
| (0.01) | ||
| Population | 26.54 *** | |
| (0.28) | ||
| January Temperature | 0.97 *** | |
| (0.01) | ||
| July Temperature | 1.46 *** | |
| (0.02) | ||
| January Precipitation | 1.02 *** | |
| (0.00) | ||
| July Precipitation | 1.27 *** | |
| (0.02) | ||
| N | 13 | 13 |
| AIC | 5048.18 | 1172.74 |
| BIC | 5048.74 | 1177.26 |
| Pseudo R2 | 0.00 | 1.00 |
| *** p < 0.001; ** p < 0.01; * p < 0.05. | ||
Residuals
##View the residuals
plot(residuals.glm(pak.nat.glm))
abline(h = 0, col = "red")
##Fitting a spline may be better for this model?plot(pak.nat.glm)
##Visualize
hist(pak.nat.glm$residuals)
pak.nat.glm.scaled2 = glm(water.num.conflicts ~ year + gdp_percent_growth +
avg_jan_celsius + avg_jul_celsius +
avg_jan_precip_mm + avg_jul_precip_mm +
offset(log(num.conflicts)),
family = poisson(link = 'log'),
data = pak.model.df.scaled)
pak.nat.glm.scaled2##
## Call: glm(formula = water.num.conflicts ~ year + gdp_percent_growth +
## avg_jan_celsius + avg_jul_celsius + avg_jan_precip_mm + avg_jul_precip_mm +
## offset(log(num.conflicts)), family = poisson(link = "log"),
## data = pak.model.df.scaled)
##
## Coefficients:
## (Intercept) year gdp_percent_growth avg_jan_celsius
## -2.85207 -0.05847 2.08433 0.02040
## avg_jul_celsius avg_jan_precip_mm avg_jul_precip_mm
## -0.08279 -0.01050 0.11766
##
## Degrees of Freedom: 12 Total (i.e. Null); 6 Residual
## Null Deviance: 390.5
## Residual Deviance: 167.1 AIC: 281pak.nat.glm.scaled3 = glm(num.conflicts ~ year + gdp_percent_growth +
avg_jan_celsius + avg_jul_celsius +
avg_jan_precip_mm + avg_jul_precip_mm +
offset(log(pak_pop)),
family = poisson(link = 'log'),
data = pak.model.df.scaled)## Warning in log(pak_pop): NaNs producedpak.nat.glm.scaled3##
## Call: glm(formula = num.conflicts ~ year + gdp_percent_growth + avg_jan_celsius +
## avg_jul_celsius + avg_jan_precip_mm + avg_jul_precip_mm +
## offset(log(pak_pop)), family = poisson(link = "log"), data = pak.model.df.scaled)
##
## Coefficients:
## (Intercept) year gdp_percent_growth avg_jan_celsius
## 10.13634 -1.27327 1.83490 -0.10573
## avg_jul_celsius avg_jan_precip_mm avg_jul_precip_mm
## 0.05949 0.00939 NA
##
## Degrees of Freedom: 5 Total (i.e. Null); 0 Residual
## (7 observations deleted due to missingness)
## Null Deviance: 10570
## Residual Deviance: -4.747e-12 AIC: 75.9##NEED TO FIX THE SCALING ISSUE TO RUN!
##Build a new data frame
pak.new.data = data.frame(year = 2030, pak_pop = 300000000, avg_jul_celsius = 30, avg_jul_precip_mm = 7, gdp = 4.5e+11, avg_jan_celsius = 10, avg_jan_precip_mm = 0.002)
##Scale the new data
#pak.new.data.scaled = as.data.frame(scale(pak.new.data))
predict.glm(pak.nat.glm,
newdata = pak.new.data,
type = "response",
se.fit = TRUE)