Помощь - Поиск - Пользователи - Календарь
Полная версия этой страницы: Определение чувствительности и специфичности комбинации маркеров
Форум врачей-аспирантов > Разделы форума > Медицинская статистика
Страницы: 1, 2
Zaycho
Здравствуйте!

Есть задача определить диагностическую ценность ряда маркеров (допустим, А, В и С) и их комбинаций в диагностике некоторого заболевания.
Я расчитал все параметры (чувствительность, специфичность, PPV, NPV, LR) и построил в SPSS ROC-кривые отдельно для А, В и С.
Как теперь сделать то же самое для комбинаций этих маркеров, допустим, А+В и В+С? Где копать?

Заранее огромное спасибо!
p2004r
Цитата(Zaycho @ 2.08.2012 - 14:20) *
Здравствуйте!

Есть задача определить диагностическую ценность ряда маркеров (допустим, А, В и С) и их комбинаций в диагностике некоторого заболевания.
Я расчитал все параметры (чувствительность, специфичность, PPV, NPV, LR) и построил в SPSS ROC-кривые отдельно для А, В и С.
Как теперь сделать то же самое для комбинаций этих маркеров, допустим, А+В и В+С? Где копать?

Заранее огромное спасибо!


можно рандомфорест применить и вклад маркеров по его результатам определить. можно логистическую регрессию для полной модели построить и для комбинаций. модели сравнивать через сравнение ROC.
Zaycho
Цитата(p2004r @ 2.08.2012 - 15:24) *
можно рандомфорест применить и вклад маркеров по его результатам определить. можно логистическую регрессию для полной модели построить и для комбинаций. модели сравнивать через сравнение ROC.


Big thx за ответ, но я не такой великий спец в статистике - если не затруднит, можно чуть подробнее? Как построить логистическую регрессию для полной модели (как сгруппировать исходные данные?) В идеале - где это все в SPSS и/или Статистике? Ну а если с примером........ smile.gif))
p2004r
Цитата(Zaycho @ 2.08.2012 - 14:37) *
Big thx за ответ, но я не такой великий спец в статистике - если не затруднит, можно чуть подробнее? Как построить логистическую регрессию для полной модели (как сгруппировать исходные данные?) В идеале - где это все в SPSS и/или Статистике? Ну а если с примером........ smile.gif))


могу в R показать smile.gif, но для полного примера нужны данные smile.gif)

как понимаю есть эксперимент или наблюдение в котором каждому случаю наличия-отсутствия заболевания соответствует строка наличия-отсутствия признаков? если случаи наблюдения-эксперимента не представляют сочетания наличия-отсутствия признаков то естественно ничего хорошего скорее всего не получится.

на такую таблицу (назовем её table) легко натравить randomForest

Код
rf <- randomForest(Болезнь ~ ., data=table,
                             importance=TRUE,
                             proximity=TRUE)


После этого смотрим вклад признаков

Код
varImpPlot(rf)


Zaycho
А можно я все-таки покажу исходные данные?
А, В и С - непрерывные величины
Выжил/умер - принимает только 2 значения smile.gif)
TheThing
Цитата(Zaycho @ 2.08.2012 - 14:20) *
Здравствуйте!

Есть задача определить диагностическую ценность ряда маркеров (допустим, А, В и С) и их комбинаций в диагностике некоторого заболевания.
Я расчитал все параметры (чувствительность, специфичность, PPV, NPV, LR) и построил в SPSS ROC-кривые отдельно для А, В и С.
Как теперь сделать то же самое для комбинаций этих маркеров, допустим, А+В и В+С? Где копать?

Заранее огромное спасибо!


Если Вы лучше знакомы с СПСС, постройте модель с помощью логистической регрессии. Заходите в Analyze-Regression-Binary. Факторы А,Б,С перемещаете в Covariates (independent variables), наличие/отсутсвие заболевания в Dependent variable. Закладка Categorical Вам не нужна, поскольку, как я понял, предикторы все у вас количественные. В закладке Save выбираете сохранить вероятности группы (group probabilities). В опциях 95% ДИ, тест Хосмера-Лемешова и все остальное по-желанию. Нажимаете ОК, строится модель, которая включает в себя набор всех предикторов А,Б,С, рассчитанные коэффициенты, р-значения, ДИ и т.д. ROC-кривая строится с помощью меню ROC-curve, где в Test-variable Вы закидываете рассчитанные групповые вероятности, а в State - указываете код, которым у Вас закодировано состояние больного человека (1 например).

Построенная логистическая модель может предсказывать аддитивный (суммационный) эффект всех факторов А,Б,С на риск развития заболевания. То есть, все факторы представлены главными эффектами (main effects) - не учитывается возможность взаимодействия фактора А например с фактором Б. Если Вас такое устраивает, можно на этом остановиться. А если Вам интересно посмотреть взаимодействует ли фактор А с фактором Б (присутствует ли синергетический эффект например), нужно построить модель, которая учитывала бы взаимодействие факторов (для этого используется кнопочка a*b в логистической регрессии). Интерпретировать результаты будет намного сложнее, могу посоветовать хорошую книгу на эту тему от Sage publishing - Interaction effects (terms) in logistic regression.

Если нравятся случайные леса, установите SPSS Clementine или используйте код, который привел p2004r smile.gif
p2004r
Цитата(Zaycho @ 2.08.2012 - 21:55) *
А можно я все-таки покажу исходные данные?
А, В и С - непрерывные величины
Выжил/умер - принимает только 2 значения smile.gif)


таблицу чуток пришлось отредактировать от пояснений

Код
rf<-randomForest(data.na[,2:4],factor(data.na$выжил.умер), proximity=TRUE)

> rf

Call:
randomForest(x = data.na[, 2:4], y = factor(data.na$выжил.умер),      proximity = TRUE)
               Type of random forest: classification
                     Number of trees: 500
No. of variables tried at each split: 1

        OOB estimate of  error rate: 31.43%
Confusion matrix:
   0  1 class.error
0 22 11   0.3333333
1 11 26   0.2972973


> round(importance(rf), 2)
  MeanDecreaseGini
A            10.93
B            12.09
C            10.97

> MDSplot(rf, factor(data.na$выжил.умер), k=2, palette=c(1, 2), pch=data.na$выжил.умер)


собственно все случаи делятся явно на две группы, и имеются два фактора "выживания". первый фактор действует в обоих группах. второй фактор собственно определяет деление на две группы в одной из которых смертность намного ниже чем во второй группе (и обусловлена действием первого фактора).
Zaycho
Огромное спасибо за разъяснения! А можете теперь ткнуть пальцем ламеру, где здесь групповые вероятности, которые можно подставить в ROC-анализ?

p2004r
Цитата(Zaycho @ 3.08.2012 - 14:27) *
Огромное спасибо за разъяснения! А можете теперь ткнуть пальцем ламеру, где здесь групповые вероятности, которые можно подставить в ROC-анализ?


для первого датасета

Код
> data.frame(votes=rf$votes, alive=data.na$выжил.умер, n=data.na$N )
      votes.0    votes.1 alive  n
