如何在R中创建模拟相关系数的for循环?

lbsnaicq  于 2023-11-14  发布在  其他
关注(0)|答案(1)|浏览(115)

我计划在现有数据集上运行10,000次模拟,并希望找到每个模拟的相关系数并将其保存到一个嵌套框架中(理论上最终会有10,000个相关系数)。我如何编写代码来运行模拟并在for循环中每次打印出系数?
这是我尝试的最新代码:

for (i in 1:10000) {
  bnsamp <- data.frame(mvrnorm(n = 48, mu = mean_combined, Sigma = sigma))
  colnames(bnsamp) <- c("Amyg_Avg", "Sys_Jus")
  cor_coeff <- cor.test(x= bnsamp$Amyg_Avg, y= bnsamp$Sys_Jus, 
                        method= "pearson", alternative= "two.sided")
   cor_coef_est <- list(cor_coeff$estimate == i)
   cor_coef_p <- list(cor_coeff$p.value == i)
   coeff_list <- data.frame(cor_coef_est, cor_coef_p)
   print(coeff_list)
}

字符串
最终,我的问题是如何循环重新计算相关性并每次打印它(以便稍后绘制成分布)。谢谢!

2wnc66cl

2wnc66cl1#

这里有两种模拟所需值的方法。
1.使用for循环:
1.关于replicate
help('replicate'),我的重点:
replicatesapply的常用 Package 器,用于重复计算表达式(通常涉及随机数生成)。

library(MASS)

mean_combined <- c(1, 2)
sigma <- matrix(c(2, 1, 3, 2), ncol = 2L)

# replicates, but not the 10,000 in the question,
# only 1,000 is enough for code tests
R <- 1000L

# with a for loop:
# 1. create the results matrix beforehand
# 2. simulate from multivariate normal
# 3. run the test
# 4. and assign the results
set.seed(2023)
coeff_mat <- matrix(nrow = R, ncol = 2,
                    dimnames = list(NULL, c("Amyg_Avg", "Sys_Jus")))

for(i in seq.int(R)) {
  bnsamp <- mvrnorm(n = 48L, mu = mean_combined, Sigma = sigma)
  cor_coeff <- cor.test(x = bnsamp[, 1L], y = bnsamp[, 2L], 
                        method= "pearson", alternative= "two.sided")
  cor_coef_est <- cor_coeff$estimate
  cor_coef_p <- cor_coeff$p.value
  coeff_mat[i, ] <- c(cor_coef_est, cor_coef_p)
}

# with help('replicate')
# 1. simulate from multivariate normal
# 2. run the test
# 3. and assign the results
# 4. the final pipe to t() puts the matrix in the wanted form
set.seed(2023)
res <- replicate(R, {
  bnsamp <- mvrnorm(n = 48L, mu = mean_combined, Sigma = sigma)
  cor_coeff <- cor.test(x = bnsamp[, 1L], y = bnsamp[, 2L], 
                        method= "pearson", alternative= "two.sided")
  c(Amyg_Avg = unname(cor_coeff$estimate), Sys_Jus = cor_coeff$p.value)
}) |> t()

identical(res, coeff_mat)
#> [1] TRUE

head(res)
#>       Amyg_Avg      Sys_Jus
#> [1,] 0.3999366 4.856605e-03
#> [2,] 0.6007393 6.353057e-06
#> [3,] 0.6641624 2.652077e-07
#> [4,] 0.6094018 4.287226e-06
#> [5,] 0.4825884 5.132400e-04
#> [6,] 0.5514744 4.854691e-05

字符串
创建于2023-11-06

相关问题