Логистическая регрессия? Опции

[mamalita]

Добрый день! Прошу помощи в анализе данных. Мы имеем 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 значений). Может быть их можно как-то объединить и логически видоизменить? Я уже просто голову сломала. Очень нужен свежий взгляд. Спасибо

[mamalita]

приложенные данные

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

 кол_во_размер_рецидив.rar ( 5,19 килобайт ) Кол-во скачиваний: 528

 

[p2004r]

приложенные данные

часть 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

[p2004r]

часть 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))
>

Теперь для способа с гистограммой надо выбрать интервалы группировки

Эскизы прикрепленных изображений

 

[p2004r]

часть 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 ...

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

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

 result.csv.gz ( 378 байт ) Кол-во скачиваний: 760

 

[p2004r]

Проведем разведочный анализ

лично я предпочитаю смотреть на трансформированные 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. видно что фактор "рецидив до года " явно неравномерно распределен в "пространстве мтс"

Эскизы прикрепленных изображений

 

[p2004r]

Построим полную модель различий

Код

> 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

[p2004r]

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

Код

mosaicplot(table(result[,1],data$рецидив.до.года, result[,5]))

вертикальная ось показывает разбиение по признаку "рецидив до года"

горизонтальная ось показывает разбиение по числу мелких мтс

внутри каждое сочетание разбито по числу крупных мтс

суть модели более чем визуализирована

Эскизы прикрепленных изображений

 

[p2004r]

Поскольку у крупных опухолей всего два уровня (есть-нет) построим графическое представление лучшей модели прямо на графике (число малых опухолей vs рецидив)

сначала склеим таблицу данных с нормальными именами (что то predict.glm() обиделась на имена со спецсимволами)

Код

model.data<-data.frame("i1"=result[,1],"i5"=result[,5],"рецидив.до.года"=data$рецидив.до.года)
model <- glm(рецидив.до.года~ i1+i5, data=model,data, family=binomial)

нанесем на график линии предсказания модели

Код

plot(model.data$i1,
       model.data$рецидив.до.года)

lines(seq(0,5, length.out=40),
       predict(model,data.frame(i1=seq(0,5, length.out=40),
                                   i5=1),
                 type="response"))

lines(seq(0,5, length.out=40),
       predict(model,data.frame(i1=seq(0,5, length.out=40),
                                   i5=0),
                  type="response"),
       col="red")

PS ну и наоборот построим тоже

Код

> plot(model.data$i5,model.data$рецидив.до.года)
> lines(seq(0,1, length.out=40), predict(model,data.frame(i5=seq(0,1, length.out=40), i1=0),type="response"))
> lines(seq(0,1, length.out=40), predict(model,data.frame(i5=seq(0,1, length.out=40), i1=1),type="response"))
> lines(seq(0,1, length.out=40), predict(model,data.frame(i5=seq(0,1, length.out=40), i1=2),type="response"))
> lines(seq(0,1, length.out=40), predict(model,data.frame(i5=seq(0,1, length.out=40), i1=3),type="response"))
> lines(seq(0,1, length.out=40), predict(model,data.frame(i5=seq(0,1, length.out=40), i1=4),type="response"))
> lines(seq(0,1, length.out=40), predict(model,data.frame(i5=seq(0,1, length.out=40), i1=5),type="response"))

тут очевидно чем больше тем хуже

на формулировку http://forum.disser.ru/index.php?showtopic...ost&p=12211 ответ получен?

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

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

Эскизы прикрепленных изображений