1  0.76243094 0.23756906     0  1
2  0.80203046 0.19796954     1  2
3  0.72067039 0.27932961     0  3
4  0.17486339 0.82513661     0  4
5  0.74054054 0.25945946     0  5
6  0.82352941 0.17647059     1  6
7  0.37055838 0.62944162     1  7
8  0.80213904 0.19786096     1  8
9  0.27472527 0.72527473     1  9
10 0.70322581 0.29677419     0 10
11 0.60215054 0.39784946     0 11
12 0.39512195 0.60487805     1 12
13 0.83815029 0.16184971     0 13
14 0.89893617 0.10106383     0 14
15 0.40101523 0.59898477     0 15
16 0.79081633 0.20918367     0 16
17 0.68333333 0.31666667     0 17
18 0.67875648 0.32124352     0 18
19 0.75126904 0.24873096     0 19
20 0.91847826 0.08152174     1 20
21 0.61931818 0.38068182     1 21
22 0.35955056 0.64044944     1 22
23 0.95480226 0.04519774     0 23
24 0.61271676 0.38728324     0 24
25 0.84574468 0.15425532     0 25
26 0.91256831 0.08743169     0 26
27 0.43930636 0.56069364     1 27
28 0.22857143 0.77142857     1 28
29 0.15346535 0.84653465     1 29
31 0.12626263 0.87373737     0 31
33 0.29444444 0.70555556     0 33
34 0.28571429 0.71428571     1 34
35 0.07692308 0.92307692     1 35
36 0.04216867 0.95783133     1 36
37 0.59259259 0.40740741     1 37
38 0.63725490 0.36274510     1 38
39 0.38596491 0.61403509     1 39
40 0.41340782 0.58659218     0 40
41 0.08074534 0.91925466     1 41
42 0.39779006 0.60220994     0 42
43 0.31428571 0.68571429     1 43
44 0.16666667 0.83333333     1 44
45 0.38636364 0.61363636     0 45
46 0.24043716 0.75956284     0 46
47 0.66883117 0.33116883     0 47
48 0.69540230 0.30459770     0 48
49 0.60679612 0.39320388     0 49
50 0.45180723 0.54819277     1 50
51 0.08152174 0.91847826     1 51
52 0.09696970 0.90303030     1 52
53 0.20329670 0.79670330     0 53
54 0.30898876 0.69101124     1 54
55 0.63783784 0.36216216     1 55
56 0.58285714 0.41714286     0 56
57 0.09826590 0.90173410     1 57
58 0.08045977 0.91954023     1 58
59 0.07608696 0.92391304     1 59
60 0.30508475 0.69491525     0 60
61 0.40000000 0.60000000     1 61
62 0.28342246 0.71657754     1 62
63 0.63068182 0.36931818     1 63
64 0.58252427 0.41747573     1 64
65 0.87730061 0.12269939     1 65
66 0.08982036 0.91017964     1 66
67 0.36898396 0.63101604     1 67
68 0.72881356 0.27118644     0 68
69 0.71022727 0.28977273     0 69
70 0.74033149 0.25966851     0 70
71 0.37158470 0.62841530     1 71
72 0.47849462 0.52150538     0 72
Zaycho
Глубокоуважаемый p2004r, правильно ли я понял, что первый столбец цифр с большим количеством знаков после запятой - это для А+В, а второй столбец - для В+С?
Если да, то можно то же самое для А+С?
СПАСИБО!
p2004r
Цитата(Zaycho @ 3.08.2012 - 16:05) *
Глубокоуважаемый p2004r, правильно ли я понял, что первый столбец цифр с большим количеством знаков после запятой - это для А+В, а второй столбец - для В+С?
Если да, то можно то же самое для А+С?
СПАСИБО!


нет, не так.

первые два столбца это вероятность для 0, и вероятность для 1 исхода (они зеркальны). нелинейная модель включает в себя все три исходных фактора.

Ваша группа неоднородна. Надо смотреть что из исходных факторов на что влияет. Иными словами --- какой вклад в два фактора про которых я написал в предыдущем посте (группообразующий и повышающий вероятность неблагоприятного исхода).
Zaycho
Я вконец запутался. А,В и С - это параметры, не влияющие на наступление исхода. А и В - это результаты биохимических тестов, С - оценка тяжести состояния по интегральной шкале. Все три параметра выше при неблагоприятном исходе. Определенный уровень каждого из них с известной чувствительностью и специфичностью является предиктором неблагоприятного исхода. Нас не устраивает их недостаточная прогностическая ценность по отдельности. Как расчитать эти показатели для комбинаций - я так и не понял...........
p2004r
Цитата(Zaycho @ 3.08.2012 - 16:38) *
Я вконец запутался. А,В и С - это параметры, не влияющие на наступление исхода. А и В - это результаты биохимических тестов, С - оценка тяжести состояния по интегральной шкале. Все три параметра выше при неблагоприятном исходе. Определенный уровень каждого из них с известной чувствительностью и специфичностью является предиктором неблагоприятного исхода. Нас не устраивает их недостаточная прогностическая ценность по отдельности. Как расчитать эти показатели для комбинаций - я так и не понял...........


B имеет несколько больший вес чем A и C. Но сама картина существенно нелинейна. Имеют место две группы --- 1я с более высоким уровнем смертности и низкой прогнозируемостью исхода, и 2я с низким уровнем смертности и визуально более предсказуемым исходом. За разбиение на эти две группы отвечает некая комбинация ABC (её надо проанализировать). В каждой из этих групп действует некий дополнительный фактор определяющий смертность, причем действует количественно разно. Что входит в этот фактор тоже нужно проанализировать.

Сначала надо построить решающее правило к какой группе относится конкретный случай, а потом построить для каждой из групп правило исхода.
TheThing
Согласно моим расчетам лишь фактор С (тяжесть состояния больного) статистически значимо ассоциирован с риском смерти (повышенный риск). Первый столбец цифр - групповые вероятности факторов А+Б+С, второй столбец - вероятности А+Б, третий столбец - А+С, четвертый - Б+С.

0,34015 0,40888 0,34989 0,39730
0,49263 0,50523 0,38566 0,55469
0,42836 0,26548 0,39794 0,62781
0,59121 0,36225 0,60615 0,70132
0,25785 0,31131 0,25945 0,36786
0,43051 0,43519 0,44927 0,47932
0,57229 0,30174 0,58240 0,72457
0,25236 0,38984 0,23751 0,31511
0,48796 0,40442 0,50717 0,56495
0,28626 0,39493 0,29287 0,34319
0,49658 0,37706 0,51123 0,59569
0,37246 0,31495 0,38091 0,50680
0,12237 0,25227 0,11330 0,21316
0,12748 0,31118 0,10785 0,19567
0,53328 0,46286 0,54719 0,57564
0,52075 0,45196 0,49164 0,58775
0,49832 0,42483 0,47445 0,58010
0,36305 0,40782 0,37485 0,42414
0,45776 0,43172 0,47040 0,51505
0,36475 0,42263 0,34469 0,43039
0,57694 0,45255 0,59446 0,62775
0,76129 0,46953 0,75743 0,81078
0,14571 0,39503 0,14859 0,16921
0,13010 0,47251 0,11741 0,13051
0,18435 0,41960 0,15968 0,21998
0,11134 0,39310 0,10980 0,12896
0,86965 0,82693 0,66838 0,83274
0,72274 0,61052 0,74657 0,68379
0,88778 0,82411 0,75763 0,84437
0,86492
0,79245 0,63566 0,81702 0,75099
0,74476
0,35380 0,73615 0,40006 0,19948
0,83108 0,69577 0,83671 0,77659
0,73978 0,64728 0,77173 0,67595
0,67489 0,61222 0,70833 0,62232
0,79319 0,67501 0,82141 0,72640
0,69415 0,70943 0,72958 0,57129
0,70840 0,70052 0,74604 0,59504
0,73655 0,79736 0,77758 0,54065
0,58925 0,63310 0,62808 0,50647
0,55043 0,71285 0,58524 0,40279
0,69639 0,67085 0,71908 0,61115
0,88188 0,67862 0,89790 0,84678
0,46330 0,67591 0,50338 0,33916
0,49648 0,53330 0,52366 0,48054
0,74687 0,58398 0,77340 0,72615
0,47468 0,57482 0,49834 0,42908
0,54642 0,51868 0,56019 0,55171
0,34062 0,60430 0,36285 0,27085
0,73011 0,55340 0,75691 0,72478
0,76456 0,50137 0,78571 0,78991
0,60758 0,46409 0,63096 0,64996
0,66616 0,59398 0,69391 0,62750
0,76595 0,57299 0,74402 0,77482
0,53226 0,56143 0,52938 0,51339
0,65210 0,57569 0,68407 0,62190
0,55687 0,54630 0,58823 0,53632
0,82779 0,56820 0,84837 0,82598
0,63646 0,52353 0,63253 0,65455
0,25158 0,48530 0,26149 0,24793
0,76566 0,53693 0,78297 0,77582
0,66544 0,55563 0,69248 0,65225
0,35473 0,50851 0,36883 0,34657
0,47536 0,54117 0,50240 0,45183
0,68043 0,58309 0,71167 0,64929
0,32747 0,47152 0,34773 0,33673
0,12440 0,56551 0,13399 0,09099
0,74684 0,58395 0,77340 0,72612
0,47667 0,58016 0,50408 0,42592
0,32795 0,59106 0,35756 0,26315
0,56675 0,64000 0,45425 0,54231

