Converting nb list containing areas with no neighbours with nb2WB()
On Wed, 25 Nov 2009, Xingli Giam wrote:
Hi Virgilio, Thank you for your clear and very quick reply. I really appreciate your help in looking at my data.
The problem is caused by the nb object encoding no-neighbour status with
an integer vector of length 1 containing 0. Please try this patched
version:
nb2WB <- function(nb)
{
if (!inherits(nb, "nb")) stop("not a neighbours list")
num <- card(nb)
if (any(num == 0)) nb[num == 0] <- NULL
adj <- unlist(nb)
weights <- rep(1, sum(num))
list(adj=adj, weights=weights, num=num)
}
It works for me so far, and num now contains the correct 0 for observation
233. I haven't tried it with WinBUGS though.
Hope this helps,
Roger
Here is my data and code, I hope it will be clear enough for you:
rm(list=ls())
library(spdep)
library(R2WinBUGS)
S <- c(66.666667, 40.000000, 50.000000, 56.666667, 250.000000,
66.666667, 183.333333, 50.000000, 7.000000, 106.666667, 90.000000,
83.333333, 100.000000,
5.500000, 146.666667, 100.000000, 160.000000, 70.000000, 160.000000,
183.333333, 70.000000, 133.333333, 130.000000, 73.333333, 93.333333,
31.666667,
183.333333, 116.666667, 116.666667, 200.000000, 110.000000, 120.000000,
50.000000, 70.000000, 133.333333, 133.333333, 106.666667, 66.666667,
106.666667,
133.333333, 9.333333, 90.000000, 33.333333, 66.666667, 200.000000,
106.666667, 90.000000, 31.666667, 110.000000, 50.000000, 100.000000,
120.000000,
110.000000, 136.666667, 31.666667, 93.333333, 70.000000, 110.000000,
60.900000, 333.333333, 133.333333, 83.333333, 83.333333, 100.000000,
36.666667,
66.666667, 86.666667, 7.333333, 53.333333, 50.000000, 96.666667,
100.000000, 60.000000, 33.333333, 60.000000, 183.333333, 133.333333,
116.666667,
126.666667, 50.000000, 166.666667, 50.000000, 12.666667, 100.000000,
50.000000, 66.666667, 143.333333, 83.333333, 126.666667, 73.333333,
45.966667,
76.666667, 93.333333, 100.000000, 233.333333, 50.000000, 76.666667,
116.666667, 93.333333, 46.666667, 90.000000, 106.666667, 36.666667,
233.333333,
26.666667, 13.333333, 216.666667, 76.666667, 110.000000, 133.333333,
46.666667, 66.666667, 73.333333, 46.666667, 53.333333, 250.000000,
116.666667,
66.666667, 53.333333, 26.666667, 133.333333, 26.666667, 23.333333,
70.000000, 66.666667, 116.666667, 100.000000, 66.666667, 90.000000,
100.000000,
150.000000, 50.000000, 166.666667, 200.000000, 200.000000, 50.000000,
300.000000, 63.333333, 150.000000, 11.666667, 136.666667, 116.666667,
126.666667,
83.333333, 300.000000, 10.666667, 123.333333, 183.333333, 110.000000,
116.666667, 216.666667, 110.000000, 2.666667, 266.666667, 266.666667,
83.333333,
126.666667, 106.666667, 266.666667, 216.666667, 93.333333, 266.666667,
166.666667, 50.000000, 200.000000, 183.333333, 150.000000, 116.666667,
106.666667,
133.333333, 150.000000, 233.333333, 216.666667, 33.333333, 216.666667,
106.666667, 56.666667, 63.333333, 23.333333, 216.666667, 100.000000,
83.333333,
216.666667, 216.666667, 60.000000, 166.666667, 166.666667, 100.000000,
53.333333, 300.000000, 91.666667, 136.666667, 233.333333, 34.466667,
116.666667,
283.333333, 283.333333, 166.666667, 96.666667, 150.000000, 70.000000,
4.000000, 216.666667, 135.000000, 150.000000, 200.000000, 90.000000,
176.666667,
53.333333, 166.666667, 53.333333, 60.000000, 38.333333, 1.666667,
4.333333, 3.666667, 48.333333, 31.166667, 3.766667, 10.600000,
16.666667,
5.833333, 1.000000, 25.833333, 20.766667, 16.000000, 4.666667,
5.000000, 1.666667, 116.666667, 200.000000)
A <- c(
233.222222, 832.888889, 312.222222, 954.000000, 19025.555556,
2971.666667, 1827.777778, 267.444444, 1.555556, 178.666667,
3879.555556,
1342.333333, 1613.666667, 4.666667, 14959.555556, 2576.000000,
3621.777778, 2141.777778, 3966.666667, 8559.333333, 11054.777778,
13593.111111,
12867.777778, 8385.777778, 464.444444, 1461.777778, 2427.888889,
8562.222222, 11489.333333, 20980.111111, 4208.777778, 45875.777778,
228.555556,
2292.111111, 5760.000000, 7244.333333, 2617.333333, 10257.222222,
20988.333333, 27526.000000, 34.444444, 3436.000000, 340.111111,
1975.444444,
12417.777778, 22108.000000, 3349.555556, 548.000000, 126.777778,
1593.666667, 7449.222222, 59007.444444, 12461.222222, 48021.111111,
114.666667,
16273.666667, 14228.888889, 22691.777778, 632.888889, 47287.000000,
12786.444444, 7463.555556, 6290.333333, 4897.111111, 4317.555556,
2259.666667,
3290.777778, 14.888889, 37796.777778, 1758.777778, 5962.888889,
3872.888889, 4236.111111, 1678.555556, 15323.333333, 73588.000000,
13239.777778,
28167.666667, 7953.888889, 919.222222, 10541.888889, 3925.666667,
60.222222, 4625.444444, 718.888889, 2013.333333, 11629.666667,
1119.888889,
15027.222222, 7370.222222, 187.555556, 5338.222222, 3425.000000,
5228.111111, 48429.888889, 1866.000000, 4667.777778, 5974.888889,
2502.444444,
2469.777778, 1586.444444, 1899.666667, 401.111111, 13882.666667,
1191.555556, 3.444444, 24858.222222, 2630.666667, 2505.000000,
5146.000000,
4067.000000, 1388.888889, 340.666667, 258.444444, 1999.444444,
28698.666667, 8068.666667, 9680.000000, 8467.222222, 1613.888889,
10770.000000,
2880.666667, 3251.333333, 29182.555556, 2907.888889, 4609.000000,
1721.000000, 451.555556, 285.555556, 3702.444444, 23963.666667,
872.222222,
12145.000000, 25471.222222, 10025.444444, 530.666667, 20373.111111,
2528.111111, 3546.111111, 10.555556, 9938.888889, 1472.333333,
639.888889,
230.777778, 8145.000000, 2.777778, 1586.888889, 7508.555556,
1211.111111, 2371.000000, 11340.222222, 336.777778, 2.111111,
16218.111111,
52904.000000, 1097.888889, 5090.444444, 12722.888889, 6484.333333,
3241.777778, 918.888889, 29827.111111, 26832.555556, 109.555556,
79631.333333,
11622.000000, 8493.333333, 9811.666667, 15737.222222, 45812.666667,
7389.555556, 27843.444444, 22325.777778, 1112.333333, 8978.444444,
841.888889,
3114.222222, 1905.888889, 854.444444, 53653.222222, 1944.666667,
2510.777778, 20662.333333, 16510.555556, 833.777778, 19640.000000,
19251.222222,
10665.111111, 529.555556, 11623.333333, 432.222222, 1246.444444,
18547.888889, 230.111111, 8350.000000, 82961.444444, 1808.222222,
37233.222222,
5630.555556, 21418.444444, 524.666667, 1.000000, 52227.555556,
12715.888889, 3252.111111, 7660.666667, 549.111111, 3772.666667,
223.777778,
29452.444444, 7720.222222, 25537.333333, 64.222222, 68.000000,
23.555556, 59.000000, 1285.777778, 746.444444, 3.777778,
119.333333,
10.666667, 51.111111, 19.777778, 342.555556, 180.000000,
104.222222, 104.777778, 15.777778, 10.333333, 29903.555556,
26650.444444)
x_long <- c(
146.97, 131.53, 135.88, 126.63, 140.71, 128.14, 146.65, 136.08,
159.09, 153.49, 151.11, 151.42, 165.57, 167.95, 143.25, 139.42,
145.73, 129.35, 155.24,
148.32, 138.58, 140.60, 121.05, 120.12, 151.03, 166.70, 133.12,
134.16, 29.68, 11.29, 10.20, 21.04, 44.25, 7.41, 8.59,
35.30, 35.83, 21.10,
-2.71, 38.95, 55.46, -12.15, 23.37, 30.48, 48.59, 47.03,
32.74, 55.55, 9.16, 6.47, 4.66, 24.46, 39.59, 14.95,
6.54, 40.02, 17.02,
-9.93, 92.81, 112.82, 113.89, 111.30, 94.31, 102.97, 100.65,
99.71, 93.96, 105.69, 81.42, 113.36, 112.43, 122.98, 82.27,
95.66, 96.06, 114.19,
98.13, 87.67, 101.55, 122.25, 122.08, 74.89, 120.18, 91.87,
98.94, 124.72, 125.83, 120.94, 92.82, 96.22, 93.32, 73.53,
73.95, 105.21, 100.91,
102.69, 100.95, 96.88, 105.81, 85.78, 119.01, 101.41, 103.31,
102.28, 105.82, 115.83, 108.12, 77.03, 76.73, 107.36, 116.59,
80.19, 80.83, 120.00,
101.54, 101.05, 100.22, 102.84, 116.12, 89.33, 99.11, 104.90,
105.10, 80.75, 108.34, 109.41, 109.65, 127.98, 120.62, 121.03,
-51.31, -50.19, -39.70,
-41.85, -66.43, -40.89, -72.04, -72.66, -76.01, -81.70, -84.67,
-89.61, -92.80, -94.19, -76.97, -87.06, -66.22, -72.12, -84.96,
-75.44, -78.26, -77.38,
-32.44, -63.22, -56.99, -53.74, -70.09, -74.28, -83.65, -82.28,
-77.28, -65.21, -66.22, -61.67, -60.03, -74.50, -75.10, -50.51,
-43.91, -52.68, -61.65,
-75.73, -66.26, -42.98, -77.11, -96.74, -61.10, -92.28, -55.58,
-52.62, -35.72, -36.34, -73.49, -90.33, -66.34, -69.45, -63.36,
-68.78, -73.89, -48.74,
-94.94, -92.28, -70.32, -80.39, -64.96, -69.17, -82.77, -54.90,
-65.46, -47.11, -61.21, -30.32, -56.78, -76.19, -70.41, -98.21,
-98.39, -79.52, -61.03,
-50.30, -89.02, -63.97, 158.22, -157.66, -159.77, 173.02, 179.31,
-155.31, -177.93, -138.98, 142.15, 134.56, -109.34, -172.48, -149.38,
-175.20, -146.36, -144.36,
-171.63, 107.97, 102.87)
y_lat <- c(
-2.10, -7.56, -0.92, -3.49, -5.21, 0.65, -6.08, -1.72, -31.56,
-11.48, -5.22, -5.48, -21.31, -29.04, -4.28, -2.79, -17.47, -3.33,
-5.94, -8.72, -7.48,
-6.23, -1.97, -2.42, -9.95, -15.14, -1.04, -2.46, -1.93, -0.75,
5.64, -1.17, -12.18, 5.53, 5.07, 0.00, -8.26, -0.30, 7.61,
13.13, -4.67, 11.23,
-33.98, -30.38, -19.27, -18.60, -26.36, -21.11, 4.23, 5.17, 6.78,
0.96, -3.94, 1.01, 0.13, -13.60, 0.43, 7.18, 12.53, 1.49,
2.07, 2.04, 26.93,
11.69, 15.11, 14.32, 21.29, -10.45, 19.93, -8.02, -7.59, 10.02,
28.07, 16.88, 20.26, 26.03, 17.76, 24.04, 18.06, 16.82, 16.61,
13.19, 13.77, 25.79,
-1.36, 7.89, 7.69, 12.68, 22.26, 18.72, 8.00, 17.17, 16.39,
17.75, 23.05, 17.75, 18.58, 26.34, 18.34, 20.23, 10.07, 4.62,
3.11, 4.08, 20.86,
10.72, 22.15, 10.54, 10.63, 15.23, 0.42, 6.78, 6.98, 5.20,
0.33, -0.13, -0.09, -0.57, 0.05, 22.68, 11.11, 11.68, 10.34,
27.65, -7.10, -6.95,
19.56, 26.48, 22.38, 23.89, -26.65, -30.05, -16.45, -16.90, -15.69,
-3.96, 0.64, 8.80, 4.84, 12.54, 14.05, 15.18, 17.09, 16.65,
5.51, 5.53, 9.99,
6.26, 10.43, 20.42, -2.04, 8.34, -3.86, 4.00, 4.89, -2.00,
19.09, -5.35, 10.22, 8.52, 18.18, -0.54, -5.88, 16.16, -8.90,
5.38, 7.75, -1.12,
-4.26, -10.02, -6.00, -1.52, 1.89, -2.75, 2.51, 17.82, 9.06,
18.16, 5.72, -23.25, -8.73, -8.79, -11.89, 17.10, 18.24, -3.99,
-6.58, 0.81, 10.74,
-25.75, 18.34, 15.15, -2.28, 25.66, -23.59, -10.48, 9.11, -6.16,
5.63, -3.04, 10.69, -20.50, -0.20, -7.38, 8.80, 21.79, 20.64,
-0.53, 14.66, -5.78,
19.67, 3.79, 6.88, 1.88, -21.23, 1.83, -16.58, 19.77, -29.27,
-9.78, 26.66, 7.53, -27.13, -13.62, -17.69, -21.17, -16.18, -27.61,
-2.82, 27.93, 26.11)
coords <- cbind(x_long, y_lat)
nb2000 <- dnearneigh(coords, 0, 2000, longlat=T)
nb2000weights <- nb2WB(nb2000)
# find zero integer value in nb2000weights adj matrix
zero.int <- (1:length(nb2000weights$adj))[nb2000weights$adj==0]
zero.int
# testing if nb2000weights adj and weights matrix is of the same length (#
Yes)
length(nb2000weights$adj)
length(nb2000weights$weights)
N=length(A)
d <- list(sp = S, area_km2 = A, N=N, adj=nb2000weights$adj,
weights=nb2000weights$weights, num=nb2000weights$num)
init1 = list(c=runif(1,0.5,1.5)*12, z=runif(1,0.5,1.5)*0.26,
prec=runif(1,0.01,100), tau=1, b=c(rep(0, N)))
init2 = list(c=runif(1,0.5,1.5)*12, z=runif(1,0.5,1.5)*0.26,
prec=runif(1,0.01,100), tau=1, b=c(rep(0, N)))
init3 = list(c=runif(1,0.5,1.5)*12, z=runif(1,0.5,1.5)*0.26,
prec=runif(1,0.01,100), tau=1, b=c(rep(0, N)))
inits = list(init1,init2,init3)
powerdat1lg.bugs <- bugs(data=d, inits,
model.file="C:/WinBUGS/CARlogcodes/powerbuglnorm.txt",
parameters = c("c", "z", "sp.rep", "b"),
n.chains =3, n.iter = 32500, n.burnin = 7500, n.sims=3000,
bugs.directory ="c:/Program Files/WinBUGS14/", debug=T)
____________________________________________
Here is my WinBUGS model file code
"C:/WinBUGS/CARlogcodes/powerbuglnorm.txt":
model
{
for (i in 1:N) # for each ecoreg in biome 1
{
S[i] <- c*(pow(area_km2[i] , z)) + b[i]
logS[i] <- log(max(0.001,S[i]))
sp[i] ~ dlnorm(logS[i],prec)
sp.rep [i] <- max(0.001, (c*pow(area_km2[i],z) + b[i]))
}
# CAR prior distribution for random effects
b[1:N] ~ car.normal(adj[], weights[], num[], tau)
c ~ dnorm(0,0.0000001)I(0,)
z ~ dnorm(0,0.0000001)I(0,)
prec ~ dgamma(0.001,0.001)
tau ~ dgamma(0.001,0.001)
}
Many thanks,
Xingli
-----Original Message-----
From: Virgilio G?mez-Rubio [mailto:Virgilio.Gomez at uclm.es]
Sent: Wednesday, 25 November, 2009 3:27 AM
To: Xingli Giam
Cc: R-sig-Geo at stat.math.ethz.ch
Subject: Re: [R-sig-Geo] Converting nb list containing areas with no
neighbours with nb2WB()
Hi,
I checked the adjacency matrix of the spatial weights list output object,
and a zero integer was placed at the position corresponding to the area
with
no neighbours.
Yes. Probably we need to fix this so that nothing is added when an area
has no neighbours.
When I re-ran the analysis after I changed my neighbourhood structure such
that all areas had neighbours, the WinBUGS simulation ran fine.
Is there a way to format the spatial weights list object via nb2WB() in a
way that allows WinBUGS to read cases with no neighbours?
It should be done automatically, but probably we forgot to consider the
case of areas with no neighbours.
Could it be possible to see the map that you are using? I'd say that
there is a mismatch between the number of neighbours and the number of
weights, but I would need to check with you data and code.
Best wishes,
Virgilio
_______________________________________________
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