R、模拟、p 值、直方图答案

【问题标题】：R, Simulation, p-value, histogramR、模拟、p 值、直方图
【发布时间】：2020-06-14 21:06:48
【问题描述】：

library(tidyverse)
library(broom)
library(dplyr)
# 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
j=exp(sim %>% select("estimate"))
OR=as.numeric(unlist(j))
mean(OR)

hist(OR,main=NULL,freq=T,breaks=10)
abline(v=mean(OR),lwd=4,col=1)

这里的问题：现在我提取所有 p

【问题讨论】：

1) 在unnest(cols = c(sum)) 之后立即停止管道； 2）创建一个simOR，就像你继续管道和一个simAll，但这次不过滤p值。

标签： r simulation

【解决方案1】：

这个解决方案重复了问题的代码，但是

在unnest(cols = c(sum))之后立即停止管道；
创建一个 simOR 就像您继续管道和一个 simAll 但这次不过滤 p 值。

首先是问题的代码。请注意，如果已加载包tidyverse，则无需加载包dplyr。
我还设置了 RNG 种子以使结果可重现。

library(tidyverse)
library(broom)
# create a tibble with an id column for each simulation and x wrapped in list()
set.seed(2020)
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))

现在创建两个要绘制的数据集。

simOR <- sim %>% 
  # drop the intercepts and every .x with a p < 0.05
  filter(term !="(Intercept)", p.value < 0.05)

j <- exp(simOR %>% select("estimate"))
OR <- as.numeric(unlist(j))
mean(OR)

所有行的数据集，只删除截距。

simAll <- sim %>% 
  filter(term !="(Intercept)")


j <- exp(simAll %>% select("estimate"))
All <- as.numeric(unlist(j))
mean(All)

现在绘制直方图（不重叠）。

op <- par(mfrow = c(2, 1))
hist(OR, main = NULL, freq = TRUE, breaks = 10)
abline(v = mean(OR), lwd = 4, col = 1)
hist(All, main = NULL, freq = TRUE, breaks = 10)
abline(v = mean(All), lwd = 4, col = 1)
par(op)

【讨论】：