Моя цель - изучить, как включение нескольких факторов риска улучшает клиническую модель, предсказывающую частоту инсульта в анализе выживаемости.Я хочу использовать NRI для сравнения двух моделей (1 базовая модель против базовой модели + новый фактор риска).Однако я бы предпочел разделить вероятности этих моделей на категориальные переменные (то есть риск 0-3%, риск 4-6%, риск> 7% ...).
Для этой цели я подготовил следующий код с R, но я не уверен, является ли он правильным или нет, и я бы предпочел уточнить его =).
Я буду использовать открытыйДанные о выживаемости «легких», чтобы показать воспроизводимый пример.ВАЖНО: Я не буду рассматривать данные, подвергнутые цензуре, в этом анализе, но их следует анализировать соответствующим образом в дальнейших анализах:
library(survival)
library(pec)##To calculate probabilities of event at a point of time.
library(Hmisc)##To calculate NRI
data(lung)
lung <- na.omit(lung)
lung$status <- lung$status-1##necessary for NRI package
stats <- lung$status
tempo <- lung$time
##Create two models, 1 baseline and 1 baseline+new predictor
model1 <- coxph(Surv(tempo,stats)~age+sex, data=lung, x=T)
model2 <- coxph(Surv(tempo,stats)~age+sex+factor(ph.ecog), data=lung,
x=T)
##Estimate the survival probability at time= 500.
lung$x <- predictSurvProb(model1, newdata=lung, times=c(500))
lung$y <- predictSurvProb(model2, newdata=lung, times=c(500))
##Calculate the probability of event (1-survival) and we categorize it.
lung$x <- cut(1-lung$x, c(0, 0.25, 0.5, 0.75, 1), include.lowest=TRUE,
right=FALSE)
lung$y <- cut(1-lung$y, c(0, 0.25, 0.5, 0.75, 1), include.lowest=TRUE,
right=FALSE)
##Confusion matrix
cases <- subset(lung, stats==1)
controls <- subset(lung, stats==0)
table(lung$x, lung$y)##All patients
table(cases$x, cases$y)##cases
table(controls$x, controls$y)##controls
##Calculate NRI
x <- as.numeric(lung$x)/4##Predictions should be ranged between 0-1
y <- as.numeric(lung$y)/4
improveProb(x, y, stats)
Мои вопросы: 1) Этот режим предназначен для расчета риска события при t = 500правильный?2) Правильно ли рассчитать NRI с помощью этой функции?Я пробовал другие пакеты (например, NRIcens) и получил те же результаты… Заранее спасибо за вашу помощь!