Hi everyone,
Thank you in advance for your advice. First I must apologize for sending the
two scrubbed messages in this thread. I thought MS Outlook Web Light sends
messages in plain text format but I was wrong. I have transferred my account
to my new email address and have set the email format to plain text. I hope
everything is fine now. I sincerely apologize for causing inconvenience by
sending multiple unreadable messages to the list.
I am trying to convert a nb object containing a single area with no
neighbours into a spatial weights list for input into WinBUGS using nb2WB().
When I tried to run the CAR models in WinBUGS using the spatial weights list
via R2WinBUGS, I encountered a WinBUGS trap. The error message was "Index
out of range".
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.
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?
In WinBUGS, if we were to enter the adjacency matrix manually it will be:
#Here is where we tell BUGS how to set up the adjacency structure. For
example, site 1 neighbors sites 19,9, and 5.
adj = c(
19, 9, 5,
10, 7,
12,
28, 20, 18,
19, 12, 1,
#Site 6 has no neighbors
17, 16, 13, 10, 2,
and so on.
)
Here are my OS details:
R version 2.10.0 (2009-10-26)
i386-pc-mingw32
locale:
[1] LC_COLLATE=English_Singapore.1252 LC_CTYPE=English_Singapore.1252
LC_MONETARY=English_Singapore.1252
[4] LC_NUMERIC=C LC_TIME=English_Singapore.1252
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] R2WinBUGS_2.1-16 spdep_0.4-52 nlme_3.1-96 coda_0.13-4
MASS_7.3-3
[6] Matrix_0.999375-32 lattice_0.17-26 spam_0.15-5 boot_1.2-41
maptools_0.7-26
[11] foreign_0.8-38 sp_0.9-47 deldir_0.0-10
loaded via a namespace (and not attached):
[1] grid_2.10.0 tools_2.10.0
Best,
Xingli Giam
Ph.D. candidate
Department of Ecology and Evolutionary Biology
Princeton University
Princeton, USA
Converting nb list containing areas with no neighbours with nb2WB()
6 messages · Xingli Giam, Virgilio Gomez Rubio, Roger Bivand
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
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
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
Hi Roger, your patch is excellent! Thank you for your help!
Your patch works fine with WinBUGS and I attempted to adapt it to the
listw2WB() function. It works well so far.
listw2WB <- function (listw)
{
if (!inherits(listw, "listw"))
stop("not listw class object")
num <- card(listw$neighbours)
if (any(num == 0)) listw$neighbours[num == 0] <- NULL
adj <- unlist(listw$neighbours)
weights <- unlist(listw$weights)
list(adj = adj, weights = weights, num = num)
}
Thanks, Roger and Virgilio, for your timely advice!
Best regards,
Xingli
-----Original Message-----
From: Roger Bivand [mailto:Roger.Bivand at nhh.no]
Sent: Wednesday, 25 November, 2009 6:20 AM
To: Xingli Giam
Cc: 'Virgilio G?mez-Rubio'; R-sig-Geo at stat.math.ethz.ch
Subject: Re: [R-sig-Geo] 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
On Thu, 26 Nov 2009, Xingli Giam wrote:
Hi Roger, your patch is excellent! Thank you for your help!
Thanks, new release on its way to CRAN. Roger
Your patch works fine with WinBUGS and I attempted to adapt it to the
listw2WB() function. It works well so far.
listw2WB <- function (listw)
{
if (!inherits(listw, "listw"))
stop("not listw class object")
num <- card(listw$neighbours)
if (any(num == 0)) listw$neighbours[num == 0] <- NULL
adj <- unlist(listw$neighbours)
weights <- unlist(listw$weights)
list(adj = adj, weights = weights, num = num)
}
Thanks, Roger and Virgilio, for your timely advice!
Best regards,
Xingli
-----Original Message-----
From: Roger Bivand [mailto:Roger.Bivand at nhh.no]
Sent: Wednesday, 25 November, 2009 6:20 AM
To: Xingli Giam
Cc: 'Virgilio G?mez-Rubio'; R-sig-Geo at stat.math.ethz.ch
Subject: Re: [R-sig-Geo] 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