Как сравнить результаты регрессии Кокса Опции

[propedevt]

Уважаемые форумчане!
Для анализа выживаемости строю модель Кокса, с разными независимыми переменными. И появилась необходимость сравнить модели.
И собственно вопрос у меня к Вам - как сравнить что одна полученная модель лучше чем другая?

[propedevt]

Никто не знает как сравнить?

[DrgLena]

У вас была проблема, вы не знали формулу, хотя в документации и к одной и к другой программе, которой вы пользуетесь, она есть, как и в любом учебнике по анализу выживаемости. Вот только вы не поняли, как использовать полученные коэффициенты, у вас, кстати, только один коэффициент значим, в том примере, который вы привели, второй не влияет на функцию выживания. Если я не права, покажите, на том же примере, как пользоваться полученными коэффициентами, используя значения предикторов для реального больного. Если вы в этом разобрались, то ознакомьте форум. Потом можно будет двигаться дальше.

[p2004r]

Уважаемые форумчане!
Для анализа выживаемости строю модель Кокса, с разными независимыми переменными. И появилась необходимость сравнить модели.
И собственно вопрос у меня к Вам - как сравнить что одна полученная модель лучше чем другая?

В R вот так:

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

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

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

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

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

anova(model1,model2)

[p2004r]

дубль

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

[propedevt]

У вас была проблема, вы не знали формулу, хотя в документации и к одной и к другой программе, которой вы пользуетесь, она есть, как и в любом учебнике по анализу выживаемости. Вот только вы не поняли, как использовать полученные коэффициенты, у вас, кстати, только один коэффициент значим, в том примере, который вы привели, второй не влияет на функцию выживания. Если я не права, покажите, на том же примере, как пользоваться полученными коэффициентами, используя значения предикторов для реального больного. Если вы в этом разобрались, то ознакомьте форум. Потом можно будет двигаться дальше.

Да Вы правы, вопрос задавал. Да и согласен второй предиктор не значим, привел данные просто для примера, чтобы посмотреть как будет выглядеть уравнение. Также Вы правы, я смотрел Халяфяна, смотрел интернет, формула действительно есть, но как применить к уравнению(формуле) коэффициенты, чтобы оно в результате давало прогнозируемые дни жизни так и не понял. До сих пор..

[propedevt]

В R вот так:

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

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

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

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

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

anova(model1,model2)

Уважаемый p2004r! Спасибо за ответ, смотрю инструмент R настолько мощный и универсальный, насколько он непонятна техника его использования новичкам в статистике:) смотрел надстройку для рисования и сравнения ROC-кривых - впечатлило! пробую разобраться в нем, если что задам еще вопросы Вам:)

[propedevt]

Уважаемый 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]

Но теперь у меня возник вопрос как (красиво) репрезентативно это показать это в виде графика или рисунка? (Сами модели Кокса в R попробовал рисовать стандартными средствами)

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


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

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

[propedevt]

Посмотрел, нашел что 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))

[p2004r]

Посмотрел, нашел что 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

[propedevt]

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

Попробовал что Вы рекомендовали Выше, и пришел к такой модели:

Код

> 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]

Вот приличный мануал

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

[propedevt]

Вот приличный мануал

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

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

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

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

Код который сейчас добавили сейчас опробую)

[DrgLena]

...

Сообщение отредактировал DrgLena - 7.11.2022 - 01:14

[propedevt]

Просто средние ковариат? тогда вот они:
gal 25,234
cys 3341,777
efmss 0,280

[DrgLena]

...

Сообщение отредактировал DrgLena - 7.11.2022 - 01:15

[propedevt]

Н0(t) это базовый риск в момент времени t, а не функция выживаемости.

Да согласен, совершенно неправильно сказал

Вот по этим средним и полученным коэффициентам и получаете базовый риск, вначале в единицу времени наблюдения

А чем является единица времени наблюдения? У меня максимальный период наблюдения 26 месяцев, а в полученных днях время жизни в днях

Сообщение отредактировал propedevt - 8.03.2012 - 23:52

[DrgLena]

Вы сами выбираете единицу измерения времени, в зависимости от того как вы определяли время наблюдения, если вычитали две даты и преобразовывали разницу в месяцы или сразу вводили месяцы наблюдения за больным, тогда единица будет месяц, в приведенной вами ссылке, по моему, год (умножение на 1 и на 5 при расчете 5 летней выживаемости). Большинство исследователей не используют кокс медели для целей оценки времени жизни. Оцениваются exp(b) для предикторов. Ограничения не столько технические, сколько нравственные.

[propedevt]

Вычитал даты и получал дни жизни (до 805 дней максимум), в месяцы не преобразовывал.

Так, банально подсчитал PI для средних, он равен (-3.75697 х 0,280) + (0.09448 х 25,234) + (0.00019 х 3341,777) = 1,96709435

Как теперь построить график функции риска? (Н0(t) подобный такому:

[propedevt]

Уважаемый p2004r, прошу меня извинить, но не сразу заметил в Вашем сообщении кроме кода, который привели, Вы указали функцию basehaz.

Написал такой простой код:
baseline <- basehaz(model4, centered=FALSE)
plot(baseline$time,baseline$hazard,xlab="Time",ylab="Hazard",col="red",type="l")

и получил в итоге график:


и для средних:
baseline1 <- basehaz(model4, centered=TRUE)
plot(baseline1$time,baseline1$hazard,xlab="Time",ylab="Hazard",col="red",type="l")

Уважаемая DrgLena, скажите пожалуйста это и есть нужный мне график H0(t) ?

Сообщение отредактировал propedevt - 9.03.2012 - 10:45

[DrgLena]

...

Сообщение отредактировал DrgLena - 7.11.2022 - 01:17

[p2004r]

Уважаемый p2004r, прошу меня извинить, но не сразу заметил в Вашем сообщении кроме кода, который привели, Вы указали функцию basehaz.

Написал такой простой код:
baseline <- basehaz(model4, centered=FALSE)
plot(baseline$time,baseline$hazard,xlab="Time",ylab="Hazard",col="red",type="l")

и получил в итоге график:


Уважаемая DrgLena, скажите пожалуйста это и есть нужный мне график H0(t) ?

Код

> basehaz
function (fit, centered = TRUE)
{
    if (!inherits(fit, "coxph"))
        stop("must be a coxph object")
    sfit <- survfit(fit)
    H <- -log(sfit$surv)
    strata <- sfit$strata
    if (!is.null(strata))
        strata <- factor(rep(names(strata), strata), levels = names(strata))
    if (!centered) {
        z0 <- fit$means
        bz0 <- sum(z0 * coef(fit))
        H <- H * exp(-bz0)
    }
    if (is.null(strata))
        return(data.frame(hazard = H, time = sfit$time))
    else return(data.frame(hazard = H, time = sfit$time, strata = strata))
}
<environment: namespace:survival>

В Вашем случае срабатывает

Код

sfit <- survfit(fit)
H <- -log(sfit$surv)
z0 <- fit$means
bz0 <- sum(z0 * coef(fit))
H <- H * exp(-bz0)

Именно hazard у одиночного случая по времени
h(t)=h0(t)*exp(sum(z_случая*coef(fit)))

h0(t) hazard для случая когда z<-0, и его считают как -log(sfit$surv)

Как показывает код функция считает именно то что Вы хотели. Подставляет средние всех ковариант.

Finally, the program lists the baseline cumulative hazard H0(t), with the cumulative hazard and survival at mean of all covariates in the model.

Вообще то для расчетов проще использовать predict(), весь этот закат солнца в ручную признаться утомителен. Если конечно это не самоцель

http://stat.ethz.ch/R-manual/R-devel/libra...dict.coxph.html

The Cox model is a relative risk model; predictions of type "linear predictor", "risk", and "terms" are all relative to the sample from which they came. By default, the reference value for each of these is the mean covariate within strata. The primary underlying reason is statistical: a Cox model only predicts relative risks between pairs of subjects within the same strata, and hence the addition of a constant to any covariate, either overall or only within a particular stratum, has no effect on the fitted results. Using the reference="strata" option causes this to be true for predictions as well.

When the results of predict are used in further calculations it may be desirable to use a fixed reference level. Use of reference="sample" will use the overall means, and agrees with the linear.predictors component of the coxph object (which uses the overall mean for backwards compatability with older code).

Predictions of type "expected" incorportate the baseline hazard and are thus absolute instead of relative; the reference option has no effect on these.

[propedevt]

Да понял, что если риск 0,8 то выживаемость 0,2
График нужен для того чтобы научному руководителю его показать, и потом когда описывать модель его нужно будет наверное показывать:)

