下面是我从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)
nested <- r %>%
group_by(rep) %>%
nest %>%
mutate(model = data %>% map(cxph_mod))
外表到第一代表看到coxph输出类似于lapply
的功能。您将看到模型对象存储在数据框的单元格中,以便于访问。
nested %>% filter(rep==1)
每个模型对象,现在用扫帚来从模型中的参数估计和预测到嵌套的数据集
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个不同的da tasets因此有100个不同的价值。
ggplot(data=ests,aes(y=p.value,x=rep))+geom_point()
维杰
太感谢了,维杰。 –