R、模拟、p 值、直方图


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 值进行比较,哪个代码可以使用它?

此解决方案重复了问题的代码,但

  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)

最新更新