Для конкретного больного например подсчитываю:
Pi (-3.75697 х 0,4) + (0.09448 х 8) + (0.00019 х 2400) = 0,290948

Дальше беру риск например на 600 день с графика и подставляю в формулу:


s(t) = exp (-0,25 x 0,290948) = 0,93 (93%)

правильно? или я чего-то не понимаю

[DrgLena]

Да понял, что если риск 0,8 то выживаемость 0,2

Нет, не все так просто, это чтобы график cum функций перевернуть

Сообщение отредактировал DrgLena - 9.03.2012 - 11:57

[DrgLena]

Закат солнца в ручную у меня не получается с такими же значениями.

[p2004r]

Закат солнца в ручную у меня не получается с такими же значениями.

Думаю надо наконец загрузить в тред конкретный датасет, зафитить конкретную модель, и пошагово сравнить все что для неё относительно этого датасета рекомендовано посчитать

[DrgLena]

...

Сообщение отредактировал DrgLena - 7.11.2022 - 01:18

[propedevt]

Что-то Вы меня совсем запутали:( перечитал что мне ответили по десять раз.. не знаю что делать дальше..

Давайте начнем с начала по порядку:
Есть формула
в нее я подставляю pi , вычисленный для конкретного больного по его данным.
Далее я должен найти h0(t) базовый риск во время t. его например беру со своего второго графика. Правильно? или Вы хотите сказать можно вычислять риск с помощью функции predict() ?

[p2004r]

Что-то Вы меня совсем запутали:( перечитал что мне ответили по десять раз.. не знаю что делать дальше..

Давайте начнем с начала по порядку:
Есть формула
в нее я подставляю pi , вычисленный для конкретного больного по его данным.
Далее я должен найти h0(t) базовый риск во время t. его например беру со своего второго графика. Правильно? или Вы хотите сказать можно вычислять риск с помощью функции predict() ?

не не не Девид Блейн (С)

давайте датасет (можно приатачить csv заархививовав его например rar)

и на какой модели Вы остановились

PS погода хорошая --- схожу в город и вернувшись непременно по шагам датасет посчитаем

Сообщение отредактировал p2004r - 9.03.2012 - 14:16

[propedevt]

Хорошо, файл прикрепил. В последней модели использовал предикторы efmss, gal, cys.
Жду Вас с нетерпением:)
только что будете делать опишите так, чтобы я понял и мог сделать тоже самое:)

Сообщение отредактировал propedevt - 9.03.2012 - 14:54

Прикрепленные файлы

 cox_m4.rar ( 1,79 килобайт ) Кол-во скачиваний: 281

 

[p2004r]

Хорошо, файл прикрепил. В последней модели использовал предикторы efmss, gal, cys.
Жду Вас с нетерпением:)
только что будете делать опишите так, чтобы я понял и мог сделать тоже самое:)

Сначала данные и подгонка модели.

Код

> library(survival)
> data<-read.csv2("cox_m4.csv") #читаем данные
> head(data) #смотрим несколько первых строк
  X day event gal proBNP  cys efmss
1 1 805     0   6    900 1800  0.30
2 2 805     0   3    400 1200  0.45
3 3 805     0   7    150 1000  0.69
4 4 805     0  13    200  900  0.57
5 5 760     1  20    900 3300  0.36
6 6 805     0   8    800 2400  0.40
> length(data[,1]) # сколько всего случаев
[1] 197
> model5<-coxph(Surv(day,event)~ efmss + gal + cys, data) # собственно выбранная модель
> summary(model5)
Call:
coxph(formula = Surv(day, event) ~ efmss + gal + cys, data = data)

  n= 197, number of events= 108

            coef  exp(coef)   se(coef)      z Pr(>|z|)    
efmss -3.757e+00  2.335e-02  1.130e+00 -3.324 0.000887 ***
gal    9.448e-02  1.099e+00  1.051e-02  8.987  < 2e-16 ***
cys    1.896e-04  1.000e+00  8.606e-05  2.203 0.027583 *  
---
Signif. codes:  0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1

      exp(coef) exp(-coef) lower .95 upper .95
efmss   0.02335    42.8186  0.002549     0.214
gal     1.09908     0.9098  1.076669     1.122
cys     1.00019     0.9998  1.000021     1.000

Rsquare= 0.654   (max possible= 0.996 )
Likelihood ratio test= 208.9  on 3 df,   p=0
Wald test            = 158.6  on 3 df,   p=0
Score (logrank) test = 215.5  on 3 df,   p=0

Прогноз который считает модель

Predicted values are available based on either the linear predictor

the risk for an individual relative to the average subject within the data set

the expected number of events for an individual over the time interval that they were
observed to be at risk, which is a component of the martingale residual,
or for individual components of the linear predictor

Код

> predict(model5, type="lp", se.fit=FALSE)
  [1] -2.18632414 -3.14706993 -3.70875790 -2.71002012 -0.80464171 -2.25929859
  [7] -3.08799757 -2.40811880 -2.41023773  0.13723479 -3.16344410 -2.52261679
