Bootstrap variables correlation in R

Question

I wrote this R function in order to test between matrix (mat1, mat2) correlation by re-sampling the second matrix numR times to obtain correlation coefficient distribution against which to test the observed value. But irrespective of the observed value I always get the distribution with SD 0.3 max. I don`t see how if the original correlation is as high as 0.9 that after 10000 re-samples it still cannot exceed 0.3. Any comments on the code or a possible test with independent data?

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)
}

Answer 1

The problem is that you're computing the correlation between the covariance matrices. When you generate the covariance matrices for the columns of mat1, mat2, and all the resampled mat2s, they will all tend to have large values on the diagonal. Thus, when you try to compute the correlation between these, there will always be some correlation.

I am not totally clear from your question, since it says "bootstrap" but all you do in the code is randomly permute the data without replacement. However, from the sumFit output argument in your code, it looks like you're trying to use permutation methods to calculate something like an empiric p-value (i.e. the proportion of the statSims determined from the resampled mat2s that are more extreme than the statObs you got using the true data). If that is the case, then here is a simplified example using the absolute value of the correlation as the test statistic:

resamplerSimAlt <- function(mat1, mat2, numR, graph = FALSE)
{  
  statObs <- abs(cor(as.vector(mat1), as.vector(mat2)))

  statSim <- replicate(numR, abs(cor(as.vector(mat1), mat2[sample(length(mat2))])))

  if (graph == TRUE)
  {
    plot(density(statSim),
      main = "resampled data density distribution",
      type = "l",
      lwd  = 2,
      xlim = c(0, max(c(statSim, statObs))))
    abline(v = statObs)
  }
  list( obs = statObs, sumFit = sum(statSim > statObs)/numR, statSim = statSim)
}

set.seed(12345)
n    <- 20
mat1 <- matrix(rnorm(n^2), n)
mat2 <- mat1 + rnorm(n^2, 0, 10)
numR <- 10000

dev.new(width=4, height=4)
results <- resamplerSimAlt(mat1, mat2, numR, graph=T)

> results[1:2]
$obs
[1] 0.1190381

$sumFit
[1] 0.0183

Also take a look at the boot package, which has a ton of related methods that can be used off-the-shelf for getting empiric confidence intervals and such.