FDR 校正 - 从 lmer() 中提取 p 值并创建用于 R 中 p.adjust 的向量
FDR correction - extracting p-values from lmer() and creating vectors for use in p.adjust in R
我正在尝试对某些感兴趣区域的神经影像数据进行 FDR 校正。我总共有 运行 18 个线性混合效应模型,并且我确保输出中系数的顺序在每个模型中都是相同的。
我已将每个模型的输出保存在以下内容中:
tidy_model1 <-tidy(model1)
tidy_model2 <-tidy(model2)
....
tidy_model18 <-tidy(model18)
我现在正在努力让我的生活更轻松,并创建一个循环,该循环遍历包含上述模型对象名称的列表,并为每个系数创建一个 p 值向量,然后我将在 p.adjust 函数检索调整后的 p 值。
所以我创建了一个列表:
model_list <- list(tidy_model1,
tidy_model2,... tidy_model18)
我试过以下循环:
for (i in 1:18) {
model_list[i] %>%
variable1_pval <- p.value[1]
}
和
for (i in 1:18) {
variable1_pval <- model_list[i]$p.value[1]
}
所以上面应该给出模型系数 1 的 p 值向量。
但是,在这两种情况下我都得到了一个空向量。
我知道我没有提供我的数据,但欢迎就这些循环可能无法正常工作的原因提出任何建议!
谢谢
我列了一个模型列表:
library(nlme)
library(broom)
models <- lapply(1:5,function(i){
idx= sample(nrow(Orthodont),replace=TRUE)
lme(distance ~ age, random=~Sex,data = Orthodont[idx,])
})
model_list <- lapply(models,tidy,effects="fixed")
在这些模型中,有用的系数是第二个:
model_list[[1]]
# A tibble: 2 x 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 15.9 1.03 15.5 7.77e-26
2 age 0.739 0.0871 8.48 9.13e-13
您可以像这样在向量中获取 p 值,例如使用 p.value1:
sapply(model_list,function(x)x$p.value[2])
跟踪模型而不用变量填充环境的更好方法是使用 purrr、dplyr(查看更多 here):
library(purrr)
library(dplyr)
models = tibble(name=1:5,models=models) %>%
mutate(tidy_res = map(models,tidy,effects="fixed"))
models
# A tibble: 5 x 3
name models tidy_res
<int> <list> <list>
1 1 <lme> <tibble [2 × 5]>
2 2 <lme> <tibble [2 × 5]>
3 3 <lme> <tibble [2 × 5]>
4 4 <lme> <tibble [2 × 5]>
5 5 <lme> <tibble [2 × 5]>
models %>% unnest(tidy_res) %>% filter(term=="age")
# A tibble: 5 x 7
name models term estimate std.error statistic p.value
<int> <list> <chr> <dbl> <dbl> <dbl> <dbl>
1 1 <lme> age 0.587 0.0601 9.77 2.44e-15
2 2 <lme> age 0.677 0.0663 10.2 3.91e-16
3 3 <lme> age 0.588 0.0603 9.74 3.05e-15
4 4 <lme> age 0.653 0.0529 12.3 2.74e-20
5 5 <lme> age 0.638 0.0623 10.2 3.34e-16
我正在尝试对某些感兴趣区域的神经影像数据进行 FDR 校正。我总共有 运行 18 个线性混合效应模型,并且我确保输出中系数的顺序在每个模型中都是相同的。
我已将每个模型的输出保存在以下内容中:
tidy_model1 <-tidy(model1)
tidy_model2 <-tidy(model2)
....
tidy_model18 <-tidy(model18)
我现在正在努力让我的生活更轻松,并创建一个循环,该循环遍历包含上述模型对象名称的列表,并为每个系数创建一个 p 值向量,然后我将在 p.adjust 函数检索调整后的 p 值。
所以我创建了一个列表:
model_list <- list(tidy_model1,
tidy_model2,... tidy_model18)
我试过以下循环:
for (i in 1:18) {
model_list[i] %>%
variable1_pval <- p.value[1]
}
和
for (i in 1:18) {
variable1_pval <- model_list[i]$p.value[1]
}
所以上面应该给出模型系数 1 的 p 值向量。
但是,在这两种情况下我都得到了一个空向量。
我知道我没有提供我的数据,但欢迎就这些循环可能无法正常工作的原因提出任何建议!
谢谢
我列了一个模型列表:
library(nlme)
library(broom)
models <- lapply(1:5,function(i){
idx= sample(nrow(Orthodont),replace=TRUE)
lme(distance ~ age, random=~Sex,data = Orthodont[idx,])
})
model_list <- lapply(models,tidy,effects="fixed")
在这些模型中,有用的系数是第二个:
model_list[[1]]
# A tibble: 2 x 5
term estimate std.error statistic p.value
<chr> <dbl> <dbl> <dbl> <dbl>
1 (Intercept) 15.9 1.03 15.5 7.77e-26
2 age 0.739 0.0871 8.48 9.13e-13
您可以像这样在向量中获取 p 值,例如使用 p.value1:
sapply(model_list,function(x)x$p.value[2])
跟踪模型而不用变量填充环境的更好方法是使用 purrr、dplyr(查看更多 here):
library(purrr)
library(dplyr)
models = tibble(name=1:5,models=models) %>%
mutate(tidy_res = map(models,tidy,effects="fixed"))
models
# A tibble: 5 x 3
name models tidy_res
<int> <list> <list>
1 1 <lme> <tibble [2 × 5]>
2 2 <lme> <tibble [2 × 5]>
3 3 <lme> <tibble [2 × 5]>
4 4 <lme> <tibble [2 × 5]>
5 5 <lme> <tibble [2 × 5]>
models %>% unnest(tidy_res) %>% filter(term=="age")
# A tibble: 5 x 7
name models term estimate std.error statistic p.value
<int> <list> <chr> <dbl> <dbl> <dbl> <dbl>
1 1 <lme> age 0.587 0.0601 9.77 2.44e-15
2 2 <lme> age 0.677 0.0663 10.2 3.91e-16
3 3 <lme> age 0.588 0.0603 9.74 3.05e-15
4 4 <lme> age 0.653 0.0529 12.3 2.74e-20
5 5 <lme> age 0.638 0.0623 10.2 3.34e-16