如何在spdep R的空间分析中输入相异矩阵
How to input dissimilarity matrix in spatial analysis in spdep R
目标： 我想在坐标对之间创建一个相异矩阵。我想使用此矩阵作为输入，使用 Moran's I (LISA) 和后者在地理加权回归 (GWR) 中计算局部空间聚类。
问题：我知道我可以使用dnearneigh{spdep}来计算距离矩阵。但是，我想使用我已经估计的多边形之间的旅行时间。实际上，我认为这就像输入一个相异矩阵，根据另一个特征告诉多边形之间的 distance/difference。我尝试将矩阵输入 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 它 returns 一个错误 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} 的使用是正确的。默认情况下，该函数保留行的名称以保持矩阵之间的对应关系。您还可以像这样指定 row.names：
listw_dissi <- mat2listw(diss_matrix_invers, row.names = row.names(diss_matrix_invers))
创建的列表将包含邻居的适当名称以及它们的距离作为权重。你可以通过查看邻居来检查这一点。
listw_dissi$neighbours[[1]][1:5]
而且您应该可以直接使用它来计算 Moran's 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")
目标： 我想在坐标对之间创建一个相异矩阵。我想使用此矩阵作为输入，使用 Moran's I (LISA) 和后者在地理加权回归 (GWR) 中计算局部空间聚类。
问题：我知道我可以使用dnearneigh{spdep}来计算距离矩阵。但是，我想使用我已经估计的多边形之间的旅行时间。实际上，我认为这就像输入一个相异矩阵，根据另一个特征告诉多边形之间的 distance/difference。我尝试将矩阵输入 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 它 returns 一个错误 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} 的使用是正确的。默认情况下，该函数保留行的名称以保持矩阵之间的对应关系。您还可以像这样指定 row.names：
listw_dissi <- mat2listw(diss_matrix_invers, row.names = row.names(diss_matrix_invers))
创建的列表将包含邻居的适当名称以及它们的距离作为权重。你可以通过查看邻居来检查这一点。
listw_dissi$neighbours[[1]][1:5]
而且您应该可以直接使用它来计算 Moran's 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")