[13] -1.74607290 -0.58176455 -2.37054908 -2.37372747 -2.45183038 -2.37266801
[19] -1.02844012 -1.37028349 -1.93966416 -2.17903068 -2.18557176 -3.29553699
[25] -2.05531376 -2.82644514 -1.94015553 -2.73267437 -0.80306783 -2.44966538
[31] -2.60104982 -1.03368376 -2.46356747 -2.10313884 -3.03069098 -2.61667161
[37] -2.59883875 -2.68933898 -2.82390395 -2.22311845 -1.54034533  0.25133355
[43] -3.59715426 -2.24247912  0.32875319 -1.31673882 -1.27992148 -1.27989844
[49] -2.54087192 -1.95878685 -0.94026614 -2.40983850 -1.84059606 -1.72958355
[55] -2.07138085 -2.50542113  1.93030325 -0.82466261 -0.76634264 -0.11003936
[61] -0.99754969  0.42022074 -0.33413724 -0.55814298 -0.69445069  0.69496973
[67] -1.21972054  1.99028925  0.64639148 -2.36953568  0.25133355  0.92986880
[73]  0.36225391 -1.39590120 -0.88269855  0.38309631  0.58988335  0.04010133
[79] -2.92235788 -1.71504426 -2.03454047 -1.71436098  0.28859619 -1.42956318
[85]  0.27215291 -1.88430764 -0.90201314 -0.63867190 -0.48286531  0.98774348
[91]  0.45748338  0.72193016 -0.88378105 -1.79018371 -2.18345283  1.76678250
[97] -2.03463262  0.55405636 -0.10584757 -0.63057998  0.40198864 -0.65904596
[103]  1.59332732  0.19256827 -0.86439734 -0.22425331 -0.52521034 -0.40772586
[109]  0.45776743  0.13723479 -1.27947618 -1.24049383 -0.19868167 -0.86263157
[115] -0.16445160  0.78097948  0.91215873 -1.86497001  0.66570607 -0.48220507
[121]  0.61098676  0.79848223  0.15555902 -1.36372016  0.06169611 -2.20566177
[127] -1.14420490  0.55335005  5.46721790  0.04452348 -0.33559594  1.19831547
[133] -0.59959741  2.74521532  4.76859263  1.93140879  3.25585374  1.72492883
[139]  0.65241815  1.35099735  2.23344833  1.16027740  2.63031810  1.91035145
[145]  1.92964301  0.70362896  3.32882819  2.00690139  2.23307213  1.12124899
[151]  4.65282785  3.76430413  1.51548542  4.82310462  3.44325707  3.25580767
[157]  4.36803005  3.02857747  2.00831402  3.01034537  2.91506983  1.32697656
[163]  2.68616600  2.97098683  0.89077128  2.36754496  2.68859201  1.44105228
[169]  0.58767227  0.51403758  2.49848623  1.49805178  2.00586497  0.58988335
[175]  2.15793268  1.70568335  2.10405789  0.77956686  1.59186863  2.08345347
[181]  0.40013072  2.28925015  0.30645209  0.87018255  2.13963148  3.32878212
[187]  2.70536541  1.10096707  1.55305515  3.55777810  0.04421639  2.38582312
[193]  3.33024081  2.21547724  1.51548542  0.23020710  1.50539054

> predict(model5, type="risk")
  [1]   0.11232890   0.04297787   0.02450795   0.06653547   0.44724814
  [6]   0.10442370   0.04559316   0.08998441   0.08979395   1.14709744
[11]   0.04227987   0.08024934   0.17445771   0.55891127   0.09342941
[16]   0.09313293   0.08613578   0.09323165   0.35756428   0.25403493
[21]   0.14375222   0.11315116   0.11241344   0.03704815   0.12805265
[26]   0.05922301   0.14368160   0.06504510   0.44795261   0.08632247
[31]   0.07419565   0.35569425   0.08513071   0.12207266   0.04828226
[36]   0.07304558   0.07435988   0.06792582   0.05937370   0.10827094
[41]   0.21430708   1.28573887   0.02740159   0.10619491   1.38923493
[46]   0.26800790   0.27805913   0.27806554   0.07879766   0.14102941
[51]   0.39052389   0.08982980   0.15872279   0.17735826   0.12601166
[56]   0.08164121   6.89159980   0.43838287   0.46470958   0.89579887
[61]   0.36878196   1.52229755   0.71595552   0.57227080   0.49934867
[66]   2.00364842   0.29531268   7.31765007   1.90864101   0.09352414
[71]   1.28573887   2.53417668   1.43656366   0.24760979   0.41366511
[76]   1.46681930   1.80377799   1.04091624   0.05380667   0.17995576
[81]   0.13074055   0.18007876   1.33455272   0.23941348   1.31278772
[86]   0.15193422   0.40575200   0.52799318   0.61701292   2.68516850
[91]   1.58009249   2.05840242   0.41321756   0.16692950   0.11265189
[96]   5.85199422   0.13072850   1.74029800   0.89956176   0.53228300
[101]   1.49479436   0.51734466   4.92009247   1.21235927   0.42130538
[106]   0.79911269   0.59143095   0.66516120   1.58054137   1.14709744
[111]   0.27818298   0.28924135   0.81981082   0.42204997   0.84835882
[116]   2.18361003   2.48969131   0.15490086   1.94586396   0.61742043
[121]   1.84224835   2.22216563   1.16831089   0.25570773   1.06363906
[126]   0.11017759   0.31847704   1.73906924 236.80047306   1.04552952
[131]   0.71491192   3.31452879   0.54903263  15.56796569 117.75340276
[136]   6.89922294  25.94175269   5.61212163   1.92017850   3.86127465
[141]   9.33199040   3.19081830  13.87818390   6.75546260   6.88705121
[146]   2.02107380  27.90562252   7.44022727   9.32848045   3.06868455
[151] 104.88115464  43.13367980   4.55163006 124.35055369  31.28870201
[156]  25.94055751  78.88807278  20.66781100   7.45074491  20.29440783
[161]  18.45010077   3.76962888  14.67530276  19.51116438   2.43700853
[166]  10.67116195  14.71094846   4.22513952   1.79979410   1.67202854
[171]  12.16406643   4.47296628   7.43251999   1.80377799   8.65323018
[176]   5.50514631   8.19937467   2.18052759   4.91292080   8.03215988
[181]   1.49201973   9.86753570   1.35859638   2.38734662   8.49630599
[186]  27.90433687  14.95978216   3.00707267   4.72588644  35.08515480
[191]   1.04520851  10.86800471  27.94507043   9.16578236   4.55163006
[196]   1.25886070   4.50591302

> predict(model5, type="expected")
         1          2          3          4          5          6          7
0.09741257 0.03727077 0.02125350 0.05770012 0.26029088 0.09055712 0.03953877
         8          9         10         11         12         13         14
0.07803524 0.07787007 0.87572766 0.03666546 0.06959290 0.15129120 0.48469256
        15         16         17         18         19         20         21
