如何在SPDEP R中的空间分析中输入差异矩阵



aim:我想在成对的坐标之间创建一个相似的矩阵。我想将此矩阵用作输入,以使用Moran的I(Lisa(计算本地空间簇,而后者则在地理位置加权回归(GWR(中。

问题:我知道我可以使用dnearneigh{spdep}来计算距离矩阵。但是,我想使用我已经估计的多边形之间的旅行时间。在实践中,我认为这就像输入一个差异矩阵,该矩阵告诉了基于另一个特征的多边形之间的距离/差异。我已经尝试将矩阵输入到dnearneigh{spdep},但是我得到了错误Error: ncol(x) == 2 is not TRUE

dist_matrix <- dnearneigh(diss_matrix_invers, d1=0, d2=5, longlat = F, row.names=rn)

有什么建议吗?下面有一个可再现的示例:

编辑:进一步挖掘,我认为我可以使用mat2listw{spdep},但我仍然不确定它是否可以保持矩阵和多边形之间的对应关系。如果添加row.names = T,它将返回错误row.names wrong length :(

listw_dissi <- mat2listw(diss_matrix_invers)
lmoran <- localmoran(oregon.tract@data$white, listw_dissi, 
                     zero.policy=T, alternative= "two.sided")

可重复的示例

library(UScensus2000tract)
library(spdep)
library(ggplot2)
library(dplyr)
library(reshape2)
library(magrittr)
library(data.table)
library(reshape)
library(rgeos)
library(geosphere)
# load data
  data("oregon.tract")
# get centroids as a data.frame
  centroids <- as.data.frame( gCentroid(oregon.tract, byid=TRUE) )
# Convert row names into first column
  setDT(centroids, keep.rownames = TRUE)[]
# create Origin-destination pairs
  od_pairs <- expand.grid.df(centroids, centroids) %>% setDT()
  colnames(od_pairs) <- c("origi_id", "long_orig", "lat_orig", "dest_id", "long_dest", "lat_dest")     
# calculate dissimilarity between each pair. 
# For the sake of this example, let's use ellipsoid distances. In my real case I have travel-time estimates
  od_pairs[ , dist := distGeo(matrix(c(long_orig, lat_orig), ncol = 2), 
                         matrix(c(long_dest, lat_dest), ncol = 2))]
# This is the format of how my travel-time estimates are organized, it has some missing values which include pairs of origin-destination that are too far (more than 2hours apart)
  od_pairs <- od_pairs[, .(origi_id, dest_id, dist)]
  od_pairs$dist[3] <- NA
  >      origi_id    dest_id         dist
  > 1:   oregon_0   oregon_0      0.00000
  > 2:   oregon_1   oregon_0           NA
  > 3:   oregon_2   oregon_0  39874.63673
  > 4:   oregon_3   oregon_0  31259.63100
  > 5:   oregon_4   oregon_0  33047.84249
# Convert to matrix
  diss_matrix <- acast(od_pairs, origi_id~dest_id, value.var="dist") %>% as.matrix()
# get an inverse matrix of distances, make sure diagonal=0
  diss_matrix_invers <- 1/diss_matrix
  diag(diss_matrix_invers) <- 0

计算简单距离矩阵

  # get row names
    rn <- sapply(slot(oregon.tract, "polygons"), function(x) slot(x, "ID"))
  # get centroids coordinates
    coords <- coordinates(oregon.tract)
  # get distance matrix
    diss_matrix <- dnearneigh(diss_matrix_invers, d1=0, d2=5, longlat =T, row.names=rn)
class(diss_matrix)
> [1] "nb"

现在如何在这里使用我的diss_matrix_invers

您对使用matlistw {spdep}的使用是正确的。默认情况下,该函数保留了行的名称,以保持矩阵之间的对应关系。您也可以指定行类似的名称:

listw_dissi <- mat2listw(diss_matrix_invers, row.names = row.names(diss_matrix_invers))  

创建的列表将包含适合邻居的适当名称,以及它们的距离。您可以通过查看邻居来检查一下。

listw_dissi$neighbours[[1]][1:5]

您应该能够直接使用它来计算Moran的I。

dnearneigh {sdep}
您无法在Dnearneigh {spdep}中使用diss_matrix,因为此功能列出了坐标列表。

但是,如果您需要使用自己的距离矩阵(旅行时间(定义给定距离阈值(D1,D2(的一组邻居。我认为此功能可以解决问题。

dis.neigh<-function(x, d1 = 0, d2=50){
  #x must be a symmetrical distance matrix
  #create empty list
  style = "M" #for style unknown
  neighbours<-list()
  weights<-list()
  #set attributes of neighbours list
  attr(neighbours, "class")<-"nb"
  attr(neighbours, "distances")<-c(d1,d2)
  attr(neighbours, "region.id")<-colnames(x)
  #check each row for neighbors that satisfy distance threshold
  neighbour<-c()
  weight<-c()
  i<-1
  for(row in c(1:nrow(x))){
    j<-1
    for(col in c(1:ncol(x))){
      if(x[row,col]>d1 && x[row,col]<d2){
        neighbour[j]<-col
        weight[j]<-1/x[row,col] #inverse distance (dissimilarity)
        j<-1+j
      }
    }
    neighbours[i]<-list(neighbour)
    weights[i]<-list(weight)
    i<-1+i
  }
  #create neighbour and weight list
  res <- list(style = style, neighbours = neighbours, weights = weights)
  class(res) <- c("listw", "nb")
  attr(res, "region.id") <- attr(neighbours, "region.id")
  attr(res, "call") <- match.call()
  return(res)
}

并像这样使用:

nb_list<-dis.neigh(diss_matrix, d1=0, d2=10000)
lmoran <- localmoran(oregon.tract@data$white, nb_lists, alternative= "two.sided")

最新更新