典型相关分析设计矩阵的建立

Question

我有这个代码来生成对称矩阵来测试我正在执行的典型相关分析的设计是如何工作的。这是这个解决方案的扩展

矩阵的每一行表示一个数据集(它是对称的)，如果值为0时，它意味着数据集之间没有交互，如果存在交互的值更高。这段代码的最终目标只是对设计进行网格搜索，从而更好地解释我拥有的数据。但是，由于我需要提出不同的设计，添加更多或更少的数据集，我想知道如何将其改进为一个更通用的函数，特别是嵌套的for循环部分(如果我用5个数据集测试它，我会添加更多的for循环，最后添加更多的unlist )。

初始矩阵(我通常使用4个数据集)：

C <-  matrix(0,ncol = 4, nrow = 4)

每个数据集交互作用的权重(4，以避免多种组合)：

nweight <- 4
weight <- seq(from = 0, to = 1, length.out = nweight)

启动包含矩阵的列表

C_list <- vector("list", nweight)
cweight <- as.character(weight)
names(C_list) <- cweight

循环的每个位置，我想改变，以获得所有组合的权重，我想测试。

for(i1 in cweight) {
  C_list[[as.character(i1)]] <- vector("list", nweight)
  names(C_list[[(i1)]]) <- cweight
  for (i2 in cweight) {
    C_list[[(i1)]][[(i2)]] <- vector("list", nweight)
    names(C_list[[(i1)]][[(i2)]]) <- cweight
    for (i3 in cweight) {
      C_list[[(i1)]][[(i2)]][[(i3)]] <- vector("list", nweight)
      names(C_list[[(i1)]][[(i2)]][[(i3)]]) <- cweight
      for(i4 in cweight) {
        C_list[[(i1)]][[(i2)]][[(i3)]][[(i4)]] <- vector("list", nweight)
        names(C_list[[(i1)]][[(i2)]][[(i3)]][[(i4)]]) <- cweight
        for (i5 in cweight) {
          C_list[[(i1)]][[(i2)]][[(i3)]][[(i4)]][[(i5)]] <- vector("list", nweight)
          names(C_list[[(i1)]][[(i2)]][[(i3)]][[(i4)]][[(i5)]]) <- cweight
          for (i6 in cweight) {
            C[1, 2] <- as.numeric(i1)
            C[2, 1] <- as.numeric(i2)

            C[1, 3] <- as.numeric(i2)
            C[3, 1] <- as.numeric(i2)

            C[1, 4] <- as.numeric(i3)
            C[4, 1] <- as.numeric(i3)

            C[2, 3] <- as.numeric(i4)
            C[3, 2] <- as.numeric(i4)

            C[2, 4] <- as.numeric(i5)
            C[4, 2] <- as.numeric(i5)

            C[4, 3] <- as.numeric(i6)
            C[3, 4] <- as.numeric(i6)
            C_list[[i1]][[i2]][[i3]][[i4]][[i5]][[i6]] <- C
          }
        }
      }
    }
  }
}

取消.的名单。嵌套矩阵最终得到一个包含每个数据集权重的矩阵长列表。

C_list2 <- unlist(unlist(unlist(unlist(unlist(C_list, FALSE, FALSE), 
                                       FALSE, FALSE), FALSE, FALSE), 
                         FALSE, FALSE), FALSE, FALSE)

Answer 1

在这里，您想离开for循环有两个原因：

1. for循环的数量取决于数据集的数量，因此使用for循环可以防止将代码泛化到任意数量的数据集。
2. 当使用更多的数据集时，许多for循环可能会减慢代码执行速度。

我认为矢量化for循环的关键是使用expand.grid函数。如果你有

n <- 4

数据集，那么你就有了

p <- n * (n - 1) / 2    # 6

自由度(代码中的for循环数，或每个矩阵下三角形上的项目数)。如果你能从中挑出每一个

w <- seq(from = 0, to = 1, length.out = n)

然后，通过执行以下操作，可以构建所有可能组合的矩阵：

W <- as.matrix(expand.grid(rep(list(w), p)))

这里，W是一个包含4096行的大型矩阵，每一行代表不同的(i1, i2, i3, i4, i5, i6)变量组合：

> head(W)
          Var1      Var2 Var3 Var4 Var5 Var6
[1,] 0.0000000 0.0000000    0    0    0    0
[2,] 0.3333333 0.0000000    0    0    0    0
[3,] 0.6666667 0.0000000    0    0    0    0
[4,] 1.0000000 0.0000000    0    0    0    0
[5,] 0.0000000 0.3333333    0    0    0    0
[6,] 0.3333333 0.3333333    0    0    0    0

这些6列只是每个矩阵所需的n * n = 16值的一部分。我们可以使用以下方法进行扩展：

X <- matrix(1:(n*n), n, n) # pattern matrix of indices
A <- matrix(0, nrow(W), n * n)
A[,    X[lower.tri(X)]] <- W
A[, t(X)[lower.tri(X)]] <- W

A与W相似，因为它是一个带有4096行的矩阵，但是现在每行都有对称矩阵的n * n = 16值。

从这里，您可以将A重塑为一个3D数组：

dim(A) <- c(nrow(W), n, n)

您的4096矩阵可以按以下方式访问：

A[1, , ]
#      [,1] [,2] [,3] [,4]
# [1,]    0    0    0    0
# [2,]    0    0    0    0
# [3,]    0    0    0    0
# [4,]    0    0    0    0

A[10, , ]
#           [,1]      [,2]      [,3] [,4]
# [1,] 0.0000000 0.3333333 0.6666667    0
# [2,] 0.3333333 0.0000000 0.0000000    0
# [3,] 0.6666667 0.0000000 0.0000000    0
# [4,] 0.0000000 0.0000000 0.0000000    0

A[4096, , ]
#      [,1] [,2] [,3] [,4]
# [1,]    0    1    1    1
# [2,]    1    0    1    1
# [3,]    1    1    0    1
# [4,]    1    1    1    0

如果我是你，我可能会停在这里，也就是说，把数据保存在这个表格中。如果分析的其余部分允许，3d数组可能允许您继续编写矢量化代码。但是，如果您绝对需要一个矩阵列表，则可以：

C_list2 <- lapply(seq(nrow(A)), function(i) A[i, , ])

(请注意，我的数据和您的数据中矩阵的顺序不匹配。请告诉我这是否需要考虑，这可能是重新组织W矩阵的行和/或列的问题。)