Коэффициент конкордации около 65%, что в принципе совпадает с классификационной способностью модели, полученной с помощью случайного леса (100 - OOB error = 100 - 31,43 = 68,57%).
p2004r
Цитата(TheThing @ 3.08.2012 - 20:07) *
Согласно моим расчетам лишь фактор С (тяжесть состояния больного) статистически значимо ассоциирован с риском смерти (повышенный риск). Первый столбец цифр - групповые вероятности факторов А+Б+С, второй столбец - вероятности А+Б, третий столбец - А+С, четвертый - Б+С.



Код
plot(cmdscale(dist(rf$proximity), k=2),
       col=c("red","green")[factor(data.na$выжил.умер)],
       cex=c(seq(from=1,
                       to=1.8,
                       along.with=levels(factor(data.na$A))))[factor(data.na$A)])

plot(cmdscale(dist(rf$proximity), k=2),
       col=c("red","green")[factor(data.na$выжил.умер)],
       cex=c(seq(from=1,
                       to=1.8,
                       along.with=levels(factor(data.na$B))))[factor(data.na$B)])

plot(cmdscale(dist(rf$proximity), k=2),
       col=c("red","green")[factor(data.na$выжил.умер)],
       cex=c(seq(from=1,
                       to=1.8,
                       along.with=levels(factor(data.na$C))))[factor(data.na$C)])


Согласно трем графикам (по A, B, и C размер кругов) в пространстве mds полученном из proximity рандомфореста. Если и можно сомневаться в связи какого то фактора так это C smile.gif.

Ну а разделяет группы больных B. Это очевидно из графика. Его маленькое значение во второй группе наступает ступенчато и ведет к низкому риску летального исхода. И этот фактор по мнению рандом фореста самый весомый.

Два остальных фактора --- маленькие значения A в обоих группах у выздоровевших. С ведет себя аналогично A но менее выражена разница у противоположных исходов.
TheThing
значит логистическая регрессия и рэндом форест дают почему-то разные результаты. Нужно разбираться..
p2004r
Цитата(TheThing @ 3.08.2012 - 21:22) *
значит логистическая регрессия и рэндом форест дают почему-то разные результаты. Нужно разбираться..


например рандом форесту практически безразлично что B надо логарифмировать.

Код
> pairs(cbind(data.na[,c(2,4)],
              log(data.na[,3])),
        col=c("red","green")[factor(data.na$выжил.умер)])

> pairs(data.na[,2:4],
        col=c("red","green")[factor(data.na$выжил.умер)])


DrgLena
А если не лес строить, а одно дерево, то при одном ветвлении, на двух конечных узлах при условии С<=16,64 классификация произойдет с точностью 65,7%
Переменная С разделяет пациентов реальной выборки, а вот В не разделяет и вообще не имеет различий в группах выживших и нет ни по каким парным критериям. И не удивительно, такие разбросы данных (скорее всего ошибки ввода) например, у умерших больных есть значение 237 и 0,14. Вряд ли В может называться маркером. Но в сгенерированных выборках значения могут повторяться или вообще не попасть.
ПО приведенным исходным данным точки разделения, и соответственно чувствительность и специфичность с моими расчетами не сходятся, но площадь ROC сходится, это очень странно. Автор должен верифицировать данные, иначе никакие многомерные методы на таких данных применять нельзя.
Автор поста проделал часть работы и задал конкретный вопрос, как объединить результаты трех тестов. Это не задача классификации. Есть понятие положительного и отрицательного теста (найдены значения для разделения на положительный и отрицательный результат теста) и нужно узнать вероятность исхода при различных сочетаниях положительных или отрицательных тестов.
Если не придираться к конкретным данным, то решить эту задачу можно следующим образом. Создать три бинарных переменных по результатам ROC, для этого использовать оптимальные точки разделения, фактор риска имеется, код 1 и нет ? 0. Эти переменные использовать в логистической регрессии. По полученным коэффициентам посчитать вероятности наступления интересующего события. При этом возможны 8 комбинаций наличия и отсутствия трех факторов риска, от отсутствия всех трех (000) до присутствия всех трех (111). Каждая из 8 групп риска будет иметь свою расчетную вероятность наступления события. При этом экспоненциальные коэффициенты уравнения регрессии для каждого предиктора будут иметь ясную интерпретацию и показывать независимый вклад каждого фактора.
Zaycho
Да, я и не предполагал, что задача окажется не такой тривиальной, как мне показалось на первый взгляд. Небольшой комментарий - "С" - это оценка тяжести состояния по шкале APACHE II, которая является общепризнанной для определения риска летальности уже более 10 лет. Ошибок ввода я не нашел - "B" действительно имеет большой диапазон значений.
В любом случае, есть над чем подумать.
p2004r
Цитата(DrgLena @ 3.08.2012 - 23:56) *
А если не лес строить, а одно дерево, то при одном ветвлении, на двух конечных узлах при условии С<=16,64 классификация произойдет с точностью 65,7%
Переменная С разделяет пациентов реальной выборки, а вот В не разделяет и вообще не имеет различий в группах выживших и нет ни по каким парным критериям. И не удивительно, такие разбросы данных (скорее всего ошибки ввода) например, у умерших больных есть значение 237 и 0,14. Вряд ли В может называться маркером. Но в сгенерированных выборках значения могут повторяться или вообще не попасть.
ПО приведенным исходным данным точки разделения, и соответственно чувствительность и специфичность с моими расчетами не сходятся, но площадь ROC сходится, это очень странно. Автор должен верифицировать данные, иначе никакие многомерные методы на таких данных применять нельзя.


1 по мимо того что B надо логарифмировать, имеются корреляции между предикторами

Код
> cor(cbind(data.na[,c(2,4)], log(data.na[,3])))                          A         C log(data.na[, 3])
A                 1.0000000 0.2813438         0.1289035
C                 0.2813438 1.0000000         0.2625643
log(data.na[, 3]) 0.1289035 0.2625643         1.0000000

> cor.test(data.na[,4], data.na[,2])

    Pearson's product-moment correlation

data:  data.na[, 4] and data.na[, 2]
t = 2.4177, df = 68, p-value = 0.01831
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.04965218 0.48430137
sample estimates:
      cor
0.2813438

> cor.test(data.na[,4], log(data.na[,3]))

    Pearson's product-moment correlation

data:  data.na[, 4] and log(data.na[, 3])
t = 2.2439, df = 68, p-value = 0.0281
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
0.02940446 0.46862610
sample estimates:
      cor
0.2625643

> cor.test(data.na[,2], log(data.na[,3]))

    Pearson's product-moment correlation

data:  data.na[, 2] and log(data.na[, 3])
t = 1.0719, df = 68, p-value = 0.2876
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
-0.1093837  0.3531800
sample estimates:
      cor
0.1289035



Группа выделена и имеет отличный состав

Код
> table(factor(cmdscale(dist(rf1$proximity), k=2)[,2]>0.2), factor(data.na$выжил.умер))
      
         0  1
  FALSE 24 21
  TRUE   9 16


во 2й группе исходы 9 к 16, в первой 24 к 21

вклад в разделение групп оказывает не один единственный предиктор, (B и С вместе)

Кружки первая группа, треугольники вторая. (На сечении С|log(B) вообще летальные исходы сгруппировались в центре и вверху. Может отработать эту нелинейность не может при разделении групп рандомфорест. Нужен какой нибудь svm или явная трансформация данных.)

Код
> pairs(cbind(data.na[,c(2,4)],
              log(data.na[,3])),
        pch=c(1,2)[factor(cmdscale(dist(rf1$proximity), k=2)[,2]>0.2)],
        col=c("red","green")[factor(data.na$выжил.умер)])


второй вариант наоборот цветом группы, значками исходы.
p2004r
собственно да, svm получше выделяет

