Как сравнить результаты регрессии Кокса |
Здравствуйте, гость ( Вход | Регистрация )
Как сравнить результаты регрессии Кокса |
21.02.2012 - 23:27
Сообщение
#1
|
|
Группа: Пользователи Сообщений: 27 Регистрация: 5.02.2012 Пользователь №: 23464 |
Уважаемые форумчане!
Для анализа выживаемости строю модель Кокса, с разными независимыми переменными. И появилась необходимость сравнить модели. И собственно вопрос у меня к Вам - как сравнить что одна полученная модель лучше чем другая? |
|
23.02.2012 - 08:14
Сообщение
#2
|
|
Группа: Пользователи Сообщений: 27 Регистрация: 5.02.2012 Пользователь №: 23464 |
Никто не знает как сравнить?
|
|
23.02.2012 - 11:24
Сообщение
#3
|
|
Группа: Пользователи Сообщений: 1325 Регистрация: 27.11.2007 Пользователь №: 4573 |
У вас была проблема, вы не знали формулу, хотя в документации и к одной и к другой программе, которой вы пользуетесь, она есть, как и в любом учебнике по анализу выживаемости. Вот только вы не поняли, как использовать полученные коэффициенты, у вас, кстати, только один коэффициент значим, в том примере, который вы привели, второй не влияет на функцию выживания. Если я не права, покажите, на том же примере, как пользоваться полученными коэффициентами, используя значения предикторов для реального больного. Если вы в этом разобрались, то ознакомьте форум. Потом можно будет двигаться дальше.
|
|
23.02.2012 - 17:50
Сообщение
#4
|
|
Группа: Пользователи Сообщений: 1091 Регистрация: 26.08.2010 Пользователь №: 22699 |
Уважаемые форумчане! Для анализа выживаемости строю модель Кокса, с разными независимыми переменными. И появилась необходимость сравнить модели. И собственно вопрос у меня к Вам - как сравнить что одна полученная модель лучше чем другая? В R вот так: допустим первая модель model1<-coxph(Surv(death,status)~strata(cohort)*gapsize) вторая модель model2<-coxph(Surv(death,status)~strata(cohort)+gapsize) сравнение моделей anova(model1,model2) |
|
23.02.2012 - 18:01
Сообщение
#5
|
|
Группа: Пользователи Сообщений: 1091 Регистрация: 26.08.2010 Пользователь №: 22699 |
дубль
Сообщение отредактировал p2004r - 23.02.2012 - 18:02 |
|
23.02.2012 - 22:29
Сообщение
#6
|
|
Группа: Пользователи Сообщений: 27 Регистрация: 5.02.2012 Пользователь №: 23464 |
У вас была проблема, вы не знали формулу, хотя в документации и к одной и к другой программе, которой вы пользуетесь, она есть, как и в любом учебнике по анализу выживаемости. Вот только вы не поняли, как использовать полученные коэффициенты, у вас, кстати, только один коэффициент значим, в том примере, который вы привели, второй не влияет на функцию выживания. Если я не права, покажите, на том же примере, как пользоваться полученными коэффициентами, используя значения предикторов для реального больного. Если вы в этом разобрались, то ознакомьте форум. Потом можно будет двигаться дальше. Да Вы правы, вопрос задавал. Да и согласен второй предиктор не значим, привел данные просто для примера, чтобы посмотреть как будет выглядеть уравнение. Также Вы правы, я смотрел Халяфяна, смотрел интернет, формула действительно есть, но как применить к уравнению(формуле) коэффициенты, чтобы оно в результате давало прогнозируемые дни жизни так и не понял. До сих пор.. |
|
23.02.2012 - 22:33
Сообщение
#7
|
|
Группа: Пользователи Сообщений: 27 Регистрация: 5.02.2012 Пользователь №: 23464 |
В R вот так: допустим первая модель model1<-coxph(Surv(death,status)~strata(cohort)*gapsize) вторая модель model2<-coxph(Surv(death,status)~strata(cohort)+gapsize) сравнение моделей anova(model1,model2) Уважаемый p2004r! Спасибо за ответ, смотрю инструмент R настолько мощный и универсальный, насколько он непонятна техника его использования новичкам в статистике:) смотрел надстройку для рисования и сравнения ROC-кривых - впечатлило! пробую разобраться в нем, если что задам еще вопросы Вам:) |
|
3.03.2012 - 21:40
Сообщение
#8
|
|
Группа: Пользователи Сообщений: 27 Регистрация: 5.02.2012 Пользователь №: 23464 |
Уважаемый p2004r! Последовал Вашему совету и начал изучение R. Хочу сказать Вам еще раз спасибо, что направили меня в это русло:) Информации в интернете много, особенно англоязычной. Почитал статьи, некоторые pdf книги, и целом разобрался. Язык R, мне понравился, вспомнилось изучение html и php по-молодости:) Сделал две модели и сравнил их дисперсионным анализом, получив такой результат: Код > anova(model1,model2) Analysis of Deviance Table Cox model: response is Surv(day, event) Model 1: ~ proBNP + efmss Model 2: ~ proBNP + efmss + gal + cys loglik Chisq Df P(>|Chi|) 1 -483.71 2 -428.51 110.41 2 < 2.2e-16 *** --- Signif. codes: 0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1 И так понимаю что различие между моделями достоверно. Но теперь у меня возник вопрос как (красиво) репрезентативно это показать это в виде графика или рисунка? (Сами модели Кокса в R попробовал рисовать стандартными средствами) Сообщение отредактировал propedevt - 3.03.2012 - 21:41 |
|
4.03.2012 - 12:54
Сообщение
#9
|
|
Группа: Пользователи Сообщений: 1091 Регистрация: 26.08.2010 Пользователь №: 22699 |
Но теперь у меня возник вопрос как (красиво) репрезентативно это показать это в виде графика или рисунка? (Сами модели Кокса в R попробовал рисовать стандартными средствами) Приятно слышать что многочисленное международное (и пока увы небольшое русскоговорящее) сообщество пользователей это мощного средства анализа приросло Вами Наверное для демонстрации выбора между моделями можно показать график как меняется AIC от полной модели до оптимальной путем применения к полной модели функции step() которая показывает как меняется AIC при исключении наименее значимых параметров модели. Показывают что выбирают модлель с минимальным для данного датасета AIC. Итоговые модели можно нарисовать как кривые survivor, density и hazard functions с соответствующими доверительными интервалами. |
|
4.03.2012 - 14:24
Сообщение
#10
|
|
Группа: Пользователи Сообщений: 27 Регистрация: 5.02.2012 Пользователь №: 23464 |
Посмотрел, нашел что AIC ? информационный критерий Акаике (первый раз слышу про такой )
Почитал про step ( http://stat.ethz.ch/R-manual/R-devel/libra.../html/step.html ) , и не совсем понял его рекомендуют применять для линейных моделей? насколько оправданно и возможно его применение для регресии кокса? Дальше попробовал его применить - мне нужно сравнить модель1, включающую proBNP + efmss, и модель2, включающую proBNP + efmss + gal + cys. Получилось это в два действия, из 2 модели программа выбирает 3 предиктора, до двух она не доходит. Поэтому вторым действием вручную сделал модель proBNP + efmss: Код > step(coxph(Surv(day,event)~ proBNP + efmss + gal + cys, data), direction = "backward") Start: AIC=865.02 Surv(day, event) ~ proBNP + efmss + gal + cys Df AIC - proBNP 1 864.48 <none> 865.02 - cys 1 865.74 - efmss 1 869.86 - gal 1 932.18 Step: AIC=864.48 Surv(day, event) ~ efmss + gal + cys Df AIC <none> 864.48 - cys 1 867.13 - efmss 1 874.56 - gal 1 930.86 Call: coxph(formula = Surv(day, event) ~ efmss + gal + cys, data = data) coef exp(coef) se(coef) z p efmss -3.75697 0.0234 1.13e+00 -3.32 0.00089 gal 0.09448 1.0991 1.05e-02 8.99 0.00000 cys 0.00019 1.0002 8.61e-05 2.20 0.02800 Likelihood ratio test=209 on 3 df, p=0 n= 197, number of events= 108 > step(coxph(Surv(day,event)~ proBNP + efmss, data), direction = "backward") Start: AIC=971.42 Surv(day, event) ~ proBNP + efmss Df AIC <none> 971.42 - proBNP 1 982.64 - efmss 1 999.83 Call: coxph(formula = Surv(day, event) ~ proBNP + efmss, data = data) coef exp(coef) se(coef) z p proBNP 0.00031 1.00031 8.15e-05 3.81 1.4e-04 efmss -6.32391 0.00179 1.23e+00 -5.13 2.9e-07 Likelihood ratio test=100 on 2 df, p=0 n= 197, number of events= 108 Посоветуйте пожалуйста как нарисовать AIC, тут я не понял Сами модели до этого рисовал просто - plot(survfit(model1)) |
|
4.03.2012 - 18:31
Сообщение
#11
|
|
Группа: Пользователи Сообщений: 1091 Регистрация: 26.08.2010 Пользователь №: 22699 |
Посмотрел, нашел что AIC ? информационный критерий Акаике (первый раз слышу про такой ) Почитал про step ( http://stat.ethz.ch/R-manual/R-devel/libra.../html/step.html ) , и не совсем понял его рекомендуют применять для линейных моделей? насколько оправданно и возможно его применение для регресии кокса? Дальше попробовал его применить - мне нужно сравнить модель1, включающую proBNP + efmss, и модель2, включающую proBNP + efmss + gal + cys. Получилось это в два действия, из 2 модели программа выбирает 3 предиктора, до двух она не доходит. Поэтому вторым действием вручную сделал модель proBNP + efmss: Код > step(coxph(Surv(day,event)~ proBNP + efmss + gal + cys, data), direction = "backward") Start: AIC=865.02 Surv(day, event) ~ proBNP + efmss + gal + cys Df AIC - proBNP 1 864.48 <none> 865.02 - cys 1 865.74 - efmss 1 869.86 - gal 1 932.18 Step: AIC=864.48 Surv(day, event) ~ efmss + gal + cys Df AIC <none> 864.48 - cys 1 867.13 - efmss 1 874.56 - gal 1 930.86 Call: coxph(formula = Surv(day, event) ~ efmss + gal + cys, data = data) coef exp(coef) se(coef) z p efmss -3.75697 0.0234 1.13e+00 -3.32 0.00089 gal 0.09448 1.0991 1.05e-02 8.99 0.00000 cys 0.00019 1.0002 8.61e-05 2.20 0.02800 Likelihood ratio test=209 on 3 df, p=0 n= 197, number of events= 108 > step(coxph(Surv(day,event)~ proBNP + efmss, data), direction = "backward") Start: AIC=971.42 Surv(day, event) ~ proBNP + efmss Df AIC <none> 971.42 - proBNP 1 982.64 - efmss 1 999.83 Call: coxph(formula = Surv(day, event) ~ proBNP + efmss, data = data) coef exp(coef) se(coef) z p proBNP 0.00031 1.00031 8.15e-05 3.81 1.4e-04 efmss -6.32391 0.00179 1.23e+00 -5.13 2.9e-07 Likelihood ratio test=100 on 2 df, p=0 n= 197, number of events= 108 Посоветуйте пожалуйста как нарисовать AIC, тут я не понял Сами модели до этого рисовал просто - plot(survfit(model1)) AIC нужен чтобы выбрать правильно модель которая сочетает оптимально и сложность своей структуры и точность подгонки к данным. Такая модель имеет наибольшую прогностическую ценность для новых данных. step() можно применять, модель чем бы она не была описана сводится к своему предсказанию для экспериментальной точки. По остаткам и проводится сравнение моделей. Вот степ пишет как уменьшается AIC упрощаемой модели, достигает минимума и потом растет. Модель с минимальным AIC должна обладать наилучшими способностями прогноза. Да, для трех параметров пожалуй достаточно перечислить в порядке убывания. Обычно параметров очень много и получается последовательность которую лучше увидеть графически. Правильно все рисовали, только вот такое замечание есть ?survival::survfit.coxph Users are strongly advised to use the newdata argument. Note that this data set needs to contain values for the main effects but not for any interaction terms. Еще можно протестировать выполнены ли предположения о пропорциональности cox.zph Вот еще графики диагностические для таких моделей http://stats.stackexchange.com/questions/2...xph-and-rms-cph Вот для визуализации самой модели plot.coxph {flexCrossHaz} Plot the estimated time varying log hazard ratio along with pointwise confidence intervals |
|
8.03.2012 - 11:49
Сообщение
#12
|
|
Группа: Пользователи Сообщений: 27 Регистрация: 5.02.2012 Пользователь №: 23464 |
Доброго времени суток, настали выходные и добрался до статистики:) Всех девушек и женщин форума с праздником!
Попробовал что Вы рекомендовали Выше, и пришел к такой модели: Код > model5<-coxph(Surv(day,event)~ efmss + gal + cys, data) > model5 Call: coxph(formula = Surv(day, event) ~ efmss + gal + cys, data = data) coef exp(coef) se(coef) z p efmss -3.75697 0.0234 1.13e+00 -3.32 0.00089 gal 0.09448 1.0991 1.05e-02 8.99 0.00000 cys 0.00019 1.0002 8.61e-05 2.20 0.02800 Likelihood ratio test=209 on 3 df, p=0 n= 197, number of events= 108 И теперь настала пора вернуться к самому первому вопросу:) как описать модель в виде функции. Нашел пока следующую информацию: http://www.medcalc.org/manual/cox_proportional_hazards.php Где указана понятная формула: pi прогностический индекс - но я не пойму где брать H0(t) в эту формулу? Если посмотреть существующие модели: 1)нева75: http://www.almazovcentre.ru/node/357 то там Н0(t) это базовый риск в момент времени t, и берут его с графика 2)сиэтлская модель: http://www.medmir.com/content/view/1727/61/ (англ вариант) то там Н0(t) = годы умноженные на коэффициент, равный 0,0405. Собственно вопрос пока таков: как мне получить базовую кумулятивную выживаемость Н0(t) исходя из собственных данных? |
|
8.03.2012 - 14:44
Сообщение
#13
|
|
Группа: Пользователи Сообщений: 1091 Регистрация: 26.08.2010 Пользователь №: 22699 |
Вот приличный мануал
http://cran.r-project.org/doc/contrib/Fox-...-regression.pdf по моему Вы пишете параметрическую модель, см. конец второй - начало третьей страницы. а так hazard (t) = density (t) / survivor (t) ха нашел еще ?basehaz This function exists primarily because users will look for the phrase 'baseline hazard' (often SAS converts looking for familiar keywords.) еще есть вот это Код > Similarly, when I do plot(zph), B(t) is fairly non-constant. > > This isn't inherently a problem for me. I don't need a hard single number > > to characterize the shape of the excess risk. However, I'd like to be > > able to say > > something qualitative about the shape of the excess risk for the predictor. > > E.g., is it linear, monotonically increasing, monotonially decreasing, etc. > > Is it safe to use the coxph diagnostic plot for this purpose? Basically - yes you can. There are a few caveats: 1. As a computational shortcut cox.zph assumes that var(X) is approximately constant over time, where X is the matrix of covariates. (Improving this has been on my to do list for some time). I have found this to be almost always true, but if you have a data set where e.g. everyone in treatment 1 is crossed over at 6 months, then you can get odd results for that covariate. I've run across 2-3 such data sets in 10+ years. 2. The spline curve on the plot is "for the eye". You can certainly use other smoothings, fit a line, etc. Often you can find a simpler fit. zpfit <- cox.zph(mycoxfit, transform='identity') plot(zpfit$x, zpfit$y[,1], xlab='Time') #look at variable 1 lines(lowess(zpfit$x, zpfit$y[,1]), col=2) abline( lm(zpfit$y[,1] ~zpfit$x), col=3) plot(zpfit$x, zpfit$y[,1], log='x') #same as transform=log etc. Sometimes the regression spline fit, the default for cox.zph, puts an extra "hook" on the end of the curve, somewhat like polynomials will. Terry T. Код -- begin included message ---
I am hoping for some advice regarding obtaining the values for the hazard function in a cox regression that I have undertaken. I have a model in the following form, analysed with the package survival (v. 2.34-1) and a log-log plot obtained using Design (v. 2.1-2). For two variables, the lines in the survival curves crossed. The statistician I been obtaining advice from (who does not use R) asked me to obtain the hazard function values. I am wanting to confirm whether basehaz is the correct command to obtain such values, to better understand what is occurring in the log plots. --- end inclusion ---- The best tool for understanding what is happening is cox.zph. cox.V <- coxph(Surv(intDaysUntilFVPO, Event_v) ~ intAgeAtMHCIndex + PRE + group + MHC + strGender, data = recidivismv) zpfit <- cox.zph(cox.V, transform='identity') plot( zpfitp[1]) #plot for the first variable plot( zpfitp[2]) #plot for the second You will get a plot of beta(t) versus time, a horizontal line corresponds to a constant hazard ratio. The Cox model coefficient is an "average" of the curve, e.g., the best horizontal line fit to it. Other than this, you can get the predicted cumulative hazard for any particular choice of covariates, the derivative of this = hazard and requires some form of smoothing. Terry Therneau Сообщение отредактировал p2004r - 8.03.2012 - 15:16 |
|
8.03.2012 - 15:24
Сообщение
#14
|
|
Группа: Пользователи Сообщений: 27 Регистрация: 5.02.2012 Пользователь №: 23464 |
Вот приличный мануал http://cran.r-project.org/doc/contrib/Fox-...-regression.pdf по моему Вы пишете параметрическую модель, см. конец второй - начало третьей страницы. а так hazard (t) = density (t) / survivor (t) Этот файлик смотрел в самом начале, когда начал кокса строить в R.. сейчас снова посмотрел, и все равно не понял как вычислить h0(t) Код который сейчас добавили сейчас опробую) |
|
8.03.2012 - 18:13
Сообщение
#15
|
|
Группа: Пользователи Сообщений: 1325 Регистрация: 27.11.2007 Пользователь №: 4573 |
...
Сообщение отредактировал DrgLena - 7.11.2022 - 01:14 |
|