R：通过Permute或Sample函数进行蒙特卡洛过程以生成空分布

问题描述

根据该数据集，我进行了聚类分析分配的所有患者样品（共69行），并且聚类标记为第3列“ Cluster.assigned”，共8个聚类，每个聚类的大小为UNEQUAL。其他列包含变量，我想测试其中的数字变量（例如Age），以查看随机变量与随机变量相比是否富集。

由于我的编码能力，我现在遇到了障碍。但是我的想法是将实际数据视为 Observed ，然后通过使用样本或置换函数（如Monte Carlo模拟）将簇的标签改组为1000次，并将模拟分布称为预期。

以“年龄”列为例：

#minimum dummy 30-row data
Patient.ID <-c("S3077497","S1041120","S162465","S563275","S2911623","S3117192","S2859024","S2088278","S3306185","S190789","S12146451","S2170842","S115594","S2024203","S1063872","S2914138","S303984","S570813","S2176683","S820460","S1235729","S3009401","S2590229","S629309","S120256","S2572773","S3180483","S3032079","S3217608","S5566943")

Cluster.assigned <- c("cluster1","cluster1","cluster2","cluster3","cluster4","cluster4")

Age <- c(61,80,78,69,57,70,60,59,72,82,66,68,62,67,77,74,64,54,73,58,87)

CLL_3S <-cbind(Patient.ID,Cluster.assigned,Age)

要查看是否有某个集群可以使患者在特定年龄富裕，则零假设是各个集群之间的年龄分布没有差异。 现在，我应该将患者标签改组或将年龄数据改组为1000次，然后应该有一个模拟的数据框，从中可以计算出模拟的（期望值）的均值和标准差

#I image to use shuffle to permute 1000 times
#And combine the simulated into a massive dataframe
 shuffled <- numeric(length=1000)
 N <-nrows(CLL_3S)

 set.seed(123)
  for (i in seq_len(length(shuffled) -1)) {
      perm <- shuffle(N)
      .........

下一步是，我将使用每个聚类中患者年龄的实际观察值，通过Z分数来计算富集度。说obs（值-预期均值）/ SD。

一旦该过程自动化，那么我可以将其应用于其他感兴趣的列和具有不同数量簇的其他数据集。我已经阅读了有关sample（）和shuffle（）的内容，但这并不能真正帮助我解决这个特定问题...

解决方法

我不确定下面的代码是否满足您的目标。如果我正确理解了您的问题，那么我应该做的是只对群集分配进行洗牌，然后添加一个按群集标签分组的z得分新列。

• sample进行随机混洗
• scale用于计算z得分
• ave有助于按聚类标签计算z得分
• replicate将多次运行模拟

replicate(1000,within(
    transform(CLL_3S,Cluster.assigned = Cluster.assigned[sample(1:nrow(CLL_3S))]
    ),zscore <- ave(Age,Cluster.assigned,FUN = scale)
  ),simplify = FALSE
)

---

更新

如果您只想对1000个模拟的平均值和标准差进行平均，则可以尝试以下代码

n <- 1000
res <- Reduce(
  `+`,replicate(n,with(
      CLL_3S,do.call(rbind,tapply(Age,Cluster.assigned[sample(1:nrow(CLL_3S))],FUN = function(x) c(Mean = mean(x),Var = var(x))))
    ),simplify = FALSE
  )
) / n
res <- within(as.data.frame(res),SD <- sqrt(Var))

给出

> res
             Mean      Var       SD
cluster1 70.21086 68.99152 8.306114
cluster2 70.06915 71.93188 8.481267
cluster3 70.03571 70.19276 8.378112
cluster4 70.12500 68.98867 8.305942

数据

> dput(CLL_3S)
structure(list(Patient.ID = c("S3077497","S1041120","S162465","S563275","S2911623","S3117192","S2859024","S2088278","S3306185","S190789","S12146451","S2170842","S115594","S2024203","S1063872","S2914138","S303984","S570813","S2176683","S820460","S1235729","S3009401","S2590229","S629309","S120256","S2572773","S3180483","S3032079","S3217608","S5566943"),Cluster.assigned = c("cluster1","cluster1","cluster2","cluster3","cluster4","cluster4"),Age = c(61,80,78,69,57,70,60,59,72,82,66,68,62,67,77,74,64,54,73,58,87)),class = "data.frame",row.names = c(NA,-30L))