Код
> model <- svm( выжил.умер ~ A+B+C, data = data.na, type= "C-classification")
> table(model$decision.values<=0,data.na$выжил.умер)
      
         0  1
  FALSE 25 12
  TRUE   8 25
> table(model$decision.values<=-0.1,data.na$выжил.умер)
      
         0  1
  FALSE 27 13
  TRUE   6 24

## логарифмирование имеет значение!!
> model <- svm( выжил.умер ~ A+log(B)+C, data = data.na, type= "C-classification")
> table(model$decision.values<=-0.1,data.na$выжил.умер)
      
         0  1
  FALSE 27  9
  TRUE   6 28


вот собственно он их сумел обвести

Код
> pairs(cbind(data.na[,c(2,4)],
              log(data.na[,3])),
        pch=c(1,2)[factor(model$decision.values<=-0.1)],
        col=c("red","green")[factor(data.na$выжил.умер)])


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

PS пробабилити

Код
> cbind(model$decision.values,data.na$выжил.умер)
           0/1  
1   0.98331246 0
2   0.61530576 1
3   1.00012653 0
4  -0.21684733 0
5   0.99955666 0
6   0.16078293 1
7  -1.00012726 1
8   1.14854667 1
9  -0.99985686 1
10  0.99966273 0
11  0.37811021 0
12 -0.42755238 1
13  1.00021827 0
14  1.00007588 0
15  0.76925049 0
16  0.88688175 0
17  1.01823453 0
18  0.86095764 0
19  1.01908459 0
20  1.18428922 1
21  0.45588013 1
22 -0.73663584 1
23  0.89897670 0
24  1.00026140 0
25  1.00020960 0
26  1.12620025 0
27 -0.71744634 1
28 -0.93837583 1
29 -0.99987350 1
31 -1.20454616 0
33  0.62211341 0
34 -1.37866569 1
35 -1.07486185 1
36 -1.09031422 1
37 -1.06053528 1
38 -0.50515338 1
39 -0.75995557 1
40  0.16392607 0
41 -1.00022084 1
42  0.34876086 0
43 -0.64127520 1
44 -1.15716373 1
45  0.11486255 0
46 -0.13282404 0
47 -1.04569130 0
48  0.37514348 0
49  0.56057748 0
50  0.49439228 1
51 -1.11937289 1
52 -0.89372931 1
53 -0.49529302 0
54 -0.65597539 1
55 -0.92934884 1
56  0.62306627 0
57 -1.15534228 1
58 -1.02889926 1
59 -1.00000111 1
60  0.11401070 0
61  0.68242740 1
62 -1.00003176 1
63 -0.65191720 1
64  0.74228345 1
65 -0.04860989 1
66 -1.13318704 1
67 -0.73920762 1
68  0.73537165 0
69 -1.04658105 0
70  0.16253365 0
71 -0.44523330 1
72  0.41346047 0
p2004r
линеаризирующее преобразование (!после отрезания всех явно тяжелых!):

1) берем точку среднее log(B) -- среднее C (или лучше среднее по шкалам mds) для летальных исходов

2) строим вектора с началом в этой точке и концом в точках наблюдений

3) новое пространство --- угол вектора vs его длинна

в таком пространстве должны работать линейные методы

Код
> data.leit<- data.na[data.na$C<25,]
> nrow(data.leit)
[1] 62
> model1 <- svm( выжил.умер ~ A+log(B)+C, data = data.leit, type= "C-classification")
> table(model1$decision.values<=0,data.leit$выжил.умер)
      
         0  1
  FALSE 19  5
  TRUE   6 32


вот по группе без явно безнадежных (критерий C>25) svm вырезает выживших до 32.

вот так решение графически выглядит

Код
> plot(cmdscale(dist(cbind(data.leit[,c(2,4)],
                                 log(data.leit[,3]))),
                      k=2),
             col=c("red","green")[factor(data.leit$выжил.умер)],
             pch=c(1,2)[factor(model1$decision.values<=0)])


Трансформация

Код
> tapply(data.leit[,3], data.leit$выжил.умер, mean)
       0        1
12.67121 20.19569
> tapply(data.leit[,4], data.leit$выжил.умер, mean)
       0        1
13.48000 11.40541
> tapply(log(data.leit[,3]), data.leit$выжил.умер, mean)
       0        1
1.775605 1.189165
> ((data.leit$C-12.67121)^2+(log(data.leit$B)-1.775605)^2)^0.5
[1]  5.3311377  5.1087745  2.1709325  6.0049577  6.3416986  2.5971127
[7]  7.6488885  9.4283268  2.9613371  7.3290071  2.0874944  2.6590571
[13]  0.7077666  1.7968509  1.6400452  4.3481172  1.3288040  5.5428207
[19]  2.6856595  9.7595860  3.7072333  4.6724333  4.3954970  7.6958694
[25]  7.8144577  5.1120372  3.4003874  6.7942236  0.7232215  2.2810260
[31]  0.6200533  3.0480755  5.3432953  1.8203928 12.6712147  7.3712229
[37]  2.4614863  6.8202322  4.3332743  0.6811193 10.3318129  7.2930514
[43]  9.8974529  4.0868944  2.6899526  6.9645308  2.7174036  3.7413277
[49]  2.2007804 12.0539703  3.0041289 11.3308822  8.6715569  3.7608580
[55]  7.3350855  3.3930430  4.4240841  7.5608856  6.8246717  4.3312637
[61] 10.4496927  6.0334436
> l<-((data.leit$C-12.67121)^2+(log(data.leit$B)-1.775605)^2)^0.5
> atan2(data.leit$C-12.67121, log(data.leit$B)-1.775605)
[1]  1.60047511  1.01092706 -0.87854661 -1.90576287  1.63461177  2.02937132
[7] -2.08195084  1.42536019 -2.91294791  1.57849395 -2.21332718  2.61831456
[13] -1.24799125  0.18401803  0.20184394  1.66511774  1.56620199  1.29199482
[19] -1.67457571 -1.43611893  0.08880546 -1.54791341 -0.65318030 -1.65087067
[25] -1.37902877 -1.98913076 -2.23799124 -1.76137710 -1.95235967 -2.31937471
[31]  2.58268440  2.69050581  1.49709553 -1.16313148 -1.57165444  1.67814707
[37]  1.90064623 -1.78022440  1.52529843  0.50375801  1.54660568 -1.98674087
[43] -1.78502201 -2.02573395 -1.68891279 -1.27953913  1.02940171 -2.34643526
[49]  2.99163425 -1.82347545 -1.09554895  1.55157912 -1.56185158 -1.78957673
[55]  1.52936211  1.76571555 -2.16280092  1.81921224 -1.78326322  1.60459531
[61]  1.72306160  1.08266063
>  a<-atan2(data.leit$C-12.67121, log(data.leit$B)-1.775605)
> plot(a,
       l,
       col=c("red","green")[factor(data.leit$выжил.умер)],
       pch=c(1,2)[factor(model1$decision.values<=0)])


и в пространстве трансформированных шкал mds

Код
> data.mds<-cmdscale(dist(cbind(data.leit[,c(2,4)],log(data.leit[,3]))))
> tapply(data.mds[,1], data.leit$выжил.умер, mean)
         0          1
-1.2605935  0.8517523
> tapply(data.mds[,2], data.leit$выжил.умер, mean)
         0          1
-0.2579902  0.1743177
> l.mds<-((data.mds[,1]+1.2605935)^2+(data.mds[,2]+0.2579902)^2)^0.5
> a.mds<-atan2(data.leit$C-12.67121, log(data.leit$B)-1.775605)
> plot(a.mds,
       l.mds,
       col=c("red","green")[factor(data.leit$выжил.умер)],
       pch=c(1,2)[factor(model1$decision.values<=0)])


похожую трансформацию, судя по внешнему сходству, рандом_форест и пытался сделать.

PS забыл отнормировать данные прежде чем считать длину вектора и угол frown.gif

Код
> data.mds.scale<-scale(data.mds, center=c(-1.2605935, -0.2579902))
> l.mds.scale<-((data.mds.scale[,1])^2+(data.mds.scale[,2])^2)^0.5
> a.mds.scale<-atan2(data.mds.scale[,1], data.mds.scale[,2])
> plot(a.mds.scale,l.mds.scale, col=c("red","green")[factor(data.leit$выжил.умер)], , pch=c(1,2)[factor(model1$decision.values<=0)])


