Forming Spatial Weights Matrix with Large Datasets (>150, 000 observations)

Dear List

I am keen to receive advice from listers on the best strategies in R for 
forming and using a spatial weights matrix for a large (ca. 150,000) 
dataset. Is it fastest and most memory-efficient to do the entire 
analysis in R, or is it better to use GeoDa to form a gal or gwt file 
and then import that into R and then procedue with commands from spdep?

The background is that I am using R.2.12.0 (64 bit) on a 16-Gb dual core 
Windows XP setup to run spatial diagnostics (Moran and LM tests) on OLS 
residuals and estimate spatial error and spatial lag models on an 
environmental dataset, with 150,000 observations (1 km pixels). To start 
with more manageable samples I have grid-sampled pixels at the 
intersection of every 4th, 8th and 16th row and column, giving samples 
of ca. 9000, 2300 and 600. To ensure no "islands" when using the 
smallest samples, the neighbourhood for the weights matrix is set as 34 
km, and the weights are inverse-distance. In case anyone asks "why are 
you even bothering to do this" it is to examine the impact of the 
spatial sampling, since we have previously based our analysis (of 
impacts of roads on forest cover) on a 1-in-25 sample and now want to 
see how sensitive inferences might be to this sort of spatial sampling.

Back to my current problem: while R breezes through the small samples 
(less than a minute for a sample of 600, about 20 minutes for a sample 
of 2300) it has been running for 3 days now on a sample of 9000 
observations, using around about 50% of the CPU and about 5.6 Gb of the 
16 Gb available on the machine. I have access to a cluster computer with 
128-CPUs to run the analysis for the full population of 150,000 
observations, but based on the slow performance with 9000 on a 16-Gb 
machine I am not optimistic about the ability of the cluster computer to 
run the data at full 1km scale.
If you ran the steps separately, you would see where the bottleneck is. My 
guess is that method="eigen" in the first spatial regression is the 
underlying problem, and that replacing this with sparse matrix techniques 
will resolve the problem: method="Matrix". See also comments inline in 
code below, as you seem to create dense matrices when sparse 
representations are all you need.

Roger
Hopefully it won't colour any answers I receive, but I am primarily a 
Stata user, and with 64-bit multi-processor stata using our own update 
of stata's spatwmat program (updated to take advantage of Mata - the 
matrix programming language that was unavailable when spatwmat was first 
written) for the sample of 9000 it takes about 2 hours to form the 
spatial weights, do the Moran and LM tests on OLS and run the spatial 
lag and spatial error models. So rather faster than we are managing in R 
(but probably due to our novice R programming, so I am certainly not 
complaining about R here).  But Stata is not a candidate for running the 
analysis on the 128-cluster computer because of licencing issues, while 
R can be installed on each CPU with no licensing issues. So I need to 
know the most memory-efficient and fastest way to form and use a spatial 
weights matrix for ca. 150,000 observations.

If people haven't completely lost interest in the question by now, I am 
pasting a fragment of the code we are using, noting that typically we 
have found the closest command to what we knew how to do in Stata, and 
so this may not be the most efficient way in R of forming the weights 
matrix and then using it. So I would be grateful for pointers, before I 
end up burning through a lot of computing time on the cluster computer 
with an inefficient set up..

Thanks,

John Gibson

#  qui nspatwmat, name(W_34km_ds) xcoord(x) ycoord(y) band(0 34000) standardize eigenval(E_34km_ds)

coords <- cbind(rawd$x,rawd$y)
W_34km_nb <- dnearneigh(coords, 0, 34000, row.names = NULL, longlat = NULL)
What does:

print(W_34km_nb)

report as the average number of neighbours? Is it too large for your 
hypothesised processes?
# Eigen value matrix for the weight matrix - slightly different from Stata results
#E_34km_ds <- eigenw(W_34km_dsls, quiet=NULL)
No, unnecessary, created on the fly, and is dense.
# Get the Moran stats for -fa00-, call it A
# nspatcorr fa00, bands(0 34000 34000) xcoord(x) ycoord(y) cumul

You may save repeated calls to nb2listw() by doing it once and storing the 
listw object for use:

lw <- nb2listw(W_34km_nb,style = "W")
fa00.moran <- moran.test(rawd$fa00, lw)

and so on ...
fa00.moran <- moran.test(rawd$fa00,nb2listw(W_34km_nb,style = "W"))
print(fa00.moran)

# Moran stats for -uhat-, B

No, use lm.morantest() passing the lm object:

uhat.moran <- moran.test(rawd.lm, lw)
uhat.moran <- moran.test(rawd.lm$residuals,nb2listw(W_34km_nb,style = "W"))
print(uhat.moran)

# Weight matrix, inverse distance, standardise
temp <- nbdists(W_34km_nb, coords, longlat = NULL)
W_34km_ivds <-lapply(temp, function(x) 1/x)
Why create the dense matrix?

W_34km_dsls <- nb2listw(W_34km_nb, glist=W_34km_ivds, style="W", ...
W_34km_dsmat <- nb2mat(W_34km_nb, glist=W_34km_ivds, style="W", zero.policy=TRUE)

#W_34km_dsmat[1:5,1:5] - Same as the Stata results

W_34km_dsls <- mat2listw(W_34km_dsmat)

# Spatial diagnose, weight matrix (W_34km_ds)
rawd.lm.diag<-lm.LMtests(rawd.lm, W_34km_dsls, test=c("LMerr", "LMlag", "RLMerr","RLMlag"))
print(rawd.lm.diag)

rm(coords,W_34km_nb, temp,W_34km_ivds,W_34km_dsmat)
rm(rawd.lm)
As stated above, don't use method="eigen" for larger N, use 
method="Matrix", or another exact sparse method, or an approximation, such 
as Monte Carlo or Chebyshev. Running N=25000 on a 1GB laptop isn't a 
problem, with sparse matrix techniques or approximations, everything 
becomes easier.

Note also that you need to use impacts() on the output of lagsarlm() 
because the coefficients should not be interpreted like OLS coefficients - 
see the references in ?impacts.
# spatial lag, weights(W_34km_ds), same RHS as the LM
 rawd.lag <- lagsarlm(fa00 ~ 
d2roadf+d2nr_ld51_00f+s_ld10_00_10+protectf+non_adjacent+d2pvcapf+slo+dem+pa+ta+soil_n+soil_ap+soil_ph+s_pop_00_10+s_gdp_00_10, 
data=rawd, W_34km_dsls, method="eigen")

_______________________________________________
R-sig-Geo mailing list
R-sig-Geo at stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo

Roger Bivand
Economic Geography Section, Department of Economics, Norwegian School of
Economics and Business Administration, Helleveien 30, N-5045 Bergen,
Norway. voice: +47 55 95 93 55; fax +47 55 95 95 43
e-mail: Roger.Bivand at nhh.no

Forming Spatial Weights Matrix with Large Datasets (>150, 000 observations)

Thread (2 messages)