1

我正在尝试针对我的用例优化 R 的 spdep 函数，因为它对于大型数据库来说非常慢。我做得很好，但我坚持了一点，我试图找到我的权重矩阵的踪迹以进行 LM 错误测试。我认为公式是 tr[(W' + W) W] （Anselin, L., Bera, AK, Florax, R. 和 Yoon, MJ 1996 的第 82 页空间依赖性的简单诊断测试。区域科学和城市经济学， 26, 77–104.) W 是一个方形权重矩阵，保存每个观测值与另一个观测值的空间关系。tr() 运算是对角线的总和。

在我的例子中，权重矩阵是对称的，对角线为零。所以，我认为公式 tr[(W' + W) W] 等于 2*sumsq(W)，速度非常快。但显然我在某个地方弄错了，因为结果与 spdep 库的结果不匹配，这很可能是正确的。

spdep 库的相关部分在这里。谁能帮助我以下函数的结果与 2*sumsq(W) 有何不同或如何使其更快？这个函数是 lm.LMtests 函数被大型数据集阻塞的地方。

tracew <- function (listw) {
    dlmtr <- 0
    n <- length(listw$neighbours)
    if (n < 1) stop("non-positive n")
    ndij <- card(listw$neighbours)
    dlmtr <- 0
    for (i in 1:n) {
        dij <- listw$neighbours[[i]]
        wdij <- listw$weights[[i]]
        for (j in seq(length=ndij[i])) {
            k <- dij[j]
# Luc Anselin 2006-11-11 problem with asymmetric listw
                dk <- which(listw$neighbours[[k]] == i)
                if (length(dk) > 0L && dk > 0L &&
                dk <= length(listw$neighbours[[k]]))
                wdk <- listw$weights[[k]][dk]
                else wdk <- 0
                dlmtr <- dlmtr + (wdij[j]*wdij[j]) + (wdij[j]*wdk)
        }
    }
    dlmtr
}

对那些不熟悉 R 的 spdep 库的人的补充说明：函数的输入 listw 包含具有两个列表列表的权重矩阵的“图形”实现。listw$neighbors 是一个列表，其中每个列表项是与观察相关的观察索引的列表。listw$weights 具有邻居的相同结构的列表，除了它保存关系的权重。

提前感谢您的任何评论和指示。

# example code

# initiliaze
library(spdep)
library(multiway)
# load the tracew function above
data(columbus)
columbus = columbus[rep(row.names(columbus), 20), ] # the difference becomes dramatic when n is high. try not replicating at first to see the results.

# manual calculation, using sumsq
w = distm(cbind(columbus$X, columbus$Y))
w[w > 1000000] = Inf # remove some relations acc. to pre-defined rule
w = 1/(1+w)
diag(w) = 0
w = w / (sum(w) / length(columbus$X)) #"C style" standardization
2*sumsq(w)

# spdep calculation
neighs.band = dnearneigh(cbind(columbus$X, columbus$Y), 0, 1000, longlat = TRUE)
w.spdep = lapply(nbdists(neighs.band, cbind(columbus$X, columbus$Y), longlat = TRUE), function(x) 1/(0.001+x))
my.listw = nb2listw(neighs.band, glist = w.spdep, style="C")
tracew(my.listw)
4

0 回答 0