на картинке я зону svm приблизительно обвел, горизонтальная ось замкнута (это радианы угла вектора)

поскольку все читать утомительно я резюмирую:

решатель состоящий из правила -- все кто тяжелее 25 считать безнадежными и svm-решателя для остальных случаев дает
11 ошибок на 70 ~ 0.16
Zaycho
Коллеги, а если вернуться к посту #14, то как определить оптимальные точки разделения (cut-off) каждого из маркеров в комбинации?
p2004r
Цитата(Zaycho @ 5.08.2012 - 11:23) *
Коллеги, а если вернуться к посту #14, то как определить оптимальные точки разделения (cut-off) каждого из маркеров в комбинации?


точка оптимальная считается для некого индекса по комбинации маркеров, рассчитанная отдельно по маркеру эта точка в функцию развернется. (фактически зона svm решения это и есть "по маркерам" построенная "точка разделения")

то есть --- сначала считать индекс, потом для индекса точку отсечения. если это имелось в виду то в ROC пакет обычно входит соответствующая функция. например в pROC это coords

Код
If ?x="best"?, the ?best.method? argument controls how the optimal
     threshold is determined.

     ?youden? Youden's J statistic (Youden, 1950) is employed. The
          optimal cut-off is the threshold that maximizes the distance
          to the identity (diagonal) line. Can be shortened to ?y?.

          The optimality criterion is:

                        max(sensitivities + specificities)              
          
     ?closest.topleft? The optimal threshold is the point closest to
          the top-left part of the plot with perfect sensitivity or
          specificity. Can be shortened to ?c? or ?t?.

          The optimality criterion is:

                min((1 - sensitivities)^2 + (1- specificities)^2)      
          
     In addition, weights can be supplied if false positive and false
     negative predictions are not equivalent: a numeric vector of
     length 2 to the ?best.weights? argument. The indices define

       1. the cost of of a false negative classification

       2. the prevalence, or the proportion of cases in the total
          population (n.cases/(n.controls+n.cases)).
:
The optimality criteria are modified as proposed by Perkins and
     Schisterman:

     ?youden?

                   max(sensitivities + r \times specificities)          
          
     ?closest.topleft?

            min((1 - sensitivities)^2 + r \times (1- specificities)^2)  
          
     with

                 r = (1 - prevalence) / (cost * prevalence)            
    
     By default, prevalence is 0.5 and cost is 1 so that no weight is
     applied in effect.

     Note that several thresholds might be equally optimal.
TheThing
Решил сделать рэндом форест, а также построить простое CART, чтобы понять почему результаты расходятся с логистической регрессией и с результатами p2004r. Как выяснилось, все со всем сходится smile.gif

Сначала рендом форест с 3 исходными предикторами:

randomForest(formula = Group ~ .,
data = crs$dataset[, c(crs$input, crs$target)],
ntree = 500, mtry = 1, importance = TRUE, replace = FALSE, na.action = na.roughfix)

Type of random forest: classification
Number of trees: 500
No. of variables tried at each split: 1

OOB estimate of error rate: 27.78%
Confusion matrix:
alive dead class.error
alive 30 9 0.2307692
dead 11 22 0.3333333

Analysis of the Area Under the Curve (AUC)
==========================================

Call:
roc.default(response = crs$rf$y, predictor = crs$rf$votes)

Data: crs$rf$votes in 39 controls (crs$rf$y alive) > 33 cases (crs$rf$y dead).
Area under the curve: 0.7584

95% CI: 0.6462-0.8706 (DeLong)

Variable Importance
===================

alive dead MeanDecreaseAccuracy
C 2.87 3.59 2.67
B 1.90 3.03 2.14
A 1.89 1.65 1.59

Ошибка OOB составляет 27,78%, что в принципе говорит, что модель неплохая (AUC = 75%)/ На первом месте по важности предикторов выступает фактор С, затем В и лишь в конце болтается А.

Решил убрать фактор А из модели и посмотреть насколько изменится классификационная способность. Вот что получилось:

OOB estimate of error rate: 31.94%
Confusion matrix:
alive dead class.error
alive 28 11 0.2820513
dead 12 21 0.3636364

Analysis of the Area Under the Curve (AUC)
==========================================

Call:
roc.default(response = crs$rf$y, predictor = crs$rf$votes)

Data: crs$rf$votes in 39 controls (crs$rf$y alive) > 33 cases (crs$rf$y dead).
Area under the curve: 0.7599

95% CI: 0.6499-0.8699 (DeLong)

Variable Importance
===================

alive dead MeanDecreaseAccuracy
C 3.51 4.77 3.26
B 1.69 3.40 2.24

Видим, что OOB увеличилась на пару процентов, а AUC остался прежним, фактор С опять на 1 месте по важности. Поэтому я бы предиктор А вообще-бы не включал в модель, а рассматривал его как ненужный "шум".

Затем построил простое CART и получил результат, который полностью совпал с методом опорных векторов (svm) в подсчетах p2004r, логистической регрессией и random forest.

Хорошо то, что хорошо заканчивается laugh.gif

Нажмите для просмотра прикрепленного файла

p2004r
и все таки, после применения правила С>25, остается смесь из 62 случаев, которую svm пилит с вероятностью ошибки (5+6)/(19+5+6+32) ~ 0.1774194

у выживших разброс в двумерном распределении C - log(B) больше чем у летальных исходов.
TheThing
Кстати, вот хорошая статья, где проводился сравнительный анализ support vector machines, logistic regression, random forest, CART, neural netss, boosted trees, bagged trees, naive Bayes и др. Согласно результатам, на первом месте - бустид-деревья, второе - рэндом форест, третье - опорные вектора.

Нажмите для просмотра прикрепленного файла
p2004r
Цитата(TheThing @ 5.08.2012 - 16:40) *
Согласно результатам, на первом месте - бустид-деревья, второе - рэндом форест, третье - опорные вектора.


Спасибо!

для интерполяции (ну и вырезания "хитрых областей", например "вложенная спираль" или вообще структурные данные типа языковых конструкций) svm хорош, правда надо вид ядра подбирать и за регуляризацией следить что бы экстраполяция не страдала.
Диана
Раз зашла речь про ROC-анализ в программе SPSS подскажите как определить Точки разделения для факторов риска чтоб выявить риск развития осложнений в зависимости от них. прилагаю результаты логистической регрессии. Есть ли программа SPSS в бесплатном доступе или лучше воспользоваться MEDCALC. у мнея еще вопрос. для фактора риска ИБС очень широкий ДИ, как мне это интерпретировать
TheThing
Цитата(Диана @ 16.08.2012 - 09:25) *
Раз зашла речь про ROC-анализ в программе SPSS подскажите как определить Точки разделения для факторов риска чтоб выявить риск развития осложнений в зависимости от них. прилагаю результаты логистической регрессии. Есть ли программа SPSS в бесплатном доступе или лучше воспользоваться MEDCALC. у мнея еще вопрос. для фактора риска ИБС очень широкий ДИ, как мне это интерпретировать


В бесплатном доступе SPSS нет, есть в пиратском.

Логистическая регрессия как линейный параметрический классификатор очень плохо справляется с задачами, где наблюдается большое количество предикторов (больше 10, у Вас их больше 20), поскольку наблюдается "curse of dimentionality" - проклятие размерности, особенно при таких небольших выборках (относительно к количеству независимых переменных). Все это выливается в biased (искаженные) коэффициенты регрессии, большие стандартный ошибки и огромные ДИ.

Расчет необходимой выборки для логит-регрессии - задача оч. непростая. Peduzzi et al. (1996) предложили упрощенный вариант:
N = 10 k / p

k - количество предикторов
р - наименьшая пропорция негативных/позитивных случаев

У Вас к равняется 23, допустим р равняется 0.3 (то есть 30% людей в выборке больных, а 70% здоровых), следовательно:

