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

[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