我如何循环进行此模拟1000次，然后将其中一列合并为矩阵？

问题描述

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,.x)),k  = map2(x,~ 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

这里的问题是这里的估计列是x的值，我想有1000个类似simALL的表（例如将整个模拟重复1000次），然后我将有1000 * 1000 x，我想将它们设为矩阵（1000 * 1000），该怎么办？

解决方法

您可以为要重复的代码编写一个函数，并仅返回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,.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())