N = 10 * 23 / 0,3 = 766 человек. Это минимальное количество людей, если Вы хотите анализировать 23 предиктора. Если это значение получилось бы меньше 100, его следует увеличить до 100, как советует Long (1997).

Выходов несколько:

1) увеличивать выборку
2) уменьшать количество предикторов
3) применять соответствующие методы (а лучше их совокупность), которые лишены проблем логит-регрессии.

P.S. p-value в логит-регрессии - это не стат. значимость отношения шансов.
Диана
Цитата(TheThing @ 16.08.2012 - 16:38) *
В бесплатном доступе SPSS нет, есть в пиратском.

Логистическая регрессия как линейный параметрический классификатор очень плохо справляется с задачами, где наблюдается большое количество предикторов (больше 10, у Вас их больше 20), поскольку наблюдается "curse of dimentionality" - проклятие размерности, особенно при таких небольших выборках (относительно к количеству независимых переменных). Все это выливается в biased (искаженные) коэффициенты регрессии, большие стандартный ошибки и огромные ДИ.

Расчет необходимой выборки для логит-регрессии - задача оч. непростая. Peduzzi et al. (1996) предложили упрощенный вариант:
N = 10 k / p

k - количество предикторов
р - наименьшая пропорция негативных/позитивных случаев

У Вас к равняется 23, допустим р равняется 0.3 (то есть 30% людей в выборке больных, а 70% здоровых), следовательно:

N = 10 * 23 / 0,3 = 766 человек. Это минимальное количество людей, если Вы хотите анализировать 23 предиктора. Если это значение получилось бы меньше 100, его следует увеличить до 100, как советует Long (1997).

Выходов несколько:

1) увеличивать выборку
2) уменьшать количество предикторов
3) применять соответствующие методы (а лучше их совокупность), которые лишены проблем логит-регрессии.

P.S. p-value в логит-регрессии - это не стат. значимость отношения шансов.

Спасибо, уже для меня начало что то проясняться. но, к сожалению, не смогу увеличить объем выборки, тогда мне остается уменьшить количество предикторов, только не знаю как- либо отсеять логически или может нужно их как то проанализировать стат.методом. Расскажите , пожалуйста, подробней про 3 вывод
p2004r
Цитата(Диана @ 16.08.2012 - 20:10) *
Спасибо, уже для меня начало что то проясняться. но, к сожалению, не смогу увеличить объем выборки, тогда мне остается уменьшить количество предикторов, только не знаю как- либо отсеять логически или может нужно их как то проанализировать стат.методом. Расскажите , пожалуйста, подробней про 3 вывод


так а где сами данные?
TheThing
Цитата(Диана @ 16.08.2012 - 20:10) *
Спасибо, уже для меня начало что то проясняться. но, к сожалению, не смогу увеличить объем выборки, тогда мне остается уменьшить количество предикторов, только не знаю как- либо отсеять логически или может нужно их как то проанализировать стат.методом. Расскажите , пожалуйста, подробней про 3 вывод


Выводов у меня не было, были возможные выХоды из ситуации smile.gif
Данные действительно бы пригодились..
Диана
КФР-контролируемый фактор риска, НФР-неконтролируемый ФР, ЗАДАЧА- выявить наиболее вероятные факторы риска для прогнозирования развития кардиальных осложнений в послеоперационном периоде
p2004r
Цитата(Диана @ 17.08.2012 - 09:02) *
КФР-контролируемый фактор риска, НФР-неконтролируемый ФР, ЗАДАЧА- выявить наиболее вероятные факторы риска для прогнозирования развития кардиальных осложнений в послеоперационном периоде


много пропущенных данных... будет интересно smile.gif

у меня после импорта данных (там иногда попадаются точки в ячейках, нельзя ли сохранить таблицу с простым заголовком в старом эксель формате? а то я не уверен в корректности импорта) получились вот такие переменные:

Код
> names(read.csv2("mp3gl.csv"))
  [1] "ФИО.шифр"                        "пол"                            
  [3] "возраст"                         "ИМТ..кг.м2"                    
  [5] "ИМТ"                             "S.тела.по.Дюб"                  
  [7] "к.во.дней.до.опер"               "к.во.дней.после.опер"          
  [9] "время.уст..я.д.за.анев"          "кардио.жалобы"                  
[11] "длит.кардиожалоб"                "хирург.жалобы"                  
[13] "группы.по.лечению"               "п.о.кард.осложн"                
[15] "Кард.осл.я"                      "летальность"                    
[17] "стресс.ЭхоКГ"                    "к.во.зон.гипо.и.акинезии"      
[19] "X"                               "X.1"                            
[21] "ГБ"                              "X.2"                            
[23] "ГБ.ГЛЖ"                          "ИБС..стенокардия"              
[25] "постинф.кардиоскл"               "НК..фк"                        
[27] "ДЛП"                             "размер.АБА"                    
[29] "синдромность"                    "прех.ишем.послеопер"            
[31] "СД"                              "доступ"                        
[33] "курение"                         "диаст..Дисфункция00"            
[35] "САД1"                            "ДАД1"                          
[37] "ЧСС1"                            "САД5"                          
[39] "ДАД5"                            "ЧСС5"                          
[41] "САД.д.о"                         "ДАД.д.о"                        
[43] "ЧСС.д.о"                         "САД.п.о"                        
[45] "ДАД.п.о"                         "ЧСС.п.о"                        
[47] "бета.блокеры"                    "ББдоза"                        
[49] "ББкол.во.суток.назн"             "антаг.кальция"                  
[51] "АКдоза.суточная"                 "АКкол.во.сут.назн"              
[53] "антиагреганты"                   "инг.АПФ"                        
[55] "ИАПФсут.доза"                    "ИАПФкол.во.сут.назн"            
[57] "нитраты"                         "дигоксин"                      
[59] "диуретики"                       "кордарон"                      
[61] "гепарин"                         "трад.тер"                      
[63] "операция1"                       "доступ1"                        
[65] "время.операции1"                 "Время.опер.11"                  
[67] "Время.опер.21"                   "Время.опер.31"                  
[69] "время.переж.аорты1"              "Вр.переж.аорты11"              
[71] "Вр.переж.аорты21"                "кровопотеря1"                  
[73] "кровопотеря11"                   "наруш.ритма00"                  
[75] "рубц.измен00"                    "ЧСС00"                          
[77] "ритм00"                          "наруш.пров00"                  
[79] "PQ"                              "QRS"                            
[81] "QT"                              "ST"                            
[83] "з.Т"                             "Соколов.Лайон"                  
[85] "корн..Индекс"                    "корн..Произв.е"                
[87] "минЧСС"                          "максЧСС"                        
[89] "средЧСС"                         "ЖЭС"                            
[91] "ЖЭС..к.во"                       "ЖТ"                            
[93] "НЖЭС..к.во"                      "НЖТ..к.во"                      
[95] "депрессия.ST"                    "АВ.блокады"                    
[97] "пароксизмы.ФП"                   "ГЛЖ00"                          
[99] "X.3"                             "зоны.гипо.и.акин00"            
[101] "КДР..см00"                       "КСР..см00"                      
[103] "КДО..мл00"                       "индекс.КДО..00"                
[105] "КСО00"                           "инд.КСО00"                      
[107] "ФВ00"                            "УО00"                          
[109] "ТМЖП00"                          "ТЗСТ.ЛЖ00"                      
[111] "X.4"                             "ММ.ЛЖ.формула.Devereux"        
[113] "ИММ.ЛЖ"                          "ПЗРВТ.00"                      
[115] "КДР.ПЖ00"                        "ЛП00"                          
[117] "d.аорты00"                       "ГДС.СДЛА00"                    
[119] "диаст..Дисфункция00.1"           "кальциноз.аорты00"              
[121] "наруш.ритма11"                   "ЧСС11"                          
[123] "PQ11"                            "QRS11"                          
[125] "QT11"                            "депрессия.ST11"                
[127] "з.Т11"                           "длит.ть.прех.ишемии11"          
[129] "очаг.измен11"                    "мин.ЧСС11"                      
[131] "макс.ЧСС11"                      "средняя.ЧСС11"                  
[133] "ЖЭС.по.Лауну11"                  "ЖЭС..к.во11"                    
[135] "парные.ЖЭС11"                    "ЖТ11"                          
[137] "НЖЭС.к.во11"                     "НЖТ.к.во11"                    
[139] "депрессия.SТ11"                  "нар.е.пров11"                  
[141] "пароксизмы.ФП11"                 "п.о.гипокинезия11"              
[143] "КДР11"                           "КСР11"                          
[145] "КДО11"                           "инд.КДО11"                      
[147] "КСО11"                           "инд.КСО11"                      
[149] "ФВ11"                            "УО11"                          
[151] "ТМЖП11"                          "ТЗСТ.ЛЖ11"                      
[153] "ММ.ЛЖ.формула.Devereux.1"        "ИММ.ЛЖ.1"                      
[155] "ПЗРВТ.11"                        "КДР.ПЖ11"                      
[157] "ЛП11"                            "d.аорты11"                      
[159] "ГДС.СДЛА11"                      "бета.блокеры11"                
[161] "доза11ББ"                        "ББкол.во.суток.назначения.пре11"
[163] "антагонисты.кальция11"           "АГдоза.суточная11"              
[165] "АГкол.во.суток.назн11"           "антиагреганты11"                
[167] "инг.АПФ11"                       "ИАПФсуточная.доза11"            
[169] "ИАПФкол.во.суток.назн11"         "нитраты11"                      
[171] "дигоксин11"                      "диуретики11"                    
[173] "кордарон11"                      "гепарин11"                      
[175] "трад.терапия11"


