2016-12-04 119 views
0

提取对数值(分数)与p值的测试结果我有100个重复的coxph模型拟合循环。我试图用数据框或列表中的每个重复项的p值提取出日志排名分数测试结果。我正在使用以下内容。但是,它只给出了我的日志排名分数,而不是p值。任何帮助将非常感激。Coxph模型

我可以共享数据集,但不知道如何在此处添加。

感谢, Krina

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))) 

回答

0

下面是我从coxph功能的帮助页面了一个模拟例子。我只是将数据集复制了100次以创建您的场景。我强烈建议开始使用tidyverse包来完成这样的工作。 broomdplyrtidyr是一个很好的补充。

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() 

维杰

+0

太感谢了,维杰。 –