r - 是否可以加快我创建相关矩阵的功能？

Question

我编写了以下函数来使用所谓的 Cramér's V 估计多项变量的成对相关性。 .我使用 vcd为此目的打包，但据我所知，没有现有的函数可以从矩阵或 data.frame 创建 V 的对称相关矩阵类似于 cor .

功能是:

require(vcd)
get.V<-function(y){
  col.y<-ncol(y)
  V<-matrix(ncol=col.y,nrow=col.y)
  for(i in 1:col.y){
    for(j in 1:col.y){
      V[i,j]<-assocstats(table(y[,i],y[,j]))$cramer
    }
  }
  return(V)
}

但是，对于大量变量，它变得相对较慢。

no.var<-5
y<-matrix(ncol=no.var,sample(1:5,100*no.var,TRUE))
get.V(y)

随着你的增加 no.var计算时间可能会爆炸。因为我需要将此应用于 data.frame长度为 100 或更高，我的问题是，是否有可能通过更优雅的编程来“加速”我的功能。谢谢你。

Answer 1

以及减少执行的测试数量，或以其他方式
优化整个函数的运行，我们也许可以使assocstats快点。我们将从建立一个测试用例开始
确保我们不会意外地创建一个不正确的更快的函数。

x <- vcd::Arthritis$Improved
y <- vcd::Arthritis$Treatment
correct <- vcd::assocstats(table(x, y))$cramer
correct

## [1] 0.3942

is_ok <- function(x) stopifnot(all.equal(x, correct))

我们将从制作 assocstats 的版本开始那非常接近
原来的。

cramer1 <- function (x, y) {
  mat <- table(x, y)

  tab <- summary(MASS::loglm(~1 + 2, mat))$tests

  phi <- sqrt(tab[2, 1] / sum(mat))
  cont <- sqrt(phi ^ 2 / (1 + phi ^ 2))

  sqrt(phi ^ 2 / min(dim(mat) - 1))
}
is_ok(cramer1(x, y))

这里最慢的操作将是 loglm ，所以在我们尝试之前
使其更快，值得寻找替代方法。一种
小谷歌搜索 a useful blogpost .
让我们也试试:

cramer2 <- function(x, y) {
  chi <- chisq.test(x, y, correct=FALSE)$statistic[[1]]

  ulength_x <- length(unique(x))
  ulength_y <- length(unique(y))

  sqrt(chi / (length(x) * (min(ulength_x, ulength_y) - 1)))
}
is_ok(cramer2(x, y))

性能如何叠加:

library(microbenchmark)

microbenchmark(
  cramer1(x, y),
  cramer2(x, y)
)

## Unit: microseconds
##           expr    min     lq median     uq  max neval
##  cramer1(x, y) 1080.0 1149.3 1182.0 1222.1 2598   100
##  cramer2(x, y)  800.7  850.6  881.9  934.6 1866   100

cramer2()是比较快的。 chisq.test()很可能是瓶颈，所以
让我们看看我们是否可以通过减少操作来使该功能更快: chisq.test()做的不仅仅是计算测试统计量，所以它是
很可能我们可以让它更快。几分钟的细心工作减少
功能:

chisq_test <- function (x, y) {
  O <- table(x, y)
  n <- sum(O)

  E <- outer(rowSums(O), colSums(O), "*")/n

  sum((abs(O - E))^2 / E)
}

然后我们可以创建一个新的 cramer3()使用 chisq.test() .

cramer3 <- function(x, y) {
  chi <- chisq_test(x, y)

  ulength_x <- length(unique(x))
  ulength_y <- length(unique(y))

  sqrt(chi / (length(x) * (min(ulength_x, ulength_y) - 1)))
}
is_ok(cramer3(x, y))
microbenchmark(
  cramer1(x, y),
  cramer2(x, y),
  cramer3(x, y)
)

## Unit: microseconds
##           expr    min     lq median     uq  max neval
##  cramer1(x, y) 1088.6 1138.9 1169.6 1221.5 2534   100
##  cramer2(x, y)  796.1  840.6  865.0  906.6 1893   100
##  cramer3(x, y)  334.6  358.7  373.5  390.4 1409   100

现在我们有了自己的简单版本 chisq.test()我们可以
使用 table() 的结果可以提高一点速度想出
找出 x 中唯一元素的数量和 y :

cramer4 <- function(x, y) {
  O <- table(x, y)
  n <- length(x)
  E <- outer(rowSums(O), colSums(O), "*")/n

  chi <- sum((abs(O - E))^2 / E)
  sqrt(chi / (length(x) * (min(dim(O)) - 1)))
}
is_ok(cramer4(x, y))
microbenchmark(
  cramer1(x, y),
  cramer2(x, y),
  cramer3(x, y),
  cramer4(x, y)
)

## Unit: microseconds
##           expr    min     lq median     uq  max neval
##  cramer1(x, y) 1097.6 1145.8 1183.3 1233.3 2318   100
##  cramer2(x, y)  800.7  840.5  860.7  895.5 2079   100
##  cramer3(x, y)  334.4  353.1  365.7  384.1 1654   100
##  cramer4(x, y)  248.0  263.3  273.2  283.5 1342   100

不错 - 我们仅使用 R 代码就使其速度提高了 4 倍。从这里，你
可以尝试通过以下方式获得更快的速度:

使用 tcrossprod()而不是 outer()

制作更快的版本 table()对于这个特殊的 (2d) 案例

使用 Rcpp 从表格数据中计算检验统计量