[R-sig-Geo] Converting nb list containing areas with no neighbours with nb2WB()

Roger Bivand Roger.Bivand at nhh.no
Wed Nov 25 12:19:33 CET 2009


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