我怎样才能循环这个模拟 1000 次并使其中一列一起成为一个矩阵？

Question

library(tidyverse)
library(broom)

# create a tibble with an id column for each simulation and x wrapped in list()
sim <- tibble(id = 1:1000,
               x = list(rbinom(1000,1,0.5))) %>% 
# to generate z, pr, y, k use map and map2 from the purrr package to loop over the list column x
# `~ ... ` is similar to `function(.x) {...}`
# `.x` represents the variable you are using map on
          mutate(z  = map(x, ~ log(1.3) * .x), 
                 pr = map(z, ~ 1 / (1 + exp(-.x))),
                 y  = map(pr, ~ rbinom(1000, 1, .x)),
                 k  = map2(x, y, ~ glm(.y ~ .x, family="binomial")),
# use broom::tidy to get the model summary in form of a tibble
                 sum = map(k, broom::tidy)) %>% 
# select id and sum and unnest the tibbles
  select(id, sum) %>% 
  unnest(cols = c(sum)) %>% 
simAll <- sim %>% 
  filter(term !="(Intercept)")

simAll 是这样的：

   id  term  estimate       std.error   statistic   p.value
1   1   .x  0.4058039189    0.1275272   3.182096892 1.462129e-03
2   2   .x  0.2515178701    0.1276719   1.970033693 4.883451e-02
3   3   .x  0.2464097082    0.1274321   1.933654251 5.315565e-02
4   4   .x  0.2308803864    0.1273598   1.812819663 6.985964e-02
5   5   .x  0.3029238760    0.1271623   2.382182816 1.721035e-02
6   6   .x  0.2452264719    0.1270829   1.929657417 5.364930e-02
7   7   .x  0.2390919312    0.1270123   1.882430831 5.977754e-02
8   8   .x  0.2437134055    0.1271373   1.916930426 5.524677e-02
9   9   .x  0.4372744410    0.1274612   3.430646232 6.021453e-04
10  10  .x  0.2915176118    0.1272609   2.290708545 2.198028e-02
11  11  .x  0.3373491310    0.1271283   2.653612132 7.963531e-03
12  12  .x  0.1991820874    0.1269380   1.569128570 1.166180e-01
13  13  .x  0.3437529981    0.1272502   2.701394595 6.904936e-03
14  14  .x  0.2229632179    0.1269851   1.755822253 7.911876e-02
15  15  .x  0.2606269011    0.1271385   2.049944533 4.036984e-02
Showing 1 to 15 of 1,000 entries, 6 total columns

这里的问题是，这里的estimate列是x的值，我想有1000个类似simALL的表（比如重复整个模拟1000次），然后我会有1000 * 1000 x，我要让它们成为一个矩阵（1000 * 1000），我该怎么办？

Answer 1

您可以为要重复的代码编写一个函数，并仅返回 estimate 列，因为这是您需要的唯一信息。

library(tidyverse)
run_fun <- function() {

sim <- tibble(id = 1:1000,
              x = list(rbinom(1000,1,0.5))) %>% 
          mutate(z  = map(x, ~ log(1.3) * .x), 
                 pr = map(z, ~ 1 / (1 + exp(-.x))),
                 y  = map(pr, ~ rbinom(1000, 1, .x)),
                 k  = map2(x, y, ~ glm(.y ~ .x, family="binomial")),
                 sum = map(k, broom::tidy)) %>%
           select(id, sum) %>% 
           unnest(cols = c(sum)) %>%
           filter(term !="(Intercept)") %>% 
           pull(estimate)
  return(sim)
}

你可以调用这个函数n次：

data <- replicate(1000, run_fun())