labten
Serena Madsen
11/14/2021
library(sf)## Warning: package 'sf' was built under R version 3.6.2## Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0library(spdep)## Warning: package 'spdep' was built under R version 3.6.2## Loading required package: sp## Warning: package 'sp' was built under R version 3.6.2## Loading required package: spData## Warning: package 'spData' was built under R version 3.6.2## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`library(spatialreg)## Warning: package 'spatialreg' was built under R version 3.6.2## Loading required package: Matrix## Warning: package 'Matrix' was built under R version 3.6.2##
## Attaching package: 'spatialreg'## The following objects are masked from 'package:spdep':
##
## aple, aple.mc, aple.plot, as_dgRMatrix_listw, as_dsCMatrix_I,
## as_dsCMatrix_IrW, as_dsTMatrix_listw, as.spam.listw, can.be.simmed,
## cheb_setup, create_WX, do_ldet, eigen_pre_setup, eigen_setup,
## eigenw, errorsarlm, get.ClusterOption, get.coresOption,
## get.mcOption, get.VerboseOption, get.ZeroPolicyOption,
## GMargminImage, GMerrorsar, griffith_sone, gstsls, Hausman.test,
## impacts, intImpacts, invIrM, invIrW, Jacobian_W, jacobianSetup,
## l_max, lagmess, lagsarlm, lextrB, lextrS, lextrW, lmSLX, localAple,
## LU_prepermutate_setup, LU_setup, Matrix_J_setup, Matrix_setup,
## mcdet_setup, MCMCsamp, ME, mom_calc, mom_calc_int2, moments_setup,
## powerWeights, sacsarlm, SE_classic_setup, SE_interp_setup,
## SE_whichMin_setup, set.ClusterOption, set.coresOption,
## set.mcOption, set.VerboseOption, set.ZeroPolicyOption,
## similar.listw, spam_setup, spam_update_setup, SpatialFiltering,
## spautolm, spBreg_err, spBreg_lag, spBreg_sac, stsls,
## subgraph_eigenw, trWlibrary(spgwr)## Warning: package 'spgwr' was built under R version 3.6.2## NOTE: This package does not constitute approval of GWR
## as a method of spatial analysis; see example(gwr)library(tmap)## Registered S3 methods overwritten by 'stars':
## method from
## st_bbox.SpatRaster sf
## st_crs.SpatRaster sfsetwd("~/Documents/geog6000/lab10")
bos <- st_read("../datafiles/boston.tr/boston.shp", quiet = TRUE)
bos$logCMEDV <- log10(bos$CMEDV)tmap_mode("view")## tmap mode set to interactive viewingtm_shape(bos) +
tm_fill(col = "logCMEDV", palette = "-magma", main = "Housing Prices in Boston", main = "Housing Prices in Boston")Basic Linear Regression:
boslm <- lm(logCMEDV ~ B+CRIM+NOX, data = bos)
summary(boslm)##
## Call:
## lm(formula = logCMEDV ~ B + CRIM + NOX, data = bos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.46704 -0.08406 -0.02183 0.06848 0.47025
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.499e+00 4.826e-02 31.072 < 2e-16 ***
## B 3.031e-04 7.519e-05 4.032 6.40e-05 ***
## CRIM -6.983e-03 8.137e-04 -8.582 < 2e-16 ***
## NOX -4.767e-01 6.027e-02 -7.909 1.67e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1374 on 502 degrees of freedom
## Multiple R-squared: 0.4027, Adjusted R-squared: 0.3991
## F-statistic: 112.8 on 3 and 502 DF, p-value: < 2.2e-16A test of spatial autocorrelation in the residuals of this model
bos.nbq <- poly2nb(bos)
bos.nbq## Neighbour list object:
## Number of regions: 506
## Number of nonzero links: 2910
## Percentage nonzero weights: 1.136559
## Average number of links: 5.750988bos.listw <- nb2listw(bos.nbq)
bos.geom <- st_geometry(bos)
plot(bos.geom,
col = "pink",
border = "blue")
coords <- st_coordinates(st_centroid(bos))## Warning in st_centroid.sf(bos): st_centroid assumes attributes are constant over
## geometries of xplot(bos.nbq, coords,
add = TRUE,
col = 1,
lwd = 1)
##I wanted to see how ridiculous I could make this map look, hence the chaotic color scheme. moran.test(residuals(boslm), bos.listw)##
## Moran I test under randomisation
##
## data: residuals(boslm)
## weights: bos.listw
##
## Moran I statistic standard deviate = 21.89, p-value < 2.2e-16
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.5874745944 -0.0019801980 0.0007250964Because the autocorrelation is significant, (Moran I = 0.58747, p-value < 2.2e-16), I am going to build a spatial regression model to try to account for autocorrelation.
Legrange Multiplier Test:
summary(lm.LMtests(boslm, bos.listw, test=c("LMerr","LMlag")))## Lagrange multiplier diagnostics for spatial dependence
## data:
## model: lm(formula = logCMEDV ~ B + CRIM + NOX, data = bos)
## weights: bos.listw
##
## statistic parameter p.value
## LMerr 469.39 1 < 2.2e-16 ***
## LMlag 450.05 1 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Robust Legrrange Multiplier Test: I chose the Legrange Multiplier Test because it is an approach that I have done before and it is an easily interpreted way to check for autocorrelation.
summary(lm.LMtests(boslm, bos.listw, test=c("RLMerr","RLMlag")))## Lagrange multiplier diagnostics for spatial dependence
## data:
## model: lm(formula = logCMEDV ~ B + CRIM + NOX, data = bos)
## weights: bos.listw
##
## statistic parameter p.value
## RLMerr 31.238 1 2.283e-08 ***
## RLMlag 11.900 1 0.0005613 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Spatial Filtering Model Because Why Not?!?!?! Because the RLMlag was not significant and RLMerr was not, I am going to run a spatial filtering model.
##Build Spatial Filter:
sf.err.bs <- SpatialFiltering(boslm,
nb = bos.nbq,
style = "C",
alpha = 0.25,
ExactEV = FALSE,
data = bos)
sf.err.bs## Step SelEvec Eval MinMi ZMinMi Pr(ZI) R2
## 0 0 0 0.0000000 0.575073712 22.510428 3.280912e-112 0.4026776
## 1 1 1 1.1190461 0.484361894 19.188836 4.588019e-82 0.4880495
## 2 2 6 1.0366953 0.444651796 17.805074 6.454667e-71 0.5223875
## 3 3 19 0.9030854 0.408750227 16.522676 2.519656e-61 0.5570746
## 4 4 15 0.9262230 0.384431738 15.695865 1.614274e-55 0.5769555
## 5 5 2 1.1085811 0.368912557 15.263472 1.339293e-52 0.5858315
## 6 6 35 0.7617034 0.352772432 14.712858 5.330728e-49 0.6021783
## 7 7 4 1.0675499 0.336265531 14.223411 6.557668e-46 0.6111581
## 8 8 22 0.8631893 0.320793258 13.718900 7.824305e-43 0.6222502
## 9 9 18 0.9095766 0.308185849 13.338740 1.377682e-40 0.6301692
## 10 10 17 0.9152685 0.295315431 12.947770 2.418791e-38 0.6378470
## 11 11 46 0.6786180 0.281939685 12.482015 9.358759e-36 0.6500586
## 12 12 33 0.7693769 0.271674266 12.160331 5.055365e-34 0.6572763
## 13 13 26 0.8257173 0.261477683 11.853063 2.074659e-32 0.6634698
## 14 14 51 0.6540746 0.251463111 11.516837 1.085224e-30 0.6718407
## 15 15 39 0.7262778 0.241953344 11.215131 3.436789e-29 0.6782841
## 16 16 10 0.9868171 0.232554899 10.975322 5.022678e-28 0.6822929
## 17 17 56 0.6064970 0.223017006 10.648000 1.781522e-26 0.6901948
## 18 18 53 0.6309415 0.213058581 10.307761 6.499765e-25 0.6975777
## 19 19 25 0.8387077 0.202996077 10.003926 1.464721e-23 0.7023647
## 20 20 34 0.7668197 0.193996906 9.728096 2.288358e-22 0.7070406
## 21 21 24 0.8424111 0.184860224 9.461359 3.039669e-21 0.7111113
## 22 22 31 0.7911629 0.175700580 9.182165 4.225084e-20 0.7154107
## 23 23 43 0.6953790 0.167004230 8.902894 5.440897e-19 0.7200946
## 24 24 57 0.6046878 0.158195099 8.602206 7.819836e-18 0.7256171
## 25 25 8 1.0129017 0.149526942 8.387996 4.945002e-17 0.7283718
## 26 26 28 0.8154846 0.140938679 8.133988 4.153940e-16 0.7318302
## 27 27 45 0.6858331 0.132132414 7.845284 4.319756e-15 0.7360952
## 28 28 11 0.9861625 0.123135587 7.607149 2.802083e-14 0.7388464
## 29 29 5 1.0548118 0.113965581 7.375037 1.642995e-13 0.7413917
## 30 30 14 0.9477277 0.104602164 7.109289 1.166425e-12 0.7442637
## 31 31 42 0.6982104 0.095170452 6.790213 1.119678e-11 0.7482635
## 32 32 82 0.4551575 0.086948302 6.483022 8.990368e-11 0.7538848
## 33 33 95 0.3882970 0.078496092 6.156537 7.435291e-10 0.7605995
## 34 34 12 0.9644707 0.069894045 5.920428 3.211052e-09 0.7629015
## 35 35 84 0.4442061 0.062164629 5.632128 1.779997e-08 0.7676985
## 36 36 59 0.5895399 0.054750800 5.379338 7.476004e-08 0.7709189
## 37 37 77 0.4848752 0.047211433 5.105169 3.304993e-07 0.7748651
## 38 38 65 0.5493895 0.039763174 4.844321 1.270452e-06 0.7781555
## 39 39 32 0.7852379 0.032365081 4.622074 3.799233e-06 0.7803355
## 40 40 78 0.4778136 0.024865798 4.347750 1.375412e-05 0.7839724
## 41 41 23 0.8576474 0.017220422 4.124336 3.718062e-05 0.7859376
## 42 42 41 0.7034361 0.010887488 3.932966 8.390407e-05 0.7878950
## 43 43 62 0.5711341 0.004441583 3.716449 2.020422e-04 0.7903076
## 44 44 13 0.9627533 -0.001452761 3.584569 3.376351e-04 0.7915895
## 45 45 61 0.5820736 -0.007184953 3.399930 6.740319e-04 0.7936169
## 46 46 88 0.4208333 -0.012765691 3.200048 1.374049e-03 0.7962732
## 47 47 118 0.2824181 -0.018243709 2.987799 2.809946e-03 0.7999851
## 48 48 83 0.4534988 -0.022865306 2.835583 4.574211e-03 0.8019256
## 49 49 89 0.4194785 -0.027182514 2.692828 7.084889e-03 0.8038401
## 50 50 149 0.1565796 -0.031561643 2.516504 1.185256e-02 0.8084059
## 51 51 101 0.3556568 -0.036054346 2.358675 1.834033e-02 0.8106033
## 52 52 9 0.9946101 -0.040181644 2.303364 2.125837e-02 0.8113587
## 53 53 73 0.5007505 -0.044188716 2.184720 2.890942e-02 0.8127459
## 54 54 29 0.8049140 -0.048177736 2.107463 3.507750e-02 0.8136215
## 55 55 133 0.2211072 -0.052177500 1.955150 5.056537e-02 0.8163493
## 56 56 138 0.2035006 -0.056219881 1.799417 7.195273e-02 0.8192077
## 57 57 30 0.7974901 -0.060141987 1.722919 8.490313e-02 0.8200345
## 58 58 72 0.5046516 -0.064051647 1.608838 1.076518e-01 0.8212717
## 59 59 120 0.2692792 -0.067978973 1.465989 1.426512e-01 0.8233529
## 60 60 156 0.1315408 -0.072020957 1.302924 1.926008e-01 0.8268605
## 61 61 7 1.0246787 -0.076057590 1.249292 2.115581e-01 0.8274954
## 62 62 153 0.1386707 -0.080176805 1.082706 2.789387e-01 0.8307424
## gamma
## 0 0.0000000
## 1 -1.1642332
## 2 0.7383634
## 3 0.7421066
## 4 0.5618237
## 5 0.3753975
## 6 0.5094475
## 7 0.3775864
## 8 0.4196506
## 9 0.3545840
## 10 -0.3491411
## 11 -0.4403199
## 12 0.3385194
## 13 0.3135812
## 14 -0.3645588
## 15 0.3198465
## 16 -0.2522816
## 17 0.3542019
## 18 0.3423693
## 19 0.2756845
## 20 -0.2724682
## 21 -0.2542236
## 22 -0.2612681
## 23 -0.2727026
## 24 0.2961062
## 25 0.2091341
## 26 -0.2343240
## 27 -0.2602231
## 28 -0.2089963
## 29 0.2010277
## 30 -0.2135378
## 31 -0.2520002
## 32 0.2987452
## 33 0.3265092
## 34 0.1911777
## 35 -0.2759720
## 36 0.2261198
## 37 0.2503077
## 38 -0.2285627
## 39 -0.1860398
## 40 -0.2402969
## 41 -0.1766390
## 42 -0.1762910
## 43 -0.1957159
## 44 -0.1426614
## 45 -0.1794113
## 46 0.2053626
## 47 -0.2427611
## 48 0.1755255
## 49 0.1743449
## 50 0.2692399
## 51 -0.1867861
## 52 0.1095154
## 53 -0.1484022
## 54 -0.1179053
## 55 -0.2081086
## 56 0.2130318
## 57 0.1145728
## 58 0.1401535
## 59 0.1817801
## 60 -0.2359853
## 61 -0.1004034
## 62 -0.2270488Visualizing the first selected eigenvector (1)
bos$ME1 = sf.err.bs$dataset[,"vec1"]
tm_shape(bos) +
tm_fill("ME1", palette = "BuGn") +
tm_borders() +
tm_layout(main.title = "Moran eigenvector (1)", main.title.size = 2)fittederr <- fitted(sf.err.bs)
lm.sf.bs <- lm(logCMEDV ~ B + CRIM + NOX + fittederr, data = bos)
summary(boslm)##
## Call:
## lm(formula = logCMEDV ~ B + CRIM + NOX, data = bos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.46704 -0.08406 -0.02183 0.06848 0.47025
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.499e+00 4.826e-02 31.072 < 2e-16 ***
## B 3.031e-04 7.519e-05 4.032 6.40e-05 ***
## CRIM -6.983e-03 8.137e-04 -8.582 < 2e-16 ***
## NOX -4.767e-01 6.027e-02 -7.909 1.67e-14 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.1374 on 502 degrees of freedom
## Multiple R-squared: 0.4027, Adjusted R-squared: 0.3991
## F-statistic: 112.8 on 3 and 502 DF, p-value: < 2.2e-16anova(boslm, lm.sf.bs)## Analysis of Variance Table
##
## Model 1: logCMEDV ~ B + CRIM + NOX
## Model 2: logCMEDV ~ B + CRIM + NOX + fittederr
## Res.Df RSS Df Sum of Sq F Pr(>F)
## 1 502 9.4836
## 2 440 2.6873 62 6.7963 17.948 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Goodness of fit:
After running the Anova test on the two models, I think the new model is a good fit because the p-value is so low.
Below, we see that the adjusted R-squared for the fitted model is 0.8057, the F-statistic is 33.22, & the p-value is < 2.2e-16
summary(lm.sf.bs)##
## Call:
## lm(formula = logCMEDV ~ B + CRIM + NOX + fittederr, data = bos)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.239194 -0.050177 0.001663 0.044416 0.214597
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.499e+00 2.744e-02 54.649 < 2e-16 ***
## B 3.031e-04 4.275e-05 7.091 5.36e-12 ***
## CRIM -6.983e-03 4.627e-04 -15.093 < 2e-16 ***
## NOX -4.767e-01 3.427e-02 -13.910 < 2e-16 ***
## fittederrvec1 -1.164e+00 7.815e-02 -14.897 < 2e-16 ***
## fittederrvec6 7.384e-01 7.815e-02 9.448 < 2e-16 ***
## fittederrvec19 7.421e-01 7.815e-02 9.496 < 2e-16 ***
## fittederrvec15 5.618e-01 7.815e-02 7.189 2.82e-12 ***
## fittederrvec2 3.754e-01 7.815e-02 4.804 2.14e-06 ***
## fittederrvec35 5.094e-01 7.815e-02 6.519 1.95e-10 ***
## fittederrvec4 3.776e-01 7.815e-02 4.832 1.87e-06 ***
## fittederrvec22 4.197e-01 7.815e-02 5.370 1.28e-07 ***
## fittederrvec18 3.546e-01 7.815e-02 4.537 7.36e-06 ***
## fittederrvec17 -3.491e-01 7.815e-02 -4.468 1.01e-05 ***
## fittederrvec46 -4.403e-01 7.815e-02 -5.634 3.14e-08 ***
## fittederrvec33 3.385e-01 7.815e-02 4.332 1.83e-05 ***
## fittederrvec26 3.136e-01 7.815e-02 4.013 7.06e-05 ***
## fittederrvec51 -3.646e-01 7.815e-02 -4.665 4.10e-06 ***
## fittederrvec39 3.198e-01 7.815e-02 4.093 5.07e-05 ***
## fittederrvec10 -2.523e-01 7.815e-02 -3.228 0.001339 **
## fittederrvec56 3.542e-01 7.815e-02 4.532 7.53e-06 ***
## fittederrvec53 3.424e-01 7.815e-02 4.381 1.48e-05 ***
## fittederrvec25 2.757e-01 7.815e-02 3.528 0.000463 ***
## fittederrvec34 -2.725e-01 7.815e-02 -3.486 0.000539 ***
## fittederrvec24 -2.542e-01 7.815e-02 -3.253 0.001230 **
## fittederrvec31 -2.613e-01 7.815e-02 -3.343 0.000899 ***
## fittederrvec43 -2.727e-01 7.815e-02 -3.489 0.000533 ***
## fittederrvec57 2.961e-01 7.815e-02 3.789 0.000172 ***
## fittederrvec8 2.091e-01 7.815e-02 2.676 0.007728 **
## fittederrvec28 -2.343e-01 7.815e-02 -2.998 0.002868 **
## fittederrvec45 -2.602e-01 7.815e-02 -3.330 0.000942 ***
## fittederrvec11 -2.090e-01 7.815e-02 -2.674 0.007768 **
## fittederrvec5 2.010e-01 7.815e-02 2.572 0.010429 *
## fittederrvec14 -2.135e-01 7.815e-02 -2.732 0.006541 **
## fittederrvec42 -2.520e-01 7.815e-02 -3.225 0.001356 **
## fittederrvec82 2.987e-01 7.815e-02 3.823 0.000151 ***
## fittederrvec95 3.265e-01 7.815e-02 4.178 3.55e-05 ***
## fittederrvec12 1.912e-01 7.815e-02 2.446 0.014824 *
## fittederrvec84 -2.760e-01 7.815e-02 -3.531 0.000457 ***
## fittederrvec59 2.261e-01 7.815e-02 2.893 0.004000 **
## fittederrvec77 2.503e-01 7.815e-02 3.203 0.001459 **
## fittederrvec65 -2.286e-01 7.815e-02 -2.925 0.003627 **
## fittederrvec32 -1.860e-01 7.815e-02 -2.381 0.017712 *
## fittederrvec78 -2.403e-01 7.815e-02 -3.075 0.002238 **
## fittederrvec23 -1.766e-01 7.815e-02 -2.260 0.024294 *
## fittederrvec41 -1.763e-01 7.815e-02 -2.256 0.024574 *
## fittederrvec62 -1.957e-01 7.815e-02 -2.504 0.012628 *
## fittederrvec13 -1.427e-01 7.815e-02 -1.825 0.068607 .
## fittederrvec61 -1.794e-01 7.815e-02 -2.296 0.022161 *
## fittederrvec88 2.054e-01 7.815e-02 2.628 0.008894 **
## fittederrvec118 -2.428e-01 7.815e-02 -3.106 0.002017 **
## fittederrvec83 1.755e-01 7.815e-02 2.246 0.025199 *
## fittederrvec89 1.743e-01 7.815e-02 2.231 0.026191 *
## fittederrvec149 2.692e-01 7.815e-02 3.445 0.000625 ***
## fittederrvec101 -1.868e-01 7.815e-02 -2.390 0.017264 *
## fittederrvec9 1.095e-01 7.815e-02 1.401 0.161816
## fittederrvec73 -1.484e-01 7.815e-02 -1.899 0.058227 .
## fittederrvec29 -1.179e-01 7.815e-02 -1.509 0.132093
## fittederrvec133 -2.081e-01 7.815e-02 -2.663 0.008030 **
## fittederrvec138 2.130e-01 7.815e-02 2.726 0.006668 **
## fittederrvec30 1.146e-01 7.815e-02 1.466 0.143347
## fittederrvec72 1.402e-01 7.815e-02 1.793 0.073598 .
## fittederrvec120 1.818e-01 7.815e-02 2.326 0.020470 *
## fittederrvec156 -2.360e-01 7.815e-02 -3.020 0.002678 **
## fittederrvec7 -1.004e-01 7.815e-02 -1.285 0.199556
## fittederrvec153 -2.270e-01 7.815e-02 -2.905 0.003854 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.07815 on 440 degrees of freedom
## Multiple R-squared: 0.8307, Adjusted R-squared: 0.8057
## F-statistic: 33.22 on 65 and 440 DF, p-value: < 2.2e-16Which variables influence house prices in Boston? (direction, strength, significance)
B 3.031e-04 4.275e-05 7.091 5.36e-12 ***
CRIM -6.983e-03 4.627e-04 -15.093 < 2e-16 ***
NOX -4.767e-01 3.427e-02 -13.910 < 2e-16 ***
B (proportion African-American) significantly increases housing prices in Boston.
CRIM (per capita crime) significantly decreases house prices.
NOX (nitric oxides concentration (parts per million)) significantly decreases house prices but not as much as per capita crime.
Testing for remaining spatial autocorrelation:
moran.test(residuals(lm.sf.bs), bos.listw)##
## Moran I test under randomisation
##
## data: residuals(lm.sf.bs)
## weights: bos.listw
##
## Moran I statistic standard deviate = -2.2402, p-value = 0.9875
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## -0.0623654367 -0.0019801980 0.0007265842With a high p-value, we can finally fail to reject the null hypothesis. This means we can accept that we have accounted for the spatial autocorrelation. After running these tests, I believe that this is an adequate model.