Здравствуйте, гость ( Вход | Регистрация )
Форум врачей-аспирантов » Все сообщения пользователя
| Отправлено: 10.03.2012 - 13:43 | |
|
Группа: Пользователи Сообщений: 27 Регистрация: 5.02.2012 Пользователь №: 23464 
|
Т.к. убрал из модели 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) |
| Форум: Медицинская статистика · Просмотр сообщения: #13174 · Ответов: 48 · Просмотров: 114142 |
| Отправлено: 9.03.2012 - 23:23 | |
|
Группа: Пользователи Сообщений: 27 Регистрация: 5.02.2012 Пользователь №: 23464
|
Итак хочу написать в математических операциях. начальная формула: 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))) Почему формулы разные? |
| Форум: Медицинская статистика · Просмотр сообщения: #13172 · Ответов: 48 · Просмотров: 114142 |
| Отправлено: 9.03.2012 - 22:39 | |
|
Группа: Пользователи Сообщений: 27 Регистрация: 5.02.2012 Пользователь №: 23464
|
Так понятно, но я вижу predict не используется же в расчетах, и потом я хочу вне R написать на javascript или php в виде интерфейса, подставил значения - получил прогноз. |
| Форум: Медицинская статистика · Просмотр сообщения: #13171 · Ответов: 48 · Просмотров: 114142 |
| Отправлено: 9.03.2012 - 22:01 | |
|
Группа: Пользователи Сообщений: 27 Регистрация: 5.02.2012 Пользователь №: 23464
|
давайте уберем ее из модели, построим только на gal и cys |
| Форум: Медицинская статистика · Просмотр сообщения: #13169 · Ответов: 48 · Просмотров: 114142 |
| Отправлено: 9.03.2012 - 21:14 | |
|
Группа: Пользователи Сообщений: 27 Регистрация: 5.02.2012 Пользователь №: 23464
|
http://pubhealth.spb.ru/COPC/STATSH/odds.htm Вот в конце статьи про доверительный интервал как его подсчитывать указано. Но думаю проще логистическую регрессию применить в statistica, которую Вы используете. Углубленные методы анализа - нелинейное оценивание - логит регрессия |
| Форум: Медицинская статистика · Просмотр сообщения: #13167 · Ответов: 4 · Просмотров: 7848 |
| Отправлено: 9.03.2012 - 20:35 | |
|
Группа: Пользователи Сообщений: 27 Регистрация: 5.02.2012 Пользователь №: 23464
|
Цитата(p2004r @ 9.03.2012 - 19:42)  Конкретно в этой модели первая переменная похоже не выполняет условия пропорциональности риска. Хм, при проверке cox.zph она немного больше 0,05 выдала |
| Форум: Медицинская статистика · Просмотр сообщения: #13166 · Ответов: 48 · Просмотров: 114142 |
| Отправлено: 9.03.2012 - 20:26 | |
|
Группа: Пользователи Сообщений: 27 Регистрация: 5.02.2012 Пользователь №: 23464
|
Цитата(p2004r @ 9.03.2012 - 19:19) Линейный предиктор он считает в момент фита ... Все рассчитывается за вычетом средней по выборке. Как это теперь подсчитать для произвольных значений предикторов? |
| Форум: Медицинская статистика · Просмотр сообщения: #13165 · Ответов: 48 · Просмотров: 114142 |
| Отправлено: 9.03.2012 - 18:49 | |
|
Группа: Пользователи Сообщений: 27 Регистрация: 5.02.2012 Пользователь №: 23464
|
Цитата(p2004r @ 9.03.2012 - 18:36) Код > exp(model5$linear.predictors) [1] 0.11232890 0.04297787 0.02450795 0.06653547 0.44724814 ... > риск считает как экспоненту от линейного предиктора Вижу predict(model5, type="risk") тоже самое выдает |
| Форум: Медицинская статистика · Просмотр сообщения: #13162 · Ответов: 48 · Просмотров: 114142 |
| Отправлено: 9.03.2012 - 17:58 | |
|
Группа: Пользователи Сообщений: 27 Регистрация: 5.02.2012 Пользователь №: 23464
|
Действия математические понятны. Но каков смысл predict и basehazard? они не используются пока, какие умозаключения мы сделали посчитав их? .. Кажется понял это делали для того чтобы увидеть, что basehaz(model5, centered=TRUE) равна -log((survfit(model5))$surv) |
| Форум: Медицинская статистика · Просмотр сообщения: #13160 · Ответов: 48 · Просмотров: 114142 |
| Отправлено: 9.03.2012 - 17:18 | |
|
Группа: Пользователи Сообщений: 27 Регистрация: 5.02.2012 Пользователь №: 23464
|
Да, до predict все делал точно также |
| Форум: Медицинская статистика · Просмотр сообщения: #13158 · Ответов: 48 · Просмотров: 114142 |
| Отправлено: 9.03.2012 - 14:40 | |
|
Группа: Пользователи Сообщений: 27 Регистрация: 5.02.2012 Пользователь №: 23464
|
Хорошо, файл прикрепил. В последней модели использовал предикторы efmss, gal, cys. Жду Вас с нетерпением:) только что будете делать опишите так, чтобы я понял и мог сделать тоже самое:)
Прикрепленные файлы
|
| Форум: Медицинская статистика · Просмотр сообщения: #13156 · Ответов: 48 · Просмотров: 114142 |
| Отправлено: 9.03.2012 - 14:04 | |
|
Группа: Пользователи Сообщений: 27 Регистрация: 5.02.2012 Пользователь №: 23464
|
Что-то Вы меня совсем запутали:( перечитал что мне ответили по десять раз.. не знаю что делать дальше.. Давайте начнем с начала по порядку: Есть формула в нее я подставляю pi , вычисленный для конкретного больного по его данным. Далее я должен найти h0(t) базовый риск во время t. его например беру со своего второго графика. Правильно? или Вы хотите сказать можно вычислять риск с помощью функции predict() ? |
| Форум: Медицинская статистика · Просмотр сообщения: #13154 · Ответов: 48 · Просмотров: 114142 |
| Отправлено: 9.03.2012 - 11:29 | |
|
Группа: Пользователи Сообщений: 27 Регистрация: 5.02.2012 Пользователь №: 23464
|
Да понял, что если риск 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%) правильно? или я чего-то не понимаю |
| Форум: Медицинская статистика · Просмотр сообщения: #13147 · Ответов: 48 · Просмотров: 114142 |
| Отправлено: 9.03.2012 - 11:00 | |
|
Группа: Пользователи Сообщений: 27 Регистрация: 5.02.2012 Пользователь №: 23464
|
Отношение шансов между группами думаю? В простейшем варианте и для понимания сути лучше построить четырехпольную таблицу. Вот например что гугл находит - http://donbas-socproject.blogspot.com/2009...og-post_25.html посмотрите |
| Форум: Медицинская статистика · Просмотр сообщения: #13145 · Ответов: 4 · Просмотров: 7848 |
| Отправлено: 9.03.2012 - 10:04 | |
|
Группа: Пользователи Сообщений: 27 Регистрация: 5.02.2012 Пользователь №: 23464
|
Уважаемый 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) ? |
| Форум: Медицинская статистика · Просмотр сообщения: #13143 · Ответов: 48 · Просмотров: 114142 |
| Отправлено: 9.03.2012 - 08:15 | |
|
Группа: Пользователи Сообщений: 27 Регистрация: 5.02.2012 Пользователь №: 23464
|
Вычитал даты и получал дни жизни (до 805 дней максимум), в месяцы не преобразовывал. Так, банально подсчитал PI для средних, он равен (-3.75697 х 0,280) + (0.09448 х 25,234) + (0.00019 х 3341,777) = 1,96709435 Как теперь построить график функции риска? (Н0(t) подобный такому:  |
| Форум: Медицинская статистика · Просмотр сообщения: #13141 · Ответов: 48 · Просмотров: 114142 |
| Отправлено: 8.03.2012 - 23:09 | |
|
Группа: Пользователи Сообщений: 27 Регистрация: 5.02.2012 Пользователь №: 23464
|
Цитата(DrgLena @ 8.03.2012 - 22:27) Н0(t) это базовый риск в момент времени t, а не функция выживаемости. Да согласен, совершенно неправильно сказал Цитата(DrgLena @ 8.03.2012 - 22:27) Вот по этим средним и полученным коэффициентам и получаете базовый риск, вначале в единицу времени наблюдения А чем является единица времени наблюдения? У меня максимальный период наблюдения 26 месяцев, а в полученных днях время жизни в днях |
| Форум: Медицинская статистика · Просмотр сообщения: #13139 · Ответов: 48 · Просмотров: 114142 |
| Отправлено: 8.03.2012 - 21:10 | |
|
Группа: Пользователи Сообщений: 27 Регистрация: 5.02.2012 Пользователь №: 23464
|
Просто средние ковариат? тогда вот они: gal 25,234 cys 3341,777 efmss 0,280 |
| Форум: Медицинская статистика · Просмотр сообщения: #13135 · Ответов: 48 · Просмотров: 114142 |
| Отправлено: 8.03.2012 - 15:24 | |
|
Группа: Пользователи Сообщений: 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) Код который сейчас добавили сейчас опробую) |
| Форум: Медицинская статистика · Просмотр сообщения: #13131 · Ответов: 48 · Просмотров: 114142 |
| Отправлено: 8.03.2012 - 11:49 | |
|
Группа: Пользователи Сообщений: 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) исходя из собственных данных? |
| Форум: Медицинская статистика · Просмотр сообщения: #13129 · Ответов: 48 · Просмотров: 114142 |
| Отправлено: 4.03.2012 - 14:24 | |
|
Группа: Пользователи Сообщений: 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)) |
| Форум: Медицинская статистика · Просмотр сообщения: #13101 · Ответов: 48 · Просмотров: 114142 |
| Отправлено: 3.03.2012 - 21:40 | |
|
Группа: Пользователи Сообщений: 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 попробовал рисовать стандартными средствами) |
| Форум: Медицинская статистика · Просмотр сообщения: #13097 · Ответов: 48 · Просмотров: 114142 |
| Отправлено: 23.02.2012 - 22:33 | |
|
Группа: Пользователи Сообщений: 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-кривых - впечатлило! пробую разобраться в нем, если что задам еще вопросы Вам:) |
| Форум: Медицинская статистика · Просмотр сообщения: #13011 · Ответов: 48 · Просмотров: 114142 |
| Отправлено: 23.02.2012 - 22:29 | |
|
Группа: Пользователи Сообщений: 27 Регистрация: 5.02.2012 Пользователь №: 23464
|
Цитата(DrgLena @ 23.02.2012 - 11:24) У вас была проблема, вы не знали формулу, хотя в документации и к одной и к другой программе, которой вы пользуетесь, она есть, как и в любом учебнике по анализу выживаемости. Вот только вы не поняли, как использовать полученные коэффициенты, у вас, кстати, только один коэффициент значим, в том примере, который вы привели, второй не влияет на функцию выживания. Если я не права, покажите, на том же примере, как пользоваться полученными коэффициентами, используя значения предикторов для реального больного. Если вы в этом разобрались, то ознакомьте форум. Потом можно будет двигаться дальше. Да Вы правы, вопрос задавал. Да и согласен второй предиктор не значим, привел данные просто для примера, чтобы посмотреть как будет выглядеть уравнение. Также Вы правы, я смотрел Халяфяна, смотрел интернет, формула действительно есть, но как применить к уравнению(формуле) коэффициенты, чтобы оно в результате давало прогнозируемые дни жизни так и не понял. До сих пор.. |
| Форум: Медицинская статистика · Просмотр сообщения: #13010 · Ответов: 48 · Просмотров: 114142 |
| Отправлено: 23.02.2012 - 08:14 | |
|
Группа: Пользователи Сообщений: 27 Регистрация: 5.02.2012 Пользователь №: 23464
|
Никто не знает как сравнить? |
| Форум: Медицинская статистика · Просмотр сообщения: #13001 · Ответов: 48 · Просмотров: 114142 |
 Открытая тема (есть новые ответы) Открытая тема (нет новых ответов) Горячая тема (есть новые ответы) Горячая тема (нет новых ответов) |
 Опрос (есть новые голоса) Опрос (нет новых голосов) Закрытая тема Тема перемещена |