Skip to content

Car.proper C[] matrix

6 messages · Jason Gasper, Uwe Ligges, Roger Bivand +1 more

#
I am hoping someone can help translate some WinBUGS code into R code. I 
would like to use R to create the C[] matrix required for a car.proper 
model in WinBUGS, but I am having a difficult time negotiating the 
coding. The C matrix provides normalized weights for each pair of 
spatial areas.   So the WinBUGS example is as follows:

# of the weight matrix with elements Cij. The first J1 elements of the 
C[] vector contain the
# weights for the J1 neighbours of area i=1; the (J1+1) to J2 elements 
of the C[] vector contain
# the weights for the J2 neighbours of area i=2; etc.
# To set up this vector, we need to define a variable cumsum, which 
gives the values of J1,
# J2, etc.; we then set up an index matrix pick[,] with N columns 
corresponding to the
# i=1,...,N areas, and with the same number of rows as there are 
elements in the C[] vector
# (i.e. sumNumNeigh). The elements C[ (cumsum[i]+1):cumsum[i+1] ] 
correspond to
# the set of weights Cij associated with area i, and so we set up ith 
column of the matrix pick[,]
# to have a 1 in all the rows k for which cumsum[i] < k <= cumsum[i+1], 
and 0's elsewhere.
# For example, let N=4 and cumsum=c(0,3,5,6,8), so area i=1 has 3 
neighbours, area i=2 has 2
# neighbours, area i=3 has 1 neighbour and area i=4 has 2 neighbours. 
The the matrix pick[,] is:
# pick
# 1, 0, 0, 0,
# 1, 0, 0, 0,
# 1, 0, 0, 0,
# 0, 1, 0, 0,
# 0, 1, 0, 0,
# 0, 0, 1, 0,
# 0, 0, 0, 1,
# 0, 0, 0, 1,
#
# We can then use the inner product (inprod(,)) function in WinBUGS and 
the kth row of pick to
# select which area corresponds to the kth element in the vector C[]; 
likewise, we can use inprod(,)
# and the ith column of pick to select the elements of C[] which 
correspond to area i.

Basically I want to do this in R to speed things up a little. Has anyone 
written a function for this conversion?
#
If you want that people help to translate *code*, you have to specify it ...

Uwe Ligges
Jason Gasper wrote:
#
Here is the WinBUGS code

model {
for(i in 1:N) {m[i] <- 1/n[ind[i]] }
cumsum[1] <- 0
for(i in 2:(N+1)) {cumsum[i] <- sum(num[1:(i-1)]) }
for(k in 1:sumNumNeigh) {
  for(i in 1:N) {
# #pick[k,i] = 1 if cumsum[i] < k <= cumsum[i=1]; otherwise, pick[k,i] = 0
##step(e) 1 if e >= 0; 0 otherwise
    pick[k,i]<-step(k-cumsum[i]-epsilon)*step(cumsum[i+1]-k) }

  C[k]<-1/ inprod(num[], pick[k,]) }
epsilon <- 0.0001
Uwe Ligges wrote:

  
    
#
Jason Gasper <Jason.Gasper <at> noaa.gov> writes:
Please look at the nb2WB() function in the spdep package, I think that you'll 
find that it provides what you need.

Roger Bivand
#
I have been using the nb2WB() package for the car.normal function in WinBUGS,
but it will not create the C[] matrix; it only creates adj[], num[], and
weights[]. I was planning on using this function to create the C[] matrix
(by using the num matix) required for the car.proper, but I got slipped up
on the R coding. My data set is large enough that WinBUGS chokes on the code
I provided below.
Roger Bivand wrote:

  
    
#
jgasper <Jason.Gasper <at> noaa.gov> writes:
OK, I understand. Then I think that in the first case on:

http://mathstat.helsinki.fi/openbugs/Manuals/GeoBUGS/Manual.html#Proper

you would use style="W" in nb2listw(my_nb) - row standardisation - and 
rename weights to C, taking M as 1/card(my_nb), and use the glist= 
argument to nb2listw() to set up the second case. The weights[] will then 
be the C[] sparse object indexed by adj[]. Something like this is 
equivalent to example 2:

glist <- vector(mode="list", length=n)
for (i in seq(along=my_nb) glist[[i]] <- sqrt(E[my_nb[[i]]]/E[i])
lw_obj <- nb2listw(my_nb, glist=glist, style="B")

run listw2WB(), and rename weights to C, taking M as 1/E.

Perhaps R-sig-geo would be a more appropriate list - maybe someone there 
has done this already?

Roger Bivand