Spatial Regression Discussion
Prabish Khadka Chhetri
2022-12-02
setwd("C:/Spatial Statistics")
suppressPackageStartupMessages(library(maptools))
library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(ggplot2)
library(sp)
library(spgwr)## 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')`## NOTE: This package does not constitute approval of GWR
## as a method of spatial analysis; see example(gwr)library(RColorBrewer)
ppdf = readShapeSpatial("C:/Spatial Statistics/pov_shp/pov_shp1.shp")## Warning: shapelib support is provided by GDAL through the sf and terra packages
## among others## Warning: shapelib support is provided by GDAL through the sf and terra paackages
## among others## Warning: shapelib support is provided by GDAL through the sf and terra packages
## among othersslotNames(ppdf)## [1] "data" "polygons" "plotOrder" "bbox" "proj4string"str(ppdf, max.level = 2)## Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
## ..@ data :'data.frame': 3070 obs. of 21 variables:
## .. ..- attr(*, "data_types")= chr [1:21] "N" "C" "N" "N" ...
## ..@ polygons :List of 3070
## ..@ plotOrder : int [1:3070] 2893 1357 2998 1360 3002 19 871 3003 882 1352 ...
## ..@ bbox : num [1:2, 1:2] -124.7 24.5 -66.9 49.4
## .. ..- attr(*, "dimnames")=List of 2
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
## ..$ comment: chr "FALSE"head(slot(ppdf, "data"))## FIPS_N OBJECTID Shape_Leng Shape_Area OID_ ag STATEFP COUNTYFP COUNTYNS
## 0 56001 35 5.382173 1.2062731 3160 3.13 56 001 01605066
## 1 56015 36 3.297442 0.6290952 3167 12.03 56 015 01605073
## 2 56045 37 3.408379 0.6956611 3182 5.70 56 045 01605086
## 3 56017 38 4.663761 0.5804015 3168 8.03 56 017 01605074
## 4 56039 39 5.118661 1.2240465 3179 2.87 56 039 01605083
## 5 56003 40 4.781177 0.9262142 3161 10.79 56 003 01605067
## GEOID NAME INTPTLAT INTPTLON black feemp foreign hisp hsch
## 0 56001 Albany +41.6655141 -105.7218826 1.11 62.74 3.80 7.49 21.90
## 1 56015 Goshen +42.0895800 -104.3534821 0.20 52.08 1.89 8.83 33.43
## 2 56045 Weston +43.8462133 -104.5700202 0.12 50.06 0.77 2.06 40.21
## 3 56017 Hot Springs +43.7208707 -108.4356515 0.35 59.48 1.31 2.38 36.13
## 4 56039 Teton +44.0493205 -110.5881023 0.15 71.88 5.88 6.49 18.91
## 5 56003 Big Horn +44.5256543 -107.9933205 0.11 48.56 2.21 6.17 34.13
## manu povty retail
## 0 4.15 23.64 9.94
## 1 4.42 13.23 10.94
## 2 7.10 7.89 10.19
## 3 3.13 10.77 6.63
## 4 1.97 5.88 10.95
## 5 5.77 12.11 10.83plot(ppdf)
Regression:
modelppdf = lm(povty ~ ag + manu + retail + foreign + feemp + hsch + black + hisp, data = ppdf)
summary(modelppdf)##
## Call:
## lm(formula = povty ~ ag + manu + retail + foreign + feemp + hsch +
## black + hisp, data = ppdf)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.9991 -2.2388 -0.3421 1.7029 29.0205
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 46.462367 1.019860 45.558 < 2e-16 ***
## ag 0.104110 0.012701 8.197 3.58e-16 ***
## manu -0.028851 0.009740 -2.962 0.00308 **
## retail -0.064335 0.039627 -1.624 0.10458
## foreign -0.171651 0.021274 -8.069 1.01e-15 ***
## feemp -0.526357 0.011166 -47.141 < 2e-16 ***
## hsch -0.188037 0.014749 -12.749 < 2e-16 ***
## black 0.081289 0.005466 14.872 < 2e-16 ***
## hisp 0.061865 0.008176 7.567 5.03e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.809 on 3061 degrees of freedom
## Multiple R-squared: 0.5861, Adjusted R-squared: 0.585
## F-statistic: 541.7 on 8 and 3061 DF, p-value: < 2.2e-16step(modelppdf)## Start: AIC=8219.82
## povty ~ ag + manu + retail + foreign + feemp + hsch + black +
## hisp
##
## Df Sum of Sq RSS AIC
## <none> 44402 8219.8
## - retail 1 38 44440 8220.5
## - manu 1 127 44529 8226.6
## - hisp 1 831 45232 8274.7
## - foreign 1 944 45346 8282.4
## - ag 1 975 45376 8284.5
## - hsch 1 2358 46759 8376.6
## - black 1 3208 47610 8432.0
## - feemp 1 32236 76638 9893.4##
## Call:
## lm(formula = povty ~ ag + manu + retail + foreign + feemp + hsch +
## black + hisp, data = ppdf)
##
## Coefficients:
## (Intercept) ag manu retail foreign feemp
## 46.46237 0.10411 -0.02885 -0.06434 -0.17165 -0.52636
## hsch black hisp
## -0.18804 0.08129 0.06186ppdf$predpov = predict(modelppdf)
range(ppdf$predpov)## [1] -0.06992971 30.39441443- All the variables are significant.
Geographic Regression:
bw = gwr.sel(povty ~ ag + manu + retail + foreign + feemp + hsch +
black + hisp, data = ppdf)## Bandwidth: 23.39464 CV score: 44181.4
## Bandwidth: 37.81553 CV score: 44532.96
## Bandwidth: 14.48204 CV score: 43438.8
## Bandwidth: 8.973749 CV score: 42079.29
## Bandwidth: 5.569439 CV score: 39784.72
## Bandwidth: 3.465459 CV score: 36558.46
## Bandwidth: 2.165128 CV score: 33847.52
## Bandwidth: 1.36148 CV score: 33293.31
## Bandwidth: 1.243318 CV score: 33726.2
## Bandwidth: 1.690287 CV score: 33118.18
## Bandwidth: 1.87166 CV score: 33331.42
## Bandwidth: 1.605416 CV score: 33071.19
## Bandwidth: 1.585686 CV score: 33067.02
## Bandwidth: 1.563302 CV score: 33065.9
## Bandwidth: 1.567936 CV score: 33065.81
## Bandwidth: 1.568232 CV score: 33065.8
## Bandwidth: 1.568177 CV score: 33065.8
## Bandwidth: 1.568136 CV score: 33065.8
## Bandwidth: 1.568177 CV score: 33065.8ppdfmodel.gwr = gwr(povty ~ ag + manu + retail + foreign + feemp + hsch +
black + hisp, data = ppdf, bandwidth = bw)
range(ppdfmodel.gwr$SDF$pred)## [1] 0.3623721 45.2962760Chorolopeth Map:
brks = seq(0,50,5)
cr = brewer.pal(10, "Blues")## Warning in brewer.pal(10, "Blues"): n too large, allowed maximum for palette Blues is 9
## Returning the palette you asked for with that many colorsspplot(ppdfmodel.gwr$SDF, "pred",
col.regions = cr, at = brks,
colorkey = list(space = "bottom"),
sub = " Prediction of Poverty")
- The model best predicts on counties of texas, kentucky, south dakota.
For best prediction of create Choropleth map of LocalR2 column in the
spatial polygon data frame.
spplot(ppdfmodel.gwr$SDF, "localR2",
col.regions = cr, at = brks,
colorkey = list(space = "bottom"),
sub = " Local Regression Fit")
x