我想通过计算每对中所有(多维)点集之间距离的平均值来量化组相似度。
我可以手动为每对组手动轻松完成此操作,如下所示:
library(dplyr)
library(tibble)
library(proxy)
# dummy data
set.seed(123)
df1 <- data.frame(x = rnorm(100,0,4),
y = rnorm(100,1,5),
z = rbinom(100, 1, 0.1))
df2 <- data.frame(x = rnorm(100,-1,3),
y = rnorm(100,0,6),
z = rbinom(100, 1, 0.1))
df3 <- data.frame(x = rnorm(100,-30,4),
y = rnorm(100,10,2),
z = rbinom(100, 1, 0.9))
# compute distance (unscaled, uncentred data)
dist(df1, df2, method = "gower") %>% mean
dist(df1, df3, method = "gower") %>% mean
dist(df2, df3, method = "gower") %>% mean
但我想以某种方式将其矢量化,因为我的实际数据有30多个组。一个简单的for循环可以实现这样:
# combine data and scale, centre
df <- rbind(df1, df2, df3) %>%
mutate(id = rep(1:3, each = 100))
df <- df %>%
select(-id) %>%
transmute_all(scale) %>%
add_column(id = df$id)
# create empty matrix for comparisons
n <- df$id %>% unique %>% length
m <- matrix(nrow = n, ncol = n)
# loop through each pair once
for(i in 1:n) {
for(j in 1:i) { #omit top right corner
if(i == j) {
m[i,j] <- NA #omit diagonal
} else {
m[i,j] <- dist(df[df$id == i,1:3], df[df$id == j,1:3], method = "gower") %>% mean
}
}
}
m
[,1] [,2] [,3]
[1,] NA NA NA
[2,] 0.2217443 NA NA
[3,] 0.8446070 0.8233932 NA
然而,这种方法可预测地严重缩放;快速基准测试表明,我的实际数据需要90多个小时,其中包含30多个组,每组1000多行。
任何人都可以建议一个更有效的解决方案,或者可能是一种根本不同的方式来解决我所缺少的问题吗?
您可以获取每对组,将它们连接起来,然后只计算该组中的相异矩阵。显然,这意味着您在一定程度上将一个组与自身进行比较,但它可能仍适用于您的用例,而使用daisy
,您的数据大小相当快。
library(cluster)
n <- 30
groups <- vector("list", 30)
# dummy data
set.seed(123)
for(i in 1:30) {
groups[[i]] = data.frame(x = rnorm(1000,ceiling(runif(1, -10, 10)),ceiling(runif(1, 2, 4))),
y = rnorm(1000,ceiling(runif(1, -10, 10)),ceiling(runif(1, 2, 4))),
z = rbinom(1000,1,runif(1,0.1,0.9)))
}
m <- matrix(nrow = n, ncol = n)
# loop through each pair once
for(i in 1:n) {
for(j in 1:i) { #omit top right corner
if(i == j) {
m[i,j] <- NA #omit diagonal
} else {
# concatenate groups
dat <- rbind(df_list[[i]], df_list[[j]])
# compute all distances (between groups and within groups), return matrix
mm <- dat %>%
daisy(metric = "gower") %>%
as.matrix
# retain only distances between groups
mm <- mm[(nrow(df_list[[i]])+1):nrow(dat) , 1:nrow(df_list[[i]])]
# write mean distance to global comparison matrix
m[i,j] <- mean(mm)
}
}
}
我不确定这是否会好起来,但这是另一种方法。你使用ls
获得矩阵的名称,combn
生成两个对,然后使用get
来获得计算dist
的矩阵
do.call(rbind,
combn(ls(pattern = "df\\d+"), 2, FUN = function(x)
data.frame(pair = toString(x),
dist = mean(dist(get(x[1]), get(x[2]), method = "gower")),
stringsAsFactors = FALSE),
simplify = FALSE
))
# pair dist
#1 df1, df2 0.2139304
#2 df1, df3 0.8315169
#3 df2, df3 0.8320911
proxy
可以使用矩阵列表作为输入,您只需要定义一个执行您想要的包装函数:
nested_gower <- function(x, y, ...) {
mean(proxy::dist(x, y, ..., method = "gower"))
}
proxy::pr_DB$set_entry(
FUN = nested_gower,
names = c("ngower"),
distance = TRUE,
loop = TRUE
)
df_list <- list(df1, df2, df3)
proxy::dist(df_list, df_list, method = "ngower")
[,1] [,2] [,3]
[1,] 0.1978306 0.2139304 0.8315169
[2,] 0.2139304 0.2245903 0.8320911
[3,] 0.8315169 0.8320911 0.2139049
这仍然会很慢,但它应该比普通R中的for
循环更快(proxy
在背景中使用C)。
重要提示:请注意,生成的交叉距离矩阵的对角线不具有零。如果你像dist
那样打电话给proxy::dist(df_list, method = "ngower")
,proxy
将假设distance(x, y) = distance(y, x)
(对称性)和distance(x, x) = 0
,后者在这种情况下不是真的。将两个参数传递给dist
可以防止这种假设。如果你真的不关心对角线,只通过一个参数来避免计算上三角形来节省一些额外的时间。或者,如果您关心对角线但仍想避免计算上三角形,请先用一个参数调用dist
,然后调用proxy::dist(df_list, df_list, method = "ngower", pairwise = TRUE)
。
旁注:如果您想使用gower
包模仿此行为(如d.b所示),您可以将包装函数定义为:
nested_gower <- function(x, y, ...) {
distmat <- sapply(seq_len(nrow(y)), function(y_row) {
gower::gower_dist(x, y[y_row, , drop = FALSE], ...)
})
mean(distmat)
}
但是,返回的值似乎会根据传递给函数的记录数而发生变化,因此很难说出最佳方法是什么。
*如果要在proxy::pr_DB$delete_entry("ngower")
中重新定义函数,请先使用proxy
。
如果您更喜欢proxy
版本的Gower跨距离矩阵,我会发现您可以利用我的dtwclust
包的一些功能来并行执行计算:
library(dtwclust)
library(doParallel)
custom_dist <- new("tsclustFamily", dist = "ngower", control = list(symmetric = TRUE))@dist
workers <- makeCluster(detectCores())
registerDoParallel(workers)
distmat <- custom_dist(df_list)
stopCluster(workers); registerDoSEQ()
对于您的实际用例而言,这可能会更快(对于此处的小样本数据而言,这不是很多)。关于对角线的相同警告(所以使用custom_dist(df_list, df_list)
或custom_dist(df_list, pairwise = TRUE)
)。如果您想了解更多信息,请参阅第3.2节here和tsclustFamily
的文档。