编写算法来识别症状对

Question

由于我是新的 R 用户，我正在努力编写循环或类似的方法来解决以下问题。

所以问题： 我有一个包含 200 名患者和 30 个症状的列表，其中信息症状存在是 (1) 或否 (0)。我想创建一个类似于相关矩阵的矩阵，每对症状都会告诉我这两种症状在同一个人身上出现的百分比。

所以：

ID /sym 1 / sym 2 / sym 3 /....
Pat 1 / 1 / 1 / 0 / ...    
Pat 2 / 1 / 1 / 1 /
Pat 3 / 1 / 1 / 0 / 
...

然后我想使用一个看起来像标准相关矩阵的矩阵，但对于每个症状对都显示两者的百分比，因为我认为关联这些分类值没有任何意义。

到目前为止，我已经开始为每个组合编写一个新变量，告诉我两者是否都存在，然后用它来计算百分比并将其填充到我的矩阵中。

不幸的是，考虑到可能的组合数量，这需要 FOREVER 来编写，但我不知道如何迭代它。也许你能帮忙？它可能非常简单，但我还不足以让程序员想到它。

Answer 1

设 x 为 200 x 30 的数据数组。根据矩阵乘法的规则，x' * x 是一个 30 x 30 的数组，其 (i,j) 项计算列 i 和 j 等于 1 的位置数。将这些计数除以 200 得到比例并将其乘以 100将其转换为所需的百分比。

这里是单行实现以及用于测试和说明它的代码。对于比您的问题大 10,000 倍的问题（20,000 名患者和 3000 个症状），在这台机器上执行大约需要 10 秒。由于它几乎与大小呈线性关系，因此您的问题的计算应该需要大约一毫秒。

#
# For columns i and j, f(x)[i,j] is the percentage of rows in which
# both columns of `x` are TRUE (optionally: nonzero).
#
f <- function(x) (t(x) %*% x) * (100 / dim(x)[1])
#
# Slow version to demonstrate `f` is correct.
#
f.direct <- function(x) {
  m <- dim(x)[1]
  n <- dim(x)[2]
  #
  # Test all elements of `x` to create a logical array.
  #
  x.indicator <- x != 0
  #
  # Initialize the result.
  #
  y <- matrix(NA_real_, n, n)
  #
  # Loop over pairs of columns.
  #
  for (i in 1:n) {
    for (j in 1:n) {
      # Compare column `i` to column `j` by averaging the times their
      # indicators are equal.  Multiply by 100 to give a percentage.
      y[i,j] <- 100 * mean(x.indicator[,i] & x.indicator[,j])
    }
  }
  return(y)
}
#------------------------------------------------------------------------------#
#
# Create some data and test `f` on them.
#
m <- 200  # Number of rows
n <- 30   # Number of columns
p <- 0.1  # Expected proportion of 1's
x <- matrix(runif(m*n) < p, m, n)

system.time(y <- f(x))                # Almost instantaneous
system.time(y.direct <- f.direct(x))  # A thousand times slower (but not bad)
#
# Display the results.
#
par(mfrow=c(1,2))
image(y, main="Matrix Result")
image(y.direct, main="Direct Result")
par(mfrow=c(1,1))
#
# Compare them and report the outcome.
#
if(all.equal(y, y.direct)) cat("Results are equal.") else cat("There's a difference!")