0.08102277 0.08076566 0.07469767 0.08085127 0.31008276 0.22030123 0.06662098
        22         23         24         25         26         27         28
0.09812564 0.09748588 0.03212846 0.11104834 0.05135869 0.12460190 0.05640766
        29         30         31         32         33         34         35
0.38846828 0.07485957 0.06434309 0.30371060 0.07382607 0.06590057 0.04187079
        36         37         38         39         40         41         42
0.06334574 0.06448551 0.05890585 0.05148937 0.09389348 0.18584891 0.92114753
        43         44         45         46         47         48         49
0.02376289 0.09209312 1.20475624 0.23241871 0.24113522 0.24114078 0.06833400
        50         51         52         53         54         55         56
0.12230189 0.33866561 0.07790116 0.13764574 0.15380657 0.10927837 0.07079994
        57         58         59         60         61         62         63
1.54319697 0.18171314 0.40300006 0.77684433 0.31981083 0.06927014 0.62088266
        64         65         66         67         68         69         70
0.48863513 0.43303938 0.97227876 0.03250515 4.78395309 0.16537057 0.08110492
        71         72         73         74         75         76         77
0.14694725 0.90498987 0.45504521 0.21472930 0.35873387 1.27203806 0.93033237
        78         79         80         81         82         83         84
0.83406021 0.04666160 0.15605915 0.11337930 0.15616582 0.07790488 0.20762139
        85         86         87         88         89         90         91
1.13846058 0.13175863 0.35187155 0.15575776 0.02807636 0.57789685 1.37026952
        92         93         94         95         96         97         98
0.54191702 0.35834575 0.10913106 0.09769267 0.26628726 0.11336886 0.87068667
        99        100        101        102        103        104        105
0.22317478 0.03107217 0.65040068 0.30108588 0.01077213 1.05136817 0.36535957
       106        107        108        109        110        111        112
0.69299726 0.51289390 0.57683340 1.37065880 0.21043827 0.24124262 0.08223406
       113        114        115        116        117        118        119
0.30781553 0.36600529 0.60779342 0.66763945 0.68130989 0.13433133 0.90179729
       120        121        122        123        124        125        126
0.40364193 0.65789260 0.74146453 0.08105875 0.22175190 0.65599925 0.09554693
       127        128        129        130        131        132        133
0.03371261 0.68596307 0.69465773 0.83775672 0.06455349 1.78933885 0.47612572
       134        135        136        137        138        139        140
0.34570904 1.65642691 0.04566908 0.01860585 0.73234819 1.66519498 3.34852992
       141        142        143        144        145        146        147
1.39505890 0.01814069 0.23856851 0.93426907 1.60596606 1.75269223 3.37110028
       148        149        150        151        152        153        154
0.11621245 1.21731071 0.37070821 1.21628368 0.40837624 3.26094425 1.44206602
       155        156        157        158        159        160        161
1.01897196 0.57604734 1.75182298 0.17576492 1.15893131 2.06286055 0.51877603
       162        163        164        165        166        167        168
0.04371556 0.56900548 0.02825276 0.42917254 1.00264005 0.47908808 0.71357947
       169        170        171        172        173        174        175
0.44651592 1.44999724 1.00980925 0.26111063 0.28818108 0.60186215 0.50513469
       176        177        178        179        180        181        182
0.02364640 0.61003936 0.17344290 0.03714014 0.15097160 0.86833034 1.96093102
       183        184        185        186        187        188        189
1.17818623 0.17762030 1.22144772 0.90875412 2.42176769 0.59758204 1.86409115
       190        191        192        193        194        195        196
2.04810554 0.90641363 0.63442276 0.12003320 0.23883016 0.90452498 0.54774348
       197
0.44053876

> predict(model5, type="terms")
           efmss         gal          cys
