Логистическая регрессия? |
Здравствуйте, гость ( Вход | Регистрация )
Логистическая регрессия? |
9.11.2011 - 11:07
Сообщение
#1
|
|
Группа: Пользователи Сообщений: 49 Регистрация: 7.04.2010 Пользователь №: 15366 |
Добрый день! Прошу помощи в анализе данных. Мы имеем 100 человека больных с метастазами в печени, лечили их хирургическим путем и наблюдали их в течение 6 лет и диагностировали у них новые метастазы и рецидивы. Суть работы заключается в том, чтобы доказать, что реже всего на 1 и 2 годах наблюдения новые МТС возникают у пациентов с количеством МТС 2-3 (у нас были варианты количества от 2-6), и размер их должен быть 2-3 см. То есть кривая частоты прогресии имеет форму колокола обращенного вершиной вниз - 1 см - часто возникает прогрессия, больше 3 - тоже. Наиболее оптимальным является размер очага для хирургического лечения 2-3 см. Вопрос как представить эти данные и их анализировать: средний и суммарный размер не учитывают разницы: то ли у больного было 3 очага по 2 см то ли 1 и 6 см что совсем не благоприятно. Если брать каждый метастаз как отдельную переменную то у разных людей будет разное количество переменных (от 2 до 6штук), но этот вариант наиболее приемлем в соответствии с поставленной задачей. Теперь вопрос каким методом воспользоваться, чтобы доказать что идеальным для лечения является количество МТС 2-3 при размере 2-3см. Еще момент : размер МТС имеет мини манимальное округление до 0, 5 разброс от 1 до 6 см (т.е. всего 12 значений). Может быть их можно как-то объединить и логически видоизменить? Я уже просто голову сломала. Очень нужен свежий взгляд. Спасибо
|
|
9.11.2011 - 11:55
Сообщение
#2
|
|
Группа: Пользователи Сообщений: 1091 Регистрация: 26.08.2010 Пользователь №: 22699 |
Может быть их можно как-то объединить и логически видоизменить? заменяем размеры плотностью их распределением plot(density(c(первый,второй, остальные если есть), from=0,to="максимального из всех ста")) density()$y это набор фиксированных точек отражающий плотность распределения размера. Вне зависимости от числа метастазов это фиксированный набор переменных. Ну или гистограмму строим, выбираем интервалы группировки и строим набор переменных соответсвующий интервалам группировки (сколько образований в каждый попало). |
|
10.11.2011 - 12:35
Сообщение
#3
|
|
Группа: Пользователи Сообщений: 49 Регистрация: 7.04.2010 Пользователь №: 15366 |
Познания у меня не очень глубокие, поэтому не совсем Вас поняла. Плотность распределения размеров МТС у каждого больного? Или вообще плотность распределения размеров на всю выборку, тогда как выявить зависимость появления новых МТС.У меня программа Statistica 6, attestat. Спасибо за помощь!
|
|
10.11.2011 - 15:50
Сообщение
#4
|
|
Группа: Пользователи Сообщений: 1114 Регистрация: 10.04.2007 Пользователь №: 4040 |
Добрый день! Прошу помощи в анализе данных. Мы имеем 100 человека больных с метастазами в печени, лечили их хирургическим путем и наблюдали их в течение 6 лет и диагностировали у них новые метастазы и рецидивы. Суть работы заключается в том, чтобы доказать, что реже всего на 1 и 2 годах наблюдения новые МТС возникают у пациентов с количеством МТС 2-3 (у нас были варианты количества от 2-6), и размер их должен быть 2-3 см. То есть кривая частоты прогресии имеет форму колокола обращенного вершиной вниз - 1 см - часто возникает прогрессия, больше 3 - тоже. Наиболее оптимальным является размер очага для хирургического лечения 2-3 см. Вопрос как представить эти данные и их анализировать: средний и суммарный размер не учитывают разницы: то ли у больного было 3 очага по 2 см то ли 1 и 6 см что совсем не благоприятно. Если брать каждый метастаз как отдельную переменную то у разных людей будет разное количество переменных (от 2 до 6штук), но этот вариант наиболее приемлем в соответствии с поставленной задачей. Теперь вопрос каким методом воспользоваться, чтобы доказать что идеальным для лечения является количество МТС 2-3 при размере 2-3см. Еще момент : размер МТС имеет мини манимальное округление до 0, 5 разброс от 1 до 6 см (т.е. всего 12 значений). Может быть их можно как-то объединить и логически видоизменить? Я уже просто голову сломала. Очень нужен свежий взгляд. Спасибо Несколько сумбурно, но, думаю, в процессе работы задача уточнится. 1. Обучение модели. Можно попробовать логистическую регрессию. Для начала взять 2 переменные: количество МТС и размер. Т.о., для каждого пациента будет одна строка из трех чисел: количество, размер, наличие рецидива (0 или 1). Эти данные ввести в программу и посчитать регрессию. Оценить ее качество путем построения ROC кривой и выбора оптимального (важного для исследователя) порога. 2. Распознавание. При вводе новых данных (размер, количество) вычисленный результат регрессии можно будет интерпретировать как вероятность появления рецидива. С учетом порога, выбранного выше, результат может быть огрублен до бинарной переменной (появится симптом или не появится). Можно посчитать и другие параметры (шансы и т.д.). При дальнейшей работе, возможно, модель уточнится путем ввода дополнительных параметров, которые исследователь сочтет важными для результата. Сообщение отредактировал Игорь - 10.11.2011 - 15:54 Ebsignasnan prei wissant Deiws ainat! As gijwans! Sta ast stas arwis!
|
|
10.11.2011 - 20:40
Сообщение
#5
|
|
Группа: Пользователи Сообщений: 49 Регистрация: 7.04.2010 Пользователь №: 15366 |
В том то и проблема - какой размер брать. Когда у одного больного 2 метастаза: 1 см и 5 см, у второго 5 метастазов по 1 см, у третьего 4 метастаза: 1см, 2см, 4см, 3см и т.д., т.е. среднее не будет полноценно отражать проблему, т.к. нужно доказать что максимально эффективно излечиваются очаги размером 2-3 см, а наличие мелких или наоборот больших образований в разы ухудшает прогноз и повышает частоту развития рецидивов. Возможно опять объяснила сумбурно, обязательно напишите, постараюсь исправиться.
|
|
10.11.2011 - 22:49
Сообщение
#6
|
|
Группа: Пользователи Сообщений: 1325 Регистрация: 27.11.2007 Пользователь №: 4573 |
Трудность понятна, есть не просто число, но и размер к каждому числу.
Если вы хотите доказать имеющуюся у вас клиническую гипотезу о лучшем исходе для определенной комбинации признаков, то и создайте те группы сравнения, которые вам необходимы. Есть два исходных параметра - число МТС и их величина. По числу, выделите две группы "0" 2-3, и альтернативная вторая группа "1" >3-х МТС, также и по величине "0" 2-3 см и "1" >3см. Т.о. вы имеете уже 4 группы по сочетанию двух созданных вами групп сравнения 00 (лучший исход), 01,10,11(худщая комбинация). А теперь, посмотрите их различия по выживаемости, длительности без рецидивного или без метастатического течения, начните с кривых Каплан- Майера, если конечно для этого есть данные, т.е. есть время операции и время диагностики нового рецидива. Если различия будут, то можно эти две бинарные переменные использовать в дальнейшем моделировании, где будут и другие предикторы, например, гистологический клеточный тип (в виде бинарной переменной - благоприятный или нет) или возраст больного. Сообщение отредактировал DrgLena - 10.11.2011 - 23:28 |
|
10.11.2011 - 22:56
Сообщение
#7
|
|
Группа: Пользователи Сообщений: 1091 Регистрация: 26.08.2010 Пользователь №: 22699 |
Познания у меня не очень глубокие, поэтому не совсем Вас поняла. Плотность распределения размеров МТС у каждого больного? Или вообще плотность распределения размеров на всю выборку, тогда как выявить зависимость появления новых МТС.У меня программа Statistica 6, attestat. Спасибо за помощь! Ну R наверняка заработает на Вашем компьютере 1 плотность считать у каждого больного, на выходе получается замена "2 и более переменных" на фиксированное число переменных отражающих плотность распределения. 2 интервалы группировки находим в гистограмме всей выборки. Затем применяем к каждому больному. Эффект тот же "2 и более чисел" заменяется числом попаданий в интервалы группировки. Выкладывайте пример данных попробуем вместе написать код. |
|
13.11.2011 - 11:29
Сообщение
#8
|
|
Группа: Пользователи Сообщений: 49 Регистрация: 7.04.2010 Пользователь №: 15366 |
Есть два исходных параметра - число МТС и их величина. По числу, выделите две группы "0" 2-3, и альтернативная вторая группа "1" >3-х МТС, также и по величине "0" 2-3 см и "1" >3см. Т.о. вы имеете уже 4 группы по сочетанию двух созданных вами групп сравнения 00 (лучший исход), 01,10,11(худщая комбинация). А теперь, посмотрите их различия по выживаемости, длительности без рецидивного или без метастатического течения, начните с кривых Каплан- Майера, если конечно для этого есть данные, т.е. есть время операции и время диагностики нового рецидива. Если различия будут, то можно эти две бинарные переменные использовать в дальнейшем моделировании, где будут и другие предикторы, например, гистологический клеточный тип (в виде бинарной переменной - благоприятный или нет) или возраст больного. В целом Вашу мысль я поняла, но как группировать тех у кого 1 мтс - 1 см, второй - 3 см, третий 5см. Пример данных прилагаю. Как загрузить R и какое-нибудь руководство к ней. |
|
13.11.2011 - 11:34
Сообщение
#9
|
|
Группа: Пользователи Сообщений: 49 Регистрация: 7.04.2010 Пользователь №: 15366 |
приложенные данные
Прикрепленные файлы
|
|
13.11.2011 - 11:42
Сообщение
#10
|
|
Группа: Пользователи Сообщений: 1091 Регистрация: 26.08.2010 Пользователь №: 22699 |
Пример данных прилагаю. Как загрузить R и какое-нибудь руководство к ней. руководство по R на русском языке http://m7876.wiki.zoho.com/Introduction-to-R.html сам R берем с крана http://cran.r-project.org/ Для виндовс http://cran.r-project.org/bin/windows/base/ |
|
13.11.2011 - 11:56
Сообщение
#11
|
|
Группа: Пользователи Сообщений: 1091 Регистрация: 26.08.2010 Пользователь №: 22699 |
приложенные данные часть 1 Импорт данных сохраняем файл экселя как csv выбирая разделителем полей ";" ввиду печальной особенности русского экселя представлять десятичную точку как "," директория в которую мы сохранили файл с данными будет нашей рабочей директорией (в линуксе достаточно стоя в ней запустить R в виндовсе возможно надо её отметить из меню оболочки), или надо указать полный путь к имени файла. после это импортируем в эксель файл функцией read.csv2() и сразу смотрим все ли пошло нормально с помощью функции str() Код > str(read.csv2("кол-во-размер-рецидив.csv")) 'data.frame': 93 obs. of 13 variables: $ номер : int 43 91 3 4 10 34 35 39 42 93 ... $ кол.во : int 5 5 3 2 2 4 5 6 6 6 ... $ средний.размер : num 1.6 1.4 2.33 3 3 2.25 2 1.17 1.5 1.25 ... $ мтс1 : num 1 0.5 1 1 1 1 0.5 1 1 0.5 ... $ мтс2 : num 1 0.5 1 5 5 1 1 1 1 1 ... $ мтс3 : num 1 1 5 0 0 2 1.5 1 1 1 ... $ мтс4 : num 1 1 0 0 0 5 2 1 2 1 ... $ мтс5 : num 4 4 0 0 0 0 5 1.5 2 4 ... $ мтс6 : num 0 0 0 0 0 0 0 1.5 2 0 ... $ суммарный.размер: num 8 7 7 6 6 9 10 7 9 7.5 ... $ рецидив : int 1 1 1 1 1 1 1 1 1 1 ... $ срок : int 1 1 1 3 3 3 3 3 3 3 ... $ рецидив.до.года : int 1 1 1 1 1 1 1 1 1 1 ... > как видим всё в порядке и мы можем смело сохранять данные для последующей обработки Код > data <- read.csv2("кол-во-размер-рецидив.csv") > ЗЫ ну а посмотреть на первые несколько строчек что у нас в data можно с помощью функции head() Код > head(data)
номер кол.во средний.размер мтс1 мтс2 мтс3 мтс4 мтс5 мтс6 суммарный.размер 1 43 5 1.60 1.0 1.0 1 1 4 0 8 2 91 5 1.40 0.5 0.5 1 1 4 0 7 3 3 3 2.33 1.0 1.0 5 0 0 0 7 4 4 2 3.00 1.0 5.0 0 0 0 0 6 5 10 2 3.00 1.0 5.0 0 0 0 0 6 6 34 4 2.25 1.0 1.0 2 5 0 0 9 рецидив срок рецидив.до.года 1 1 1 1 2 1 1 1 3 1 1 1 4 1 3 1 5 1 3 1 6 1 3 1 Сообщение отредактировал p2004r - 13.11.2011 - 12:24 |
|
13.11.2011 - 13:16
Сообщение
#12
|
|
Группа: Пользователи Сообщений: 1091 Регистрация: 26.08.2010 Пользователь №: 22699 |
часть 2 Выбираем интервал группировки
Нас интересует как распределены размеры "мтс" Код > names(data) [1] "номер" "кол.во" "средний.размер" "мтс1" [5] "мтс2" "мтс3" "мтс4" "мтс5" [9] "мтс6" "суммарный.размер" "рецидив" "срок" [13] "рецидив.до.года" это переменные с 4 по 9, у них отсутствующие значения закодированы нолем (а надо NA) выкусываем часть таблицы данных, заменяем 0 на NA: Код > # прицеливаемся :) > head(data[,4:9]) мтс1 мтс2 мтс3 мтс4 мтс5 мтс6 1 1.0 1.0 1 1 4 0 2 0.5 0.5 1 1 4 0 3 1.0 1.0 5 0 0 0 4 1.0 5.0 0 0 0 0 5 1.0 5.0 0 0 0 0 6 1.0 1.0 2 5 0 0 > mts<-as.vector(data[,4:9]) > # заменяем 0 на NA > mts[mts==0]<-NA и склеиваем в одну переменную: 1й способ Код c(mts[,1],mts[,2],mts[,3],mts[,4],mts[,5],mts[,6]) 2й способ Код mts.new<-as.matrix(mts) dim(mts.new)<-6*93 mts.new<-as.vector(mts.new) Очистка от NA производится функцией na.omit() Строим плотность распределения размера мтс Код > plot(density(na.omit(mts.new))) > Ну и голую гистограмму для сравнения Код > hist(na.omit(mts.new)) > Теперь для способа с гистограммой надо выбрать интервалы группировки |
|
13.11.2011 - 14:14
Сообщение
#13
|
|
Группа: Пользователи Сообщений: 1091 Регистрация: 26.08.2010 Пользователь №: 22699 |
часть 3 Применяем интервал группировки Допустим что мы выбрали группировку 0.5-1, 1-2, 2-3, 3-4, 4-5 функция hist( кроме картинки ) возвращает много нужных нам параметров: сколько куда попало и позволят задать произвольные интервалы группировки Код > str(hist(na.omit(mts.new), plot=FALSE)) List of 7 $ breaks : num [1:10] 0.5 1 1.5 2 2.5 3 3.5 4 4.5 5 $ counts : int [1:9] 147 40 78 12 36 2 25 9 24 $ intensities: num [1:9] 0.7882 0.2145 0.4182 0.0643 0.193 ... $ density : num [1:9] 0.7882 0.2145 0.4182 0.0643 0.193 ... $ mids : num [1:9] 0.75 1.25 1.75 2.25 2.75 3.25 3.75 4.25 4.75 $ xname :8322456 "na.omit(mts.new)" $ equidist : logi TRUE - attr(*, "class")=8322456 "histogram" > hist(na.omit(mts.new), plot=FALSE)$breaks [1] 0.5 1.0 1.5 2.0 2.5 3.0 3.5 4.0 4.5 5.0 > hist(na.omit(mts.new), plot=FALSE)$counts [1] 147 40 78 12 36 2 25 9 24 > таким образом если мы скажем Код > hist(na.omit(mts.new), plot=FALSE, breaks=c(0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0))$counts [1] 147 40 78 12 36 2 25 9 24 > мы получим разбиение na.omit(mts.new) по границам c(0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0). Поменяв разбиение на нужное нам мы получаем число попаданий в интересующие нас интервалы: Код > hist(na.omit(mts.new), plot=FALSE, breaks=c(0.5, 1.0, 2.0, 3.0, 4.0, 5.0))$counts [1] 147 118 48 27 33 Осталось применить эту функцию подсчета попаданий в выбранные интервалы индивидуально к каждому случаю. Код вот например как последовательно трансформируются мтс данные первого случая > mts[1,] мтс1 мтс2 мтс3 мтс4 мтс5 мтс6 1 1 1 1 1 4 NA > as.matrix(mts[1,]) мтс1 мтс2 мтс3 мтс4 мтс5 мтс6 1 1 1 1 1 4 NA > as.vector(as.matrix(mts[1,])) [1] 1 1 1 1 4 NA > na.omit(as.vector(as.matrix(mts[1,]))) [1] 1 1 1 1 4 attr(,"na.action") [1] 6 attr(,"class") [1] "omit" > ну подсчет группировки для него Код > hist(as.vector(as.matrix(mts[1,])), plot=FALSE, breaks=c(0.5, 1.0, 2.0, 3.0, 4.0, 5.0))$counts [1] 4 0 0 1 0 > теперь применим данную операцию ко всем случаям, это проще потому что функция apply() применяет написанную нами функцию сразу к строке таблицы (или столбцу если написать apply(данные,2,функция)). t() поворачивает таблицу "набок", транспонирует. Код > result<-t(apply(mts,1,function (x) hist(x, plot=FALSE, breaks=c(0.5, 1.0, 2.0, 3.0, 4.0, 5.0))$counts)) > head(result) [,1] [,2] [,3] [,4] [,5] [1,] 4 0 0 1 0 [2,] 4 0 0 1 0 [3,] 2 0 0 0 1 [4,] 1 0 0 0 1 [5,] 1 0 0 0 1 [6,] 2 1 0 0 1 > str(result) int [1:93, 1:5] 4 4 2 1 1 2 2 4 3 4 ... вот мы и получили трансформированные данные...
Прикрепленные файлы
|
|
13.11.2011 - 15:21
Сообщение
#14
|
|
Группа: Пользователи Сообщений: 1091 Регистрация: 26.08.2010 Пользователь №: 22699 |
Проведем разведочный анализ
лично я предпочитаю смотреть на трансформированные PCA данные, фактически на "голую ковариацию" (соответствующие картинки присоединены) Код > result.pca<-prcomp(result, scale.=FALSE, center=FALSE) > plot(result.pca) > biplot(result.pca) > biplot(result.pca,choices=3:2) посмотрим какие у нас есть исходы Код > summary(as.factor(data[,11])) 0 1 3 90 > summary(as.factor(data[,12])) 1 3 6 9 12 15 18 21 24 27 30 36 3 12 18 18 16 7 5 4 4 2 3 1 > summary(as.factor(data[,13])) 0 1 27 66 > names(data[,11:13]) [1] "рецидив" "срок" "рецидив.до.года" > теперь отобразим в пространство PCA исходы Код > plot(result.pca$x[,1], result.pca$x[,2], bg=c("grey50", "white")[as.factor(data$рецидив)], pch=c(21,22)[as.factor(data$рецидив.до.года)], cex=c(1.8,1)[as.factor(data$рецидив.до.года)]) > plot(result.pca$x[,3], result.pca$x[,2], bg=c("grey50", "white")[as.factor(data$рецидив)], pch=c(21,22)[as.factor(data$рецидив.до.года)], cex=c(1.8,1)[as.factor(data$рецидив.до.года)]) 1. очевидно что "рецидив" более менее равномерно расположился в "пространстве размеров мтс" (даже удивительно для всего 3х случаев). весь он попал в число "рецидив до года" 2. видно что фактор "рецидив до года " явно неравномерно распределен в "пространстве мтс" |
|
13.11.2011 - 16:31
Сообщение
#15
|
|
Группа: Пользователи Сообщений: 1091 Регистрация: 26.08.2010 Пользователь №: 22699 |
Построим полную модель различий
Код > summary(glm(data$рецидив.до.года~ result, family=binomial)) Call: glm(formula = data$рецидив.до.года ~ result, family = binomial) Deviance Residuals: Min 1Q Median 3Q Max -3.0450 -0.4758 0.3491 0.6423 1.8174 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -0.7664 1.1000 -0.697 0.485961 result1 1.2934 0.3929 3.292 0.000994 *** result2 -0.4517 0.3396 -1.330 0.183483 result3 0.2313 0.4565 0.507 0.612424 result4 0.1051 0.8383 0.125 0.900231 result5 1.4076 0.8880 1.585 0.112934 --- Signif. codes: 0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 112.05 on 92 degrees of freedom Residual deviance: 79.59 on 87 degrees of freedom AIC: 91.59 Number of Fisher Scoring iterations: 5 > достоверно влияние первого диапазона размеров сравним модели с разным числом параметров Код model.1 <- glm(data$рецидив.до.года~ result[,1], family=binomial) model.15 <- glm(data$рецидив.до.года~ result[,c(1,5)], family=binomial) model.125 <- glm(data$рецидив.до.года~ result[,c(1,2,5)], family=binomial) model.full <- glm(data$рецидив.до.года~ result, family=binomial) > anova(model.1,model.15,model.125,model.full) Analysis of Deviance Table Model 1: data$рецидив.до.года ~ result[, 1] Model 2: data$рецидив.до.года ~ result[, c(1, 5)] Model 3: data$рецидив.до.года ~ result[, c(1, 2, 5)] Model 4: data$рецидив.до.года ~ result Resid. Df Resid. Dev Df Deviance 1 91 90.709 2 90 83.080 1 7.6290 3 89 79.866 1 3.2139 4 87 79.590 2 0.2764 выберем лучшую для прогноза по минимуму AIC Код > AIC(model.1,model.15,model.125,model.full) df AIC model.1 2 94.70903 model.15 3 89.08005 model.125 4 87.86617 model.full 6 91.58975 и добьем кси квадратом Код > anova(model.1,model.15,model.125,model.full,test="Chi") Analysis of Deviance Table Model 1: data$рецидив.до.года ~ result[, 1] Model 2: data$рецидив.до.года ~ result[, c(1, 5)] Model 3: data$рецидив.до.года ~ result[, c(1, 2, 5)] Model 4: data$рецидив.до.года ~ result Resid. Df Resid. Dev Df Deviance P(>|Chi|) 1 91 90.709 2 90 83.080 1 7.6290 0.005744 ** 3 89 79.866 1 3.2139 0.073016 . 4 87 79.590 2 0.2764 0.870918 --- Signif. codes: 0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1 отобразим модель Код > summary(model.15)
Call: glm(formula = data$рецидив.до.года ~ result[, c(1, 5)], family = binomial) Deviance Residuals: Min 1Q Median 3Q Max -2.8694 -0.7329 0.3294 0.7098 1.7005 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -1.1773 0.4975 -2.367 0.017954 * result[, c(1, 5)]1 1.2137 0.3245 3.741 0.000183 *** result[, c(1, 5)]2 1.6365 0.6453 2.536 0.011214 * --- Signif. codes: 0 ?***? 0.001 ?**? 0.01 ?*? 0.05 ?.? 0.1 ? ? 1 (Dispersion parameter for binomial family taken to be 1) Null deviance: 112.05 on 92 degrees of freedom Residual deviance: 83.08 on 90 degrees of freedom AIC: 89.08 Number of Fisher Scoring iterations: 5 Сообщение отредактировал p2004r - 13.11.2011 - 16:39 |
|