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

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