1   -0.076855836 -1.81712524 -0.292343064
2   -0.640401728 -2.10055649 -0.406111708
3   -1.542075156 -1.72264816 -0.444034589
4   -1.091238442 -1.15578565 -0.462996030
5   -0.302274193 -0.49444606 -0.007921455
6   -0.452553097 -1.62817107 -0.178574420
7   -0.978529263 -1.81712524 -0.292343064
8   -0.828250359 -1.43921690 -0.140651539
9   -0.602832002 -1.43921690 -0.368188827
10   0.110992795  0.16689353 -0.140651539
11  -1.091238442 -1.62817107 -0.444034589
12  -0.715541180 -1.53369399 -0.273381624
13  -0.903389811 -0.77787731 -0.064805777
14   0.298841426 -0.77787731 -0.102728658
15  -0.790680633 -1.43921690 -0.140651539
16  -0.452553097 -1.43921690 -0.481957471
17  -0.001716383 -2.00607941 -0.444034589
18  -0.565262276 -1.43921690 -0.368188827
19  -1.016098989  0.35584770 -0.368188827
20  -0.677971454 -0.39996898 -0.292343064
21  -0.339843919 -1.15578565 -0.444034589
22  -1.203947620 -0.87235440 -0.102728658
23  -0.227134740 -1.62817107 -0.330265945
24  -1.053668716 -1.91160233 -0.330265945
25   0.035853343 -1.62817107 -0.462996030
26  -0.527692550 -1.91160233 -0.387150267
27  -0.076855836 -1.72264816 -0.140651539
28  -0.452553097 -1.81712524 -0.462996030
29  -0.715541180  0.16689353 -0.254420183
30  -0.302274193 -1.81712524 -0.330265945
31  -0.264704466 -2.00607941 -0.330265945
32  -0.001716383 -0.87235440 -0.159612980
33  -1.016098989 -1.34473982 -0.102728658
34  -1.279087073 -0.68340023 -0.140651539
35  -1.128808168 -1.53369399 -0.368188827
36  -0.865820085 -1.34473982 -0.406111708
37  -0.640401728 -1.62817107 -0.330265945
38  -1.203947620 -1.34473982 -0.140651539
39  -0.903389811 -1.62817107 -0.292343064
40  -0.302274193 -1.53369399 -0.387150267
41  -0.565262276 -0.87235440 -0.102728658
42   0.110992795  0.07241645  0.067924308
43  -1.241517346 -1.91160233 -0.444034589
44  -0.189565014 -1.72264816 -0.330265945
45  -0.302274193  0.73375604 -0.102728658
46  -0.076855836 -1.06130857 -0.178574420
47   0.110992795 -1.25026274 -0.140651539
48   0.073423069 -1.15578565 -0.197535861
49  -0.790680633 -1.53369399 -0.216497302
50  -0.076855836 -1.81712524 -0.064805777
51   0.110992795 -0.77787731 -0.273381624
52  -0.715541180 -1.25026274 -0.444034589
53  -0.828250359 -1.06130857  0.048962867
54  -0.640401728 -0.77787731 -0.311304505
55  -0.377413645 -1.34473982 -0.349227386
56  -0.527692550 -1.53369399 -0.444034589
57   0.674538687  1.01718729  0.238577274
58  -0.189565014 -0.49444606 -0.140651539
59  -0.377413645 -0.39996898  0.011039986
60   0.298841426 -0.11653772 -0.292343064
61   0.223701973 -0.96683148 -0.254420183
62   0.298841426  0.07241645  0.048962867
63   0.073423069 -0.49444606  0.086885749
64  -0.302274193 -0.49444606  0.238577274
65   0.148562521 -0.68340023 -0.159612980
66   0.073423069  0.73375604 -0.112209378
67  -0.452553097 -0.68340023 -0.083767217
68  -0.076855836  2.15091230 -0.083767217
69   0.373980878  0.26137061  0.011039986
70  -0.828250359 -1.62817107  0.086885749
71   0.110992795  0.07241645  0.067924308
72   0.298841426  0.73375604 -0.102728658
73   0.449120330 -0.02206064 -0.064805777
74   0.486690057 -1.62817107 -0.254420183
75   0.073423069 -0.87235440 -0.083767217
76   0.073423069  0.45032478 -0.140651539
77   0.298841426  0.35584770 -0.064805777
78   0.674538687 -0.68340023  0.048962867
79  -0.339843919 -2.10055649 -0.481957471
80   0.110992795 -1.53369399 -0.292343064
81  -0.227134740 -1.43921690 -0.368188827
82   0.073423069 -1.62817107 -0.159612980
83   0.110992795  0.26137061 -0.083767217
84  -0.001716383 -1.53369399  0.105847189
85  -0.227134740  0.45032478  0.048962867
86  -0.001716383 -1.62817107 -0.254420183
87   0.110992795 -0.87235440 -0.140651539
88   0.336411152 -0.87235440 -0.102728658
89  -0.452553097  0.07241645 -0.102728658
90   0.298841426  0.45032478  0.238577274
91   0.298841426  0.26137061 -0.102728658
92   0.336411152  0.45032478 -0.064805777
93   0.223701973 -0.96683148 -0.140651539
94   0.035853343 -1.53369399 -0.292343064
95  -0.452553097 -1.62817107 -0.102728658
96  -0.189565014  1.30061854  0.655728968
97  -0.076855836 -1.81712524 -0.140651539
98   0.110992795  0.26137061  0.181692952
99  -0.076855836 -0.30549189  0.276500155
100 -1.016098989  0.45032478 -0.064805777
101  0.186132247  0.16689353  0.048962867
102  0.486690057 -0.87235440 -0.273381624
103  0.073423069  1.39509563  0.124808630
104  0.486690057 -0.39996898  0.105847189
105  0.073423069 -0.68340023 -0.254420183
106  0.486690057 -0.49444606 -0.216497302
107  0.298841426 -0.68340023 -0.140651539
108 -0.377413645  0.07241645 -0.102728658
109  0.373980878 -0.02206064  0.105847189
110  0.110992795  0.16689353 -0.140651539
111 -0.076855836 -0.87235440 -0.330265945
112 -0.189565014 -0.87235440 -0.178574420
113 -0.602832002  0.54480187 -0.140651539
114 -0.114425562 -0.68340023 -0.064805777
115  0.035853343 -0.11653772 -0.083767217
116  0.035853343  0.63927895  0.105847189
117  0.411550604  0.07241645  0.428191680
118 -0.076855836 -1.53369399 -0.254420183
119  0.336411152  0.26137061  0.067924308
120 -0.452553097 -0.11653772  0.086885749
121  0.035853343  0.45032478  0.124808630
122  0.261271700  0.45032478  0.086885749
123  0.073423069  0.45032478 -0.368188827
124 -0.227134740 -0.77787731 -0.358708106
125  0.148562521 -0.02206064 -0.064805777
126 -0.001716383 -1.91160233 -0.292343064
127 -0.452553097 -0.58892314 -0.102728658
128  0.186132247  0.26137061  0.105847189
129  0.073423069  5.17417900  0.219615833
130 -0.076855836  0.07241645  0.048962867
131  0.298841426 -0.68340023  0.048962867
132 -0.264704466  1.39509563  0.067924308
133  0.073423069 -0.49444606 -0.178574420
134  0.298841426  2.15091230  0.295461596
135  0.110992795  4.22940816  0.428191680
136  0.486690057  1.20614146  0.238577274
137  0.298841426  2.52882064  0.428191680
138  0.298841426  1.11166437  0.314423036
139 -0.302274193  0.35584770  0.598844646
140 -0.264704466  1.11166437  0.504037443
141  0.524259783  1.48957271  0.219615833
142 -0.076855836  0.92271020  0.314423036
143  0.524259783  1.86748105  0.238577274
144  0.674538687  1.30061854 -0.064805777
145  0.674538687  1.20614146  0.048962867
146  0.336411152  0.26137061  0.105847189
147  0.674538687  2.33986647  0.314423036
148  0.524259783  1.20614146  0.276500155
149  0.599399235  1.39509563  0.238577274
150  0.110992795  0.73375604  0.276500155
151  0.674538687  3.28463731  0.693651849
152  0.636968961  2.62329772  0.504037443
153  0.636968961  0.45032478  0.428191680
154  0.749678140  3.19016023  0.883266255
155  0.674538687  2.15091230  0.617806086
156  0.373980878  2.33986647  0.541960324
157  0.749678140  3.19016023  0.428191680
158  0.411550604  2.15091230  0.466114561
159  0.373980878  1.20614146  0.428191680
160  0.298841426  2.24538938  0.466114561
161  0.524259783  1.77300396  0.617806086
162  0.449120330  0.63927895  0.238577274
163  0.599399235  1.96195813  0.124808630
164  0.486690057  2.15091230  0.333384477
165  0.599399235  0.26137061  0.030001427
166  0.449120330  1.30061854  0.617806086
167  0.411550604  1.77300396  0.504037443
168  0.486690057  0.45032478  0.504037443
169  0.674538687 -0.02206064 -0.064805777
170  0.298841426  0.35584770 -0.140651539
171  0.674538687  1.20614146  0.617806086
172  0.298841426  0.92271020  0.276500155
173  0.599399235  1.30061854  0.105847189
174  0.298841426  0.35584770 -0.064805777
175  0.524259783  1.39509563  0.238577274
176  0.223701973  1.39509563  0.086885749
177 -0.076855836  2.15091230  0.030001427
178  0.186132247  0.63927895 -0.045844336
179  0.298841426  1.20614146  0.086885749
180  0.449120330  1.20614146  0.428191680
181  0.524259783 -0.21101481  0.086885749
182  0.674538687  1.39509563  0.219615833
183  0.298841426  0.07241645 -0.064805777
184  0.298841426  0.45032478  0.121016342
185  0.524259783  1.20614146  0.409230239
186  0.749678140  2.15091230  0.428191680
187  0.749678140  1.48957271  0.466114561
188  0.110992795  1.11166437 -0.121690098
189  0.674538687  0.45032478  0.428191680
190  0.524259783  2.33986647  0.693651849
191 -0.114425562  0.26137061 -0.102728658
192  0.486690057  1.39509563  0.504037443
193  0.524259783  2.33986647  0.466114561
194  0.524259783  1.20614146  0.485076002
195  0.636968961  0.45032478  0.428191680
196  0.411550604 -0.11653772 -0.064805777
197  0.561829509  0.82823312  0.115327909
attr(,"constant")
[1] 1.967401

