r studio:模拟我的代码 1000 次并选择 p 值<0.05 的东西

r studio: simulate my code 1000 times and pick the things which p value<0.05

这是我的原始代码:

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 的那个?

我不会为你做这个,但这些是你可能想要完成的步骤:

  1. 将您的代码编写为函数,returns 您感兴趣的值(大概是 t)
  2. 多次使用类似replicate到运行此功能并记录所有答案
  3. 使用类似 quantile 的方法提取您感兴趣的百分位数

这是 tidyverse 及其列表列的完美应用程序。请参阅内联评论中的解释。

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

reprex package (v0.3.0)

于 2020-05-18 创建