如何解决我如何循环进行此模拟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())
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。