#r #predict #survival-analysis #cox-regression #hazard
#r #прогнозировать #анализ выживаемости #cox-регрессия #опасность
Вопрос:
Моя цель — предсказать (предсказать кумулятивную опасность для нового наблюдения из приведенной ниже модели) значение кумулятивной опасности от временной шкалы 0 до времени начала из подобранной модели.
Я установил модель Кокса, используя 2 раза (время начала, которое не равно 0, и время окончания). Таким образом, я могу найти разницу между кумулятивной опасностью в конечное время (т. Е. кумулятивной опасностью от 0 до конечного времени, которую я уже рассчитал, используя ту же подобранную модель) и кумулятивной опасностью во время начала (т. Е. кумулятивной опасностью от 0 до конечного времени, которую я хочу рассчитать здесь), что в конечном итоге даст суммарную опасность между временем начала и окончания каждого наблюдения.
Итак, для получения ожидаемого количества событий я использовал predict(coxph(), newdata, type= "expected") .
Данные, которые я использовал, следующие:
N <- 10^4 # population
H <- within(data.frame(start_time=runif(N, 0, 50), x1=rnorm(N, 2, 1), x2=rnorm(N, -2, 1)), {
lp <- 0.05*x1 0.2*x2
Tm <- qweibull(runif(N,pweibull(start_time,shape = 7.5, scale = 84*exp(-lp/7.5)),1), shape=7.5, scale=84*exp(-lp/7.5))
Cens1 <- 100
event_time <- pmin(Tm,Cens1)
status <- as.numeric(event_time == Tm)})
и код для прогнозирования является:
H$X <- rep(1,nrow(H))
D = coxph(Surv(start_time, event_time, status) ~ X, data = H, x = TRUE )
pred2 <- predict(D, newdata = data.frame(start_time = rep(0,nrow(H)),event_time = H$start_time, status = rep(0,nrow(H)), X = rep(1, nrow(H))), type = "expected")
Но pred2 приводит только к значениям «NA». Может кто-нибудь указать, есть ли какая-либо ошибка в моей идее или в коде
Пожалуйста, дайте мне знать, если потребуются дополнительные разъяснения.
Комментарии:
1. Ваш код не воспроизводим. Пожалуйста, предоставьте
d.2. Извините за это, я исправил код, и я думаю, что теперь он воспроизводим. Можете ли вы попробовать это сейчас. Дайте мне знать, требуется ли еще исправление
3. Смысл модели Кокса заключается в оценке <выделено жирным шрифтом>относительных</bold> опасностей. Это всегда относительно гипотетических случаев со средним значением ковариат. Вы просите его оценить относительные опасности, с которыми не с чем сравнивать риски, поскольку ковариат нет. Если все, что вам нужно, это километровая кривая, то
cphэто неправильная функция.
Ответ №1:
Есть две проблемы. Во-первых, вы сталкиваетесь с проблемой, потому что при указании ~1 , что означает подгонку модели только для перехвата без коэффициентов. итак, все ваши прогнозы дадут вам одно значение?
library(survival)
D <- coxph(Surv(H$start_time, H$event_time, H$status) ~ 1, x = TRUE )
names(D)
[1] "loglik" "linear.predictors" "method"
[4] "residuals" "n" "nevent"
[7] "terms" "assign" "concordance"
[10] "x" "y" "timefix"
[13] "formula" "call"
table(predict(D))
0
10000
Я не думаю, что это имеет большой смысл, и, следовательно, вы сталкиваетесь со всеми ошибками. Итак, вам нужно прогнозировать с независимыми переменными, которые используются, например, в регрессии:
D <- coxph(Surv(start_time,event_time,status) ~ x1 x2, data=H )
pred2 <- predict(D, newdata = data.frame(t_0 = rep(0,nrow(H)),time = H$start_time, event_M = rep(0,nrow(H)), X = rep(1, nrow(H))), type = "expected")
predict(D,newdata=data.frame(x1=runif(10,0,1),x2=runif(10,-1,1)))
1 2 3 4 5 6 7 8
0.3033206 0.4213120 0.3952827 0.3879701 0.4798670 0.2170032 0.3385253 0.4141698
9 10
0.3690579 0.4128084
Когда вы подбираете модель со всеми значениями X = 1, это дает вам все NAs, потому что уже есть перехват, который делает эту переменную избыточной. Вы можете проверить:
H$X = 1
D <- coxph(Surv(start_time, event_time, status) ~ X,data=H)
Call:
coxph(formula = Surv(start_time, event_time, status) ~ X, data = H)
coef exp(coef) se(coef) z p
X NA NA 0 NA NA
Это работает только тогда, когда X является фактической переменной в подогнанных данных, поэтому я использую пример с 2 ковариатами:
H$X = runif(nrow(H))
D <- coxph(Surv(start_time, event_time, status) ~ X x1,data=H)
И вы можете предсказать, например, установив X равным 1 и изменив x1:
predict(D,newdata=data.frame(X=1,x1=c(0.1,0.2,0.3)))
1 2 3
-0.1132548 -0.1084592 -0.1036637
или X в точке 2:
predict(D,newdata=data.frame(X=2,x1=c(0.1,0.2,0.3)))
1 2 3
-0.1579480 -0.1531524 -0.1483568
Комментарии:
1. Как и в вопросе, я хочу оценить базовую кумулятивную опасность без влияния ковариат, и когда я добавил еще одну переменную
H$X = rep(1, nrow(H))и установил модель с использованиемD = coxph(Surv(start_time, event_time, status) ~ X, data = H, x = TRUE ), код выполняется без каких-либо ошибок, но в результате я получаю значения «NA». Я не уверен, верна ли моя концепция. PS: я только что изменил код, и я думаю, что вы не столкнетесь с промежуточной ошибкой (я также объяснил, что я хочу в нем более подробно в вопросе)2. В данных, которые вы вводите в модель, X — это все одно. Чем это отличается от подгонки модели
coxph(Surv(start_time, event_time, status) ~ 1? Теперь вы сталкиваетесь с ошибкой, потому что у вас тоже есть перехват, что делает X избыточным3. Можете ли вы на самом деле подогнать модель под разумные данные? Я обновил свой ответ. Похоже, что это больше не вопрос программирования. Возможно, вы можете сначала проверить это перед публикацией? socialsciences.mcmaster.ca/jfox/Books/Companion/appendices /…
Ответ №2:
Я сам нашел ответ, это просто быстрый трюк, который, я не уверен, будет работать всегда. Если я использую следующую строку перед predict() функцией:
D$coefficients["X"] <- 0
Но я получаю правильные значения, которые проверяются с помощью nelsonaalen() функции, которая не принимает время запуска (или две переменные одновременно)
Дайте мне знать, есть ли какой-либо другой правильный способ решить эту проблему.