Кластерный анализ в R

[womak]

Возникла необходимость сделать кластерный анализ для довольно большой выборки (5200 значений).
Применил следующий набор команд статистического пакета R ?

Код: Выделить всё

> setwd("C:/data/R Project/earthqwake")                             # установил рабочую директорию
> x <- read.table('earth.csv', sep=',', header=TRUE)                # прочитал данные из файла в матрицу Х
> x <- x[-1]                                                       # удалил первый столбец с номером записи
> d <-  dist(x, method = "euclidean", diag = FALSE, upper = FALSE)  # расчет матрицы расстояний 
> hc <- hclust(d)                                                  # построение дерева кластеров 
> plot(hc)                                                         # отрисовал результат 

[/size]
Могу ли я в этом пакете получить каким-нибудь образом матрицу объединения наподобие той, которую выдает STATISTICA? Кстати, STATISTICA с поставленной задачей не справилась: вылетала молча на этапе объединения

Код: Выделить всё

Связь  Объект 1 Объект 2 Объект 3 Объект 4 ... ... ... ... ... ... ... ... Объект N
0.0001 c_1 c_2
0.0001 c_5 c_345 c_515
...........................
0.2324 c_3 c_6 c_453 c_345 c_123 c_543 c_231
...........................
2.6543 c_3 c_6 c_453 c_345 c_123 c_543 c_231 ... ... ... ... ... ... ... ... ...  c_N

где с_i - порядковый номер записи в файле данных i [1, N] 

Или нечто подобное, но в другом формате.
Стоит "обратная" задача: полученные кластеры с определенным уровнем связи "вынести на местность" и заполнить значениями.
Задача является продолжением и расширением темы

earth.csv

Исходные данные

(92.53 КБ) 1219 скачиваний

Кластеры .rar

Результаты расчетов

(17.96 КБ) 797 скачиваний

[Максим Дубинин]

матрицу объединения наподобие той, которую выдает STATISTICA

покажете матрицу?

[gamm]

1) по умолчанию в hclust метод "complete", что не всегда хорошо - "ward" дает более компактные классы, но все зависит от задач. И лучше, если данные были нормированы (scale) перед вычисление матрицы (иначе разный масштаб сильно повлияет), или вообще использовать расстояние Махалонобиса, которое уберет и масштаб, и корреляцию (особенно, если много переменных)

2) Примерно понятно, чего нужно - можно использовать cutree() - но проблема, что она перенумеровывает классы. Единственное решение, которое я нашел (давно это было ...), это получить первичную нумерацию классов, а потом отслеживать перекодировку. Кусок кода приведен ниже (мне была нужна структура последних need.classes.num, "бахрома" не интересна; final.cls содержит историю номеров с сохранением номера, если не было разделения). В коде остались следы еще чего-то, но разберетесь, кусок вроде рабочий - только что проверял. Исходные данные в таблице t, максимальное число классов - в need.classes.num

Код: Выделить всё

  p.dist<-dist(t)
  p.clust<-hclust(p.dist,method="ward",members=NULL)
  p.nrow<-length(t[,1])
  p.nclass<-need.classes.num
  p.classes<-matrix(ncol=p.nclass-1,nrow=p.nrow)
  N<-nrow(t)
  m<-matrix(ncol=p.nclass+1,nrow=N)
  m[,1]<-seq(1,N)
  m[,2:p.nclass]<-cutree(p.clust,k=2:p.nclass)

  recode<-matrix(ncol=p.nclass,nrow=p.nclass)
  code<-recode
  code[2,1:2]<-1:2
  for(p.cur_cls in (p.nclass-1):2) {
    t.tab<-table(m[,p.cur_cls],m[,p.cur_cls+1])
    for(i in 1:p.cur_cls) {
      for(j in 1:(p.cur_cls+1)) {
        if(t.tab[i,j] > 0) {
          recode[p.cur_cls+1,j]<-i
        }
      }
    }
  }

  for(p.cur_cls in 2:(p.nclass-1)) {
    used<-rep(FALSE,p.cur_cls+1)
    for(j in 1:(p.cur_cls+1)) {
      if(used[recode[p.cur_cls+1,j]]) {
        code[p.cur_cls+1,j] <- p.cur_cls+1
        j0<-p.cur_cls+1
        j1<-code[p.cur_cls,recode[p.cur_cls+1,j]]
        cat(sprintf("Split %d => %d + %d\n",j1,j1,j0))
      } else {
        code[p.cur_cls+1,j] <- code[p.cur_cls,recode[p.cur_cls+1,j]]
      }
      used[recode[p.cur_cls+1,j]]<-TRUE
    }  
  }

  m2<-m
  m2[,1]<-1
  final.cls<-data.frame(num=1:nrow(t))
  for(i_cls in 2:p.nclass) {
    m2[,i_cls]<-code[i_cls,m2[,i_cls]]
  }
  m2.copy<-m2

  final.cls<-data.frame(num=(1:nrow(t)))
  for(i_cls in 2:p.nclass) {
    p.cls2<-m2[,i_cls]
    tmp.names<-names(final.cls)
    final.cls<-data.frame(final.cls,foo=p.cls2)
    names(final.cls)<-c(tmp.names,sprintf("Class%d",i_cls))
  }

