Roger,
thanks very much for your reply!
I tried to update spdep and downloaded the new version. Installation worked
fine and I'm now working with spdep, version 0.5-11, 2010-05-31.
The problem is, however, that I now can no longer use the read.dat2listw
function, which I use to obtain network information from an external file.
When I run
Network <- read.dat2listw("C:/Network.txt")
I get the error message:
Error in `[.data.frame`(sn, , 3) : undefined columns selected
I therefore had to move back to spdep, version 0.4-56, 2009-12-14.
I also used moran.test under debug() as you suggested but if I understood
the output correctly, the error message only comes up at the end.
I have attached the full debug report below.
lag.listw(Network,x) seems to work fine although the values are very, very
small:
y<-lag.listw(Network,x)
summary(y)
Min. 1st Qu. Median Mean 3rd Qu. Max.
1.497e-310 1.200e-300 1.698e-300 1.094e-226 2.414e-300 1.994e-223
This is a bit surprising to me as the values of x are several orders of
magnitude larger:
summary(x)
Min. 1st Qu. Median Mean 3rd Qu. Max.
1082 1543 2119 2250 2886 3965
Thanks,
Michael
> debug(moran.test)
> moran.test(x,Network)
debugging in: moran.test(x, Network)
debug: {
alternative <- match.arg(alternative, c("greater", "less",
"two.sided"))
if (!inherits(listw, "listw"))
stop(paste(deparse(substitute(listw)), "is not a listw object"))
if (!is.numeric(x))
stop(paste(deparse(substitute(x)), "is not a numeric vector"))
if (is.null(spChk))
spChk <- get.spChkOption()
if (spChk && !chkIDs(x, listw))
stop("Check of data and weights ID integrity failed")
xname <- deparse(substitute(x))
wname <- deparse(substitute(listw))
NAOK <- deparse(substitute(na.action)) == "na.pass"
x <- na.action(x)
na.act <- attr(x, "na.action")
if (!is.null(na.act)) {
subset <- !(1:length(listw$neighbours) %in% na.act)
listw <- subset(listw, subset, zero.policy = zero.policy)
}
n <- length(listw$neighbours)
if (n != length(x))
stop("objects of different length")
wc <- spweights.constants(listw, zero.policy = zero.policy,
adjust.n = adjust.n)
S02 <- wc$S0 * wc$S0
res <- moran(x, listw, wc$n, wc$S0, zero.policy = zero.policy,
NAOK = NAOK)
I <- res$I
K <- res$K
if (rank)
K <- (3 * (3 * wc$n^2 - 7))/(5 * (wc$n^2 - 1))
EI <- (-1)/wc$n1
if (randomisation) {
VI <- wc$n * (wc$S1 * (wc$nn - 3 * wc$n + 3) - wc$n *
wc$S2 + 3 * S02)
tmp <- K * (wc$S1 * (wc$nn - wc$n) - 2 * wc$n * wc$S2 +
6 * S02)
VI <- (VI - tmp)/(wc$n1 * wc$n2 * wc$n3 * S02)
VI <- VI - EI^2
}
else {
VI <- (wc$nn * wc$S1 - wc$n * wc$S2 + 3 * S02)/(S02 *
(wc$nn - 1))
VI <- VI - EI^2
}
ZI <- (I - EI)/sqrt(VI)
statistic <- ZI
names(statistic) <- "Moran I statistic standard deviate"
if (alternative == "two.sided")
PrI <- 2 * pnorm(abs(ZI), lower.tail = FALSE)
else if (alternative == "greater")
PrI <- pnorm(ZI, lower.tail = FALSE)
else PrI <- pnorm(ZI)
if (!is.finite(PrI) || PrI < 0 || PrI > 1)
warning("Out-of-range p-value: reconsider test arguments")
vec <- c(I, EI, VI)
names(vec) <- c("Moran I statistic", "Expectation", "Variance")
method <- paste("Moran's I test under", ifelse(randomisation,
"randomisation", "normality"))
data.name <- paste(xname, ifelse(rank, "using rank correction",
""), "\nweights:", wname, ifelse(is.null(na.act), "",
paste("\nomitted:", paste(na.act, collapse = ", "))),
"\n")
res <- list(statistic = statistic, p.value = PrI, estimate = vec,
alternative = alternative, method = method, data.name = data.name)
if (!is.null(na.act))
attr(res, "na.action") <- na.act
class(res) <- "htest"
res
}
Browse[1]>
debug: alternative <- match.arg(alternative, c("greater", "less",
"two.sided"))
Browse[1]>
debug: if (!inherits(listw, "listw"))
stop(paste(deparse(substitute(listw)), "is not a listw object"))
Browse[1]>
debug: if (!is.numeric(x)) stop(paste(deparse(substitute(x)), "is not a
numeric vector"))
Browse[1]>
debug: if (is.null(spChk)) spChk <- get.spChkOption()
Browse[1]>
debug: if (spChk && !chkIDs(x, listw)) stop("Check of data and weights ID
integrity failed")
Browse[1]>
debug: xname <- deparse(substitute(x))
Browse[1]>
debug: wname <- deparse(substitute(listw))
Browse[1]>
debug: NAOK <- deparse(substitute(na.action)) == "na.pass"
Browse[1]>
debug: x <- na.action(x)
Browse[1]>
debug: na.act <- attr(x, "na.action")
Browse[1]>
debug: if (!is.null(na.act)) {
subset <- !(1:length(listw$neighbours) %in% na.act)
listw <- subset(listw, subset, zero.policy = zero.policy)
}
Browse[1]>
debug: n <- length(listw$neighbours)
Browse[1]>
debug: if (n != length(x)) stop("objects of different length")
Browse[1]>
debug: wc <- spweights.constants(listw, zero.policy = zero.policy, adjust.n
= adjust.n)
Browse[1]>
debug: S02 <- wc$S0 * wc$S0
Browse[1]>
debug: res <- moran(x, listw, wc$n, wc$S0, zero.policy = zero.policy, NAOK =
NAOK)
Browse[1]>
debug: I <- res$I
Browse[1]>
debug: K <- res$K
Browse[1]>
debug: if (rank) K <- (3 * (3 * wc$n^2 - 7))/(5 * (wc$n^2 - 1))
Browse[1]>
debug: EI <- (-1)/wc$n1
Browse[1]>
debug: if (randomisation) {
VI <- wc$n * (wc$S1 * (wc$nn - 3 * wc$n + 3) - wc$n * wc$S2 +
3 * S02)
tmp <- K * (wc$S1 * (wc$nn - wc$n) - 2 * wc$n * wc$S2 + 6 *
S02)
VI <- (VI - tmp)/(wc$n1 * wc$n2 * wc$n3 * S02)
VI <- VI - EI^2
} else {
VI <- (wc$nn * wc$S1 - wc$n * wc$S2 + 3 * S02)/(S02 * (wc$nn -
1))
VI <- VI - EI^2
}
Browse[1]>
debug: VI <- wc$n * (wc$S1 * (wc$nn - 3 * wc$n + 3) - wc$n * wc$S2 + 3 *
S02)
Browse[1]>
debug: tmp <- K * (wc$S1 * (wc$nn - wc$n) - 2 * wc$n * wc$S2 + 6 * S02)
Browse[1]>
debug: VI <- (VI - tmp)/(wc$n1 * wc$n2 * wc$n3 * S02)
Browse[1]>
debug: VI <- VI - EI^2
Browse[1]>
debug: ZI <- (I - EI)/sqrt(VI)
Browse[1]>
debug: statistic <- ZI
Browse[1]>
debug: names(statistic) <- "Moran I statistic standard deviate"
Browse[1]>
debug: if (alternative == "two.sided") PrI <- 2 * pnorm(abs(ZI), lower.tail
= FALSE) else if (alternative ==
"greater") PrI <- pnorm(ZI, lower.tail = FALSE) else PrI <- pnorm(ZI)
Browse[1]>
debug: if (!is.finite(PrI) || PrI < 0 || PrI > 1) warning("Out-of-range
p-value: reconsider test arguments")
Browse[1]>
debug: vec <- c(I, EI, VI)
Browse[1]>
debug: names(vec) <- c("Moran I statistic", "Expectation", "Variance")
Browse[1]>
debug: method <- paste("Moran's I test under",
ifelse(randomisation, "randomisation", "normality"))
Browse[1]>
debug: data.name <- paste(xname, ifelse(rank, "using rank correction",
""), "\nweights:", wname, ifelse(is.null(na.act), "",
paste("\nomitted:",
paste(na.act, collapse = ", "))), "\n")
Browse[1]>
debug: res <- list(statistic = statistic, p.value = PrI, estimate = vec,
alternative = alternative, method = method, data.name = data.name)
Browse[1]>
debug: if (!is.null(na.act)) attr(res, "na.action") <- na.act
Browse[1]>
debug: class(res) <- "htest"
Browse[1]>
debug: res
Browse[1]>
exiting from: moran.test(x, Network)
Moran's I test under randomisation
data: x
weights: Network
Moran I statistic standard deviate = NA, p-value = NA
alternative hypothesis: greater
sample estimates:
Moran I statistic Expectation Variance
-Inf -0.0002926544 NA
Warning message:
In moran.test(x, Network) : Out-of-range p-value: reconsider test arguments
>
-----Original Message-----
From: r-sig-geo-boun...@stat.math.ethz.ch [mailto:
r-sig-geo-boun...@stat.math.ethz.ch] On Behalf Of Roger Bivand
Sent: Wednesday, June 16, 2010 15:08
To: Michael Haenlein
Cc: r-sig-geo@stat.math.ethz.ch
Subject: Re: [R-sig-Geo] Problem with moran.test function
Well, Moran's I is -Inf, and the analytical variance is NA, so something is
not right. The problem could lie in x, Network, or the lag of x (when x and
Network are OK but their combination is unhappy). Can you run moran.test()
under debug() and check which values lead the value of I to go to -Inf? Is
this what is going on:
data(columbus)
set.seed(1)
x <- log(rpois(n=49, 2))
x
moran.test(x, nb2listw(col.gal.nb))
where the current spdep release fails reporting:
Error in lag.listw(listw, z, zero.policy = zero.policy, NAOK =
NAOK): Variable contains non-finite values
which was a fix introduced four weeks ago, changed to a test on |Inf| from a
test on NA:
https://r-forge.r-project.org/scm/viewvc.php/pkg/src/lagw.c?root=spdep&r1=244&r2=282
If you update spdep, you'll pick up the improvement (made thanks to a bug
report by Matias Mayor Fernandez), and if this is the case, the problem is
in the x.
Hope this helps,
Roger
--
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.biv...@nhh.no
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
[[alternative HTML version deleted]]
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@stat.math.ethz.ch
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
