[R-sig-Geo] Converting nb list containing areas with no neighbours with nb2WB()
Roger Bivand
Roger.Bivand at nhh.no
Thu Nov 26 08:45:13 CET 2009
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
More information about the R-sig-Geo
mailing list