Как создать бинарный предитор из многомерной модели glmnet (coxnet)?

Question

Давайте использовать следующий пример:

генерировать данные о выживании (1000 образцов с 30 переменными)

library(glmnet)
library(survival)
set.seed(10101)
N=1000;p=30
nzc=p/3
x=matrix(rnorm(N*p),N,p)
beta=rnorm(nzc)
fx=x[,seq(nzc)]%*%beta/3
hx=exp(fx)
ty=rexp(N,hx)
tcens=rbinom(n=N,prob=.3,size=1) 
y=cbind(time=ty,status=1-tcens)  

используйте glmnet для определения переменных, связанных с выживанием

fit=glmnet(x,y,family="cox")
cvfit <- cv.glmnet(x, y, family="cox")
plot(cvfit)
coefficients <- coef(fit, s = cvfit$lambda.min)
active_coefficients <- coefficients[,1] != 0

подмножество матрицы и сохранить только те параметры (n=17), которые были определены как соответствующие glmnet

x_selected <- x[,active_coefficients]

генерировать модель Кокса с соответствующими параметрами (n=17)

summary(coxph(Surv(y[,1],y[,2])~x_selected))

Вопрос, который сейчас возникает у меня, заключается в том, можно ли и как мне собрать информацию из параметров n = 17, чтобы получить одну (идеально двоичную) переменную предиктора для создания графика Каплана-Мейера, который иллюстрирует прогностические характеристики этого 17-параметрического подпись. Я мог бы использовать PCA и преобразовать в двоичную форму основной компонент (и затем использовать это для графика Каплана-Мейера), но я уверен, что должен быть более элегантный способ, поскольку в основном недавно был выполнен идентичный анализ, который я хотел бы выполнить другими (см. http://ascopubs.org/doi/pdf/10.1200/JCO.2012.45.5626 & http://ascopubs.org/doi/suppl/10.1200/jco.2012.45.5626/suppl_file/DS2_JCO.2012.45.5626.pdf -> авторы использовали glmnet и определили 20 генов, которые имеют отношение к прогнозированию выживания (пока мой код идентичен), но затем они также показывают графики Каплана-Мейера, где они объединяют эту "20 генную сигнатуру" в одну переменная с 3 уровнями ["низкий","средний","высокий"] - посмотрите на рисунок 1 C & D. Я не уверен, как я могу воспроизвести это на моем примере. Есть идеи?

Спасибо!

Answer 1

Уже найдено решение - продолжить анализ следующим образом:

cox_model <- coxph(Surv(y)~x_selected)

#generate a linear predictor from my cox_model
linear_predictor <- predict(cox_model, type="lp")

#check the linear predictor
coxph(Surv(y) ~ linear_predictor)

#stone-beran estimate of survival curve
df <- cbind.data.frame(y,linear_predictor)
s <- prodlim(Surv(time,status) ~ linear_predictor, data=df)

#plot survival curve
xl <- c(0,60)
plot(s, xlab="Time (months)", ylab="Survival rate",
     col=c("green","blue","red"), automar=TRUE, axes=FALSE, atrisk=FALSE,
     confint=FALSE, legend=TRUE,
     legend.title="Coxnet signature", legend.legend=c("low levels", "medium
levels","high levels"), legend.x="bottomright", legend.cex=0.8, xlim=xl)
axis(side=1, at=seq(0,240,12))
axis(side=2, at=seq(0,1,.2))