R中的引导变量相关性

Question

我的目的是编写几个函数，旨在找出两个协方差矩阵之间的总体相似性，或者通过将它们与随机向量相乘并关联响应向量，或者通过引导其中一个矩阵来获得可用于的相关系数分布比较。但在这两种情况下，我都得到了错误的结果。观察到的矩阵间相关性高达 0.93，但分布范围最高仅达到 0.2。这是函数的代码：

resamplerSimAlt <- function(mat1, mat2, numR, graph = FALSE)
{
  statSim <- numeric(numR)
  mat1vcv <- cov(mat1)
  mat2vcvT <- cov(mat2)
  ltM1 <- mat1vcv[col(mat1vcv) <= row(mat1vcv)]
  ltM2T <- mat2vcvT[col(mat2vcvT) <= row(mat2vcvT)]
  statObs <- cor(ltM1, ltM2T)                           
  indice <- c(1:length(mat2))
  resamplesIndices <- lapply(1:numR, function(i) sample(indice, replace = F))
  for (i in 1:numR)
  {
    ss <- mat2[sample(resamplesIndices[[i]])]
    ss <- matrix(ss, nrow = dim(mat2)[[1]], ncol = dim(mat2)[[2]])
    mat2ss <- cov(ss)
    ltM2ss <- mat2ss[col(mat2ss) <= row(mat2ss)]
    statSim[i] <- cor(ltM1, ltM2ss)
  }
  if (graph == TRUE)
  {
    plot(1, main = "resampled data density distribution", xlim = c(0, statObs+0.1), ylim = c(0,14))
    points(density(statSim), type="l", lwd=2)
    abline(v = statObs)
    text(10, 10, "observed corelation = ")
  }
  list( obs = statObs , sumFit = sum(statSim > statObs)/numR)
}  

事实上，我很难相信两个原始矩阵之间的相关系数很高，并且第一个原始矩阵和第二个重新采样的矩阵之间的相关系数在 10000 次引导重复后最大为 0.2。

任何关于代码有效性的cmets？

Answer 1

抱歉，我没有受过足够的教育，无法实现您关于检查两个协方差矩阵之间的相关效率的目标，但我试图理解您的代码本身。

如果我是对的，您将通过重新排列所有单元格，并重新计算 mat1 的协方差矩阵与重采样的协方差矩阵之间的相关性，从同一个矩阵 (mat2) 组成 10.000 个不同的矩阵大批。这些存储在statSim 变量中。

你说原来的关联效率很高（statObs），但是statSim的最大值很低，这很奇怪。我认为问题在于您的结果列表：

list( obs = statObs , sumFit = sum(statSim > statObs)/numR)

返回原始相关系数 (obs)，而不是 sumFit 的书面最大值。在那里你可能会使用例如。 max(statSim)。我看到返回 sumFit 以检查重采样是否对相关效率有任何改进的意义，但是根据您的代码，我认为该理论没有问题。

使用模拟相关系数max 更新函数：

resamplerSimAlt <- function(mat1, mat2, numR, graph = FALSE)
{
  statSim <- numeric(numR)
  mat1vcv <- cov(mat1)
  mat2vcvT <- cov(mat2)
  ltM1 <- mat1vcv[col(mat1vcv) <= row(mat1vcv)]
  ltM2T <- mat2vcvT[col(mat2vcvT) <= row(mat2vcvT)]
  statObs <- cor(ltM1, ltM2T)                           
  indice <- c(1:length(mat2))
  resamplesIndices <- lapply(1:numR, function(i) sample(indice, replace = F))
  for (i in 1:numR)
  {
    ss <- mat2[sample(resamplesIndices[[i]])]
    ss <- matrix(ss, nrow = dim(mat2)[[1]], ncol = dim(mat2)[[2]])
    mat2ss <- cov(ss)
    ltM2ss <- mat2ss[col(mat2ss) <= row(mat2ss)]
    statSim[i] <- cor(ltM1, ltM2ss)
  }
  if (graph == TRUE)
  {
    plot(1, main = "resampled data density distribution", xlim = c(0, statObs+0.1), ylim = c(0,14))
    points(density(statSim), type="l", lwd=2)
    abline(v = statObs)
    text(10, 10, "observed corelation = ")
  }
  list( obs = statObs , sumFit = sum(statSim > statObs)/numR, max=max(statSim))
}

我跑过什么：

> mat1 <- matrix(runif(25),5,5)
> mat2 <- mat1+0.2
> resamplerSimAlt(mat1, mat2, 10000)
$obs
[1] 1

$sumFit
[1] 0

$max
[1] 0.94463

随机mat2:

> mat2 <- matrix(runif(25),5,5)
> resamplerSimAlt(mat1, mat2, 10000)
$obs
[1] 0.31144

$sumFit
[1] 0.9124

$max
[1] 0.9231

我的答案可能不是真正的答案。如果是这种情况，请提供有关问题的更多详细信息:)