[womak]

матрицу объединения наподобие той, которую выдает STATISTICA

покажете матрицу?

Для 4000 значений это выглядит так (и правильно называется список объединения или каталог объединения)

maxrix.rar

Матрица объединения

(36.69 КБ) 736 скачиваний

Первая строка - подписи объектов (лишняя в нашем деле вещь)
Первый столбец - уровень связи, каждому уровню соответствуют свои записи. По этим записям можно восстановить данные, принадлежащие к тому или иному кластеру.

gamm, спасибо за программу. Для моего уровня сложновато, но попытаюсь разобраться...

На сколько я знаю, в MatLab-e можно открыть интересующую функцию и посмотреть содержимое, а в R можно сделать также?

[gamm]

gamm, спасибо за программу. Для моего уровня сложновато, но попытаюсь разобраться...

На сколько я знаю, в MatLab-e можно открыть интересующую функцию и посмотреть содержимое, а в R можно сделать также?

1) программа работает много лет, выдрал я кусок вроде аккуратно и проверил. Выдает столбцы номеров классов, сохраняя нумерацию, и говорит (на экране) кто на что разделился, можно запоминать - будет аналог Статистики

Код: Выделить всё

Split 2 => 2 + 3
Split 1 => 1 + 4
Split 1 => 1 + 5

2) можно взять исходники с сайта R, либо использовать debug("имя функции"), тогда при вызове функции попадете в отладчик (точки останова, к сожалению, ставить нельзя). Если набрать имя функции без скобок, то для старых функций выдается текст. Учтите, что язык замысловатый, для обычного программиста непривычный.

[womak]

gamm, запустил вашу программу, теперь пытаюсь осмыслить результат...
При need.classes.num = 10, я получил такой результат

Спойлер

Split 1 => 1 + 3
Split 2 => 2 + 4
Split 1 => 1 + 5
Split 3 => 3 + 6
Split 2 => 2 + 7
Split 4 => 4 + 8
Split 3 => 3 + 9
Split 1 => 1 + 10

При need.classes.num = 50, я получил такой результат

Спойлер

Split 1 => 1 + 3
Split 2 => 2 + 4
Split 1 => 1 + 5
Split 3 => 3 + 6
Split 2 => 2 + 7
Split 4 => 4 + 8
Split 3 => 3 + 9
Split 1 => 1 + 10
Split 6 => 6 + 11
Split 7 => 7 + 12
Split 8 => 8 + 13
Split 5 => 5 + 14
Split 2 => 2 + 15
Split 4 => 4 + 16
Split 2 => 2 + 17
Split 3 => 3 + 18
Split 9 => 9 + 19
Split 11 => 11 + 20
Split 4 => 4 + 21
Split 14 => 14 + 22
Split 6 => 6 + 23
Split 1 => 1 + 24
Split 4 => 4 + 25
Split 10 => 10 + 26
Split 8 => 8 + 27
Split 16 => 16 + 28
Split 7 => 7 + 29
Split 24 => 24 + 30
Split 14 => 14 + 31
Split 5 => 5 + 32
Split 15 => 15 + 33
Split 10 => 10 + 34
Split 3 => 3 + 35
Split 20 => 20 + 36
Split 2 => 2 + 37
Split 1 => 1 + 38
Split 12 => 12 + 39
Split 13 => 13 + 40
Split 34 => 34 + 41
Split 20 => 20 + 42
Split 7 => 7 + 43
Split 21 => 21 + 44
Split 22 => 22 + 45
Split 19 => 19 + 46
Split 39 => 39 + 47
Split 6 => 6 + 48
Split 32 => 32 + 49
Split 25 => 25 + 50

