Извлечь результаты теста ранг (балл) с p-значением для модели Кокша
У меня есть 100 копий модели Coxph, встроенной в петлю. Я пытаюсь извлечь результаты теста оценки лог-ранга с 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)))
1 ответ
Вот пример, который я взял со страницы помощи функции 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))
Ниже я повторяю набор данных 100 раз, используя replicate
функция.
r <- replicate(test,n = 100,simplify = FALSE) %>% bind_rows %>%
mutate(rep = rep(seq(1,100,1),each=7))
Я настроил модель Кокса как небольшую функцию, которую я могу передать каждой копии кадра данных.
cxph_mod <- function(df) {
coxph(Surv(time, status) ~ x + strata(sex), df)
}
Ниже приведен пошаговый процесс подгонки модели и извлечения значений.
tidyr::nest
кадр данных purrr::map
модель в каждое гнездоnest
функция в library(tidyr)
map
функция похожа на lapply
в library(purrr)
nested <- r %>%
group_by(rep) %>%
nest %>%
mutate(model = data %>% map(cxph_mod))
загляните в первый повтор, чтобы увидеть вывод coxph. Вы увидите объект модели, хранящийся в ячейках кадра данных, что облегчает доступ.
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 раз, значение будет одинаковым, но в вашем случае у вас будет 100 различных наборов данных и, следовательно, 100 различных значений p.value.
ggplot(data=ests,aes(y=p.value,x=rep))+geom_point()
Виджай