[R-sig-Geo] Converting nb list containing areas with no neighbours with nb2WB()
Xingli Giam
xgiam at princeton.edu
Wed Nov 25 11:11:04 CET 2009
Hi Virgilio,
Thank you for your clear and very quick reply. I really appreciate your help
in looking at my data.
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
More information about the R-sig-Geo
mailing list