При need.classes.num = 200, я получил такой результат

Спойлер

Split 1 => 1 + 3
Split 2 => 2 + 4
Split 1 => 1 + 5
Split 3 => 3 + 6
Split 2 => 2 + 7
Split 4 => 4 + 8
Split 3 => 3 + 9
Split 1 => 1 + 10
Split 6 => 6 + 11
Split 7 => 7 + 12
Split 8 => 8 + 13
Split 5 => 5 + 14
Split 2 => 2 + 15
Split 4 => 4 + 16
Split 2 => 2 + 17
Split 3 => 3 + 18
Split 9 => 9 + 19
Split 11 => 11 + 20
Split 4 => 4 + 21
Split 14 => 14 + 22
Split 6 => 6 + 23
Split 1 => 1 + 24
Split 4 => 4 + 25
Split 10 => 10 + 26
Split 8 => 8 + 27
Split 16 => 16 + 28
Split 7 => 7 + 29
Split 24 => 24 + 30
Split 14 => 14 + 31
Split 5 => 5 + 32
Split 15 => 15 + 33
Split 10 => 10 + 34
Split 3 => 3 + 35
Split 20 => 20 + 36
Split 2 => 2 + 37
Split 1 => 1 + 38
Split 12 => 12 + 39
Split 13 => 13 + 40
Split 34 => 34 + 41
Split 20 => 20 + 42
Split 7 => 7 + 43
Split 21 => 21 + 44
Split 22 => 22 + 45
Split 19 => 19 + 46
Split 39 => 39 + 47
Split 6 => 6 + 48
Split 32 => 32 + 49
Split 25 => 25 + 50
Split 15 => 15 + 51
Split 29 => 29 + 52
Split 18 => 18 + 53
Split 24 => 24 + 54
Split 3 => 3 + 55
Split 11 => 11 + 56
Split 14 => 14 + 57
Split 23 => 23 + 58
Split 2 => 2 + 59
Split 9 => 9 + 60
Split 27 => 27 + 61
Split 8 => 8 + 62
Split 31 => 31 + 63
Split 33 => 33 + 64
Split 21 => 21 + 65
Split 26 => 26 + 66
Split 4 => 4 + 67
Split 58 => 58 + 68
Split 41 => 41 + 69
Split 16 => 16 + 70
Split 35 => 35 + 71
Split 11 => 11 + 72
Split 37 => 37 + 73
Split 46 => 46 + 74
Split 32 => 32 + 75
Split 19 => 19 + 76
Split 6 => 6 + 77
Split 44 => 44 + 78
Split 17 => 17 + 79
Split 45 => 45 + 80
Split 40 => 40 + 81
Split 4 => 4 + 82
Split 13 => 13 + 83
Split 10 => 10 + 84
Split 30 => 30 + 85
Split 7 => 7 + 86
Split 46 => 46 + 87
Split 36 => 36 + 88
Split 72 => 72 + 89
Split 43 => 43 + 90
Split 42 => 42 + 91
Split 62 => 62 + 92
Split 25 => 25 + 93
Split 33 => 33 + 94
Split 38 => 38 + 95
Split 39 => 39 + 96
Split 31 => 31 + 97
Split 53 => 53 + 98
Split 6 => 6 + 99
Split 12 => 12 + 100
Split 37 => 37 + 101
Split 45 => 45 + 102
Split 2 => 2 + 103
Split 28 => 28 + 104
Split 68 => 68 + 105
Split 17 => 17 + 106
Split 18 => 18 + 107
Split 92 => 92 + 108
Split 60 => 60 + 109
Split 40 => 40 + 110
Split 75 => 75 + 111
Split 47 => 47 + 112
Split 27 => 27 + 113
Split 49 => 49 + 114
Split 16 => 16 + 115
Split 33 => 33 + 116
Split 67 => 67 + 117
Split 50 => 50 + 118
Split 9 => 9 + 119
Split 11 => 11 + 120
Split 35 => 35 + 121
Split 90 => 90 + 122
Split 48 => 48 + 123
Split 34 => 34 + 124
Split 20 => 20 + 125
Split 71 => 71 + 126
Split 24 => 24 + 127
Split 51 => 51 + 128
Split 32 => 32 + 129
Split 5 => 5 + 130
Split 44 => 44 + 131
Split 55 => 55 + 132
Split 25 => 25 + 133
Split 4 => 4 + 134
Split 91 => 91 + 135
Split 98 => 98 + 136
Split 84 => 84 + 137
Split 56 => 56 + 138
Split 29 => 29 + 139
Split 80 => 80 + 140
Split 1 => 1 + 141
Split 96 => 96 + 142
Split 70 => 70 + 143
Split 22 => 22 + 144
Split 14 => 14 + 145
Split 23 => 23 + 146
Split 35 => 35 + 147
Split 3 => 3 + 148
Split 13 => 13 + 149
Split 52 => 52 + 150
Split 65 => 65 + 151
Split 17 => 17 + 152
Split 68 => 68 + 153
Split 61 => 61 + 154
Split 28 => 28 + 155
Split 100 => 100 + 156
Split 48 => 48 + 157
Split 74 => 74 + 158
Split 31 => 31 + 159
Split 76 => 76 + 160
Split 41 => 41 + 161
Split 38 => 38 + 162
Split 55 => 55 + 163
Split 79 => 79 + 164
Split 26 => 26 + 165
Split 3 => 3 + 166
Split 90 => 90 + 167
Split 15 => 15 + 168
Split 59 => 59 + 169
Split 5 => 5 + 170
Split 75 => 75 + 171
Split 142 => 142 + 172
Split 87 => 87 + 173
Split 21 => 21 + 174
Split 10 => 10 + 175
Split 46 => 46 + 176
Split 63 => 63 + 177
Split 37 => 37 + 178
Split 124 => 124 + 179
Split 73 => 73 + 180
Split 24 => 24 + 181
Split 2 => 2 + 182
Split 139 => 139 + 183
Split 7 => 7 + 184
Split 110 => 110 + 185
Split 137 => 137 + 186
Split 36 => 36 + 187
Split 77 => 77 + 188
Split 130 => 130 + 189
Split 58 => 58 + 190
Split 106 => 106 + 191
Split 66 => 66 + 192
Split 82 => 82 + 193
Split 53 => 53 + 194
Split 30 => 30 + 195
Split 64 => 64 + 196
Split 60 => 60 + 197
Split 57 => 57 + 198
Split 153 => 153 + 199
Split 8 => 8 + 200