[propedevt]

Да, до predict все делал точно также

[p2004r]

Теперь basehazard и счет на прямую. Центрированная версия равна случаю, когда h0(t) hazard считается для случая когда z<-0. Первый вариант считает когда значения предикторов заменены их средними.

Код

> basehaz(model5, centered=FALSE)
         hazard time
1  0.0001002811   16
2  0.0002024633   45
3  0.0003061236   48
4  0.0004101635   55
5  0.0006646868   60
6  0.0007949146   62
7  0.0009255306   69
8  0.0010569939   75
9  0.0011890672   76
10 0.0013237703   83
11 0.0017833064   90
12 0.0019668331   91
13 0.0021839121   92
14 0.0024035279   99
15 0.0026280378  100
16 0.0033604084  120
17 0.0036432445  126
18 0.0039314237  140
19 0.0048821851  150
20 0.0056043508  160
21 0.0067440470  165
22 0.0092375477  180
23 0.0097008584  210
24 0.0106394178  220
25 0.0111215067  240
26 0.0116072477  250
27 0.0121144210  254
28 0.0126251301  267
29 0.0131371764  270
30 0.0136700471  278
31 0.0142122284  280
32 0.0148007218  300
33 0.0153900052  307
34 0.0159800228  356
35 0.0172082186  380
36 0.0186000766  390
37 0.0193368429  400
38 0.0201008041  420
39 0.0209019571  430
40 0.0217483690  440
41 0.0226347604  450
42 0.0236140231  475
43 0.0246231476  480
44 0.0256503391  488
45 0.0289022541  500
46 0.0300917523  505
47 0.0313090582  540
48 0.0326040644  550
49 0.0353892636  560
50 0.0368104561  570
51 0.0382620189  580
52 0.0397520960  590
53 0.0412467805  605
54 0.0427499492  620
55 0.0442892539  650
56 0.0474533644  660
57 0.0507720453  680
58 0.0524983401  690
59 0.0560595003  700
60 0.0579564065  705
61 0.0618094617  710
62 0.0658048238  720
63 0.0678481443  735
64 0.0699531004  740
65 0.0721146898  745
66 0.0766243540  750
67 0.0837802983  760
68 0.0862338454  766
69 0.0941286515  770
70 0.1033435078  780
71 0.1067425188  785
72 0.1138292088  790
73 0.1212529131  800
74 0.1212529131  805
> basehaz(model5)
         hazard time
1  0.0007172164   16
2  0.0014480302   45
3  0.0021894155   48
4  0.0029335149   55
5  0.0047538818   60
6  0.0056852791   62
7  0.0066194528   69
8  0.0075596863   75
9  0.0085042832   76
10 0.0094676883   83
11 0.0127543187   90
12 0.0140669133   91
13 0.0156194759   92
14 0.0171901823   99
15 0.0187958908  100
16 0.0240338513  120
17 0.0260567126  126
18 0.0281177885  140
19 0.0349176933  150
20 0.0400826682  160
21 0.0482338464  165
22 0.0660675194  180
23 0.0693811469  210
24 0.0760937821  220
25 0.0795417127  240
26 0.0830157622  250
27 0.0866430977  254
28 0.0902957213  267
29 0.0939579080  270
30 0.0977690327  278
31 0.1016467478  280
32 0.1058556895  300
33 0.1100702806  307
34 0.1142901239  356
35 0.1230742570  380
36 0.1330289128  390
37 0.1382983114  400
38 0.1437622096  420
39 0.1494921065  430
40 0.1555456970  440
41 0.1618852240  450
42 0.1688889721  475
43 0.1761062938  480
44 0.1834528314  488
45 0.2067107316  500
46 0.2152180967  505
47 0.2239243447  540
48 0.2331863105  550
49 0.2531062296  560
50 0.2632706871  570
51 0.2736523556  580
52 0.2843094803  590
53 0.2949995571  605
54 0.3057503141  620
55 0.3167595173  650
56 0.3393894339  660
57 0.3631248475  680
58 0.3754714163  690
59 0.4009410571  700
60 0.4145078485  705
61 0.4420651411  710
62 0.4706402210  720
63 0.4852541769  735
64 0.5003089551  740
65 0.5157687777  745
66 0.5480221780  750
67 0.5992019402  760
68 0.6167498625  766
69 0.6732140103  770
70 0.7391192393  780
71 0.7634291793  785
72 0.8141136296  790
73 0.8672084276  800
74 0.8672084276  805
> -log((survfit(model5))$surv)
[1] 0.0007172164 0.0014480302 0.0021894155 0.0029335149 0.0047538818
[6] 0.0056852791 0.0066194528 0.0075596863 0.0085042832 0.0094676883
[11] 0.0127543187 0.0140669133 0.0156194759 0.0171901823 0.0187958908
[16] 0.0240338513 0.0260567126 0.0281177885 0.0349176933 0.0400826682
[21] 0.0482338464 0.0660675194 0.0693811469 0.0760937821 0.0795417127
[26] 0.0830157622 0.0866430977 0.0902957213 0.0939579080 0.0977690327
[31] 0.1016467478 0.1058556895 0.1100702806 0.1142901239 0.1230742570
[36] 0.1330289128 0.1382983114 0.1437622096 0.1494921065 0.1555456970
[41] 0.1618852240 0.1688889721 0.1761062938 0.1834528314 0.2067107316
[46] 0.2152180967 0.2239243447 0.2331863105 0.2531062296 0.2632706871
[51] 0.2736523556 0.2843094803 0.2949995571 0.3057503141 0.3167595173
[56] 0.3393894339 0.3631248475 0.3754714163 0.4009410571 0.4145078485
[61] 0.4420651411 0.4706402210 0.4852541769 0.5003089551 0.5157687777
[66] 0.5480221780 0.5992019402 0.6167498625 0.6732140103 0.7391192393
[71] 0.7634291793 0.8141136296 0.8672084276 0.8672084276

hazard у одиночного случая по времени
h(t)=h0(t)*exp(sum(z_случая*coef(fit)))

считаем формулу по шагам для первого случая выборки

Код

> coef(model5)
        efmss           gal           cys
-3.7569726150  0.0944770843  0.0001896144
> data[1,c(7,4,6)]
  efmss gal  cys
1   0.3   6 1800
> data[1,c(7,4,6)]*coef(model5)
      efmss       gal       cys
