我进行了模拟,其中创建了 2 组,每组 100 名患者。数据框包含患者、干预类型(R 用于参考,T 用于测试)和恢复状态(1 或 0)。然后我需要做一个实验:为概率差异制定置信区间。 我需要重复实验 1000 次(每次实验都会创建一个新的患者样本组)并在每次重复中保存结果(df_prop2)。最后我应该得到一个包含 1000 行的数据框(df_prop2)(7 个变量的 1000 个观测值)。我不明白为什么使用我的代码我得到的数据框有 1 行和 1007 列。请帮忙
这是我的代码:
library(dplyr)
library(tidyr)
set.seed(123)
p_R <- 0.1 #true probability for reference
p_T <- 0.2 #true probabiluty for Test
sample_size <- 100 #I have 100 patients in 2 groups
n_rep <- 1000 #number of repeats
df_all_repeats <- data.frame(
n_exp = rep(1:n_rep, each = 2*sample_size),
arm = rep(c('R', 'T'), each = sample_size),
patient_ID = rep(1:sample_size, n_rep))
df_all_repeats <- df_all_repeats %>%
group_by(n_exp) %>%
mutate(recovery_status = c(sample(c(1,0), sample_size, replace = TRUE, prob = c(p_R, 1 - p_R)),
sample(c(1,0), sample_size, replace = TRUE, prob = c(p_T, 1 - p_T))))
#Experiment
df_prop2 <- df_all_repeats %>%
group_by(arm) %>%
group_by(n_exp) %>%
dplyr::summarise(x = sum(recovery_status),
n = n()) %>%
ungroup() %>%
dplyr::summarise(X = list(x), N = list(n)) %>%
rowwise() %>%
mutate(tst = list(broom::tidy(prop.test(X, N)))) %>%
unnest(tst)
print(df_prop2)
这是一个也简化了数据生成的解决方案:
library(dplyr)
library(tidyr)
library(broom)
p_R <- 0.1 # true probability for reference
p_T <- 0.2 # true probabiluty for Test
sample_size <- 100 # number of patients
n_rep <- 1000 # number of repeats
set.seed(123)
df_all_repeats <- expand.grid(
patient_id = 1:sample_size,
arm = c("R", "T"),
n_exp = 1:n_rep,
stringsAsFactors = FALSE) %>%
as_tibble() %>% ## data is sorted n_exp, arm, patient_id
mutate(recovery_status = c(sample(1:0, sample_size, replace = TRUE,
prob = c(p_R, 1 - p_R)),
sample(1:0, sample_size, replace = TRUE,
prob = c(p_T, 1 - p_T))),
.by = n_exp)
df_all_repeats %>%
group_by(n_exp, arm) %>%
summarise(x = sum(recovery_status),
n = n()) %>%
summarise(broom::tidy(prop.test(x, n)))
# # A tibble: 1,000 × 10
# n_exp estimate1 estimate2 statistic p.value parameter conf.low conf.high method alternative
# <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
# 1 1 0.07 0.19 5.35 0.0207 1 -0.222 -0.0183 2-sample test for equality of proportions with continuity corre… two.sided
# 2 2 0.1 0.19 2.58 0.108 1 -0.197 0.0168 2-sample test for equality of proportions with continuity corre… two.sided
# 3 3 0.09 0.27 9.79 0.00175 1 -0.294 -0.0665 2-sample test for equality of proportions with continuity corre… two.sided
# 4 4 0.08 0.21 5.81 0.0160 1 -0.236 -0.0241 2-sample test for equality of proportions with continuity corre… two.sided
# 5 5 0.07 0.21 7.02 0.00807 1 -0.244 -0.0358 2-sample test for equality of proportions with continuity corre… two.sided
# 6 6 0.03 0.18 10.4 0.00124 1 -0.242 -0.0576 2-sample test for equality of proportions with continuity corre… two.sided
# 7 7 0.08 0.16 2.32 0.128 1 -0.179 0.0194 2-sample test for equality of proportions with continuity corre… two.sided
# 8 8 0.12 0.15 0.171 0.679 1 -0.135 0.0746 2-sample test for equality of proportions with continuity corre… two.sided
# 9 9 0.08 0.28 12.2 0.000471 1 -0.313 -0.0872 2-sample test for equality of proportions with continuity corre… two.sided
# 10 10 0.11 0.13 0.0473 0.828 1 -0.120 0.0800 2-sample test for equality of proportions with continuity corre… two.sided
# # ℹ 990 more rows
# # ℹ Use `print(n = ...)` to see more rows