Alexandre, bom dia!

Se ainda lhe for útil, segue uma sugestão...

### <code r>
var1 <- rnorm(200,5,0.25)
var2 <- 1:500
var3 <- rnorm(100,5,0.25)
var4 <- 500:1

cor.prob <- function(X){ 
      pair.SampSize <- crossprod(!is.na(X)) 
      above1 <- row(pair.SampSize) < col(pair.SampSize) 
      pair.df <- pair.SampSize[above1] - 2 

      R <- cor(X, use="pair") 
      above2 <- row(R) < col(R) 
      r2 <- R[above2]^2 
      Fstat <- (r2 * pair.df)/(1 - r2) 
      R[above2] <- 1 - pf(Fstat, 1, pair.df) 
      R 
      } 

tam   <- 100 # amostra 
n_sim <-  99 # simulações 

RES <- NULL 
for(i in 1:n_sim){
     dados   <- sapply(c("var1", "var2", "var3", "var4"), function(x) sample(get(x), 100))
     correla <- round(cor.prob(dados),4) 
     res     <- cbind(i, correla) 
     RES     <- rbind(RES, res) 

RES[c(1:5,392:396),]
#       i    var1    var2    var3   var4
# var1  1  1.0000  0.8025  0.8451 0.9178
# var2  1 -0.0253  1.0000  0.6609 0.0106
# var3  1 -0.0198  0.0444  1.0000 0.4850
# var4  1 -0.0104 -0.2545 -0.0706 1.0000
# var1  2  1.0000  0.5200  0.1704 0.8461
# var4 98  0.1309  0.0332 -0.1361 1.0000
# var1 99  1.0000  0.1433  0.2047 0.3260
# var2 99 -0.1474  1.0000  0.5391 0.9524
# var3 99 -0.1279 -0.0621  1.0000 0.2885
# var4 99 -0.0992 -0.0060  0.1072 1.0000
### </code>



================================================
Éder Comunello
Agronomist (UEM), MSc in Environ. Sciences (UEM)
DSc in Agricultural Systems Engineering (USP/Esalq)
Brazilian Agricultural Research Corporation (Embrapa)
Dourados, MS, Brazil |<O>|
================================================
GEO, -22.2752, -54.8182, 408m
UTC-04:00 / DST: UTC-03:00