R studio：模拟我的代码 1000 次，然后选择 p 值<0.05 的东西

Question

这是我的原始代码：

x = rbinom(1000,1,0.5)
z = log(1.3)*x       
pr = 1/(1+exp(-z)) 
y = rbinom(1000,1,pr)     
k=glm(y~x,family="binomial")$coef
t=exp(k)

如何模拟 1000 次并选择 p 值为 <0.05 的那个？

Answer 1

这是整洁宇宙的完美应用程序，它是列表列。请参阅内联注释中的说明。

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)) %>% 
# drop the intercepts and every .x with a p < 0.05
filter(term !="(Intercept)",
p.value < 0.05)
sim  
#> # A tibble: 545 x 6
#>       id term  estimate std.error statistic  p.value
#>    <int> <chr>    <dbl>     <dbl>     <dbl>    <dbl>
#>  1     3 .x       0.301     0.127      2.37 0.0176  
#>  2     7 .x       0.263     0.127      2.06 0.0392  
#>  3     8 .x       0.293     0.127      2.31 0.0211  
#>  4    11 .x       0.377     0.128      2.96 0.00312 
#>  5    12 .x       0.265     0.127      2.08 0.0373  
#>  6    13 .x       0.366     0.127      2.88 0.00403 
#>  7    14 .x       0.461     0.128      3.61 0.000305
#>  8    17 .x       0.274     0.127      2.16 0.0309  
#>  9    18 .x       0.394     0.127      3.09 0.00200 
#> 10    19 .x       0.371     0.127      2.92 0.00354 
#> # … with 535 more rows

^{创建于 2020-05-18 由 reprex 软件包 (v0.3.0(}

Answer 2

我不会为你做这个，但这些是你可能想要经历的步骤：

1. 将代码编写为返回您感兴趣的值的函数(大概是 t(
2. 使用类似replicate的东西多次运行此功能并记录所有答案
3. 使用类似quantile的内容来提取您感兴趣的百分位数