Load packages
library(tidyr)
library(dplyr)
library(ggplot2)
library(sf)
library(skimr)
library(ggrepel)
Attaching package: ‘dplyr’
The following objects are masked from ‘package:stats’:
filter, lag
The following objects are masked from ‘package:base’:
intersect, setdiff, setequal, union
Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
Load data
sh <- st_read("vnm-final.gpkg")
df <- read.csv("vnm-long.csv")
Reading layer `vnm-append' from data source `/work/vnm-final.gpkg' using driver `GPKG'
Simple feature collection with 63 features and 55 fields
geometry type: MULTIPOLYGON
dimension: XY
bbox: xmin: 102.1421 ymin: 6.953306 xmax: 116.9473 ymax: 23.3939
geographic CRS: WGS 84
Plot Median Local Moran's I (2018)
palette <- c('#f2f2f2', '#ca0020', '#0571b0', '#92c5de', '#f4a582')
labelss <- c("Not Significant", "High-High (4)", "Low-Low (8)", "Low-High (2)", "High-Low (0)")
#map
ggplot() +
geom_sf(data = sh, aes(fill = factor(SAT_CL)), lwd=0.05) +
labs(title = "Average night-time lights (Radiance), 2018") +
scale_fill_manual(name ="Median LISAs clusters", values = palette, labels = labelss) +
theme_set(theme_bw()) +
theme(plot.title=element_text(hjust=0.5, lineheight=1.2),
axis.title.x = element_blank(),
axis.title.y = element_blank())
#save
ggsave("LISAs SAT.png")
Saving 6.67 x 6.67 in image
palette <- c('#f2f2f2', '#ca0020', '#0571b0', '#92c5de', '#f4a582')
labelss <- c("Not Significant", "High-High (3)", "Low-Low (5)", "Low-High (1)", "High-Low (0)")
#map
ggplot() +
geom_sf(data = sh, aes(fill = factor(GRDP_CL)), lwd=0.05) +
labs(title = "Gross Regional Domestic Product (GRDP), 2018") +
scale_fill_manual(name ="Median LISAs clusters", values = palette, labels = labelss) +
theme_set(theme_bw()) +
theme( plot.title=element_text(hjust=0.5, lineheight=1.2),
axis.title.x = element_blank(),
axis.title.y = element_blank())
#save
ggsave("LISAs GRDP.png", width = 6.67 , height = 6.67 )
Plot relationship between GDP and Light
Proportion of variation in GRDP that is predictable by night-time lights
ggplot(data= df, aes(x = sat, y = grdp, col = factor(year))) +
geom_point() +
geom_smooth(method=lm, se=FALSE) +
geom_label(aes(label = province), alpha=0.5) +
scale_color_manual(name ="Year (R²)",
values = c("#edf8fb","#ccece6","#99d8c9","#66c2a4","#41ae76","#238b45","#005824"),
labels = c("2012 (0.08)","2013 (0.16)", "2014 (0.23)", "2015 (0.15)", "2016 (0.18)", "2017 (0.17)", "2018 (0.10)")) +
labs(title = "Proportion of variation in GRDP predictable by lights", x ="Average night-time lights intensity (Radiance)", y = "Gross Regional Domestic Product (GRDP), VND") +
theme_set(theme_bw()) +
theme_update(plot.title = element_text(hjust = 0.5))
ggsave("rel.png")
`geom_smooth()` using formula 'y ~ x'
Saving 6.67 x 6.67 in image
`geom_smooth()` using formula 'y ~ x'
Compute explanatory power (R2) across time
lm2012 <- lm(sh$grdp_2012 ~ sh$sat_2012)
summary(lm2012)
lm2013 <- lm(sh$grdp_2013 ~ sh$sat_2013)
summary(lm2013)
lm2014 <- lm(sh$grdp_2014 ~ sh$sat_2014)
summary(lm2014)
lm2015 <- lm(sh$grdp_2015 ~ sh$sat_2015)
summary(lm2015)
lm2016 <- lm(sh$grdp_2016 ~ sh$sat_2016)
summary(lm2016)
lm2017 <- lm(sh$grdp_2017 ~ sh$sat_2017)
summary(lm2017)
lm2018 <- lm(sh$grdp_2018 ~ sh$sat_2018)
summary(lm2018)
Plot explanatory power (R2) across space
note: value is median of 7 years
GRDP
#map
ggplot() +
geom_sf(data = sh, aes(fill = grdp_r2), lwd=0.05) +
labs(title = "Gross regional domestic product (GRDP)") +
theme_set(theme_bw()) +
scale_fill_gradient(name = "R²", low = "white", high = "darkgreen") +
theme(plot.title=element_text(hjust=0.5, lineheight=1.2),
axis.title.x = element_blank(),
axis.title.y = element_blank())
#save
ggsave("grdp_r2.png")
Saving 6.67 x 6.67 in image
Industrial
#map
ggplot() +
geom_sf(data = sh, aes(fill = indus_r2), lwd=0.05) +
labs(title = "Industry and construction") +
theme_set(theme_bw()) +
scale_fill_gradient(name = "R²", low = "white", high = "darkgreen") +
theme(plot.title=element_text(hjust=0.5, lineheight=1.2),
axis.title.x = element_blank(),
axis.title.y = element_blank())
#save
ggsave("indus_r2.png")
Saving 6.67 x 6.67 in image
Service
#map
ggplot() +
geom_sf(data = sh, aes(fill = service_r2), lwd=0.05) +
labs(title = "Services") +
scale_fill_gradient(name = "R²", low = "white", high = "darkgreen") +
theme_set(theme_bw()) +
theme(plot.title=element_text(hjust=0.5, lineheight=1.2),
axis.title.x = element_blank(),
axis.title.y = element_blank())
#save
ggsave("service_r2.png")
Saving 6.67 x 6.67 in image
Agriculture
#map
ggplot() +
geom_sf(data = sh, aes(fill = agro_r2), lwd=0.05) +
labs(title = "Agriculture, forestry and fishing") +
scale_fill_gradient(name = "R²", low = "white", high = "darkgreen") +
theme_set(theme_bw()) +
theme(plot.title=element_text(hjust=0.5, lineheight=1.2),
axis.title.x = element_blank(),
axis.title.y = element_blank())
#save
ggsave("agro_r2.png")
Saving 6.67 x 6.67 in image
Plot p-values
In null hypothesis significance testing, the p-value is the probability of obtaining test results at least as extreme as the results actually observed, under the assumption that the null hypothesis is correct.
GRDP
#map
ggplot() +
geom_sf(data = sh, aes(fill = grdp_p), lwd=0.05) +
labs(title = "Gross regional domestic product (GRDP)") +
theme_set(theme_bw()) +
scale_fill_gradient(name = "P-value") +
theme(plot.title=element_text(hjust=0.5, lineheight=1.2),
axis.title.x = element_blank(),
axis.title.y = element_blank())
#save
ggsave("grdp_p.png")
Saving 6.67 x 6.67 in image
#map
ggplot() +
geom_sf(data = sh, aes(fill = indus_p), lwd=0.05) +
labs(title = "Industry and construction") +
scale_fill_gradient(name = "P-value") +
theme_set(theme_bw()) +
theme(plot.title=element_text(hjust=0.5, lineheight=1.2),
axis.title.x = element_blank(),
axis.title.y = element_blank())
#save
ggsave("indus_p.png")
Saving 6.67 x 6.67 in image
#map
ggplot() +
geom_sf(data = sh, aes(fill = service_p), lwd=0.05) +
labs(title = "Services") +
scale_fill_gradient(name = "P-value") +
theme_set(theme_bw()) +
theme(plot.title=element_text(hjust=0.5, lineheight=1.2),
axis.title.x = element_blank(),
axis.title.y = element_blank())
#save
ggsave("service_p.png")
Saving 6.67 x 6.67 in image
#map
ggplot() +
geom_sf(data = sh, aes(fill = agro_p), lwd=0.05) +
labs(title = "Agriculture, forestry and fishing") +
scale_fill_gradient(name = "P-value") +
theme_set(theme_bw()) +
theme(plot.title=element_text(hjust=0.5, lineheight=1.2),
axis.title.x = element_blank(),
axis.title.y = element_blank())
#save
ggsave("agro_p.png")
Saving 6.67 x 6.67 in image