r - 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。

对代码的有效性有什么意见吗？

Answer 1

对不起，我没有足够的教育来达到你检查两个协方差矩阵之间相关效率的目标，但我试图理解你的代码本身。

如果我是对的，您正在通过重新排列所有单元格并重新计算 mat2 的协方差矩阵之间的相关性，从同一矩阵 ( mat1 ) 中组合出 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

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