1 -1.127092 0.5668625 0.3413059
> sum(data[1,c(7,4,6)]*coef(model5))
[1] -0.2189233
> exp(sum(data[1,c(7,4,6)]*coef(model5)))
[1] 0.8033833
>  -log((survfit(model5))$surv)*exp(sum(data[1,c(7,4,6)]*coef(model5)))
[1] 0.0005761997 0.0011633233 0.0017589398 0.0023567369 0.0038191892
[6] 0.0045674583 0.0053179578 0.0060733257 0.0068321990 0.0076061827
[11] 0.0102466066 0.0113011232 0.0125484260 0.0138103053 0.0151003048
[16] 0.0193083947 0.0209335277 0.0225893616 0.0280522916 0.0322017461
[21] 0.0387502666 0.0530775416 0.0557396545 0.0611324736 0.0639024834
[26] 0.0666934767 0.0696076175 0.0725420743 0.0754842140 0.0785460079
[31] 0.0816612994 0.0850426929 0.0884286250 0.0918187766 0.0988758024
[36] 0.1068732067 0.1111065535 0.1154961580 0.1200994615 0.1249628150
[41] 0.1300558851 0.1356825793 0.1414808550 0.1473829406 0.1660679492
[46] 0.1729026242 0.1798970785 0.1873379871 0.2033413173 0.2115072727
[51] 0.2198477318 0.2284094878 0.2369977169 0.2456346956 0.2544793055
[56] 0.2726598026 0.2917284374 0.3016474645 0.3221093486 0.3330086822
[61] 0.3551477507 0.3781044927 0.3898451007 0.4019398581 0.4143600213
[66] 0.4402718645 0.4813888306 0.4954865382 0.5408488915 0.5937960517
[71] 0.6133262514 0.6540452922 0.6967007661 0.6967007661

[propedevt]

Действия математические понятны. Но каков смысл predict и basehazard? они не используются пока, какие умозаключения мы сделали посчитав их?

..
Кажется понял это делали для того чтобы увидеть, что basehaz(model5, centered=TRUE) равна -log((survfit(model5))$surv)

Сообщение отредактировал propedevt - 9.03.2012 - 18:04

[p2004r]

Код

> exp(model5$linear.predictors)
  [1]   0.11232890   0.04297787   0.02450795   0.06653547   0.44724814
  [6]   0.10442370   0.04559316   0.08998441   0.08979395   1.14709744
[11]   0.04227987   0.08024934   0.17445771   0.55891127   0.09342941
[16]   0.09313293   0.08613578   0.09323165   0.35756428   0.25403493
[21]   0.14375222   0.11315116   0.11241344   0.03704815   0.12805265
[26]   0.05922301   0.14368160   0.06504510   0.44795261   0.08632247
[31]   0.07419565   0.35569425   0.08513071   0.12207266   0.04828226
[36]   0.07304558   0.07435988   0.06792582   0.05937370   0.10827094
[41]   0.21430708   1.28573887   0.02740159   0.10619491   1.38923493
[46]   0.26800790   0.27805913   0.27806554   0.07879766   0.14102941
[51]   0.39052389   0.08982980   0.15872279   0.17735826   0.12601166
[56]   0.08164121   6.89159980   0.43838287   0.46470958   0.89579887
[61]   0.36878196   1.52229755   0.71595552   0.57227080   0.49934867
[66]   2.00364842   0.29531268   7.31765007   1.90864101   0.09352414
[71]   1.28573887   2.53417668   1.43656366   0.24760979   0.41366511
[76]   1.46681930   1.80377799   1.04091624   0.05380667   0.17995576
[81]   0.13074055   0.18007876   1.33455272   0.23941348   1.31278772
[86]   0.15193422   0.40575200   0.52799318   0.61701292   2.68516850
[91]   1.58009249   2.05840242   0.41321756   0.16692950   0.11265189
[96]   5.85199422   0.13072850   1.74029800   0.89956176   0.53228300
[101]   1.49479436   0.51734466   4.92009247   1.21235927   0.42130538
[106]   0.79911269   0.59143095   0.66516120   1.58054137   1.14709744
[111]   0.27818298   0.28924135   0.81981082   0.42204997   0.84835882
[116]   2.18361003   2.48969131   0.15490086   1.94586396   0.61742043
[121]   1.84224835   2.22216563   1.16831089   0.25570773   1.06363906
[126]   0.11017759   0.31847704   1.73906924 236.80047306   1.04552952
[131]   0.71491192   3.31452879   0.54903263  15.56796569 117.75340276
[136]   6.89922294  25.94175269   5.61212163   1.92017850   3.86127465
[141]   9.33199040   3.19081830  13.87818390   6.75546260   6.88705121
[146]   2.02107380  27.90562252   7.44022727   9.32848045   3.06868455
[151] 104.88115464  43.13367980   4.55163006 124.35055369  31.28870201
[156]  25.94055751  78.88807278  20.66781100   7.45074491  20.29440783
[161]  18.45010077   3.76962888  14.67530276  19.51116438   2.43700853
[166]  10.67116195  14.71094846   4.22513952   1.79979410   1.67202854
[171]  12.16406643   4.47296628   7.43251999   1.80377799   8.65323018
[176]   5.50514631   8.19937467   2.18052759   4.91292080   8.03215988
[181]   1.49201973   9.86753570   1.35859638   2.38734662   8.49630599
[186]  27.90433687  14.95978216   3.00707267   4.72588644  35.08515480
[191]   1.04520851  10.86800471  27.94507043   9.16578236   4.55163006
[196]   1.25886070   4.50591302
>

риск считает как экспоненту от линейного предиктора

[propedevt]

Код

> exp(model5$linear.predictors)
  [1]   0.11232890   0.04297787   0.02450795   0.06653547   0.44724814
  ...
>

риск считает как экспоненту от линейного предиктора

Вижу predict(model5, type="risk") тоже самое выдает

[p2004r]

Линейный предиктор он считает в момент фита


names(coef) <- dimnames(x)[[2]]
lp <- c(x %*% coef) + offset - sum(coef*coxfit$means)

Код

# напоминаю что было вот так
195  0.636968961  0.45032478  0.428191680
196  0.411550604 -0.11653772 -0.064805777
197  0.561829509  0.82823312  0.115327909
attr(,"constant")
[1] 1.967401
# константа получается вот так
> model5$mean
[1]    0.2795431   25.2335025 3341.7766497
> model5$mean*model5$coefficients
    efmss       gal       cys
-1.050236  2.383988  0.633649
> sum(model5$mean*model5$coefficients)
[1] 1.967401

теперь линейный предиктор для 197го случая

Код

# так он в принципе и все сразу посчитать может, а не только один случай
>  c(as.numeric(data[197,c(7,4,6)])%*%as.numeric(model5$coefficients)) - sum(model5$mean*model5$coefficients)
[1] 1.505391

# вот оно же "по простому", для одиночного случая

