在 for 循环中附加 t.test 结果的问题

如何解决在 for 循环中附加 t.test 结果的问题

让我用模拟数据集来解释：

我有数据集 dt 和 dt1

# dataset 1 `dt`
set.seed(12)
dt <- rnorm(5000,mean=10,sd=1)
dt <- data.frame(dt)
dt$group <- c("case","control")
colnames(dt) <- c("severity","group")
head(dt)

severity   group
1  8.519432    case
2 11.577169 control
3  9.043256    case
4  9.079995 control
5  8.002358    case
6  9.727704 control


# dataset 2 `dt2`
set.seed(12)
dt2 <- rnorm(200,mean=12,sd=1)
dt2 <- data.frame(dt2)
dt2$group <- c("case2","control2")
colnames(dt2) <- c("severity","group")
head(dt2)

  severity    group
1 10.51943    case2
2 13.57717 control2
3 11.04326    case2
4 11.07999 control2
5 10.00236    case2
6 11.72770 control2

我正在构建一个 1000 次迭代的 for 循环来执行以下步骤：

• 从 dt 中随机抽取 500 行并保存为 dt_sub
• 将 dt_sub 与 dt2 绑定并另存为 bd
• 仅从 bd 数据集中选择 group 为 case2 或 control 的行（只关心这两个组之间的差异）
• t. 对案例 2 和对照组之间可变严重程度的测试
• 将 t.tests 结果输出到 t
• 使用 for 循环重复 1000 次
• 迭代地将所有 t.test 结果附加到数据框 results

以下是我在r中构建的代码

library(broom)
library(dplyr)

iter <- 1000
t <- data.frame()
for (i in 1:iter) {
  dt_sub <- dt[sample(nrow(dt),500),]
  bd <- rbind(dt_sub,dt2)
  compare <- filter(bd,group %in% c("case2","control"))
  compare %>% group_by(group) %>% do(tidy(t.test(severity ~ group,data = compare))) -> t
  t$iter <- i
}

  results <- do.call(rbind,t)

我的问题是，这段代码在 iter=1 时运行良好，但是我应该如何设置 compare %>% group_by(group) %>% do(tidy(t.test(severity ~ group,data = compare))) -> t 行以确保当 iter ≥ 1 时不会覆盖每次运行的 t.test 结果？我试过 t[i] 但失败了，有人可以建议吗？

谢谢。

解决方法

创建一个运行该进程一次的函数。

library(broom)
library(dplyr)

t_test_function <- function() {
  dt_sub <- dt[sample(nrow(dt),500),]
  bd <- rbind(dt_sub,dt2)
  compare <- filter(bd,group %in% c("case2","control"))
  compare %>% 
    group_by(group) %>% 
    do(tidy(t.test(severity ~ group,data = compare))) %>%
    ungroup
}
t_test_function()

#  group estimate estimate1 estimate2 statistic  p.value parameter conf.low
#  <chr>    <dbl>     <dbl>     <dbl>     <dbl>    <dbl>     <dbl>    <dbl>
#1 case2     1.94      11.9      9.99      17.4 9.40e-42      199.     #1.72
#2 cont…     1.94      11.9      9.99      17.4 9.40e-42      199.     1.72
# … with 3 more variables: conf.high <dbl>,method <chr>,#   alternative <chr>

现在您可以使用 iter 调用 replicate 次并合并数据集。

iter <- 5
results <- bind_rows(replicate(iter,t_test_function(),simplify = FALSE),.id = 'iter')

# A tibble: 10 x 12
#   iter  group estimate estimate1 estimate2 statistic  p.value parameter
#   <chr> <chr>    <dbl>     <dbl>     <dbl>     <dbl>    <dbl>     <dbl>
# 1 1     case2     1.88      11.9     10.1       17.3 1.05e-40      189.
# 2 1     cont…     1.88      11.9     10.1       17.3 1.05e-40      189.
# 3 2     case2     1.96      11.9      9.97      17.8 9.88e-43      194.
# 4 2     cont…     1.96      11.9      9.97      17.8 9.88e-43      194.
# 5 3     case2     1.94      11.9      9.99      17.9 3.76e-42      184.
# 6 3     cont…     1.94      11.9      9.99      17.9 3.76e-42      184.
# 7 4     case2     2.03      11.9      9.90      18.6 1.82e-44      189.
# 8 4     cont…     2.03      11.9      9.90      18.6 1.82e-44      189.
# 9 5     case2     1.96      11.9      9.97      18.1 7.05e-43      187.
#10 5     cont…     1.96      11.9      9.97      18.1 7.05e-43      187.
# … with 4 more variables: conf.low <dbl>,conf.high <dbl>,#   alternative <chr>