Теперь хотелось бы понять, что получилось... Как я понял, программа выдаёт список последних N - need.classes.num кластеров. То есть "Split 106 => 106 + 191 означает, что 106-ая запись объединилась с 191-ой. И подразумевается, что 106-ая запись до этого также объединялась с различными записями, но 106-ая там присутствует точно. Тогда, чтобы охватить все облако точек принадлежащих кластеру, мне необходимо присвоить переменной need.classes.num <- N-1 ?

[womak]

Продолжаю пляски с бубном вокруг R...
Ниже представлена функция hclust

Код: Выделить всё

> hclust
function (d, method = "complete", members = NULL) 
{
    METHODS <- c("ward", "single", "complete", "average", "mcquitty", 
        "median", "centroid")
    method <- pmatch(method, METHODS)
    if (is.na(method)) 
        stop("invalid clustering method")
    if (method == -1) 
        stop("ambiguous clustering method")
    n <- as.integer(attr(d, "Size"))
    if (is.null(n)) 
        stop("invalid dissimilarities")
    if (n < 2) 
        stop("must have n >= 2 objects to cluster")
    len <- as.integer(n * (n - 1)/2)
    if (length(d) != len) 
        (if (length(d) < len) 
            stop
        else warning)("dissimilarities of improper length")
    if (is.null(members)) 
        members <- rep(1, n)
    else if (length(members) != n) 
        stop("invalid length of members")
    hcl <- .Fortran(C_hclust, n = n, len = len, method = as.integer(method), 
        ia = integer(n), ib = integer(n), crit = double(n), members = as.double(members), 
        nn = integer(n), disnn = double(n), flag = logical(n), 
        diss = as.double(d), PACKAGE = "stats")
    hcass <- .Fortran(C_hcass2, n = as.integer(n), ia = as.integer(hcl$ia), 
        ib = as.integer(hcl$ib), order = integer(n), iia = integer(n), 
        iib = integer(n), PACKAGE = "stats")
    tree <- list(merge = cbind(hcass$iia[1L:(n - 1)], hcass$iib[1L:(n - 
        1)]), height = hcl$crit[1L:(n - 1)], order = hcass$order, 
        labels = attr(d, "Labels"), method = METHODS[method], 
        call = match.call(), dist.method = attr(d, "method"))
    class(tree) <- "hclust"
    tree
}