какие номера являются осложнениями? или какие являются предикторами?

вот например я предполагаю что осложнение это:

Код
> (read.csv2("mp3gl.csv")$"постинф.кардиоскл")
[1] 1 1 1 1 0 1 1 1 0 0 0 1 1 0 1 0 1 0 0 1 1 0 1 0 0 1 0 0 0 1 0 0 0 0 0 0 0 0
[39] 1 0 1 1 0 0 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 1 1 0 0 0 0 0 0 0 1
[77] 0 0 0 0 0 0 0 1 0 1 1 0 0 1 1 0 0 0 0


а предикторы это со 2 по 18ю

Код
> names(data[,2:18])
[1] "пол"                      "возраст"                
[3] "ИМТ..кг.м2"               "ИМТ"                    
[5] "S.тела.по.Дюб"            "к.во.дней.до.опер"      
[7] "к.во.дней.после.опер"     "время.уст..я.д.за.анев"  
[9] "кардио.жалобы"            "длит.кардиожалоб"        
[11] "хирург.жалобы"            "группы.по.лечению"      
[13] "п.о.кард.осложн"          "Кард.осл.я"              
[15] "летальность"              "стресс.ЭхоКГ"            
[17] "к.во.зон.гипо.и.акинезии"


поскольку среди предикторов много пропусков перед нами два пути:

1) исключить случаи в которых есть пропуски

Код
> data.na<-na.omit(data[,c(2:18,25)])

> table(data.na$постинф.кардиоскл)

0  1
8 18


> rf<-randomForest(data.na[,-18],factor(data.na$"постинф.кардиоскл"), proximity=TRUE)
> rf

Call:
randomForest(x = data.na[, -18], y = factor(data.na$постинф.кардиоскл),      proximity = TRUE)
               Type of random forest: classification
                     Number of trees: 500
No. of variables tried at each split: 4

        OOB estimate of  error rate: 38.46%
Confusion matrix:
  0  1 class.error
0 0  8   1.0000000
1 2 16   0.1111111


> MDSplot(rf, factor(data.na$"постинф.кардиоскл"), k=3)
> varImpPlot(rf)


как видим ничего в таком варианте модель не разделяет ничего.

2) заполним пропущенные случаи прибегнув к множественной импутации
Код
> data.2.18.impute<-rfImpute(data[,2:18],factor(data$"постинф.кардиоскл"), iter=5)
ntree      OOB      1      2
  300:  26.32% 10.94% 58.06%
ntree      OOB      1      2
  300:  27.37% 14.06% 54.84%
ntree      OOB      1      2
  300:  27.37% 15.62% 51.61%
ntree      OOB      1      2
  300:  26.32% 14.06% 51.61%
ntree      OOB      1      2
  300:  29.47% 17.19% 54.84%
> rf.impute<-randomForest(data.2.18.impute[,-1],data.2.18.impute[,1], proximity=TRUE)
> rf.impute

Call:
randomForest(x = data.2.18.impute[, -1], y = data.2.18.impute[,      1], proximity = TRUE)
               Type of random forest: classification
                     Number of trees: 500
No. of variables tried at each split: 4

        OOB estimate of  error rate: 29.47%
Confusion matrix:
   0  1 class.error
0 54 10   0.1562500
1 18 13   0.5806452

> MDSplot(rf.impute, data.2.18.impute[,1], k=3)
> varImpPlot(rf.impute)
> plot(importance(rf.impute),varUsed(rf.impute))
> text(importance(rf.impute),varUsed(rf.impute), labels= names(data.2.18.impute[,-1]))


получается имеют влияние предикторы: "ИМТ..кг.м2", "к.во.зон.гипо.и.акинезии", "длит.кардиожалоб", "возраст", "к.во.дней.после.опер".

меньше влияет "к.во.дней.до.опер".

и совсем слабо (но еще влияет) "S.тела.по.Дюб", "кардио.жалобы", "стресс.ЭхоКГ", "хирург.жалобы", "группы.по.лечению".

Естественно эту информацию можно накопить относительно различных осложнений и построить общую картину.

Так какие показатели являются осложнениями? или какие являются предикторами к осложнениям? Скорее всего "группы.по.лечению" я зря в предикторы включал?
Диана
я немножко сократила таблицу, убрала малоинтересующие данные. желтым цветом выделила неконтролируемые факторы риска, синим- контролируемые. красный цвет- осложнений. цель- узнать как влияет кардиопротекторная терапия- то есть группы лечения-серый цвет.[ в группе 1 и 2 под влиянием терапии происходит урежение ЧсС(ЧСС д/о) и меньше случаев тахикардии (чсс больше 100) в послеоперационном периоде( считается, что тахикардия индуцирует ишемию миокарда и, в свою очередь, инфаркт миокарда)] на развитие кардиальных осложнений
p2004r
Цитата(Диана @ 17.08.2012 - 18:05) *
я немножко сократила таблицу, убрала малоинтересующие данные. желтым цветом выделила неконтролируемые факторы риска, синим- контролируемые. красный цвет- осложнений. цель- узнать как влияет кардиопротекторная терапия- то есть группы лечения-серый цвет.[ в группе 1 и 2 под влиянием терапии происходит урежение ЧсС(ЧСС д/о) и меньше случаев тахикардии (чсс больше 100) в послеоперационном периоде( считается, что тахикардия индуцирует ишемию миокарда и, в свою очередь, инфаркт миокарда)] на развитие кардиальных осложнений


я попробую сказать конкретнее smile.gif

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

давайте обсуждать список переменных датасета который я уже один раз почистил.
мне нужно знать какие номера переменных являются предикторами, а какие "факторами риска".

это всего два списка smile.gif
Диана
Факторы риска:
2
3
5
17
20
22 - 34
64
68
70
73
колонки 44-46 закодировала 0/1(44А, 45А,46А)
предикторы
44А
45А
46А
118
извините, что не сохранила файл в старом формате smile.gif
Заранее спасибо
p2004r
Цитата(Диана @ 18.08.2012 - 08:58) *
Факторы риска:
2
3
5
17
20
22 - 34
64
68
70
73
колонки 44-46 закодировала 0/1(44А, 45А,46А)
предикторы
44А
45А
46А
118
извините, что не сохранила файл в старом формате smile.gif
Заранее спасибо


вот это --- "колонки 44-46 закодировала 0/1(44А, 45А,46А)" было сделано зря. В таком виде никакого разделения практически нет. надо вернуть прежнюю шкалу и тогда (как мне кажется) удастся построить с помощью случайного леса регрессию.

118 дал разделение

Код
> data2.varset1.impute.X118<-rfImpute(data2.varset1,factor(data2$X118), iter=5)
ntree      OOB      1      2
  300:  16.84% 10.91% 25.00%
