Alguém pode me explicar meu modelo de Cox em inglês simples?
Eu ajustei o seguinte modelo de regressão de Cox em todos os meus dados usando a cph
função Meus dados são salvos em um objeto chamado Data
. As variáveis w
, x
e y
são contínuos; z
é um fator de dois níveis. O tempo é medido em meses. Alguns dos meus pacientes estão faltando dados para a variável z
( NB : observei devidamente a sugestão do Dr. Harrell, abaixo, de atribuir esses valores a fim de evitar distorções no meu modelo, e o farei no futuro).
> fit <- cph(formula = Surv(time, event) ~ w + x + y + z, data = Data, x = T, y = T, surv = T, time.inc = 12)
Cox Proportional Hazards Model
Frequencies of Missing Values Due to Each Variable
Surv(time, event) w x y z
0 0 0 0 14
Model Tests Discrimination
Indexes
Obs 152 LR chi2 8.33 R2 0.054
Events 64 d.f. 4 g 0.437
Center 0.7261 Pr(> chi2) 0.0803 gr 1.548
Score chi2 8.07
Pr(> chi2) 0.0891
Coef S.E. Wald Z Pr(>|Z|)
w -0.0133 0.0503 -0.26 0.7914
x -0.0388 0.0351 -1.11 0.2679
y -0.0363 0.0491 -0.74 0.4600
z=1 0.3208 0.2540 1.26 0.2067
Também tentei testar a suposição de riscos proporcionais usando o cox.zph
comando abaixo, mas não sei como interpretar seus resultados. A colocação plot()
do comando emite uma mensagem de erro.
cox.zph(fit, transform="km", global=TRUE)
rho chisq p
w -0.1125 1.312 0.2520
x 0.0402 0.179 0.6725
y 0.2349 4.527 0.0334
z=1 0.0906 0.512 0.4742
GLOBAL NA 5.558 0.2347
Primeiro Problema
- Alguém pode me explicar os resultados da saída acima em inglês simples? Tenho formação médica e nenhum treinamento formal em estatística.
Segundo Problema
Conforme sugerido pelo Dr. Harrell, eu gostaria de validar internamente meu modelo executando 100 iterações de validação cruzada 10 vezes usando o
rms
pacote (pelo que entendi, isso implicaria a construção de100 * 10 = 1000
modelos diferentes e, em seguida, solicitando que eles previssem os tempos de sobrevivência de pacientes que nunca haviam visto).Eu tentei usar a
validate
função, como mostrado.> v1 <- validate(fit, method="crossvalidation", B = 10, dxy=T) > v1 index.orig training test optimism index.corrected n Dxy -0.2542 -0.2578 -0.1356 -0.1223 -0.1320 10 R2 0.0543 0.0565 0.1372 -0.0806 0.1350 10 Slope 1.0000 1.0000 0.9107 0.0893 0.9107 10 D 0.0122 0.0128 0.0404 -0.0276 0.0397 10 U -0.0033 -0.0038 0.0873 -0.0911 0.0878 10 Q 0.0155 0.0166 -0.0470 0.0636 -0.0481 10 g 0.4369 0.4424 0.6754 -0.2331 0.6700 10
Como você realiza a reamostragem 100x? Eu acho que meu código acima só executa a validação cruzada uma vez.
Eu queria saber o quão bom meu modelo era na previsão. Eu tentei o seguinte:
> c_index <- abs(v1[1,5])/2 + 0.5 > c_index [1] 0.565984
Isso significa que meu modelo é apenas um pouco melhor do que lançar uma moeda?
Terceiro problema
O Dr. Harrell ressalta que eu assumi linearidade para os efeitos covariáveis e que o número de eventos em minha amostra é apenas grande o suficiente para caber em um modelo confiável se todos os efeitos covariáveis forem lineares.
- Isso significa que devo incluir algum tipo de termo de interação no meu modelo? Se sim, algum conselho sobre o que colocar?
fonte
cph
saída acima em inglês simples ou me indicasse uma referência que o faria. Dr. Harrell, muito obrigado por sua ajuda até agora!Respostas:
Para começar, considere algumas coisas. Primeiro, você está excluindo muitas observações com dados ausentes e isso causará um viés. Considere imputação múltipla. Segundo, existe um método de plotagem para2 × 2
cox.zph
que é útil na avaliação de riscos proporcionais. Terceiro, você assumiu linearidade para os efeitos covariáveis. Quarto, o número de eventos em sua amostra de treinamento é apenas grande o suficiente para caber em um modelo confiável, se todos os efeitos covariáveis forem lineares (o que é raro). E sua amostra de teste precisaria ter talvez 400 eventos antes de produzir uma avaliação confiável da precisão da previsão. Não está claro se você tinha dados suficientes para dividi-los em duas partes. A validação de reamostragem (100 repetições de 10 vezes a validação cruzada ou o bootstrap) é uma solução melhor. Ambos os seus validação externa original (funçõesrcorr.cens
eval.surv
) e remontando validação interna (funçõesvalidate
,calibrate
) são implementados no Rrms
pacote. Estudos de caso para orms
fonte
plot(cox.zph(fit[[1]], transform="km", global=TRUE))
, no entanto, isso rendeuError in plot.cox.zph(cox.zph(fit[[1]], transform = "km", global = TRUE)) : Spline fit is singular, try a smaller degrees of freedom
. Estou chamando essa função incorretamente?plot(cox.zph(...), df=2)
. Veja os estudos de caso nas notas do curso pararms
exemplos ou instale o pacote (que também precisa doHmisc
pacote) e digite estes comandos para exibir os arquivos de ajuda:?cph
?validate.cph
?calibrate.cph
A saída da função R CPH, com base em um exemplo relevante, é explicado neste easy-to-follow papel por J. Fox.
Eu recomendaria fortemente a leitura deste artigo, se você ainda não o fez.
fonte
cph
saída?