r - 在一个级别中置换一列,对 2 列执行测试,并保存 pvalues
问题描述
我有一个数据框
> dput(df)
structure(list(id = c(1, 2, 3, 4, 1, 2, 3, 4), level = structure(c(1L,
1L, 1L, 1L, 2L, 2L, 2L, 2L), .Label = c("g01", "g02"), class = "factor"),
m_col = c(1, 2, 3, 4, 11, 22, 33, 44), u_col = c(11, 12,
13, 14, 21, 22, 23, 24), group = c(0, 0, 1, 1, 0, 0, 1, 1
)), row.names = c(NA, -8L), class = "data.frame")
看起来像这样
id level m_col u_col group
1 1 g01 1 11 0
2 2 g01 2 12 0
3 3 g01 3 13 1
4 4 g01 4 14 1
5 1 g02 11 21 0
6 2 g02 22 22 0
7 3 g02 33 23 1
8 4 g02 44 24 1
我想对每个“级别”执行二项式加权测试(本质上,我需要比较每个 id 的 u_col 和 m_col)......所以使用tidyverse
并且broom
我可以执行以下操作:
res <- df %>%
group_by(level) %>%
do(tidy(glm(cbind(.$m_col,.$u_col) ~ .$group, family="binomial"))) %>%
filter(term == ".$group")
这给了我每个级别的一些 p 值:
> res
# A tibble: 2 x 6
# Groups: level [2]
level term estimate std.error statistic p.value
<fct> <chr> <dbl> <dbl> <dbl> <dbl>
1 g01 .$group 0.687 0.746 0.921 0.357
2 g02 .$group 0.758 0.296 2.56 0.0105
然后我可以问有多少 p<0.05
length(which(res$p.value < 0.05)
我现在想置换数据,重复二项式检验,询问有多少 p < 0.05,然后存储该值,然后再重复 999 次。
但是,排列需要对每个“级别”内的“组”列进行洗牌。我正在努力寻找一种方法来做到这一点,所以例如一个排列看起来像这样
id level m_col u_col group
1 1 g01 1 11 1
2 2 g01 2 12 0
3 3 g01 3 13 1
4 4 g01 4 14 0
5 1 g02 11 21 1
6 2 g02 22 22 0
7 3 g02 33 23 1
8 4 g02 44 24 0
一秒钟看起来像
id level m_col u_col group
1 1 g01 1 11 0
2 2 g01 2 12 1
3 3 g01 3 13 1
4 4 g01 4 14 0
5 1 g02 11 21 0
6 2 g02 22 22 1
7 3 g02 33 23 1
8 4 g02 44 24 0
ETC
让测试依赖于 2 列限制了随机播放选项,我很难过。我会很感激任何建议。
解决方案
如果你想要一个数据框,你可以试试这个:
library(tidyverse)
map_dfr(1:1000, ~ df %>%
group_by(level) %>%
mutate(group = group[sample(row_number())]) %>% # permutation shuffle the 'group' column within each 'level'.
do(tidy(glm(cbind(.$m_col,.$u_col) ~ .$group, family="binomial"))) %>%
filter(term == ".$group") %>%
ungroup() %>%
summarise(sum(p.value < 0.05))) # ask how many p<0.05
如果你想要一个向量:
map_dbl(1:1000, ~ df %>%
group_by(level) %>%
mutate(group = group[sample(row_number())]) %>% # permutation shuffle the 'group' column within each 'level'.
do(tidy(glm(cbind(.$m_col,.$u_col) ~ .$group, family="binomial"))) %>%
filter(term == ".$group") %>%
ungroup() %>%
summarise(sum(p.value < 0.05)) %>% # ask how many p<0.05
pull())
推荐阅读
- android - 是否可以在 android.defaultConfig 中为 versionCod 设置“字符串”占位符?
- caching - 如何实现动态索引?
- file - 通过 telnet 从 IP 摄像机下载文件
- php - 如何正确地将字符串形式的 post 变量转换为对象或数组
- java - 谷歌地图:当用户足够近时,如何使标记仅可点击?
- java - 按集合中的条件过滤实体
- angular - 地图在强制刷新或鼠标点击之前不会更新
- python - 将 PCollection 分配回全局窗口
- list - SwiftUI:防止视图在 VStack 内展开
- python - Python如何为时间范围设置多个计划并以设定的时间间隔执行