efficient code/function for rectangular SP weight Matrix and gwr

On Tue, 22 May 2007, Sam Field wrote:

  

In the event that others are following this thread,  I have included 
some r-code for constructing a rectangular weight matrix from two sets 
of coordinates (in feet). The first set of coords (nhood_pts[,1] and 
nhood_pts[,2]) come from a regular grid, the second set refer to the 
location of patients ( combo2$X and combo2$Y).  I don't provide the data 
only the code below.  The original data consists of 1000 grids and 
30000+ patients. It takes 30 minutes on my fairly speedy computer to 
complete the first step- extracting pair-wise distances less then a 
given bandwidth (D). I store the neighbor ids and counts in separate 
lists.   The rest of the code calculates kernel weights based on these 
distances. The code is not efficient from any number of standpoints, but 
storing these distances in a list rather then a matrix avoided running 
into memory allocation problems. Here is the code.
    

Would spDistsN1() in sp have helped to speed things up a little? Looping
over the 1k grid points and passing out the 30k patient points shouldn't
be very demanding (34 seconds on a oldish box, more if you need Great
Circle distances). For heavyweight things, quadtrees would be a better
choice, I have an unpublished package for approximate nearest neighbours.

Roger

  

# define bandwidth
D <- 2*5280


#This creates a list of distances < 1 mile
w_raw <- vector("list",length(nhood_pts[,1]))
for (i in 1:length(nhood_pts[,1])){
w_raw[[i]] <- rep(NA,length(combo2$X))}

system.time(for (i in 1:length(nhood_pts[,1])){for(j in 1:length(combo2$X)){
if(sqrt((nhood_pts[,1][i]-combo2$X[j])^2 + 
(nhood_pts[,2][i]-combo2$Y[j])^2) < D) w_raw[[i]][j] <-  
sqrt((nhood_pts[,1][i]-combo2$X[j])^2 + 
(nhood_pts[,2][i]-combo2$Y[j])^2)}} )



w_raw.id <-   vector("list",length(nhood_pts[,1]))
for (i in 1:length(nhood_pts[,1])){
w_raw.id[[i]] <-    which(!is.na(w_raw[[i]]))}

for (i in 1:length(nhood_pts[,1])){
w_raw[[i]] <-  w_raw[[i]][w_raw.id[[i]]]}


w_raw.count <-   rep(1,length(nhood_pts[,1]))
for (i in 1:length(nhood_pts[,1])){
w_raw.count[i] <-  length(w_raw[[i]])}

w_weights <-  vector("list",length(nhood_pts[,1]))
for (i in 1:length(nhood_pts[,1])){
w_weights[[i]] <- (1-(w_raw[[i]]/D)^2)^2}

w <-  vector("list",length(nhood_pts[,1]))
for (i in 1:length(nhood_pts[,1])){
w[[i]] <- w_weights[[i]]/sum(w_weights[[i]])}





Roger Bivand wrote:
    

On Fri, 11 May 2007, St?phane Dray wrote:

  
      

Hi Sam,

I think that this question is quite general and could interest other 
people, including me, with very different aims. I have developed a 
method to look for the relationships between two data sets that have 
been sampled on the same area but for different locations. In my 
example, the two samples are two polygons layers. In this approach, I 
compute a rectangular weighting matrix where each weight correspond to 
the area of intersection between polygons of each layer. I have used 
also the matrix form to store these weights (my data set was very small 
compared to you). I remember that Roger was also interested by these 
rectangular weights in another context. Here we  have different  problems:
- how to compute these kind  of weights
- how to store them.

For the first point, I think that for each method/application, the 
solution  is different. We could develop/extend classical tools for 
square weights (one set of spatial units) to rectangular weights (two 
sets of spatial units).
For the second one, It would be probably interesting to define a class 
of object in spdep. nb objects are lists, and I think that it would be 
the solution for rectangular neighborhood.

If I consider two sets of spatial units (A and B) where the number of 
units is equal to na and nb.  We could store the neighbors in a list of 
length 2. The first element of this list is a list of length na. In this 
list, the j-th element is a vector of the neighbors of the j-th unit of 
the layer A. These neighbors are spatial units of the layer B.  The 
second element of the global list is a list of length nb where each 
element is a vector of neighbors.

I think that we have to think to a class of object that could be useful 
for everybody dealing with this kind of rectangular weights. If this 
class is properly defined (second point), we could then develop tools to 
construct this kind of neighborhoods (first point). The eventual 
extension to more than two data sets could also be taken into account in 
this reflexion.

    
        

I would welcome input on this. I'm looking at an alternative weights 
representation through classes in the Matrix package, which is evolving 
fast, and which seems to be promising. If the dimnames slot is used to 
hold the region.id values, it might be possible to make progress.

Best wishes,

Roger

  
      

Cheers,


Sam Field wrote:
    
        

List,

I need to create a rectangular spatial weight matrix for a set of n and 
m objects. I quickly run in to memory allocation problems when 
constructing the full matrix in a single pass. I am looking for a more 
efficient way of doing this. There appears to be efficient procedures in 
spdep for constructing SQUARE spatial weight matrices (e.g. 
dnearneigh()). Are there analogous procedures for constructing distance 
based weights between two different point patterns? I am doing this in 
preparation for implementing an approximate geographically weighted 
logistic regression procedure. I was thinking about using re sampling 
procedure as an inferential frame- perhaps I might get some feedback. 
This is what I was going to do.

I have a point pattern of 30,000 diabetic people based on where they 
lived during a 2 year period. During that period, approximately 4% of 
them developed diabetes. I am interested in isolating the impact of 
ecological factors on the geographic variation" of the disease, so it is 
necessary to control for the spatial clustering of individual level risk 
factors associated with the disease (diabetes).

Step 1: Estimate a logistic regression using the full sample and predict 
incidence diabetes using individual level covariates (i.e. who developed 
diabetes over the two year period).

Step 2. Estimate a weighted logit model at each location (grid). The 
observations would be the people (not the geographic units) and the 
weights would be kernel weights based on distance. The model would only 
contain a single freely estimated parameter, the intercept, but it would 
also contain an offset term. For each patient, the offset term would 
simply be an evaluation of the linear predictor of the global model 
estimated above (based on the observed covariate values), but without 
the intercept. This would effectively fix the estimates of the patient 
level coefficients to their global values, requiring only a local 
estimate of the intercept. My hope is that I could interpret geographic 
variability in the intercept as evidence for a "location effect" net of 
the patient composition or "risk profile" at a particular location. It 
would probably make sense to center the X variables so that the 
intercept was interpretable and estimated in a region of the response 
plane where their is plenty of data. I would let the other covariates 
vary as well, but I doubt the model could be estimated in large portions 
of the study area because of sparse data.

Step 3. If I were going to do inference on the location specific 
intercepts, I would generate a sampling distribution at each location by 
re sampling from the global model, and repeat Step 2 for each randomly 
drawn sample. This would give me a local sampling distribution of 
intercept estimates at each location and I could compare it to the the 
single one generated from the observed data. The global model represents 
a kind of null because the intercept is fixed to its global value and 
geographic variability is driven entirely by the spatial clustering of 
patient level factors.


thanks!

Sam