ntree      OOB      1      2
  300:  22.11% 14.55% 32.50%
ntree      OOB      1      2
  300:  21.05% 12.73% 32.50%
ntree      OOB      1      2
  300:  20.00% 12.73% 30.00%
ntree      OOB      1      2
  300:  20.00% 12.73% 30.00%
> rf.data2.varset1.impute.X118.best<-tuneRF(data2.varset1.impute.X118[,-1],factor(data2$X118), proximity=TRUE, doBest=TRUE)
mtry = 4  OOB error = 24.21%
Searching left ...
mtry = 2     OOB error = 23.16%
0.04347826 0.05
Searching right ...
mtry = 8     OOB error = 22.11%
0.08695652 0.05
mtry = 16     OOB error = 20%
0.0952381 0.05
mtry = 22     OOB error = 23.16%
-0.1578947 0.05
> rf.data2.varset1.impute.X118.best

Call:
randomForest(x = x, y = y, mtry = res[which.min(res[, 2]), 1],      proximity = TRUE)
               Type of random forest: classification
                     Number of trees: 500
No. of variables tried at each split: 16

        OOB estimate of  error rate: 20%
Confusion matrix:
   0  1 class.error
0 46  9   0.1636364
1 10 30   0.2500000
> MDSplot(rf.data2.varset1.impute.X118.best, factor(data2$X118), k=3)
> varImpPlot(rf.data2.varset1.impute.X118.best)
> plot(importance(rf.data2.varset1.impute.X118.best),varUsed(rf.data2.varset1.impute.X118.best))
> text(importance(rf.data2.varset1.impute.X118.best),varUsed(rf.data2.varset1.impute.X118.best), labels= names(data2.varset1.impute.X118[,-1]), pos=4)


имеет значение только

[30,] "X30" "прех.ишем.послеопер"
[3,] "X3" "возраст"

что вполне подтверждается графиком

> mosaicplot(table(data2$X30, cut(data2$X3, hist(data2$X3, plot=FALSE)$breaks), data2$X118))

Оценим пригодность модели для диагностики с помощью ROC.

Код
> plot(roc(data2$X118, rf.data2.varset1.impute.X118.best$votes[,2]), print.auc=TRUE)

Call:
roc.default(response = data2$X118, predictor = rf.data2.varset1.impute.X118.best$votes[,     2])

Data: rf.data2.varset1.impute.X118.best$votes[, 2] in 55 controls (data2$X118 0) < 40 cases (data2$X118 1).
Area under the curve: 0.8111


Судя по области решения, которую надо вырезать в пространстве сформированном mds по восстановленным случайным лесом расстояниям между экспериментальными случаями, лучше с разделением справится svm. Случайный лес похоже не дотягивается до случаев со средним возрастом из основной группы.
Диана
спасибо за помощь. колонки 44-46 не имеют влияние? посоветуйте мне, пожалуйста, литературу, чтоб я поняла все это, поскольку я все еще плохо дружу со статистикой. по моим расчетам вид лечения(группы по лечению) до операции не влияет на осложнения, а мне желательно доказать обратное. я сокращала все осложнения до фатальные и нефатальные-тоже нет желаемого результата, только отмечено снижение летальности в 1 группе по сравнению с 3й.
p2004r
Цитата(Диана @ 18.08.2012 - 20:57) *
спасибо за помощь. колонки 44-46 не имеют влияние? посоветуйте мне, пожалуйста, литературу, чтоб я поняла все это, поскольку я все еще плохо дружу со статистикой. по моим расчетам вид лечения(группы по лечению) до операции не влияет на осложнения, а мне желательно доказать обратное. я сокращала все осложнения до фатальные и нефатальные-тоже нет желаемого результата, только отмечено снижение летальности в 1 группе по сравнению с 3й.


0. Прочитать любую книгу про методы многомерного анализа данных.

1. Анализ показывает только то, что колонки 44А, 45А, 46А (в виде 0-1) невозможно разделить с помощью 2 3 5 17 20 22 - 34 64 68 70 73. Колонки 44, 45, 46 в набор тех с помощью которых пытаемся разделять я не включал, поскольку это не следовало из Вашего описания. Увы, но за Вас я не могу придумать какой набор переменных какую переменную должен попытаться оценить.

2. (еще раз) То что 44А, 45А, 46А перекодировали в 0-1 плохо, это потеря информации. Случайным лесом можно строить не только модель классифицирующую, но и регрессионную. возможно что тогда зависимость станет очевидной.

3. Группы лечения это переменная 13? или какая?
Диана
smile.gif колонки 44А 45А 46А были перекодированные колонки 44 45 46 соответственно, группа лечения- колонка 13
Диана
Подскажите кто знает: возможно ли удалить выложенные данные
Диана
ничего не получается, от проводимого лечения нет достоверного влияния на снижение осложнений
p2004r
Цитата(Диана @ 26.08.2012 - 08:59) *
ничего не получается, от проводимого лечения нет достоверного влияния на снижение осложнений


Немного был занят.

1) "группы.по.лечению" это что за переменная? это она описывает проводимое лечение?

2) в конце первого файла данных были всякие "дигоксины", они имеют отношение к проводимому лечению?

Диана
Цитата(p2004r @ 27.08.2012 - 16:12) *
Немного был занят.

1) "группы.по.лечению" это что за переменная? это она описывает проводимое лечение?

2) в конце первого файла данных были всякие "дигоксины", они имеют отношение к проводимому лечению?

здравствуйте.
лечение - колонки 47-62, а колонка 13 группы по лечению- делит всех пациентов на 3 группы
задачи исследования: 1. влияет ли снижение ЧСС до операции на снижение частоты кардиальных осложнений
2. влияет ли тип операции на снижение частоты кардиальных осложнений
3. влияние факторов риска
p2004r
Цитата(Диана @ 27.08.2012 - 16:58) *
здравствуйте.
лечение - колонки 47-62, а колонка 13 группы по лечению- делит всех пациентов на 3 группы
задачи исследования: 1. влияет ли снижение ЧСС до операции на снижение частоты кардиальных осложнений
2. влияет ли тип операции на снижение частоты кардиальных осложнений
3. влияние факторов риска


Не совсем понимаю frown.gif

Хорошо, 13 колонка делит пациентов на 3 группы. Лечение в группах было разное? Можно эту колонку использовать как признак назначенного лечения? Иными словами какой принцип разбиения пациентов на эти три группы?

1) влияет ли снижение ЧСС до операции на снижение частоты кардиальных осложнений

в какой колонке находится "снижение ЧСС до операции", и в какой колонке соответственно "частота кардиальных осложнений"?

2) влияет ли тип операции на снижение частоты кардиальных осложнений

колонка "тип операции" какой номер?

3) влияние факторов риска

на что влияние? на какую колонку? ну и номера колонок "факторов риска" укажите.
Диана

Хорошо, 13 колонка делит пациентов на 3 группы. Лечение в группах было разное? Можно эту колонку использовать как признак назначенного лечения? Иными словами какой принцип разбиения пациентов на эти три группы?

можно, принцип разбиения-по проводимой терапии(1-бетаблокаторы, 2 - антагонисты кальция, 3 - ингибиторы АПФ
Диана



1) влияет ли снижение ЧСС до операции на снижение частоты кардиальных осложнений

в какой колонке находится "снижение ЧСС до операции", - колонка 43, ( может надо подсчитать на сколько снизилась от первоначального при поступлении в стационар), и в какой колонке соответственно "частота кардиальных осложнений"?-колонка 14
Диана


2) "влияет ли тип операции на снижение частоты кардиальных осложнений

колонка "тип операции" какой номер? " - колонка- 64 (два вида доступа), и обьем операции - 63,65, 69, 72

3) влияние факторов риска

"на что влияние? на какую колонку?" - 14 ( п/о кард. осложнения), ну и номера колонок "факторов риска" укажите колонки- 2, 3, 4, 17, 20 - 34, 44, 45, 46.
[/quote]
Для просмотра полной версии этой страницы, пожалуйста, пройдите по ссылке.
Форум IP.Board © 2001-2025 IPS, Inc.