[R-sig-Geo] Problem with moran.mc() in spdep package
Roger Bivand
Roger.Bivand at nhh.no
Wed Jan 7 10:29:27 CET 2015
On Wed, 7 Jan 2015, Roger Bivand wrote:
> Please provide the xy and VAR objects in a form that permits reproduction,
> for example using dput(), and/or post only plain text. The present rendering
> in most email clients does not permit copy&paste.
>
> Roger
Oh, well (copying & pasting from the list archive rendering) - comments
below:
VAR <- c(0.39, 0.38, 0.39, 0.28, 0.32, 0.29, 0.46, 0.32, 0.33, 0.29,
0.33, 0.31, 0.35, 0.36, 0.39, 0.48, 0.42, 0.47, 0.34, 0.38, 0.42,
0.42, 0.71, 0.71, 0.57, 0.55, 0.56, 0.6, 0.59, 0.58, 0.65, 0.66,
0.6, 0.58, 0.57, 0.57, 0.56, 0.65, 0.68, 0.69, 0.67, 0.66, 0.63,
0.6, 0.61, 0.63, 0.62, 0.63, 0.64, 0.64, 0.63, 0.61, 0.67, 0.66,
0.68, 0.68, 0.62, 0.65, 0.68, 0.66, 0.66, 0.6, 0.54, 0.58, 0.58,
0.59, 0.58, 0.22, 0.38, 0.38, 0.35, 0.35, 0.26, 0.33, 0.45, 0.44,
0.24, 0.29, 0.37, 0.34, 0.4, 0.47, 0.5, 0.52, 0.51, 0.42, 0.44,
0.6, 0.46, 0.36, 0.32, 0.58, 0.42, 0.27, 0.24, 0.41, 0.3, 0.43,
0.46, 0.48, 0.52, 0.54, 0.34, 0.6, 0.55, 0.63, 0.59, 0.58, 0.41,
0.29, 0.33, 0.5, 0.49, 0.29, 0.32, 0.4, 0.5, 0.56, 0.42, 0.48,
0.27, 0.27, 0.33, 0.32, 0.38, 0.36, 0.37, 0.44, 0.47, 0.33, 0.34,
0.37, 0.34, 0.41, 0.53, 0.47, 0.41, 0.5, 0.35, 0.41, 0.38, 0.42,
0.36, 0.35, 0.35, 0.36, 0.28, 0.26, 0.26, 0.3, 0.32, 0.23, 0.32,
0.45, 0.35, 0.31, 0.48, 0.38, 0.5, 0.42, 0.31, 0.38, 0.42, 0.53,
0.37, 0.44, 0.49, 0.49, 0.48, 0.3, 0.6, 0.29, 0.32, 0.3, 0.33,
0.26, 0.33, 0.43, 0.37, 0.29, 0.27, 0.32, 0.27, 0.32, 0.36, 0.37,
0.32, 0.36, 0.3, 0.39, 0.34, 0.29)
# the matrix
xy <- structure(c(381739, 383080, 396833, 397637, 379144, 379347,
379488,
379543, 380018, 380064, 380085, 380138, 380277, 380334, 380404,
380968, 381037, 399015, 399029, 398980, 399013, 399083, 399209,
399285, 399351, 399389, 399427, 399485, 399526, 399545, 399590,
399545, 399528, 399350, 399348, 399300, 399222, 399127, 399119,
399087, 399063, 399024, 399154, 399193, 399187, 399185, 399173,
399124, 399082, 399122, 399081, 399003, 398982, 398946, 398900,
398893, 398860, 398908, 398897, 398908, 398941, 398924, 398972,
399043, 399034, 399015, 399048, 402002, 402120, 401403, 399127,
398471, 398196, 397346, 396927, 396472, 394941, 411075, 410840,
410705, 410602, 410321, 410364, 410233, 410230, 410209, 410217,
410221, 410266, 410364, 410444, 410401, 410431, 410513, 410480,
410414, 410347, 410385, 410417, 410430, 410465, 410367, 410259,
410268, 410334, 410328, 410344, 410341, 410398, 410408, 410352,
410498, 410315, 410233, 410183, 410210, 410215, 410246, 410481,
410494, 390915, 395462, 378616, 378898, 378972, 379005, 379024,
379055, 379080, 379122, 379199, 379254, 379404, 379402, 404070,
403921, 403531, 403027, 414950, 414867, 414685, 414587, 414454,
414011, 413896, 413804, 413597, 413513, 413474, 413414, 413473,
413217, 413099, 412908, 412837, 410570, 378883, 379124, 379330,
384327, 384735, 385000, 385116, 387763, 385189, 385298, 385419,
387738, 387935, 393172, 394798, 377709, 377633, 377478, 377469,
377486, 377498, 377523, 377541, 377562, 377592, 377620, 377664,
377683, 377737, 377773, 377838, 377883, 378099, 378135, 378208,
378516, 6925207, 6923900, 6904017, 6903938, 6899689, 6898704,
6897820, 6897673, 6895093, 6894715, 6894596, 6894417, 6893664,
6893358, 6892984, 6891698, 6891385, 6908553, 6908485, 6908457,
6908424, 6908510, 6908558, 6908540, 6908537, 6908522, 6908562,
6908557, 6908533, 6908484, 6908587, 6908630, 6908688, 6908649,
6908583, 6908675, 6908720, 6908546, 6908592, 6908651, 6908703,
6908824, 6908684, 6908731, 6908775, 6908831, 6908882, 6908884,
6908907, 6908955, 6908994, 6909005, 6908881, 6908964, 6909071,
6909110, 6909215, 6909170, 6909087, 6908999, 6908918, 6909204,
6909196, 6909165, 6909115, 6909066, 6909028, 6918317, 6918254,
6918655, 6919937, 6920385, 6920528, 6921056, 6921354, 6921611,
6922552, 6910046, 6909714, 6909495, 6909328, 6910006, 6910111,
6910101, 6910064, 6909673, 6909630, 6909439, 6909368, 6909345,
6909339, 6909237, 6909517, 6909536, 6909604, 6909616, 6909791,
6909902, 6909985, 6910059, 6910096, 6910151, 6909876, 6910196,
6909415, 6909504, 6909554, 6909591, 6909603, 6909702, 6909858,
6910075, 6910067, 6909747, 6909774, 6909922, 6909968, 6910150,
6909847, 6909943, 6910899, 6906497, 6912776, 6912329, 6911976,
6911662, 6911381, 6911210, 6911059, 6910730, 6910494, 6910306,
6909713, 6909620, 6917160, 6917244, 6917468, 6917746, 6909363,
6909406, 6909513, 6909567, 6909647, 6910024, 6910181, 6910330,
6910584, 6910708, 6910822, 6910898, 6911012, 6911214, 6911321,
6911614, 6911690, 6913180, 6927706, 6927573, 6927499, 6922144,
6921488, 6921083, 6920895, 6916349, 6920789, 6920632, 6920453,
6916581, 6915727, 6903415, 6903553, 6918153, 6917985, 6917619,
6917449, 6917261, 6917044, 6916877, 6916742, 6916548, 6916181,
6915971, 6915724, 6915572, 6915110, 6914943, 6914747, 6914650,
6914112, 6913978, 6913765, 6912930), .Dim = c(192L, 2L))
The underlying problem is the weights object, with n=192:
table(card(lista$neighbours))
# 0 1 2
# 153 32 7
With this number of no-neighbour observations, the measure is pretty
meaningless. Taking the weights constants:
unlist(spweights.constants(lista, zero.policy=TRUE, adjust.n=TRUE))
# n n1 n2 n3 nn S0 S1 S2
# 39 38 37 36 1521 39 69 162
unlist(spweights.constants(lista, zero.policy=TRUE, adjust.n=FALSE))
# n n1 n2 n3 nn S0 S1 S2
# 192 191 190 189 36864 39 69 162
we see that S0 (sum of weights) is constant, but n changes, so:
moran(VAR, lista, n=192, S0=39, zero.policy=TRUE)$I
# [1] 1.259683
moran(VAR, lista, n=39, S0=39, zero.policy=TRUE)$I
# [1] 0.2558731
as S0 is part of the denominator and n the numerator of the measure. The
adjustment of n is done to prevent people shooting themselves in the foot
by specifying improbable weights objects.
It is, however, true that moran.mc() could take the adjust.n= argument; it
currently takes adjust.n=FALSE by definition, change committed.
Hope this clarifies,
Roger
>
> On Wed, 7 Jan 2015, Leandro Roser wrote:
>
>> Hello everyone. I'm new in this group. I will appreciate any
>> suggestion for solving the following problem.
>> I am using the package spdep, with the purpose computing bootstraped
>> values of Moran's I over a range of geographical distances. The
>> weighing matrix for each range was obtained with the package adegenet:
>>
>> library(adegenet)
>> library(spdep)
>>
>>
>
>>
>> xy<-matrix(xy,192,2)
>>
>>
>>
>> lista<-chooseCN(xy,type=5,d1=5,d2=50, result.type = "listw",plot.nb = F)
>> #listw object obtained with adegenet
>>
>>
>> moran.mc(VAR,lista,999,na.action = na.exclude, zero.policy = T)
>>
>> # Monte-Carlo simulation of Moran's I
>>
>> #data: VAR
>> #weights: lista
>> #number of simulations + 1: 1000
>> #
>> #*statistic = 1.2597*, observed rank = 1000, p-value = 0.001
>> #alternative hypothesis: greater
>>
>>
>> #How can I intepret this result (Moran >1)? Extracting the values in the
>> object "lista", I obtained:
>>
>> print(lista,zero.policy=TRUE)
>> #n nn S0 S1 S2
>> #W 39 1521 39 69 162
>>
>> #And using these parameters in the function moran():
>>
>> moran(VAR,lista,39,39, zero.policy =T)
>> #$I
>> #[1] *0.2558731*
>> #
>> #$K
>> #[1] 1.911139
>>
>> #The result of moran() is different of the obtained with moran.mc (and is
>> not higher than 1). Might be a problem in the moran.mc() function?
>>
>> Thanks in advance,
>>
>> Leandro.
>>
>>
>>
>>
>>
>>
>>
>> --
>> Lic. Leandro Gabriel Roser
>> Laboratorio de Genética
>> Dto. de Ecología, Genética y Evolución,
>> F.C.E.N., U.B.A.,
>> Ciudad Universitaria, PB II, 4to piso,
>> Nuñez, Cdad. Autónoma de Buenos Aires,
>> Argentina.
>> tel ++54 +11 4576-3300 (ext 219)
>> fax ++54 +11 4576-3384
>>
>> [[alternative HTML version deleted]]
>>
>> _______________________________________________
>> R-sig-Geo mailing list
>> R-sig-Geo at r-project.org
>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>
>
--
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 91 00
e-mail: Roger.Bivand at nhh.no
More information about the R-sig-Geo
mailing list