使用 Coxph 模型的 p 值提取对数秩(分数)测试结果
extract log rank (score) test result wiht p-value for Coxph Model
我在循环中安装了 100 个 coxph 模型副本。我正在尝试为数据框或列表中的每个副本提取具有 p 值的对数秩分数测试结果。我正在使用以下内容。但是,它只给我对数排名分数,而不是 p 值。任何帮助将不胜感激。
我可以共享数据集,但不确定如何附加到这里。
谢谢,
克里纳
Repl_List <- unique(dat3$Repl)
doLogRank = function(sel_name) {
dum <- dat3[dat3$Repl == sel_name,]
reg <- with(dum, coxph(Surv(TIME_day, STATUS) ~ Treatment, ties = "breslow"))
LogRank <- with(reg, reg$score)
}
LogRank <- t(as.data.frame(lapply(Repl_List, doLogRank)))
这是我从 coxph 函数的帮助页面中获取的模拟示例。我只是将数据集复制了 100 次来创建您的场景。我强烈建议开始使用 tidyverse
包来完成此类工作。 broom
是 dplyr
和 tidyr
的重要补充。
library(survival)
library(tidyverse)
library(broom)
test <- data.frame(time=c(4,3,1,1,2,2,3),
status=c(1,1,1,0,1,1,0),
x=c(0,2,1,1,1,0,0),
sex=c(0,0,0,0,1,1,1))
下面我使用 replicate
函数将数据集复制了 100 次。
r <- replicate(test,n = 100,simplify = FALSE) %>% bind_rows %>%
mutate(rep = rep(seq(1,100,1),each=7))
我将 cox 模型设置为一个小函数,我可以将其传递给数据帧的每个副本。
cxph_mod <- function(df) {
coxph(Surv(time, status) ~ x + strata(sex), df)
}
下面是拟合模型和提取值的分步过程。
tidyr::nest
数据框
purrr::map
将模型放入每个巢中
nest
是 library(tidyr)
中的函数
map
是library(purrr)
中类似于lapply
的函数
nested <- r %>%
group_by(rep) %>%
nest %>%
mutate(model = data %>% map(cxph_mod))
查看第一个代表以查看 coxph 输出。您将看到模型对象存储在数据框的单元格中,以便于访问。
nested %>% filter(rep==1)
对于每个模型对象,现在使用 broom 将模型的参数估计和预测放入嵌套数据集中
nested <- nested %>%
mutate(
ests = model %>% map(broom::tidy)
)
tidyr::unnest
查看您对每个重采样数据集的拟合预测
ests <- unnest(nested,ests,.drop=TRUE) %>% dplyr::select(rep,estimate:conf.high)
在这种情况下,因为我重复同一个数据集 100 次,pvalue 将是相同的,但在您的情况下,您将有 100 个不同的数据集,因此有 100 个不同的 p.values。
ggplot(data=ests,aes(y=p.value,x=rep))+geom_point()
维杰
我在循环中安装了 100 个 coxph 模型副本。我正在尝试为数据框或列表中的每个副本提取具有 p 值的对数秩分数测试结果。我正在使用以下内容。但是,它只给我对数排名分数,而不是 p 值。任何帮助将不胜感激。
我可以共享数据集,但不确定如何附加到这里。
谢谢, 克里纳
Repl_List <- unique(dat3$Repl)
doLogRank = function(sel_name) {
dum <- dat3[dat3$Repl == sel_name,]
reg <- with(dum, coxph(Surv(TIME_day, STATUS) ~ Treatment, ties = "breslow"))
LogRank <- with(reg, reg$score)
}
LogRank <- t(as.data.frame(lapply(Repl_List, doLogRank)))
这是我从 coxph 函数的帮助页面中获取的模拟示例。我只是将数据集复制了 100 次来创建您的场景。我强烈建议开始使用 tidyverse
包来完成此类工作。 broom
是 dplyr
和 tidyr
的重要补充。
library(survival)
library(tidyverse)
library(broom)
test <- data.frame(time=c(4,3,1,1,2,2,3),
status=c(1,1,1,0,1,1,0),
x=c(0,2,1,1,1,0,0),
sex=c(0,0,0,0,1,1,1))
下面我使用 replicate
函数将数据集复制了 100 次。
r <- replicate(test,n = 100,simplify = FALSE) %>% bind_rows %>%
mutate(rep = rep(seq(1,100,1),each=7))
我将 cox 模型设置为一个小函数,我可以将其传递给数据帧的每个副本。
cxph_mod <- function(df) {
coxph(Surv(time, status) ~ x + strata(sex), df)
}
下面是拟合模型和提取值的分步过程。
tidyr::nest
数据框
purrr::map
将模型放入每个巢中
nest
是 library(tidyr)
中的函数
map
是library(purrr)
lapply
的函数
nested <- r %>%
group_by(rep) %>%
nest %>%
mutate(model = data %>% map(cxph_mod))
查看第一个代表以查看 coxph 输出。您将看到模型对象存储在数据框的单元格中,以便于访问。
nested %>% filter(rep==1)
对于每个模型对象,现在使用 broom 将模型的参数估计和预测放入嵌套数据集中
nested <- nested %>%
mutate(
ests = model %>% map(broom::tidy)
)
tidyr::unnest
查看您对每个重采样数据集的拟合预测
ests <- unnest(nested,ests,.drop=TRUE) %>% dplyr::select(rep,estimate:conf.high)
在这种情况下,因为我重复同一个数据集 100 次,pvalue 将是相同的,但在您的情况下,您将有 100 个不同的数据集,因此有 100 个不同的 p.values。
ggplot(data=ests,aes(y=p.value,x=rep))+geom_point()
维杰