[R-sig-Geo] Converting nb list containing areas with no neighbours with nb2WB()
Xingli Giam
xgiam at princeton.edu
Thu Nov 26 06:35:44 CET 2009
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
More information about the R-sig-Geo
mailing list