R，模拟，p值，直方图

Question

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<0.05的所有值，现在我使用代码"hist(OR，main=NULL，freq=T，breaks=10)“来制作赔率比的直方图。我想做的新事情是让另一个直方图(比如没有任何条件的p值)重叠在原来的直方图上，然后我可以在一个图中将直方图与不同的p值进行比较，哪个代码可以使用它？

Answer 1

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

1. 紧跟在unnest(cols = c(sum))之后停止管道；
2. 创建一个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)