Como fazer análise ROC em R com um modelo de Cox

10

Eu criei alguns modelos de regressão de Cox e gostaria de ver o desempenho desses modelos e pensei que talvez uma curva ROC ou uma estatística c possa ser útil semelhante ao uso deste artigo:

JN Armitage e JH van der Meulen, “Identificando comorbidades em pacientes cirúrgicos usando dados administrativos com o Royal College of Surgeons Charlson Score”, British Journal of Surgery, vol. 97, num. 5, ss. 772-781, Maj 2010.

Armitage usou regressão logística, mas me pergunto se é possível usar um modelo do pacote de sobrevivência, o survivalROC dá uma dica disso, mas não consigo descobrir como fazer isso funcionar com uma regressão Cox regular.

Ficaria muito grato se alguém me mostrasse como fazer uma análise ROC neste exemplo:

library(survival)
data(veteran)

attach(veteran)
surv <- Surv(time, status)
fit <- coxph(surv ~ trt + age + prior, data=veteran)
summary(fit)

Se possível, eu apreciaria a saída c-statics bruta e um bom gráfico

Obrigado!

Atualizar

Muito obrigado pelas respostas. @ Dwin: Gostaria apenas de ter certeza de que entendi bem antes de selecionar sua resposta.

O cálculo como eu o entendo de acordo com a sugestão de DWin:

library(survival)
library(rms)
data(veteran)

fit.cph <- cph(surv ~ trt + age + prior, data=veteran, x=TRUE, y=TRUE, surv=TRUE)

# Summary fails!?
#summary(fit.cph)

# Get the Dxy
v <- validate(fit.cph, dxy=TRUE, B=100)
# Is this the correct value?
Dxy = v[rownames(v)=="Dxy", colnames(v)=="index.corrected"]

# The c-statistic according to the Dxy=2(c-0.5)
Dxy/2+0.5

Não estou familiarizado com a função validar e a inicialização, mas depois de analisar o prof. A resposta de Frank Harrel aqui no R-help eu achei que provavelmente é o caminho para obter o Dxy. A ajuda para validar estados:

... A correlação de classificação Dxy de Somers deve ser calculada a cada nova amostra (isso leva um pouco mais do que as estatísticas baseadas em probabilidade). Os valores correspondentes à linha Dxy são iguais a 2 * (C - 0,5), onde C é o índice C ou a probabilidade de concordância.

Acho que estou principalmente confuso com as colunas. Achei que o valor corrigido é o que devo usar, mas ainda não entendi a saída validada:

      index.orig training    test optimism index.corrected   n
Dxy      -0.0137  -0.0715 -0.0071  -0.0644          0.0507 100
R2        0.0079   0.0278  0.0037   0.0242         -0.0162 100
Slope     1.0000   1.0000  0.2939   0.7061          0.2939 100
...

Na pergunta R-help, entendi que deveria ter "surv = TRUE" no cph se tiver estratos, mas não tenho certeza sobre qual é o objetivo do parâmetro "u = 60" na função validate. Ficaria muito grato se você pudesse me ajudar a entender isso e verificar se não cometi nenhum erro.

Max Gordon
fonte
2
Eu provavelmente daria uma olhada no pacote rms e seu cph()comando.
chl
2
index.correctedé o que deve ser enfatizado. Essas são estimativas de provável desempenho futuro. u=60não é necessário, validatepois você não possui estratos. Se você tivesse estratos, as curvas de sobrevivência podem se cruzar e é necessário especificar um ponto de tempo específico para obter a área ROC generalizada.
Frank Harrell

Respostas:

2

@chl apontou para uma resposta específica à sua pergunta. A cphfunção do pacote 'rms' produzirá um Somers-D que pode ser transformado trivialmente em um índice-c. Entretanto, Harrell (que introduziu o índice c na prática bioestatística) acha que isso não é sensato como uma estratégia geral para avaliar medidas prognósticas, porque tem baixo poder de discriminação entre alternativas. Em vez de confiar na literatura cirúrgica para sua orientação metodológica, seria mais sensato procurar a sabedoria acumulada no texto de Harrell, "Estratégias de Modelagem de Regressão" ou "Modelos de Previsão Clínica" de Steyerberg.

DWin
fonte
4
Obrigado pela observação. Eu acho que e não são ruins para descrever a discriminação preditiva de um único modelo pré-especificado. Mas, como você disse, eles não têm poder para fazer mais do que isso. CDxyC
precisa
Obrigado por suas respostas, minha situação é que tenho três pontuações diferentes que quero comparar e ver como elas se saem. Não tive tempo de analisar a parte de Somers-D e voltarei assim que tiver tempo (dei uma olhada rápida e não achei nada útil). Também encomendei o livro @FrankHarrell, "Estratégias de modelagem de regressão", ISBN 13: 978-0387952321, e espero que ele me guie nas minhas escolhas.
Max Gordon
2
Como Dxy = 2 * (c-0,5), o cálculo de c dado Dxy deve ser trivial.
DWin 27/10/11
3

Dependendo das suas necessidades, a incorporação de um modelo em um modelo maior e a realização de um teste de taxa de verossimilhança para o valor agregado das variáveis ​​adicionais fornecerão um teste poderoso. Meu livro fala sobre um índice decorrente dessa abordagem (o "índice de adequação").χ2

Frank Harrell
fonte
+1 por me guiar na direção certa. Acabei de terminar a estatística C e a pontuação mais detalhada que eu estou vendo tinha uma estatística C de 0,4365081 enquanto a outra tinha 0,4414625 (acho que eu deveria contar 0,5-Dxy / 2 no meu caso). Demorei bastante para fazer o cálculo na minha amostra de 140.000; Eu tive que diminuir o bootstraps para 10 e não sei ao certo qual é o impacto disso. Estou ansioso para ler seu livro (está no correio) e espero que isso me ajude a entender melhor a metodologia e comparar a estatística C com o índice de adequação.
Max Gordon
Boa. Não é fácil dizer se 0,44 x 0,43 significa muito sem observar as distribuições dos valores previstos.
precisa
Eu entendo que é difícil comentar sobre números assim. Vou tentar olhar para a distribuição. Minha principal interpretação do resultado é que muito pouco é explicado pelo meu modelo e, embora exista uma pequena diferença, provavelmente não é muito importante. Seria interessante o que esperar em um cenário de sobrevivência - atingir um valor de 0,8, como fizeram na análise que referenciei na minha pergunta, parece muito distante ... mas, novamente, minha sobrevivência é a sobrevivência de uma prótese implantada e não sobrevida do paciente. Eles também usaram regressão logística que talvez mude a estimativa.
Max Gordon
A regressão logística não funcionaria se o tempo fosse importante ou o tempo de acompanhamento varia entre os sujeitos. De volta à pergunta original, os riscos previstos terão uma distribuição estreita se uma variação muito pequena for explicada pelo modelo.
precisa
Acabei de receber seu livro ... Eu tive um bloqueio rápido na parte de sobrevivência, mas quando eu experimento seu Estudo de caso no capítulo 20, mas recebo um erro na parte de imputação (w, sz): 'variable sz não tem um atributo names () '. Eu segui chapt. 8: carregou o dataframe com getHdata (próstata) (não conseguiu encontrar o site no livro), fez o w <- transcan (~ sz + sg + ap + sbp + dbp + age + wt + hg + ekg + pf + bm + hx, imputada = T, transformado = T, imcat = "árvore", data = próstata), mas eu não encontrar nada sobre nomear ...
Max Gordon