Раз в функции нет циклов, то, на мой взгляд дилетанта, все действия по формированию кластерного дерева совершаются в
hcl <- .Fortran(C_hclust, n = n, len = len, method = as.integer(method), ia = integer(n), ib = integer(n), crit = double(n), members = as.double(members), nn = integer(n), disnn = double(n), flag = logical(n), diss = as.double(d), PACKAGE = "stats")
hcass <- .Fortran(C_hcass2, n = as.integer(n), ia = as.integer(hcl$ia), ib = as.integer(hcl$ib), order = integer(n), iia = integer(n), iib = integer(n), PACKAGE = "stats")
Могу я где-нибудь посмотреть процедуры C_hclust и C_hcass2 или они присутствуют в R уже в скомпилированные? Если это возможно, то появляется надежда сделать вывод процесса объединения кластеров в файл.

[gamm]

Теперь хотелось бы понять, что получилось... Как я понял, программа выдаёт список последних N - need.classes.num кластеров.

программа выдает таблицу, в столбцах которой номера кластеров, которым принадлежат строки. Каждый столбец соответствует числу классов, записанному в заголовке. Строки с одинаковым номером в столбце это и есть кластер при заданном (в заголовке столбца) числе классов. Что такое "охватить все облако", я не понимаю

[gamm]

Могу я где-нибудь посмотреть процедуры C_hclust и C_hcass2 или они присутствуют в R уже в скомпилированные? Если это возможно, то появляется надежда сделать вывод процесса объединения кластеров в файл.

процесс объединения кластеров выводится на экран, причем примерно в том же виде, как в Статистике (нет только порога расстояния), сделайте sink(filename) - пойдет в файл. Исходные тексты (в данном случае на Фортране) загружаются с сайта проекта R. Зачем они вам нужны - непонятно, поскольку сам процесс неоднократно описан в литературе, начиная с книжек Айвазяна 70-х годов.

Собственно, вся информация лежит в объекте, полученном из hclust(), можете его парсить, там опять примерно то, что выдает Статистика (см. merge)

merge an n-1 by 2 matrix. Row i of merge describes the merging of clusters at step i of the clustering. If an element j in the row is negative, then observation -j was merged at this stage. If j is positive then the merge was with the cluster formed at the (earlier) stage j of the algorithm. Thus negative entries in merge indicate agglomerations of singletons, and positive entries indicate agglomerations of non-singletons.

height a set of n-1 non-decreasing real values. The clustering height: that is, the value of the criterion associated with the clustering method for the particular agglomeration.

order a vector giving the permutation of the original observations suitable for plotting, in the sense that a cluster plot using this ordering and matrix merge will not have crossings of the branches.
labels labels for each of the objects being clustered.
call the call which produced the result.
method the cluster method that has been used.

[womak]

Возможно, я не четко сформулировал задачу… Я лучше нарисую…

В конечном итоге я хочу, задав пределы изменения расстояния от min до max, получить набор кластеров (A,B,C, …..) в виде:
l1;32;33;34;30;31;55;60;58;59;57;61
l2;50;51;52;54;47;53;48;49;28;6;39;4;27;3;1;25;2;26
l3;41;44;43;45;40;46;8;13;9;15;11;35;7;14;10;12;36
l4 ……
……
Lm
Где min < l1.....lm < max


А «облако» получится тогда, когда я вынесу все записи одного кластера на плоскость по координатам, которые находятся в каждой записи исходных данных. А если рассматривать только записи (случаи), то это правильнее назвать «веткой» кластерного дерева (ветки А, В, С).

[gamm]

