Install and load R packages
install.packages(c("sf","gdverse"), dep = T)
# install.packages("devtools")
devtools::install_github("stscl/sesp",build_vignettes = TRUE,dep = TRUE)
library(sf)
## Linking to GEOS 3.12.1, GDAL 3.8.4, PROJ 9.3.1; sf_use_s2() is TRUE
library(sesp)
library(gdverse)Using the same data as the gdverse idsa vignette:
depression = system.file('extdata/Depression.csv',package = 'gdverse') |>
readr::read_csv() |>
sf::st_as_sf(coords = c('X','Y'), crs = 4326)
## Rows: 1072 Columns: 13
## ── Column specification ───────────────────────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (13): X, Y, Depression_prevelence, PopulationDensity, Population65, NoHealthInsurance, Neig...
##
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
depression
## Simple feature collection with 1072 features and 11 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -83.1795 ymin: 32.11464 xmax: -78.6023 ymax: 35.17354
## Geodetic CRS: WGS 84
## # A tibble: 1,072 × 12
## Depression_prevelence PopulationDensity Population65 NoHealthInsurance Neighbor_Disadvantage
## * <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 23.1 61.5 22.5 7.98 -0.0525
## 2 22.8 58.3 16.8 11.0 -0.254
## 3 23.2 35.9 24.5 9.31 -0.0540
## 4 21.8 76.1 21.8 13.2 0.0731
## 5 20.7 47.3 22.0 11 0.763
## 6 21.3 32.5 19.2 13.0 0.422
## 7 22 36.9 19.2 10.8 0.113
## 8 21.2 61.5 15.9 8.57 -0.154
## 9 22.7 67.2 15.7 17.8 -0.320
## 10 20.6 254. 11.3 12.7 0.457
## # ℹ 1,062 more rows
## # ℹ 7 more variables: Beer <dbl>, MentalHealthPati <dbl>, NatureParks <dbl>, Casinos <dbl>,
## # DrinkingPlaces <dbl>, X.HouseRent <dbl>, geometry <POINT [°]>SESP With Linear Regression
system.time({
g1 = sesp(Depression_prevelence ~ ., data = depression,
model = 'ols', overlay = 'intersection', cores = 8)
})
## user system elapsed
## 4.70 0.68 28.17
g1
## *** Spatially Explicit Stratified Power
##
## Q values are estimated using *Ordinary Least Square*
##
## -------------- Global Power of Determinant : ------------
## Variable Qvalue AIC BIC LogLik
## Neighbor_Disadvantage 0.296 4370.901 4370.901 -2179.450
## MentalHealthPati 0.260 4428.205 4428.205 -2206.103
## NoHealthInsurance 0.255 4435.513 4435.513 -2209.756
## X.HouseRent 0.225 4476.344 4476.344 -2231.172
## Casinos 0.205 4503.507 4503.507 -2244.754
## DrinkingPlaces 0.203 4503.855 4503.855 -2245.927
## PopulationDensity 0.202 4505.365 4505.365 -2246.682
## Beer 0.200 4507.758 4507.758 -2247.879
## NatureParks 0.181 4533.589 4533.589 -2260.794
## Population65 0.159 4561.304 4561.304 -2274.652
##
## ------------- Global Variable Interaction : ------------
## Variable Interaction
## PopulationDensity ∩ Population65 Enhance, bi-
## PopulationDensity ∩ NoHealthInsurance Enhance, bi-
## PopulationDensity ∩ Neighbor_Disadvantage Enhance, bi-
## PopulationDensity ∩ Beer Enhance, bi-
## PopulationDensity ∩ MentalHealthPati Enhance, bi-
## PopulationDensity ∩ NatureParks Enhance, bi-
## PopulationDensity ∩ Casinos Enhance, bi-
## PopulationDensity ∩ DrinkingPlaces Enhance, bi-
## PopulationDensity ∩ X.HouseRent Enhance, bi-
## Population65 ∩ NoHealthInsurance Enhance, bi-
##
## ! Only the top ten items of global scale are displayed.
## ! The others can be accessed through specific subsets.
plot(g1,slicenum = 10) +
ggplot2::theme(axis.text.x = ggplot2::element_text(
angle = 30,vjust = 0.85,hjust = 0.75)
)
Figure 1. Results of SESP With
Linear Regression
SESP With Spatial Lag Regression
system.time({
g2 = sesp(Depression_prevelence ~ .,
data = depression,
model = 'lag', cores = 8)
})
## user system elapsed
## 7.80 1.07 252.03
g2
## *** Spatially Explicit Stratified Power
##
## Q values are estimated using *Spatial Lag Model*
##
## -------------- Global Power of Determinant : ------------
## Variable Qvalue AIC BIC LogLik
## Neighbor_Disadvantage 0.289 4248.866 4248.866 -2117.433
## MentalHealthPati 0.240 4288.518 4288.518 -2135.259
## NoHealthInsurance 0.216 4267.870 4267.870 -2124.935
## X.HouseRent 0.187 4316.759 4316.759 -2150.380
## Beer 0.187 4319.936 4319.936 -2149.968
## DrinkingPlaces 0.160 4330.763 4330.763 -2158.382
## PopulationDensity 0.154 4333.170 4333.170 -2159.585
## Casinos 0.153 4346.299 4346.299 -2166.149
## NatureParks 0.138 4364.690 4364.690 -2175.345
## Population65 0.122 4370.485 4370.485 -2178.242
##
## ------------- Global Variable Interaction : ------------
## Variable Interaction
## PopulationDensity ∩ Population65 Enhance, bi-
## PopulationDensity ∩ NoHealthInsurance Enhance, bi-
## PopulationDensity ∩ Neighbor_Disadvantage Enhance, bi-
## PopulationDensity ∩ Beer Enhance, bi-
## PopulationDensity ∩ MentalHealthPati Weaken, uni-
## PopulationDensity ∩ NatureParks Weaken, uni-
## PopulationDensity ∩ Casinos Enhance, bi-
## PopulationDensity ∩ DrinkingPlaces Weaken, uni-
## PopulationDensity ∩ X.HouseRent Enhance, bi-
## Population65 ∩ NoHealthInsurance Weaken, uni-
##
## ! Only the top ten items of global scale are displayed.
## ! The others can be accessed through specific subsets.
plot(g2,slicenum = 10) +
ggplot2::theme(axis.text.x = ggplot2::element_text(
angle = 30,vjust = 0.85,hjust = 0.75)
)
Figure 2. Results of SESP With
Spatial Lag Regression
SESP With Spatial Error Regression
system.time({
g3 = sesp(Depression_prevelence ~ .,
data = depression,
model = 'error', cores = 8)
})
## user system elapsed
## 5.35 0.23 135.22
g3
## *** Spatially Explicit Stratified Power
##
## Q values are estimated using *Spatial Error Model*
##
## -------------- Global Power of Determinant : ------------
## Variable Qvalue AIC BIC LogLik
## Neighbor_Disadvantage 0.289 4238.682 4238.682 -2112.341
## NoHealthInsurance 0.253 4230.432 4230.432 -2106.216
## MentalHealthPati 0.247 4262.110 4262.110 -2122.055
## X.HouseRent 0.219 4279.959 4279.959 -2131.979
## PopulationDensity 0.199 4288.640 4288.640 -2137.320
## DrinkingPlaces 0.198 4295.621 4295.621 -2140.811
## Beer 0.193 4292.395 4292.395 -2139.197
## Casinos 0.184 4305.020 4305.020 -2145.510
## NatureParks 0.171 4325.679 4325.679 -2155.839
## Population65 0.156 4335.830 4335.830 -2160.915
##
## ------------- Global Variable Interaction : ------------
## Variable Interaction
## PopulationDensity ∩ Population65 Weaken, uni-
## PopulationDensity ∩ NoHealthInsurance Enhance, bi-
## PopulationDensity ∩ Neighbor_Disadvantage Enhance, bi-
## PopulationDensity ∩ Beer Enhance, bi-
## PopulationDensity ∩ MentalHealthPati Enhance, bi-
## PopulationDensity ∩ NatureParks Weaken, uni-
## PopulationDensity ∩ Casinos Enhance, bi-
## PopulationDensity ∩ DrinkingPlaces Enhance, bi-
## PopulationDensity ∩ X.HouseRent Enhance, bi-
## Population65 ∩ NoHealthInsurance Weaken, uni-
##
## ! Only the top ten items of global scale are displayed.
## ! The others can be accessed through specific subsets.
plot(g3,slicenum = 10) +
ggplot2::theme(axis.text.x = ggplot2::element_text(
angle = 30,vjust = 0.85,hjust = 0.75)
)
Figure 3. Results of SESP With
Spatial Error Regression
SESP With Spatial Durbin Regression
system.time({
g4 = sesp(Depression_prevelence ~ ., data = depression,
model = 'lag', durbin = TRUE, cores = 8)
})
## user system elapsed
## 8.67 1.39 248.84
g4
## *** Spatially Explicit Stratified Power
##
## Q values are estimated using *Spatial Durbin Model*
##
## -------------- Global Power of Determinant : ------------
## Variable Qvalue AIC BIC LogLik
## Neighbor_Disadvantage 0.316 4226.864 4226.864 -2102.432
## MentalHealthPati 0.296 4247.222 4247.222 -2108.611
## NoHealthInsurance 0.261 4238.839 4238.839 -2104.420
## NatureParks 0.252 4285.895 4285.895 -2125.948
## X.HouseRent 0.239 4280.458 4280.458 -2127.229
## Casinos 0.237 4291.052 4291.052 -2132.526
## PopulationDensity 0.211 4291.372 4291.372 -2134.686
## DrinkingPlaces 0.210 4294.833 4294.833 -2136.417
## Beer 0.209 4291.057 4291.057 -2134.529
## Population65 0.163 4340.250 4340.250 -2159.125
##
## ------------- Global Variable Interaction : ------------
## Variable Interaction
## PopulationDensity ∩ Population65 Enhance, bi-
## PopulationDensity ∩ NoHealthInsurance Enhance, bi-
## PopulationDensity ∩ Neighbor_Disadvantage Enhance, bi-
## PopulationDensity ∩ Beer Enhance, bi-
## PopulationDensity ∩ MentalHealthPati Enhance, bi-
## PopulationDensity ∩ NatureParks Enhance, bi-
## PopulationDensity ∩ Casinos Enhance, bi-
## PopulationDensity ∩ DrinkingPlaces Enhance, bi-
## PopulationDensity ∩ X.HouseRent Enhance, bi-
## Population65 ∩ NoHealthInsurance Enhance, bi-
##
## ! Only the top ten items of global scale are displayed.
## ! The others can be accessed through specific subsets.
plot(g4,slicenum = 10) +
ggplot2::theme(axis.text.x = ggplot2::element_text(
angle = 30,vjust = 0.85,hjust = 0.75)
)
Figure 4. Results of SESP With
Spatial Durbin Regression
SESP With Geographically Weighted Regression
system.time({
g5 = sesp(Depression_prevelence ~ .,
data = depression,
model = 'gwr', cores = 8)
})
## user system elapsed
## 8.35 0.78 147.83
g5
## *** Spatially Explicit Stratified Power
##
## Q values are estimated using *Geographically Weighted Regression*
##
## -------------- Global Power of Determinant : ------------
## Variable Qvalue AIC
## X.HouseRent 0.370 4255.989
## MentalHealthPati 0.367 4258.876
## PopulationDensity 0.366 4260.457
## Neighbor_Disadvantage 0.354 4277.147
## Beer 0.350 4288.285
## Population65 0.317 4337.745
## Casinos 0.307 4354.933
## DrinkingPlaces 0.294 4373.148
## NatureParks 0.266 4412.671
## NoHealthInsurance 0.254 4429.183
##
## ------------- Global Variable Interaction : ------------
## Variable Interaction
## PopulationDensity ∩ Population65 Enhance, bi-
## PopulationDensity ∩ NoHealthInsurance Enhance, bi-
## PopulationDensity ∩ Neighbor_Disadvantage Enhance, bi-
## PopulationDensity ∩ Beer Enhance, bi-
## PopulationDensity ∩ MentalHealthPati Enhance, bi-
## PopulationDensity ∩ NatureParks Weaken, uni-
## PopulationDensity ∩ Casinos Weaken, uni-
## PopulationDensity ∩ DrinkingPlaces Weaken, uni-
## PopulationDensity ∩ X.HouseRent Enhance, bi-
## Population65 ∩ NoHealthInsurance Weaken, uni-
##
## ! Only the top ten items of global scale are displayed.
## ! The others can be accessed through specific subsets.
plot(g5,slicenum = 10) +
ggplot2::theme(axis.text.x = ggplot2::element_text(
angle = 30,vjust = 0.85,hjust = 0.75)
)
Figure 5. Results of SESP With
Geographically Weighted Regression
Results of optimal spatial discretization
plot_optdisc = \(g){
gmap = sf::st_set_geometry(g$optdisc,sf::st_geometry(depression))
fig1 = seq_along(g$optdisc) |>
purrr::map(\(.x) ggplot2::ggplot(data = gmap) +
ggplot2::geom_sf(ggplot2::aes(color = factor(g$optdisc[,.x,drop = TRUE])),
alpha = .65, size = 0.5) +
ggplot2::labs(color = 'zones') +
ggplot2::theme_void() +
ggplot2::theme(
legend.position = "none")
) %>%
patchwork::wrap_plots(ncol = 3, byrow = TRUE) +
patchwork::plot_annotation(
tag_levels = 'a',
tag_prefix = '(',
tag_suffix = ')',
tag_sep = '',
theme = ggplot2::theme(plot.tag = ggplot2::element_text(family = "serif"),
plot.tag.position = "topleft"))
return(fig1)
}
plot_optdisc(g1)
Figure 6. Optimal spatial
discretization result with linear regression operator. Subfigures
(a)–(j) depict the optimal spatial discretizations for the variables
PopulationDensity, Population65, NoHealthInsurance,
Neighbor_Disadvantage, Beer, MentalHealthPati, NatureParks, Casinos,
DrinkingPlaces, and X.HouseRent, respectively.
plot_optdisc(g2)
Figure 7. Optimal spatial
discretization result with spatial lag regression operator. Subfigures
(a)–(j) depict the optimal spatial discretizations for the variables
PopulationDensity, Population65, NoHealthInsurance,
Neighbor_Disadvantage, Beer, MentalHealthPati, NatureParks, Casinos,
DrinkingPlaces, and X.HouseRent, respectively.
plot_optdisc(g3)
Figure 8. Optimal spatial
discretization result with spatial error regression operator. Subfigures
(a)–(j) depict the optimal spatial discretizations for the variables
PopulationDensity, Population65, NoHealthInsurance,
Neighbor_Disadvantage, Beer, MentalHealthPati, NatureParks, Casinos,
DrinkingPlaces, and X.HouseRent, respectively.
plot_optdisc(g4)
Figure 9. Optimal spatial
discretization result with spatial durbin regression operator.
Subfigures (a)–(j) depict the optimal spatial discretizations for the
variables PopulationDensity, Population65, NoHealthInsurance,
Neighbor_Disadvantage, Beer, MentalHealthPati, NatureParks, Casinos,
DrinkingPlaces, and X.HouseRent, respectively.
plot_optdisc(g5)
Figure 10. Optimal spatial
discretization result with geographically weighted regression operator.
Subfigures (a)–(j) depict the optimal spatial discretizations for the
variables PopulationDensity, Population65, NoHealthInsurance,
Neighbor_Disadvantage, Beer, MentalHealthPati, NatureParks, Casinos,
DrinkingPlaces, and X.HouseRent, respectively.