Форум врачей-аспирантов

Здравствуйте, гость ( Вход | Регистрация )

4 страниц V   1 2 3 > »   
Добавить ответ в эту темуОткрыть тему
> Как сравнить результаты регрессии Кокса
propedevt
сообщение 21.02.2012 - 23:27
Сообщение #1





Группа: Пользователи
Сообщений: 27
Регистрация: 5.02.2012
Пользователь №: 23464



Уважаемые форумчане!
Для анализа выживаемости строю модель Кокса, с разными независимыми переменными. И появилась необходимость сравнить модели.
И собственно вопрос у меня к Вам - как сравнить что одна полученная модель лучше чем другая?
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
propedevt
сообщение 23.02.2012 - 08:14
Сообщение #2





Группа: Пользователи
Сообщений: 27
Регистрация: 5.02.2012
Пользователь №: 23464



Никто не знает как сравнить? frown.gif
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
DrgLena
сообщение 23.02.2012 - 11:24
Сообщение #3





Группа: Пользователи
Сообщений: 1325
Регистрация: 27.11.2007
Пользователь №: 4573



У вас была проблема, вы не знали формулу, хотя в документации и к одной и к другой программе, которой вы пользуетесь, она есть, как и в любом учебнике по анализу выживаемости. Вот только вы не поняли, как использовать полученные коэффициенты, у вас, кстати, только один коэффициент значим, в том примере, который вы привели, второй не влияет на функцию выживания. Если я не права, покажите, на том же примере, как пользоваться полученными коэффициентами, используя значения предикторов для реального больного. Если вы в этом разобрались, то ознакомьте форум. Потом можно будет двигаться дальше.
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
p2004r
сообщение 23.02.2012 - 17:50
Сообщение #4





Группа: Пользователи
Сообщений: 1091
Регистрация: 26.08.2010
Пользователь №: 22699



Цитата(propedevt @ 21.02.2012 - 23:27) *
Уважаемые форумчане!
Для анализа выживаемости строю модель Кокса, с разными независимыми переменными. И появилась необходимость сравнить модели.
И собственно вопрос у меня к Вам - как сравнить что одна полученная модель лучше чем другая?



В R вот так:

допустим первая модель

model1<-coxph(Surv(death,status)~strata(cohort)*gapsize)

вторая модель

model2<-coxph(Surv(death,status)~strata(cohort)+gapsize)

сравнение моделей

anova(model1,model2)


Signature
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
p2004r
сообщение 23.02.2012 - 18:01
Сообщение #5





Группа: Пользователи
Сообщений: 1091
Регистрация: 26.08.2010
Пользователь №: 22699



дубль

Сообщение отредактировал p2004r - 23.02.2012 - 18:02


Signature
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
propedevt
сообщение 23.02.2012 - 22:29
Сообщение #6





Группа: Пользователи
Сообщений: 27
Регистрация: 5.02.2012
Пользователь №: 23464



Цитата(DrgLena @ 23.02.2012 - 11:24) *
У вас была проблема, вы не знали формулу, хотя в документации и к одной и к другой программе, которой вы пользуетесь, она есть, как и в любом учебнике по анализу выживаемости. Вот только вы не поняли, как использовать полученные коэффициенты, у вас, кстати, только один коэффициент значим, в том примере, который вы привели, второй не влияет на функцию выживания. Если я не права, покажите, на том же примере, как пользоваться полученными коэффициентами, используя значения предикторов для реального больного. Если вы в этом разобрались, то ознакомьте форум. Потом можно будет двигаться дальше.

Да Вы правы, вопрос задавал. Да и согласен второй предиктор не значим, привел данные просто для примера, чтобы посмотреть как будет выглядеть уравнение. Также Вы правы, я смотрел Халяфяна, смотрел интернет, формула действительно есть, но как применить к уравнению(формуле) коэффициенты, чтобы оно в результате давало прогнозируемые дни жизни так и не понял. До сих пор..
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
propedevt
сообщение 23.02.2012 - 22:33
Сообщение #7





Группа: Пользователи
Сообщений: 27
Регистрация: 5.02.2012
Пользователь №: 23464



Цитата(p2004r @ 23.02.2012 - 17:50) *
В R вот так:

допустим первая модель

model1<-coxph(Surv(death,status)~strata(cohort)*gapsize)

вторая модель

model2<-coxph(Surv(death,status)~strata(cohort)+gapsize)

сравнение моделей

anova(model1,model2)


Уважаемый p2004r! Спасибо за ответ, смотрю инструмент R настолько мощный и универсальный, насколько он непонятна техника его использования новичкам в статистике:) смотрел надстройку для рисования и сравнения ROC-кривых - впечатлило! пробую разобраться в нем, если что задам еще вопросы Вам:)
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
propedevt
сообщение 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
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
p2004r
сообщение 4.03.2012 - 12:54
Сообщение #9





Группа: Пользователи
Сообщений: 1091
Регистрация: 26.08.2010
Пользователь №: 22699



Цитата(propedevt @ 3.03.2012 - 21:40) *
Но теперь у меня возник вопрос как (красиво) репрезентативно это показать это в виде графика или рисунка? (Сами модели Кокса в R попробовал рисовать стандартными средствами)


Приятно слышать что многочисленное международное (и пока увы небольшое русскоговорящее) сообщество пользователей это мощного средства анализа приросло Вами smile.gif


Наверное для демонстрации выбора между моделями можно показать график как меняется AIC от полной модели до оптимальной путем применения к полной модели функции step() которая показывает как меняется AIC при исключении наименее значимых параметров модели. Показывают что выбирают модлель с минимальным для данного датасета AIC.

Итоговые модели можно нарисовать как кривые survivor, density и hazard functions с соответствующими доверительными интервалами.


Signature
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
propedevt
сообщение 4.03.2012 - 14:24
Сообщение #10





Группа: Пользователи
Сообщений: 27
Регистрация: 5.02.2012
Пользователь №: 23464



Посмотрел, нашел что AIC ? информационный критерий Акаике (первый раз слышу про такой insane.gif )
Почитал про 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, тут я не понял frown.gif Сами модели до этого рисовал просто - plot(survfit(model1))
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
p2004r
сообщение 4.03.2012 - 18:31
Сообщение #11





Группа: Пользователи
Сообщений: 1091
Регистрация: 26.08.2010
Пользователь №: 22699



Цитата(propedevt @ 4.03.2012 - 14:24) *
Посмотрел, нашел что AIC ? информационный критерий Акаике (первый раз слышу про такой insane.gif )
Почитал про 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, тут я не понял frown.gif Сами модели до этого рисовал просто - 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


Signature
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
propedevt
сообщение 8.03.2012 - 11:49
Сообщение #12





Группа: Пользователи
Сообщений: 27
Регистрация: 5.02.2012
Пользователь №: 23464



Доброго времени суток, настали выходные и добрался до статистики:) Всех девушек и женщин форума с праздником! smile.gif

Попробовал что Вы рекомендовали Выше, и пришел к такой модели:
Код
> 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) исходя из собственных данных?
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
p2004r
сообщение 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)

ха нашел еще smile.gif

?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


Signature
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
propedevt
сообщение 8.03.2012 - 15:24
Сообщение #14





Группа: Пользователи
Сообщений: 27
Регистрация: 5.02.2012
Пользователь №: 23464



Цитата(p2004r @ 8.03.2012 - 14:44) *
Вот приличный мануал

http://cran.r-project.org/doc/contrib/Fox-...-regression.pdf

по моему Вы пишете параметрическую модель, см. конец второй - начало третьей страницы.

а так hazard (t) = density (t) / survivor (t)


Этот файлик смотрел в самом начале, когда начал кокса строить в R.. сейчас снова посмотрел, и все равно не понял как вычислить h0(t)

Код который сейчас добавили сейчас опробую)
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 
DrgLena
сообщение 8.03.2012 - 18:13
Сообщение #15





Группа: Пользователи
Сообщений: 1325
Регистрация: 27.11.2007
Пользователь №: 4573



...

Сообщение отредактировал DrgLena - 7.11.2022 - 01:14
Вернуться в начало страницы
 
+Ответить с цитированием данного сообщения
 

4 страниц V   1 2 3 > » 
Добавить ответ в эту темуОткрыть тему