В конечном итоге я хочу, задав пределы изменения расстояния от min до max, получить набор кластеров (A,B,C, …..) в виде
l1;32;33;34;30;31;55;60;58;59;57;61
l2;50;51;52;54;47;53;48;49;28;6;39;4;27;3;1;25;2;26

не подскажете, где в строке l1;32;33;34;30;31;55;60;58;59;57;61 кластер? или вы вместо меток классов написали уровни по ошибке?

у вас же все есть - есть соответствие между числом классов и порогом (берете из height и номера строки, они упорядочены), для порога смотрите, между какими height попало ваше значение (используете approx() и floor()), получаете число классов k и используете cutree(p.clust,k=k), чтобы получить списки. И готовы ваши облака (строки с одинаковыми номерами). Порядок вывода лежит в order, если посмотрите исходники, сможете посчитать и координаты на картинке если что-то красить нужно. Вам даже перенумеровывать не нужно.

[womak]

gamm, огромнейшее Вам спасибо за разъяснения!
постараюсь разобраться с приведенными Вами функциями.

не подскажете, где в строке l1;32;33;34;30;31;55;60;58;59;57;61 кластер? или вы вместо меток классов написали уровни по ошибке?

Возможно мы по разному понимаем понятие "кластер". Я - как группа объектов с определенной степени похожести.
Строка l1;32;33;34;30;31;55;60;58;59;57;61 означает, в моем понимании, что записи из таблицы исходных данных с №№ 32;33;34;30;31;55;60;58;59;57;61 входят в один кластер A с расстоянием связи l1

можете его парсить

- слаб в терминологии

[gamm]

Возможно мы по разному понимаем понятие "кластер".

возможно, и мое понимание - общепринятое в математике. Успехов!

[womak]

Имеется вот такой код и результаты:

Код: Выделить всё

> x <- read.table('test.dat', sep=';', header= FALSE) 
>  x <- x[-1] 
>  d <-  dist(x, method = "euclidean", diag = FALSE, upper = FALSE) 
>  hc <- hclust(d) 
>  plot(hc)
> 
> hc$order
 [1] 21  4 18  5  6 14  9 15 19  8  7 26 22 10 11 12 16 20 24 23 25  2 13  1  3 17
> hc$height
 [1] 0.02251455 0.02847891 0.02847891 0.03184039 0.04151478 0.05034460 0.05134108
 [8] 0.05871077 0.05871077 0.06447187 0.08055051 0.08117737 0.08117859 0.09552116
[15] 0.10119033 0.12124554 0.16141626 0.16141626 1.37560364 1.52518881 2.03540001
[22] 2.48539311 4.36002875 5.00284447 8.52590503
> hc$merge
      [,1] [,2]
 [1,]   -3  -17
 [2,]   -7  -26
 [3,]  -16  -20
 [4,]  -15  -19
 [5,]  -12    3
 [6,]   -4  -18
 [7,]   -6  -14
 [8,]   -9    4
 [9,]   -2  -13
[10,]   -1    1
[11,]   -8    2
[12,]   -5    7
[13,]  -11    5
[14,]    8   11
[15,]    9   10
[16,]  -10   13
[17,]  -25   15
[18,]    6   12
[19,]  -21   18
[20,]  -22   16
[21,]  -23   17
[22,]   14   20
[23,]   19   22
[24,]  -24   21
[25,]   23   24

Получено такое дерево:

пример.jpeg (55.1 КБ) 28215 просмотров

Очень хотелось бы увидеть комментарии к результатам выполнения функций hc$order, hc$merge и hc$height.
hc$order - выведены номера строк так, как они показаны на рисунке.
hc$height - высоты в порядке возрастания, в том же порядке выведены данные функции hc$merge?
hc$merge - знак минус говорит о том, что запись впервые участвует в объединении. Рассматривая рисунок, видим, что объединяются -3 с -17, -7 с -26, -16 с -20, -15 с -19. А вот [5,] -12 3, указывает на то, что 12 запись соединилась с [3]-кластером[-16 -20]. Я правильно понял? А [10,] -1 1 означает, что первая запись соединилась с [1] кластером?

[gamm]

Очень хотелось бы увидеть комментарии к результатам выполнения функций hc$order, hc$merge и hc$height

во-первых, это не функции, во-вторых, описание (выдержка из help) уже было приведено выше. Успехов!