A guide to perform Bayesian Model Averaging and Geographically Weighted Regression using BMS and MGWR application
Environment: quarcslab/deepnote-r-spatial
Load libraries
#devtools::install_version("BMS", version = "0.3.4", repos = "http://cran.us.r-project.org")
library(sf)
library(tidyverse)
library(ggplot2)
library(RColorBrewer)
library(ggrepel)
library(skimr)
library(imager)
library(spdep)
library(ggrepel)
library(dplyr)
library(BMS)
Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
── Attaching packages ─────────────────────────────────────── tidyverse 1.3.0 ──
✔ ggplot2 3.3.3 ✔ purrr 0.3.4
✔ tibble 3.1.0 ✔ dplyr 1.0.4
✔ tidyr 1.1.2 ✔ stringr 1.4.0
✔ readr 1.4.0 ✔ forcats 0.5.1
── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag() masks stats::lag()
Loading required package: magrittr
Attaching package: ‘magrittr’
The following object is masked from ‘package:purrr’:
set_names
The following object is masked from ‘package:tidyr’:
extract
Attaching package: ‘imager’
The following object is masked from ‘package:magrittr’:
add
The following object is masked from ‘package:stringr’:
boundary
The following object is masked from ‘package:tidyr’:
fill
The following objects are masked from ‘package:stats’:
convolve, spectrum
The following object is masked from ‘package:graphics’:
frame
The following object is masked from ‘package:base’:
save.image
Loading required package: sp
Attaching package: ‘sp’
The following object is masked from ‘package:imager’:
bbox
Loading required package: spData
To access larger datasets in this package, install the spDataLarge
package with: `install.packages('spDataLarge',
repos='https://nowosad.github.io/drat/', type='source')`
Load shapefiles and datasets
phl_adm1 <- st_read("phl_adm1.gpkg") #for correlation between ecnomic output vs energy consumption
phl_adm2 <- st_read("phl_adm2.gpkg") #for borders
phl_adm3 <- st_read("phl_adm3.gpkg") #for energy consumption vs environmental and climatic determinants
phl_bma = read_csv("phl_bma.csv") #for running bma
phl_adm3_df <- st_drop_geometry(phl_adm3)
skim(phl_adm3_df)
Reading layer `PHL_ADM1' from data source `/work/phl_adm1.gpkg' using driver `GPKG'
Simple feature collection with 17 features and 6 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 114.2779 ymin: 4.587294 xmax: 126.605 ymax: 21.12189
Geodetic CRS: WGS 84
Reading layer `PHL_adm1' from data source `/work/phl_adm2.gpkg' using driver `GPKG'
Simple feature collection with 81 features and 12 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 116.9283 ymin: 4.58694 xmax: 126.6053 ymax: 21.07014
Geodetic CRS: WGS 84
Reading layer `PHL_ADM3' from data source `/work/phl_adm3.gpkg' using driver `GPKG'
Simple feature collection with 1647 features and 19 fields
Geometry type: MULTIPOLYGON
Dimension: XY
Bounding box: xmin: 114.2779 ymin: 4.587294 xmax: 126.605 ymax: 21.12189
Geodetic CRS: WGS 84
── Column specification ────────────────────────────────────────────────────────
cols(
`Night-time lights` = col_double(),
`Physical elevation` = col_double(),
`CO2 concentration` = col_double(),
`Land surface temperature` = col_double(),
`Vegetation index` = col_double(),
`Population density` = col_double(),
`Wind speed potential` = col_double(),
`Wind power density potential` = col_double(),
`Photovoltaic power potential` = col_double(),
`Distance to coast` = col_double(),
`Distance to water` = col_double(),
`Air temperature` = col_double(),
Precipitation = col_double(),
`Ozone concentration` = col_double(),
`Particulate matter (PM2.5)` = col_double()
)
── Data Summary ────────────────────────
Values
Name phl_adm3_df
Number of rows 1647
Number of columns 19
_______________________
Column type frequency:
character 1
numeric 18
________________________
Group variables None
── Variable type: character ────────────────────────────────────────────────────
skim_variable n_missing complete_rate min max empty n_unique whitespace
1 region 0 1 3 27 0 1424 0
── Variable type: numeric ──────────────────────────────────────────────────────
skim_variable n_missing complete_rate mean sd
1 gqid 0 1 823 476.
2 Night.time.lights 0 1 1.05 4.73
3 Physical.elevation 0 1 234. 270.
4 CO2.concentration 0 1 413. 0.986
5 Land.surface.temperature 0 1 301. 7.78
6 Vegetation.index 0 1 4918. 1739.
7 Population.density 0 1 1133. 4276.
8 Wind.speed.potential 0 1 4.16 1.10
9 Wind.power.density.potential 0 1 128. 103.
10 Photovoltaic.power.potential 0 1 3.80 0.267
11 Distance.to.coast 0 1 18447. 19141.
12 Distance.to.water 0 1 14197. 13750.
13 Air.temperature 0 1 26.4 3.65
14 Precipitation 0 1 194. 57.0
15 Ozone.concentration 0 1 37.6 9.08
16 Particulate.matter..PM2.5. 0 1 7.29 2.82
17 X 0 1 123. 1.82
18 Y 0 1 12.0 3.52
p0 p25 p50 p75 p100 hist
1 0 412. 823 1234. 1646 ▇▇▇▇▇
2 0 0.021 0.072 0.300 52.9 ▇▁▁▁▁
3 2.21 59.4 142. 298. 1902. ▇▁▁▁▁
4 410. 412. 413. 414. 416. ▁▃▇▃▁
5 0 299. 301. 302. 310. ▁▁▁▁▇
6 0 3881. 5072. 6194. 9414. ▁▃▇▇▁
7 0 168. 303. 589. 66301. ▇▁▁▁▁
8 1.70 3.32 4.10 4.89 7.82 ▃▇▇▂▁
9 8.87 59.5 98.9 160. 780. ▇▂▁▁▁
10 0 3.64 3.78 3.95 4.51 ▁▁▁▂▇
11 741. 5328. 10125. 24701. 96374. ▇▂▁▁▁
12 458. 4972. 8921. 17912. 70723. ▇▂▁▁▁
13 0 26.3 27.3 27.6 28.6 ▁▁▁▁▇
14 0 157. 189. 228. 408. ▁▅▇▂▁
15 0 33.4 36.7 45.2 54.2 ▁▁▃▇▆
16 0 5.46 7.19 9.11 30.2 ▅▇▁▁▁
17 114. 121. 123. 124. 126. ▁▁▇▇▇
18 4.75 9.00 11.7 14.8 20.8 ▅▇▇▇▂
1. Correlate GRP with Lights using OLS
lm <- lm(phl_adm1$sat2012 ~ phl_adm1$gdp2012)
summary(lm)
theme_update(plot.title = element_text(hjust = 0.5))
ggplot(phl_adm1, aes(x=gdp2012, y=sat2012)) +
geom_point() +
geom_smooth(method=lm, se=FALSE) +
geom_label(aes(label = region)) +
labs(title = "Year 2012 (R2: 0.79)", x ="Night-time lights intensity (radiance)", y = "GRP in current prices (PHP)") +
theme_set(theme_bw())
ggsave("sat-grp12.png")
`geom_smooth()` using formula 'y ~ x'
Warning message:
“Computation failed in `stat_smooth()`:
'what' must be a function or character string”
Saving 6.67 x 6.67 in image
`geom_smooth()` using formula 'y ~ x'
Warning message:
“Computation failed in `stat_smooth()`:
'what' must be a function or character string”
lm20 <- lm(phl_adm1$sat2020 ~ phl_adm1$gdp2020)
summary(lm20)
theme_update(plot.title = element_text(hjust = 0.5))
ggplot(phl_adm1, aes(x=gdp2020, y=sat2020)) +
geom_point() +
geom_smooth(method=lm, se=FALSE) +
geom_label(aes(label = region)) +
labs(title = "Year 2020 (R2: 0.83)", x ="Night-time lights intensity (radiance)", y = "GRP in current prices (PHP)") +
theme_set(theme_bw())
ggsave("sat-grp20.png")
`geom_smooth()` using formula 'y ~ x'
Warning message:
“Computation failed in `stat_smooth()`:
'what' must be a function or character string”
Saving 6.67 x 6.67 in image
`geom_smooth()` using formula 'y ~ x'
Warning message:
“Computation failed in `stat_smooth()`:
'what' must be a function or character string”
2. Explore economic activities hotspots using ESDA using GeoDa
3. Feature selection using BMA
phl_bms = bms(phl_bma, burn =100000, iter=200000, g="BRIC", mprior="uniform", mcmc="bd", nmodel=2000)
Warning message in bms(phl_bma, burn = 1e+05, iter = 2e+05, g = "BRIC", mprior = "uniform", :
“Argument 'X.data' contains NAs. The corresponding rows have not been taken into account.”
PIP Post Mean Post SD Cond.Pos.Sign
Land surface temperature 1.000000 3.919104e-01 4.726009e-02 1.00000000
Population density 1.000000 8.922475e-04 1.672902e-05 1.00000000
Physical elevation 0.996900 1.357517e-03 3.365230e-04 1.00000000
Photovoltaic power potential 0.957110 -1.141649e+00 4.226841e-01 0.00000000
Ozone concentration 0.873615 -2.271840e-02 1.104061e-02 0.00000000
Distance to coast 0.807960 -9.651360e-06 5.701190e-06 0.00000000
Vegetation index 0.548310 -5.785023e-05 5.945347e-05 0.00000000
Precipitation 0.389450 1.096783e-03 1.538110e-03 1.00000000
Wind speed potential 0.276495 3.760712e-02 6.905772e-02 1.00000000
Particulate matter (PM2.5) 0.148095 -9.713770e-03 2.630807e-02 0.00023633
Distance to water 0.096015 -1.043330e-06 3.839172e-06 0.00015623
CO2 concentration 0.089395 8.406063e-03 3.200541e-02 1.00000000
Wind power density potential 0.058460 3.127410e-05 2.744894e-04 0.83484434
Air temperature 0.025680 1.359287e-05 3.087728e-03 0.57651869
Idx
Land surface temperature 3
Population density 5
Physical elevation 1
Photovoltaic power potential 8
Ozone concentration 13
Distance to coast 9
Vegetation index 4
Precipitation 12
Wind speed potential 6
Particulate matter (PM2.5) 14
Distance to water 10
CO2 concentration 2
Wind power density potential 7
Air temperature 11
Mean no. regressors Draws Burnins Time
"7.2675" "2e+05" "1e+05" "16.32012 secs"
No. models visited Modelspace 2^K % visited % Topmodels
"64354" "16384" "393" "100"
Corr PMP No. Obs. Model Prior g-Prior
"0.9994" "1646" "uniform / 7" "BRIC"
Shrinkage-Stats
"Av=0.9994"
Time difference of 16.32012 secs
image(phl_bms)
summary(phl_bms)
4. Run GWR for the highest probability covariates/variables
im<-load.image("gwr_spec.png")
plot(im)
#Load MGWR session results
phl_gwr <- read_csv("phl_gwr.csv")
#select relevant variables
var <- c("gqid", "x_coor", "y_coor", "beta_land", "beta_pop", "beta_ele", "beta_solar", "beta_ozone", "beta_coast", "p_land", "p_pop", "p_ele", "p_solar", "p_ozone", "p_coast")
phl_gwr <- phl_gwr%>% select(var)
##coefficient is negative, meaning negative correlation
phl_gwr <- phl_gwr %>% mutate(solar_label = ifelse(beta_solar < 0, "Negative", "Positive"))
phl_gwr <- phl_gwr %>% mutate(elevation_label = ifelse(beta_ele < 0, "Negative", "Positive"))
phl_gwr <- phl_gwr %>% mutate(ozone_label = ifelse(beta_ozone < 0, "Negative", "Positive"))
phl_gwr <- phl_gwr %>% mutate(landt_label = ifelse(beta_land < 0, "Negative", "Positive"))
phl_gwr <- phl_gwr %>% mutate(coast_label = ifelse(beta_coast < 0, "Negative", "Positive"))
phl_gwr <- phl_gwr %>% mutate(pop_label = ifelse(beta_pop < 0, "Negative", "Positive"))
##0.01 significance threshold
phl_gwr <- phl_gwr %>% mutate(solar_p_label = ifelse(p_solar <= 0.01, "Significant", "Not Significant"))
phl_gwr <- phl_gwr %>% mutate(elevation_p_label = ifelse(p_ele <= 0.01, "Significant", "Not Significant"))
phl_gwr <- phl_gwr %>% mutate(ozone_p_label = ifelse(p_ozone <= 0.01, "Significant", "Not Significant"))
phl_gwr <- phl_gwr %>% mutate(landt_p_label = ifelse(p_land <= 0.01, "Significant", "Not Significant"))
phl_gwr <- phl_gwr %>% mutate(coast_p_label = ifelse(p_coast <= 0.01, "Significant", "Not Significant"))
phl_gwr <- phl_gwr %>% mutate(pop_p_label = ifelse(p_pop <= 0.01, "Significant", "Not Significant"))
##significant positive/negative regions
phl_gwr <- phl_gwr %>% mutate(sig_solar = ifelse(solar_p_label == "Significant" & solar_label == "Positive", "Positive",
ifelse(solar_p_label == "Significant" & solar_label == "Negative", "Negative", "Not Significant")))
phl_gwr <- phl_gwr %>% mutate(sig_elevation = ifelse(elevation_p_label == "Significant" & elevation_label == "Positive", "Positive",
ifelse(elevation_p_label == "Significant" & elevation_label == "Negative", "Negative", "Not Significant")))
phl_gwr <- phl_gwr %>% mutate(sig_ozone= ifelse(ozone_p_label == "Significant" & ozone_label == "Positive", "Positive",
ifelse(ozone_p_label == "Significant" & ozone_label == "Negative", "Negative", "Not Significant")))
phl_gwr <- phl_gwr %>% mutate(sig_landt = ifelse(landt_p_label == "Significant" & landt_label == "Positive", "Positive",
ifelse(landt_p_label == "Significant" & landt_label == "Negative", "Negative", "Not Significant")))
phl_gwr <- phl_gwr %>% mutate(sig_coast = ifelse(coast_p_label == "Significant" & coast_label == "Positive", "Positive",
ifelse(coast_p_label == "Significant" & coast_label == "Negative", "Negative", "Not Significant")))
phl_gwr <- phl_gwr %>% mutate(sig_pop = ifelse(pop_p_label == "Significant" & pop_label == "Positive", "Positive",
ifelse(pop_p_label == "Significant" & pop_label == "Negative", "Negative", "Not Significant")))
#Merge, save csv and display
phl_gwr <- merge(phl_adm3, phl_gwr <- phl_gwr , by.x=c("gqid"), by.y=c("gqid"), all.x = TRUE)
phl_gwr
── Column specification ────────────────────────────────────────────────────────
cols(
.default = col_double()
)
ℹ Use `spec()` for the full column specifications.
ERROR while rich displaying an object: Error in loadNamespace(name): there is no package called ‘geojsonio’
Traceback:
1. FUN(X[[i]], ...)
2. tryCatch(withCallingHandlers({
. if (!mime %in% names(repr::mime2repr))
. stop("No repr_* for mimetype ", mime, " in repr::mime2repr")
. rpr <- repr::mime2repr[[mime]](obj)
. if (is.null(rpr))
. return(NULL)
. prepare_content(is.raw(rpr), rpr)
. }, error = error_handler), error = outer_handler)
3. tryCatchList(expr, classes, parentenv, handlers)
4. tryCatchOne(expr, names, parentenv, handlers[[1L]])
5. doTryCatch(return(expr), name, parentenv, handler)
6. withCallingHandlers({
. if (!mime %in% names(repr::mime2repr))
. stop("No repr_* for mimetype ", mime, " in repr::mime2repr")
. rpr <- repr::mime2repr[[mime]](obj)
. if (is.null(rpr))
. return(NULL)
. prepare_content(is.raw(rpr), rpr)
. }, error = error_handler)
7. repr::mime2repr[[mime]](obj)
8. repr_geojson.sf(obj)
9. repr_geojson(geojsonio::geojson_list(obj), ...)
10. geojsonio::geojson_list
11. getExportedValue(pkg, name)
12. asNamespace(ns)
13. getNamespace(ns)
14. loadNamespace(name)
15. withRestarts(stop(cond), retry_loadNamespace = function() NULL)
16. withOneRestart(expr, restarts[[1L]])
17. doWithOneRestart(return(expr), restart)
ggplot() +
geom_sf(data = phl_gwr, aes(fill = sig_elevation), lwd = 0) +
geom_sf(data = phl_adm1, aes(fill = NA, alpha=0), show.legend = FALSE) +
scale_fill_manual(name = "Correlation",
values = c('#CE1126', '#FFFFFF', '#0038A8')) +
theme_set(theme_bw()) +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank()) +
labs(title="Physical Elevation (metres)") +
theme(plot.title=element_text(hjust=0.5, lineheight=1.2))
ggsave("elevation.png")
Saving 6.67 x 6.67 in image
ggplot() +
geom_sf(data = phl_gwr, aes(fill = sig_ozone), lwd = 0) +
geom_sf(data = phl_adm1, aes(fill = NA, alpha=0), show.legend = FALSE) +
scale_fill_manual(name = "Correlation",
values = c('#CE1126', '#FFFFFF', '#0038A8')) +
theme_set(theme_bw()) +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank()) +
labs(title="Ozone concentration (ug/m3)") +
theme(plot.title=element_text(hjust=0.5, lineheight=1.2))
ggsave("ozone.png")
Saving 6.67 x 6.67 in image
ggplot() +
geom_sf(data = phl_gwr, aes(fill = sig_solar), lwd = 0) +
geom_sf(data = phl_adm1, aes(fill = NA, alpha=0), show.legend = FALSE) +
scale_fill_manual(name = "Correlation",
values = c('#CE1126', '#FFFFFF', '#0038A8')) +
theme_set(theme_bw()) +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank()) +
labs(title="Photovoltaic power potential (kWh/kWp)") +
theme(plot.title=element_text(hjust=0.5, lineheight=1.2))
ggsave("solar.png")
Saving 6.67 x 6.67 in image
ggplot() +
geom_sf(data = phl_gwr, aes(fill = sig_pop), lwd = 0) +
geom_sf(data = phl_adm1, aes(fill = NA, alpha=0), show.legend = FALSE) +
scale_fill_manual(name = "Correlation",
values = c('#FFFFFF', '#0038A8')) +
theme_set(theme_bw()) +
theme(axis.title.x = element_blank(),
axis.title.y = element_blank()) +
labs(title="Population density (person/km2)") +
theme(plot.title=element_text(hjust=0.5, lineheight=1.2))
ggsave("pop.png")
Saving 6.67 x 6.67 in image
================================================================================
MGWR Version: 2.2.1
Released on: 03/20/2020
Source code is available at: https://github.com/pysal/mgwr
Development Team: Ziqi Li, Taylor Oshan, Stewart Fotheringham, Wei Kang,
Levi Wolf, Hanchen Yu, Mehak Sachdeva, and Sarah Bardin
Spatial Analysis Research Center (SPARC)
Arizona State University, Tempe, USA
================================================================================
Model type: Gaussian
Number of observations: 1647
Number of covariates: 7
Dependent variable: Night-time lights
Variable standardization: On
Total runtime: 1:27:39
Global Regression Results
--------------------------------------------------------------------------------
Residual sum of squares: 346.844
Log-likelihood: -1054.113
AIC: 2122.226
AICc: 2124.314
R2: 0.789
Adj. R2: 0.789
Variable Est. SE t(Est/SE) p-value
------------------------------------ ---------- ---------- ---------- ----------
Intercept 0.000 0.011 0.000 1.000
Land surface temperature 0.025 0.013 1.924 0.054
Population density 0.885 0.012 75.865 0.000
Physical elevation -0.007 0.013 -0.540 0.589
Photovoltaic power potential -0.003 0.014 -0.231 0.817
Ozone concentration -0.029 0.012 -2.461 0.014
Distance to coast 0.003 0.013 0.203 0.839
Geographically Weighted Regression (GWR) Results
--------------------------------------------------------------------------------
Coordinates type: Projected
Spatial kernel: Adaptive bisquare
Criterion for optimal bandwidth: AICc
Bandwidth used: 57.000
Bandwidth confidence interval (95%): (57.0, 57.0)
Diagnostic Information
--------------------------------------------------------------------------------
Residual sum of squares: 123.996
Effective number of parameters (trace(S)): 362.849
Degree of freedom (n - trace(S)): 1284.151
Sigma estimate: 0.311
Log-likelihood: -207.039
Degree of Dependency (DoD): 0.467
AIC: 1141.775
AICc: 1348.849
BIC: 3109.001
R2: 0.925
Adj. R2: 0.903
Adj. alpha (95%): 0.001
Adj. critical t value (95%): 3.307
Monte Carlo Test for Spatial Variability
--------------------------------------------------------------------------------
Variable p-value
Intercept 0.000
Land surface temperature 0.779
Population density 0.725
Physical elevation 0.000
Photovoltaic power potential 0.000
Ozone concentration 0.000
Distance to coast 0.000
Summary Statistics For GWR Parameter Estimates
--------------------------------------------------------------------------------
Variable Mean STD Min Median Max
-------------------- ---------- ---------- ---------- ---------- ----------
Intercept 0.188 1.614 -3.337 -0.049 26.031
Land surface temperature 0.074 0.464 -3.555 0.015 4.176
Population density 0.736 0.418 -0.191 0.737 2.555
Physical elevation -0.040 0.585 -11.458 0.009 1.354
Photovoltaic power potential -0.072 0.421 -4.368 -0.004 0.886
Ozone concentration -0.274 1.856 -31.683 0.004 3.379
Distance to coast -0.081 0.572 -4.825 -0.001 1.078
================================================================================
Acknowledgement:
We acknowledge the support of the National Science Foundation under Award 1758786
from the Geography and Spatial Sciences Program to A. S. Fotheringham which
enabled this software to be written and made freely available.
================================================================================
Execution error
ERROR: Error in parse(text = x, srcfile = src): <text>:1:1: unexpected '=='
1: ==
^
Execution error
KernelInterrupted: Execution interrupted by the Jupyter kernel.