Thank you for the Geographical Analysis reference, very helpful. I have
deduced that the method argument should be "Matrix" for my data rather
than "LU", since the "W" weights are symmetric. However, I am still
getting the warnings on non convergence and since you state that this
maybe a specification problem, I would like to find the error. I did
eliminate all polygons without neighbours for all my datasets before
calling errorsarlm, so this should not be the problem. I include a full
code example inline below (including documentation of the intermediate
objects and the R outputs) for one dataset with warning and one without
below in case you have time to look at it and find a possible error
source that I have not identified.
comments on crash also inline
Am 30.09.2013 12:31, schrieb Roger Bivand:
On Mon, 30 Sep 2013, Philine Gaffron wrote:
Thank you for the response, Roger. My answers below. I have also
added a question on a warning message I have received and would be
glad to find out whether this means I need to change the analysis
approach or just have to take it as additional information.
Thanks for reporting back: some points for clarification inline.
Am 25.09.2013 00:53, schrieb Roger Bivand:
On Wed, 25 Sep 2013, Philine Gaffron wrote:
Hello,
I am trying to run a spatial regression analysis using the errorsarlm
function in the sdpep package. I first ran into memory problems (on a
larger dataset than the one described below) but have adjusted the
memory size limit:
memory.limit(size = 800000)
[1] 8e+05
Which sets your memory limit under Windows to 781 Gb (size is in Mb) -
that probably explains the error exit unless you have a lot of memory.
Best avoid making any changes from the default in a fresh session
unless you know exactly what you are doing.
I did up the virtual memory of the machine to 800Gb first (it has 1Tb
in total and there is currently nothing else running on it) so I hope
that that was not the problem. Thanks for the hint, though.
However, when I try to execute the following code in R Studio, I
receive the message "R Session aborted - R encountered a fatal error
- The session was terminated"
Always rerun any code with problems outside R Studio, because it can
interfere with the reporting of warnings and errors. Certainly
errorsarlm() has no reason to cause R to crash - and this is what you
are suggesting. If you can reproduce this outside R Studio, and
without setting memory limits, I may need to see your data to try to
reproduce any error exit from the R session.
I tried this and had no more crashes - especially after I made the
changes suggested/outlined below
Does this mean that the R error exit occurred when you used R Studio,
and did not occur when R was run directly, using exactly the same data
and code? If so, please also let the R Studio people know, so that
they can continue debugging.
I will try to reproduce this once I have understood / solved my
non-convergence problem. If and when I do, I will send the code and data
to the RStudio people and cc you (?).
ElDU.slm.PM25.white <- errorsarlm(PM_PM25 ~ race_white, data =
ElDorado_urban, listw=ElDoU_listw)
Did you read the advice on ?errorsarlm about the use of alternative
methods, such as "LU" or "Matrix" for larger data sets? You are using
the default "eigen" method, but this should certainly have worked for
the modest n=2081. Can you reproduce the problem with a built-in data
set, such as elect80?
I opted for the "LU" method and changed the weights style (see
below). Is there any documentation that tells me the implications of
using method ="LU"? I only found the note in the spdep documentation
"LU_order default FALSE; used in ?LU_prepermutate?, note warnings
given for lu method" but nothing further. I did receive warnings in
the output, though (for several but not all of my datasets) and I am
not sure what they mean / what their implications are for the
results. Do I need to choose another method? The warnings in every
case were:
See Bivand, Hauke & Kossowski (2013) Geographical Analysis for some
details; LU_prepermutate emulates the implementation by finding an
initial fill-in re-ordering, but in R has little speed-up effect,
hence it is turned off by default.
1: In powerWeights(W = W, rho = lambda, order = con$pWOrder, X = B,:
not converged within order iterations
2: In powerWeights(W = t(W), rho = lambda, order = con$pWOrder, X = C,:
not converged within order iterations
These are unusual, and may signal a mis-specified setting. The default
number of powers is 250, and lambda^250 is a very small number unless
lambda is close to its (probably upper) limit. The computation in fact
only affects the accuracy of the Hausman test on the difference
between OLS and spatial error coefficients, so may be adjusted if
lambda is close to limits by setting control argument pWOrder to a
larger value than 250, see ?powerWeights. As a symptom of lambda close
to its limits, it suggests mis-specification, possibly in terms of
observation support, with neighbouring properties sharing important
drivers. This may be suggested by your data scenario described below;
scrutinising the intermediate objects in your script may also put you
on the track of things behaving in unexpected ways. I would expect
many observations to have no neighbours with a 150m threshold, but you
don't show no-neighbour control being turned off.
I removed all polygons without neighbours prior to calling errorsarlm -
see example below.
Hope the helps,
Roger
Here is the worked example of the two different datasets. No warning for
the first one, warnings received for the second one included (sorry,
this is a bit lengthy but I did not want to leave out anything that
might be relevant).
####READING IN ArcGIS SHAPEFILE###############################
Parcels2008_inc_POP <- readOGR(".", "Parcels_2008_inc_pop")
Parcels2008_UNinc_POP <- readOGR(".", "Parcels_2008_UNinc_pop")
########EXTRACT SUBSETS OF Parcels2008 SpatialPolygonsDataframe####
ElDorado_rural <- subset(Parcels2008_UNinc_POP, COUNTYFP10 == "017")
ElDorado_urban <- subset(Parcels2008_inc_POP, COUNTYFP10 == "017")
##############TO CREATE A DISTANCE BASED NEIGHBOURS LIST IN
R########################
# FIRST, calculate coords of feature centroids####
Eld_R_coords <- coordinates(ElDorado_rural)
Eld_U_coords <- coordinates(ElDorado_urban)
##SECOND, create vector with row names ####
Eld_R_IDs <- row.names(ElDorado_rural)
Eld_U_IDs <- row.names(ElDorado_urban)
###THIRD, creating the nb objects for calculating spatial weights####
##DATASET ONE: NO WARNING FROM errorsarlm (see below)####
ElDoR_nb <- dnearneigh(Eld_R_coords, d1 = 0, d2 = 984, row.names =
Eld_R_IDs)
ElDoR_nb
# Neighbour list object:
# Number of regions: 26863
# Number of nonzero links: 1591926
# Percentage nonzero weights: 0.2206042
# Average number of links: 59.26092
# 34 regions with no links:
# 0 83 872 993 2301 2420 7518 12579 14649 17136 17148 18526 18601 18786
20413
# 20427 20662 21445 21824 22497 23309 23421 23535 24366 24743 24784
25007 25009
# 25107 25184 25203 25229 25421 25465
#remove parcels/regions with no links
ElDorado_rural <- subset(ElDorado_rural, subset=card(ElDoR_nb) > 0)
Eld_R_coords <- coordinates(ElDorado_rural)
summary(Eld_R_coords)
# V1 V2
# Min. :6816265 Min. :1953502
# 1st Qu.:6845579 1st Qu.:2007471
# Median :6876165 Median :2016977
# Mean :6884229 Mean :2021703
# 3rd Qu.:6918361 3rd Qu.:2034467
# Max. :7091505 Max. :2108970
Eld_R_IDs <- row.names(ElDorado_rural)
summary(Eld_R_IDs)
# Length Class Mode
# 26829 character character
ElDoR_nb <- dnearneigh(Eld_R_coords, d1 = 0, d2 = 984, row.names =
Eld_R_IDs)
ElDoR_nb
# Neighbour list object:
# Number of regions: 26829
# Number of nonzero links: 1591926
# Percentage nonzero weights: 0.2211637
# Average number of links: 59.33602
is.symmetric.nb(ElDoR_nb)
# [1] TRUE
##DATASET TWO:NON-CONVERGENCE WARNING FROM errorsarlm (see below)########
ElDoU_nb <- dnearneigh(Eld_U_coords, d1 = 0, d2 = 492, row.names =
Eld_U_IDs)
ElDoU_nb
# Neighbour list object:
# Number of regions: 2094
# Number of nonzero links: 57172
# Percentage nonzero weights: 1.303857
# Average number of links: 27.30277
# 13 regions with no links:
# 425 468 628 929 1216 1218 1224 1252 1363 1416 1927 2007 2088
#remove parcels/regions with no links
ElDorado_urban <- subset(ElDorado_urban, subset=card(ElDoU_nb) > 0)
Eld_U_coords <- coordinates(ElDorado_urban)
summary(Eld_U_coords)
# V1 V2
# Min. :6892551 Min. :2024721
# 1st Qu.:6901757 1st Qu.:2028604
# Median :6904635 Median :2030283
# Mean :6905218 Mean :2030021
# 3rd Qu.:6908098 3rd Qu.:2031299
# Max. :6916815 Max. :2035333
Eld_U_IDs <- row.names(ElDorado_urban)
summary(Eld_U_IDs)
# Length Class Mode
# 2081 character characte
ElDoU_nb <- dnearneigh(Eld_U_coords, d1 = 0, d2 = 492, row.names =
Eld_U_IDs)
ElDoU_nb
# Neighbour list object:
# Number of regions: 2081
# Number of nonzero links: 57172
# Percentage nonzero weights: 1.320198
# Average number of links: 27.47333
is.symmetric.nb(ElDoU_nb)
# [1] TRUE
####WEIGHTS######################################
#ElDorado rural
ElDoR_listw <- nb2listw(ElDoR_nb, glist=NULL, style="W", zero.policy=NULL)
names(ElDoR_listw)
# [1] "style" "neighbours" "weights"
summary(ElDoR_listw)
# Characteristics of weights list object:
# Neighbour list object:
# Number of regions: 26829
# Number of nonzero links: 1591926
# Percentage nonzero weights: 0.2211637
# Average number of links: 59.33602
# Link number distribution:
#
# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
# 85 160 267 454 692 870 1007 1016 940 802 751 583 541 457 377 347 301
261 209 221 234 193
# 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44
# 191 207 216 183 185 185 161 184 175 209 189 163 178 155 164 138 138
151 134 135 116 125
# 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66
# 96 97 106 101 107 116 100 122 110 103 125 103 105 105 107 112 77 101
95 95 102 102
# 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88
# 69 102 110 87 90 96 135 94 107 90 111 98 99 96 135 125 98 118 110 123
94 76
# 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108
109 110
# 101 113 95 105 84 99 89 102 94 98 114 107 91 104 116 90 98 97 110 93
103 97
# 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127
128 129 130 131 132
# 107 98 114 111 93 100 77 103 101 95 109 96 79 103 99 80 78 86 83 93 84 80
# 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149
150 151 152 153 154
# 72 90 91 93 82 81 82 71 77 70 84 67 57 74 69 87 74 74 77 52 76 70
# 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171
172 173 174 175 176
# 49 60 52 42 62 57 56 47 56 45 59 58 53 53 61 47 32 26 39 37 35 39
# 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193
194 195 196 197 198
# 30 28 32 17 26 22 16 14 16 18 7 13 9 12 15 13 18 6 7 18 12 11
# 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215
216 217 218 219 220
# 13 13 13 16 13 10 17 12 12 11 8 9 1 4 8 6 2 7 6 8 5 3
# 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237
238 239 240 241 242
# 8 7 4 7 3 10 1 2 4 2 6 9 1 6 6 9 7 4 7 6 1 5
# 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259
261 262 263 264 265
# 4 2 5 3 5 4 3 2 5 4 1 2 5 2 4 5 4 2 3 4 2 3
# 266 267 268 269 270 272 273 274 275 276 278 279 280 282 283 284 285
286 287 289 290 292
# 2 2 2 4 6 2 1 3 1 1 4 1 2 6 4 8 6 5 2 3 1 3
# 293 295 297
# 7 2 1
# 85 least connected regions:
# 37 261 262 263 268 334 585 1011 1012 1035 2419 2546 2621 2659 3258
4544 4557 5501 5502 5859 10005 10006 10007 12577 12578 16301 17256 17257
18697 18788 19339 19445 20499 20658 21104 21105 21440 21441 21494 21829
21981 21982 22039 22056 22154 22208 22310 22315 22316 22396 22397 22398
22423 22816 22823 22989 22990 23121 23133 23554 23621 23771 23772 23823
23891 23903 24313 24348 24437 24464 24676 24802 24803 24886 25008 25306
25354 25466 25634 26060 26432 26433 26452 26453 26509 with 1 link
# 1 most connected region:
# 9783 with 297 links
#
# Weights style: W
# Weights constants summary:
# n nn S0 S1 S2
# W 26829 719795241 26829 3194.683 108332.1
summary(sapply(ElDoR_listw$weights, sum))
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 1 1 1 1 1 1
#El Dorado urban
ElDoU_listw <- nb2listw(ElDoU_nb, glist=NULL, style="W", zero.policy=NULL)
summary(ElDoU_listw)
# Characteristics of weights list object:
# Neighbour list object:
# Number of regions: 2081
# Number of nonzero links: 57172
# Percentage nonzero weights: 1.320198
# Average number of links: 27.47333
# Link number distribution:
#
# 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26
27 28 29 30 31 32 33 34 35 36 37
# 14 15 29 29 21 27 34 36 35 32 43 58 46 44 49 53 60 40 56 42 56 49 50
58 54 46 48 40 49 45 45 38 31 38 37 26 29
# 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
61 62 63 64 65 66
# 32 37 33 38 32 32 33 31 34 34 43 28 26 15 17 17 18 16 14 5 14 9 3 4 5
3 2 3 1
# 14 least connected regions:
# 147 148 159 378 466 467 596 629 642 1217 1705 1710 2074 2087 with 1 link
# 1 most connected region:
# 865 with 66 links
#
# Weights style: W
# Weights constants summary:
# n nn S0 S1 S2
# W 2081 4330561 2081 260.5586 8405.665
summary(sapply(ElDoU_listw$weights, sum))
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 1 1 1 1 1 1
####SPATIAL REGRESSION ANALYSIS#####
##El Dorado rural################################
ElDR.slm.PM25.white <- errorsarlm(PM_PM25 ~ race_white, data =
ElDorado_rural, listw=ElDoR_listw, method = "Matrix")
ElDR.PM25.white <- summary(ElDR.slm.PM25.white, Nagelkerke = TRUE,
digits = 4, signif.stars = TRUE)
ElDR.PM25.white
# Call:errorsarlm(formula = PM_PM25 ~ race_white, data = ElDorado_rural,
# listw = ElDoR_listw, method = "Matrix")
#
# Residuals:
# Min 1Q Median 3Q Max
# -0.37948413 -0.00236403 -0.00096381 0.00088396 0.69478272
#
# Type: error
# Coefficients: (asymptotic standard errors)
# Estimate Std. Error z value Pr(>|z|)
# (Intercept) 3.0372e-02 3.0068e-03 10.1013 < 2.2e-16
# race_white -5.7658e-05 1.4887e-05 -3.8731 0.0001074
#
# Lambda: 0.94999, LR test value: 32485, p-value: < 2.22e-16
# Approximate (numerical Hessian) standard error: 0.0022514
# z-value: 421.96, p-value: < 2.22e-16
# Wald statistic: 178050, p-value: < 2.22e-16
#
# Log likelihood: 62497.7 for error model
# ML residual variance (sigma squared): 0.00049093, (sigma: 0.022157)
# Nagelkerke pseudo-R-squared: 0.70418
# Number of observations: 26829
# Number of parameters estimated: 4
# AIC: -124990, (AIC for lm: -92505)
##El Dorado urban####################################
ElDU.slm.PM25.white <- errorsarlm(PM_PM25 ~ race_white, data =
ElDorado_urban, listw=ElDoU_listw, method = "Matrix")
# Warning messages:
# 1: In powerWeights(W = W, rho = lambda, order = con$pWOrder, X = B, :
# not converged within order iterations
# 2: In powerWeights(W = t(W), rho = lambda, order = con$pWOrder, X = C, :
# not converged within order iterations
ElDU.PM25.white <- summary(ElDU.slm.PM25.white, Nagelkerke = TRUE,
digits = 4, signif.stars = TRUE)
ElDU.PM25.white
# Call:errorsarlm(formula = PM_PM25 ~ race_white, data = ElDorado_urban,
# listw = ElDoU_listw, method = "Matrix")
#
# Residuals:
# Min 1Q Median 3Q Max
# -0.1277026 -0.0080154 -0.0023227 0.0049338 0.3288472
#
# Type: error
# Coefficients: (asymptotic standard errors)
# Estimate Std. Error z value Pr(>|z|)
# (Intercept) 5.3598e-02 1.5882e-02 3.3747 0.0007389
# race_white -2.8711e-05 8.8431e-05 -0.3247 0.7454294
#
# Lambda: 0.95004, LR test value: 2875.6, p-value: < 2.22e-16
# Approximate (numerical Hessian) standard error: 0.0070736
# z-value: 134.31, p-value: < 2.22e-16
# Wald statistic: 18039, p-value: < 2.22e-16
#
# Log likelihood: 4093.031 for error model
# ML residual variance (sigma squared): 0.00099886, (sigma: 0.031605)
# Nagelkerke pseudo-R-squared: 0.74927
# Number of observations: 2081
# Number of parameters estimated: 4
# AIC: -8178.1, (AIC for lm: -5304.5)
###Session Info######
# R Studio version 0.97.551
sessionInfo()
# R version 3.0.1 (2013-05-16)
# Platform: x86_64-w64-mingw32/x64 (64-bit)
#
# locale:
# [1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United
States.1252
# [3] LC_MONETARY=English_United States.1252 LC_NUMERIC=C
# [5] LC_TIME=English_United States.1252
#
# attached base packages:
# [1] stats graphics grDevices utils datasets methods base
#
# other attached packages:
# [1] spdep_0.5-65 Matrix_1.0-12 lattice_0.20-15 rgdal_0.8-11 sp_1.0-13
#
# loaded via a namespace (and not attached):
# [1] boot_1.3-9 coda_0.16-1 deldir_0.0-22 grid_3.0.1 LearnBayes_2.12
MASS_7.3-26
# [7] nlme_3.1-109 splines_3.0.1 tools_3.0.1
I have not tried reproducing the error with the built-in dataset but
if that is of interest to you from the point of view of spdep
development, I could give it a try.
Is 492 a well-chosen upper distance threshold? I assume that the input
data are projected, so that these are measures in the units of the
projection - metres? If Eld_U_coords are actually geographical
coordinates, but longlat=NULL, d2 is taken as a planar measure, if you
declare longlat=TRUE, it is in km. There is no easy way to trap
undeclared geographical coordinates.
The distance thresholds are in feet
(NAD_1983_StatePlane_California_II_FIPS_0402_Feet) and are related to
dispersion distances for road traffic emissions. For explanation: my
polygons are land use parcels for which I have calculated emission
loads in grams/area from traffic loads on a road network. Since I am
assuming a dispersion distance of 492 feet / 150 m of emissions away
from the network in the example above, I am further assuming, that
parcels which are within that distance of one another are affected by
the same emission sources (i.e. traffic loads on the network).
Does choosing a different style make a difference ("S" is not often
used)?
I chose "S" as a variance stabilising coding schemes following
Bivand, Pebesma et al. (2008), p. 253 but have now switched back to
"W". Apart from the warnings cited above I encountered no more issues.
Thank you for your help and for making this package available.
Philine
Roger
Is there anything inherent in my approach that could cause this
problem or would one need to look at the data? I include the session
info below.
Thanks for any help,
Philine
sessionInfo()
R version 3.0.1 (2013-05-16)
Platform: x86_64-w64-mingw32/x64 (64-bit)
locale:
[1] LC_COLLATE=English_United States.1252 LC_CTYPE=English_United
States.1252 LC_MONETARY=English_United States.1252 LC_NUMERIC=C
[5] LC_TIME=English_United States.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] spdep_0.5-65 Matrix_1.0-12 lattice_0.20-15 rgdal_0.8-11 sp_1.0-13
loaded via a namespace (and not attached):
[1] boot_1.3-9 coda_0.16-1 deldir_0.0-22 grid_3.0.1 LearnBayes_2.12
MASS_7.3-26 nlme_3.1-109 splines_3.0.1 tools_3.0.1
library(help=spdep)
Information on package ?spdep?
Description:
Package: spdep
Version: 0.5-65
Date: 2013-09-04