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

[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")