如何从数据子集中随机抽取并在 R 中引导统计测试

Posted 2023-03-24

【中文标题】如何从数据子集中随机抽取并在 R 中引导统计测试

【英文标题】：How to randomly draw from subsets of data and bootstrap a statistic test in R

【发布时间】：2013-12-06 05:04:44

【问题描述】：

我有一个包含两个变量的数据集，我希望在引导循环中统计测试它们是否相关（即使用 Spearman 的秩校正和cor.test(...)）。

我的数据集中的大多数测量值都来自独立的样本单位（我们称单位为植物），尽管有些测量值来自同一个植物。为了处理伪复制问题，我希望多次引导统计测试，在每次测试运行中只使用每个工厂的一次测量。因此，我需要编写一个引导循环，在执行相关测试之前为每个植物随机抽取一个测量值（然后重复此过程 99 次）。

我希望最终得到一个 csv 文件，其中包含 99 个测试中每一个测试的 p 值、rho 和 S 统计数据。

示例数据：

dput(df)

structure(list(Plant = c(1L, 2L, 3L, 4L, 5L, 6L, 6L, 7L, 8L, 
9L, 10L, 10L, 10L, 11L, 12L, 13L, 14L, 15L, 16L, 17L, 18L, 18L, 
19L, 20L, 21L), Length = c(170L, 232L, 123L, 190L, 112L, 207L, 
93L, 291L, 178L, 206L, 141L, 257L, 304L, 222L, 279L, 192L, 101L, 
253L, 176L, 278L, 311L, 129L, 191L, 205L, 226L), Count = c(7L, 
9L, 5L, 7L, 5L, 6L, 2L, 10L, 6L, 7L, 4L, 8L, 11L, 7L, 8L, 5L, 
5L, 9L, 7L, 6L, 9L, 4L, 5L, 7L, 6L)), .Names = c("Plant", "Length", 
"Count"), class = "data.frame", row.names = c(NA, -25L))


   Plant Length Count
1      1    170     7
2      2    232     9
3      3    123     5
4      4    190     7
5      5    112     5
6      6    207     6
7      6     93     2   
8      7    291    10  etc....

到目前为止，我已经整理了以下代码，首先为由多行表示的每个植物随机绘制一行，然后在运行统计测试之前将这些值与其余数据结合起来。但是，我现在正在努力合并一个引导功能（即boot() 或bootstrap()）来运行统计测试并多次执行循环：

# 1. create dataframe without plants with >1 measurement/row (in this example plant 6,10 & 18 have multiple rows)
df_uniq = df[ ! df$Plant %in% c(6,10,18), ]

# 2. create data subsets for each plant with >1 measurement/row
dup1 = df[6:7,]
dup2 = df[11:13,] 
dup3 = df[21:22,]

# 3. randomly draw one row for each plant with multiple measurements
d1_draw = dup1[sample(nrow(dup1), 1), ]
d2_draw = dup2[sample(nrow(dup2), 1), ]
d3_draw = dup3[sample(nrow(dup3), 1), ]

# 4. merge df_uniq with randomly drawn rows for each plant with multiple measurements
df_merge = rbind(df_uniq, d1_draw, d2_draw, d3_draw)

# 5. Test whether the two variables (length & Count) are related and write results to file
cor_res <- cor.test(df_merge$Length, df_merge$Count, method= "spearman")
write.csv(matrix(c(cor_res$statistic, cor_res$p.value, cor_res$estimate)), row.names=c("statistic", "p.value", "rho"), "test_output.csv")

我确信有一种快速而优雅的方法可以解决问题。任何帮助将不胜感激！非常感谢。

【参考方案1】：

为什么首先要提取唯一的行？如果只有一行，则对该植物进行一次抽样将导致该行保持不变，但仍会从多行的植物中随机抽样。

所以你可以这样做：

set.seed(123)
library(plyr)
ddply(df, .(Plant), function(x)  y <- x[sample(nrow(x), 1) ,] )

#   Plant Length Count howmany
#1      1    170     7       1
#2      2    232     9       1
#3      3    123     5       1
#4      4    190     7       1
#5      5    112     5       1
#6      6    207     6       2
#7      7    291    10       1
#8      8    178     6       1
#9      9    206     7       1
#10    10    257     8       3
#11    11    222     7       1
#12    12    279     8       1
#13    13    192     5       1
#14    14    101     5       1
#15    15    253     9       1
#16    16    176     7       1
#17    17    278     6       1
#18    18    311     9       2
#19    19    191     5       1
#20    20    205     7       1
#21    21    226     6       1

和你的cor.test

# first create your own function:
myrandomcors <- function(P)
ss <- ddply(P, .(Plant), function(x)  y <- x[sample(nrow(x), 1) ,] )
cor_res <- cor.test(ss$Length, ss$Count, method= "spearman")
return(c(stat = cor_res$statistic, p = cor_res$p.value, est = cor_res$estimate))


# then repeat it 5 times...
answer <- do.call( rbind, replicate(5, myrandomcors(df), simplify=FALSE ) )

#    > answer
#       stat.S            p   est.rho
#[1,] 352.4557 4.275291e-05 0.7711327
#[2,] 461.2733 4.060286e-04 0.7004719
#[3,] 340.2024 3.159626e-05 0.7790893
#[4,] 368.3967 6.227648e-05 0.7607814
#[5,] 342.4391 3.341956e-05 0.7776369

【讨论】：

@user1317221：非常感谢关于二次抽样的建议。关于如何将其与统计测试结合到引导循环中的任何想法？

@user1317221。太棒了——正是我想要的！感谢您的帮助！