>  sum(data[197,c(7,4,6)]*model5$coefficients) - sum(model5$mean*model5$coefficients)
[1] 1.505391

Все рассчитывается за вычетом средней по выборке.

[p2004r]

Конкретно в этой модели первая переменная похоже не выполняет условия пропорциональности риска.

[propedevt]

Линейный предиктор он считает в момент фита
...

Все рассчитывается за вычетом средней по выборке.

Как это теперь подсчитать для произвольных значений предикторов?

[propedevt]

Конкретно в этой модели первая переменная похоже не выполняет условия пропорциональности риска.

Хм, при проверке cox.zph она немного больше 0,05 выдала

[p2004r]

Хм, при проверке cox.zph она немного больше 0,05 выдала

совсем немного я бы сказал, и на взгляд оно как то ступенькой уж очень

[propedevt]

давайте уберем ее из модели, построим только на gal и cys

[p2004r]

Как это теперь подсчитать для произвольных значений предикторов?

так и считать, за вычетом среднего.... В predict() подставить датафрейм с новыми значениями, оно все применит само

[propedevt]

Так понятно, но я вижу predict не используется же в расчетах, и потом я хочу вне R написать на javascript или php в виде интерфейса, подставил значения - получил прогноз.

[propedevt]

Итак хочу написать в математических операциях.
начальная формула: h(t)=h0(t)*exp(sum(z_случая*coef(fit)))
в итоге она сводится к такому выражению -log((survfit(model5))$surv)*exp(sum(data[1,c(7,4,6)]*coef(model5)))

Напишу ее так
h(t)= h0(t) * exp((-3.7569726150 * efmss) + (0.0944770843 * gal) + (0.0001896144 * cys))

где h0(t) = -log((survfit(model5))$surv) для случая.

Но я ведь хочу получить h(t) в произвольное время, например 1 месяц, 2, 6, 12 месяцев. то как быть?
Может просто предрасчитать h0(t) для определенных периодов времени? для 1, 6, 12, 24 месяцев?
и написать 4 выражения
h(t)1 = h0(t) для 1 месяца * exp((-3.7569726150 * efmss) + (0.0944770843 * gal) + (0.0001896144 * cys))
...
h(t)24 = h0(t) для 24 месяцев * exp((-3.7569726150 * efmss) + (0.0944770843 * gal) + (0.0001896144 * cys))

И еще у меня сомнения почему http://www.medcalc.org/manual/cox_proportional_hazards.php - тут формула как я уже приводил -, а в источнике http://cran.r-project.org/doc/contrib/Fox-...-regression.pdf как мы и применяли h(t)=h0(t)*exp(sum(z_случая*coef(fit)))
Почему формулы разные?

[p2004r]

Итак хочу написать в математических операциях.
начальная формула: h(t)=h0(t)*exp(sum(z_случая*coef(fit)))
в итоге она сводится к такому выражению -log((survfit(model5))$surv)*exp(sum(data[1,c(7,4,6)]*coef(model5)))

Напишу ее так
h(t)= h0(t) * exp((-3.7569726150 * efmss) + (0.0944770843 * gal) + (0.0001896144 * cys))

где h0(t) = -log((survfit(model5))$surv) для случая.

Но я ведь хочу получить h(t) в произвольное время, например 1 месяц, 2, 6, 12 месяцев. то как быть?

h0(t) это вектор

Код

> -log((survfit(model5))$surv)-log((survfit(model5))$surv)
[1] 0.001434433 0.002896060 0.004378831 0.005867030 0.009507764 0.011370558
[7] 0.013238906 0.015119373 0.017008566 0.018935377 0.025508637 0.028133827
[13] 0.031238952 0.034380365 0.037591782 0.048067703 0.052113425 0.056235577
[19] 0.069835387 0.080165336 0.096467693 0.132135039 0.138762294 0.152187564
[25] 0.159083425 0.166031524 0.173286195 0.180591443 0.187915816 0.195538065
[31] 0.203293496 0.211711379 0.220140561 0.228580248 0.246148514 0.266057826
[37] 0.276596623 0.287524419 0.298984213 0.311091394 0.323770448 0.337777944
[43] 0.352212588 0.366905663 0.413421463 0.430436193 0.447848689 0.466372621
[49] 0.506212459 0.526541374 0.547304711 0.568618961 0.589999114 0.611500628
[55] 0.633519035 0.678778868 0.726249695 0.750942833 0.801882114 0.829015697
[61] 0.884130282 0.941280442 0.970508354 1.000617910 1.031537555 1.096044356
[67] 1.198403880 1.233499725 1.346428021 1.478238479 1.526858359 1.628227259
[73] 1.734416855 1.734416855

и в predict() вычитается предиктор для средних выборки по которой шло обучение.

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

[propedevt]

Т.к. убрал из модели efmss, добавил в нее возраст. Он подошел под требования пропорциональности, в итоге модель составил на gal+cys+age

Далее пошел эмпирическим путем слегка, взял формулу как в существующих моделях, и получил для себя такую формулу:


,где PI=(0,0985149930 * gal)+(0,0002328366 * cys)+(0,0328478062 *age)
H0(t) - риск в момент времени t, когда ковариаты равны нулю.
и я наверно этим уже надоел:) но опять беру его с графика полученного
baseline <- basehaz(model6, centered=FALSE)
plot(baseline$time,baseline$hazard,xlab="Time",ylab="Hazard",col="red",type="l")


И тогда у меня все получается как хотелось.

И у меня поэтому наверно последний вопрос: нельзя ли в виде математического выражения написать что выдает мне basehaz(model6, centered=FALSE) ? чтобы для любого времени от нуля до максимального периода наблюдения получать h0(t), в любой программе, на любом языке программирования.
или если это нельзя написать математически, то как в R считать h0(t) для какого-то определенного времени t чтобы выдало аналогию что выдает basehaz(model6, centered=FALSE)

Сообщение отредактировал propedevt - 10.03.2012 - 13:44

[p2004r]

И у меня поэтому наверно последний вопрос: нельзя ли в виде математического выражения написать что выдает мне basehaz(model6, centered=FALSE) ? чтобы для любого времени от нуля до максимального периода наблюдения получать h0(t), в любой программе, на любом языке программирования.
или если это нельзя написать математически, то как в R считать h0(t) для какого-то определенного времени t чтобы выдало аналогию что выдает basehaz(model6, centered=FALSE)

Код

        sfit <- survfit(fit)
        H <- -log(sfit$surv)
        z0 <- fit$means
        bz0 <- sum(z0 * coef(fit))
        H <- H * exp(-bz0)

считает именно вот этот метод
?survfit.coxph

пример из его хелпа

Код

     #fit a Cox proportional hazards model and plot the  
     #predicted survival for a 60 year old
     fit <- coxph(Surv(futime, fustat) ~ age, data = ovarian)
     plot(survfit(fit, newdata=data.frame(age=60)),
          xscale=365.25, xlab = "Years", ylab="Survival")