English
РусскийНа этой странице
Кластерный анализ — методические заметки
Решения, на которых построен модуль кластерного анализа: какие ограничения, отказы, запасные варианты и аргументы пакетов стоят за каждой карточкой, и подтверждающие их числа. Все числа получены на наборе-примере (cafe.ru.csv, спецификация лежит рядом) после трёх исправлений со страницы о данных. Разделы следуют разделам руководства.
Требования
Любое k должно быть меньше числа объектов, и максимум диапазона на шаге 2 тоже. При k, равном числу объектов, каждый объект образует отдельный кластер: cluster::silhouette() возвращает в этом случае просто NA, из-за которого не вычисляется весь блок метрик качества, а индекс Калинского — Харабаша делит на n − k. Поэтому в режиме наблюдений нужно k + 1 полных наблюдений, а в режиме переменных — k + 1 переменных и по-прежнему два полных наблюдения: dist() по транспонированной матрице нужны хотя бы две строки, чтобы между двумя переменными был какой-то разброс. Отказ выводится во всплывающем уведомлении, где названы оба числа, ещё до обращения к R.
Проверка. По оценке, чеку и числу визитов после смены меню полные данные есть у 192 участников кафе, поэтому k = 192 отклоняется с сообщением «Недостаточно полных наблюдений (192) для 192 кластеров», а наибольшее допустимое значение — k = 191.
Для бикластеризации нужны два полных наблюдения и не меньше, чем запрошено бикластеров, а для FABIA число бикластеров не может превышать числа переменных. У алгоритмов, сами выбирающих k, нет k, с которым можно сравнить число строк, а матрица из одной строки доходит до R и вызывает необработанную ошибку пакета. Поэтому два полных наблюдения — нижняя граница для всех алгоритмов. fabia() завершается ошибкой «Too many biclusters. Stopped.», как только число бикластеров превышает min(nrow, ncol), а ограничивающими здесь всегда оказываются столбцы, так что k для FABIA ограничивается числом переменных ещё до запуска. Алгоритмы Ченга — Чёрча и BiMax в таком случае просто возвращают меньше бикластеров, чем запрошено, и ограничение им не нужно.
Постоянные переменные исключаются до запуска и перечисляются. scale() превращает постоянный столбец в NaN, который kmeans() отвергает, а dist() без предупреждения отбрасывает, поэтому постоянная на полных наблюдениях переменная либо прервала бы расчёт, либо бесследно из него исчезла. Модуль исключает её заранее, выводит в строке Исключены (без дисперсии) и отказывается от запуска, если осталось меньше двух переменных.
Числовые поля читаются после исправления, а не как введены. Пустое числовое поле разбирается как NaN, а любое сравнение с NaN ложно, поэтому непроверенное значение проходит все проверки диапазона и попадает в код R буквальным NaN. Значение, которое браузер принял, но не проверил, доходит до R как недопустимый аргумент (probs = 1.5, iter.max = 0, p = 0). Поэтому каждый числовой параметр читается как значение, ограниченное объявленными в разметке min и max, а при пустом поле — как значение по умолчанию. Исключение — два поля k на шаге 3: они читаются как введены, и пустое значение или значение ниже минимума отклоняется всплывающим уведомлением. Диапазон шага 2 тоже поддерживается в рабочем виде: если исправленный максимум опустился до минимума, он поднимается на единицу выше, поскольку для сравнения нужен диапазон.
Шаг 1: метод и настройки
В режиме переменных стандартизация принудительно включена, при метрике Гауэра — выключена, а выбор пользователя тем временем сохраняется. Квадрат евклидова расстояния между двумя столбцами в z-оценках равен , так что именно стандартизация превращает расстояние в режиме переменных в корреляцию. Без стандартизации расчёт измерял бы под тем же названием нечто иное, поэтому флажок установлен и заблокирован. Коэффициент Гауэра делит разность по каждой числовой переменной на её собственный размах, а это отношение не меняется ни при каком линейном преобразовании столбца. Предварительная z-стандартизация не изменила бы ни одного коэффициента, хотя карточка сообщала бы о ней, и в этом случае флажок, наоборот, снят и заблокирован. Пока действует любая из блокировок, выбор пользователя запоминается и восстанавливается, когда блокировка снимается.
Проверка. На шести пунктах анкеты в режиме переменных дерево Уорда по евклидову расстоянию и дерево по 1 − r при разрезе на два дают одно и то же разбиение (медленно против остальных пяти), а первое слияние евклидова дерева, пунктов вкус и свежесть на высоте 13.95, равно , где 0.489 — высота первого слияния дерева по 1 − r.
Поле метрики скрыто при К-средних в режиме наблюдений и закреплено на евклидовой метрике в режиме переменных. Недоступные пункты отключаются, а не скрываются. Метод k-средних минимизирует квадрат евклидова расстояния (именно среднее — точка, которая это делает), поэтому другую метрику он учесть не может. В режиме наблюдений поле скрывается без сброса, и метрика, выбранная для К-медоидов или иерархической кластеризации, сохраняется после временного переключения на К-средних. В режиме переменных поле остаётся на экране, потому что здесь метрика и определяет, что значит «похожие переменные», а скрытие сделало бы недоступными оба корреляционных пункта. Пункты, недоступные для К-средних, отключаются, а не скрываются: скрытый пункт остаётся доступен с клавиатуры и через поиск по вводу, и поле могло бы хранить значение, которое расчёт заменяет без предупреждения. Что бы ни было выбрано в поле, расчёт использует метрику, которую может учесть: евклидову при К-средних и евклидову же, если в режиме переменных осталась метрика Гауэра (усреднение оценок по переменным теряет смысл, когда объектами становятся столбцы). Поэтому оставшееся от прежнего выбора значение никогда не оценивает решение в геометрии, которую метод не оптимизировал. Закрепление ничего не отнимает: на стандартизованных переменных евклидов метод k-средних по тождеству выше и так воспроизводит кластеризацию по 1 − r.
Метод Уорда при неевклидовой метрике сопровождается предупреждением, а каждое дерево строится в той шкале, которую предполагает метод связи. Центроидный и медианный методы — рекуррентные формулы Ланса — Уильямса для квадратов евклидовых расстояний, ward.D2 сам возводит входные расстояния в квадрат, а 1 − r и 1 − s Гауэра уже являются квадратами расстояний: обе меры сходства положительно полуопределены, поэтому их дополнения вкладываются в евклидово пространство как квадраты расстояний. Что именно получает hclust(), определяется сочетанием метрики и метода связи. Центроидный и медианный методы на обычном dist() получают d², и из их высот затем извлекается корень. Метод ward.D2 на корреляционной мере или мере Гауэра получает √d, и его высоты затем возводятся в квадрат. Остальные сочетания передаются как есть. В результате высоты любого дерева выражены в шкале самой меры различия, и по ней строятся ось дендрограммы и кофенетическая корреляция. Предупреждение о методе Уорда выводится только при неевклидовой метрике (метрика Минковского при p = 2 считается евклидовой). Предупреждение о центроидном и медианном методах выводится при любой метрике, поскольку hclust определяет оба метода для квадратов евклидовых расстояний и на обычном объекте dist() даёт немонотонные высоты слияний. Инверсия — содержательный диагностический признак, поэтому высоты рисуются как есть и подсчитываются в примечании на графике, а не сглаживаются.
Проверка. На трёх стандартизованных переменных кафе центроидный метод даёт 9 инверсий на 191 слияние, и примечание появляется. Три последних слияния метода Уорда лежат на высотах 12.48, 12.77 и 16.21 в евклидовой шкале.
Для коэффициента Гауэра важен класс столбца, и приложение порождает три из шести подтипов daisy. daisy() выбирает способ сравнения по классу столбца. Модуль передаёт категориальную переменную как обычный фактор, а числовой столбец оставляет без изменений. Порядковая переменная передаётся как упорядоченный фактор с уровнями, отсортированными как числа (иначе ordered() сортирует их как текст и ставит 10 перед 2), что соответствует ранговой трактовке Подани. Коды подтипов считываются из атрибута Types меры различия и выводятся строкой Типы переменных для Гауэра: больше нигде на карточке не видно, как оценивалась переменная. Для асимметричного бинарного, симметричного бинарного и логарифмического подтипов daisy нужен аргумент type =, который модуль не передаёт, поэтому в строке могут стоять только «интервальная», «номинальная» и «порядковая»: категориальная переменная с двумя уровнями будет номинальной, а числовая переменная 0/1 — интервальной. Статья Подани цитируется, только если в расчёте была порядковая переменная, иначе его трактовка не применялась.
Проверка. Для оценки, чека, визитов, филиал и рекомендует строка показывает «интервальная, интервальная, интервальная, номинальная, номинальная», а коэффициент между первыми двумя участниками равен 0.257 — среднему из 0.111, 0.079, 0.093, 0 и 1.
При метрике Гауэра центроидные величины не выводятся, а индекс Дэвиса — Болдина переходит на медоиды. Суммы квадратов, объяснённая дисперсия, индексы Калинского — Харабаша и Хартигана, статистика Хопкинса, критерий Дуды — Харта, статистика разрыва и график локтя — центроидные величины для числовой матрицы. У таблицы смешанных типов нет среднего, нет ограничивающего параллелепипеда для равномерных пробных точек Хопкинса, нет Je для Дуды — Харта и нет параллелепипеда главных компонент для эталона статистики разрыва. Вычислять их по евклидовой матрице, по которой расчёт не кластеризовал, значило бы оценивать решение в другой геометрии, поэтому они не выводятся, а на карточке об этом сообщает примечание. Индекс Дэвиса — Болдина пересчитывается по медоидам кластеров на мере различия: это альтернативная параметризация самого индекса, а не подмена статистики. Силуэт, индекс Данна, кофенетическая корреляция, бутстрэп-устойчивость и консенсусное разбиение работают с мерой различия напрямую и не меняются. Биплот опирается на prcomp(), которой нужны числовая матрица и по стрелке на каждую переменную кластеризации, поэтому он тоже не строится и заменяется примечанием. К-средних недоступен, поскольку ему нужны центроиды. Профили кластеров и медоиды показывают модальную категорию там, где среднее не определено. F и η² в таблице вклада переменных сравнивают средние, поэтому для категориальной переменной обе ячейки пусты, о чём сообщает примечание.
Несходимость К-средних определяется по ifault и числу итераций. kmeans() сообщает ifault 2 (нет сходимости за iter.max итераций) и 4 (превышено число шагов быстрого переноса) только для варианта Хартигана — Вонга, остальные три варианта ничего не сообщают. Поэтому расчёт для любого варианта считает сбоем и число итераций, достигшее предела. Предупреждение называет достигнутое число и два средства — поднять предел или увеличить число запусков. Случайные начальные точки задаются зерном воспроизводимости, так что при тех же данных и настройках выбираются те же начальные точки и лучший запуск даёт то же решение.
Проверка. На трёх стандартизованных переменных кафе при k = 3 лучший из 25 запусков сходится за 2 итерации, предупреждения нет.
Фон в модели Plaid — общее среднее плюс эффекты строк и столбцов, а слой сохраняется, только если превосходит три перемешанные копии. В документации biclust background = TRUE функции BCPlaid описан как фоновый слой, «постоянный для всех строк и столбцов», но при background.layer, оставленном NA, этот слой подгоняется той же моделью, что и остальные слои, y ~ m + a + b, по всей матрице, с 1 + (n − 1) + (p − 1) степенями свободы. Иначе говоря, перед поиском первого бикластера вычитается не константа, а общее среднее плюс эффект каждого наблюдения и каждой переменной. При FALSE фон не подгоняется, и первый слой ищется по исходным значениям. Сумма квадратов каждого кандидата в слои сравнивается со слоями, подогнанными на shuffle = 3 перемешанных копиях матрицы, и кандидат сохраняется, только если превосходит все три. Так модель и останавливается сама. По правилу удаления наблюдение остаётся в слое, пока остаточный средний квадрат слоя для него меньше произведения (1 − release) на его общий средний квадрат по переменным слоя. Для переменной то же правило действует по наблюдениям слоя. Поэтому чем больше значение удаления, тем жёстче отсев.
Проверка. На шести пунктах анкеты Plaid не находит ни одного слоя, со стандартизацией и без неё: ни один кандидат не превосходит перемешанные копии, поскольку структура пунктов — два фактора, затрагивающих всех участников, а не отдельные блоки.
BiMax бинаризует каждый столбец по его собственному квантилю. binarize() из biclust использует один порог для всей матрицы, а без стандартизации такой порог для каждой переменной означает разное. Поэтому расчёт режет каждый столбец по выбранному квантилю его собственных значений (выше — 1, на пороге и ниже — 0) и передаёт бинарную матрицу в BCBimax() вместе с двумя минимальными размерами.
Проверка. На шести пунктах анкеты, разрезанных по медианам, при k = 3 получаются блоки 36 участников × {цена, приветливо}, 11 × {быстро, медленно} и 36 × {вкус, цена}.
Порог δ алгоритма Ченга — Чёрча задаётся как доля собственного остатка матрицы. Найденные строки маскируются, а не удаляются. biclust принимает delta как абсолютный потолок среднего квадрата остатков, который без шкалы матрицы ничего не значит, поэтому значение пользователя умножается на остаток всей матрицы. BCCC() не передаёт аргумент rand, а по умолчанию biclust ставит его в FALSE, и тогда строки каждого найденного бикластера удаляются из матрицы, а не маскируются случайными значениями, как предписывают Ченг и Чёрч. Поэтому расчёт вызывает внутреннюю процедуру пакета с rand = TRUE. alpha — коэффициент группового удаления узлов: строки и столбцы, средний остаток которых больше остатка блока в α раз, удаляются все сразу, до поочерёдного удаления и добавления узлов.
Проверка. Средний квадрат остатков стандартизованной матрицы шести пунктов равен 0.735, так что δ = 0.2 даёт потолок 0.147. Первый найденный бикластер — 42 участника на пяти пунктах, кроме медленно, с остатком 0.141.
Принадлежность в FABIA определяется порогом, а не априорной разреженностью. Разреженность в FABIA мягкая, и в точности нулевыми нагрузки и факторы не становятся никогда, поэтому принадлежность берётся из индексов строк и столбцов, отобранных по порогу функцией extractBic(fit, thresZ). Порог для нагрузок extractBic вычисляет как √mom / thresZ, и поэтому с ростом порога в бикластеры попадает больше наблюдений и меньше переменных. Бикластер, в котором не осталось ни одного наблюдения или ни одной переменной, отбрасывается, поэтому при слишком низком пороге не находится ничего. spl и spz определяют подгонку, cyc — число её итераций. Заставку FABIA и сообщения на каждом цикле модуль подавляет, чтобы в консоли было только то, что выводит сам расчёт.
Числа групп и векторов у спектрального алгоритма ограничены размерами матрицы, а порог разброса задаётся как доля собственного разброса матрицы. BCSpectral() не передаёт n_clusters, а без него biclust запускает автоматический поиск, окно которого на узкой матрице схлопывается, и расчёт завершается ошибкой «missing value where TRUE/FALSE needed». Поэтому расчёт вызывает внутреннюю процедуру пакета с заданным числом. numberOfEigenvalues ограничивается значением min(nrow, ncol) и в R, и на панели, поскольку максимум элемента управления актуален лишь на момент последнего выбора. Число групп ограничено на панели наименьшим из трёх значений: min(nCols, nRows) − 1 (обе оси группируются методом Хартигана — Вонга, которому нужно строго меньше центров, чем точек), nCols − minc и nRows − minr. Дело в том, что при g группах самая широкая группа содержит m − g + 1 столбцов, а бикластер сохраняется, только если он шире minc, так что любое g больше m − minc не даёт ничего при сколь угодно мягком пороге разброса. withinVar в biclust — не дисперсия, а среднее расстояние от строки до центроида во входной шкале, поэтому доля, заданная пользователем, умножается на среднее расстояние от строки до центроида для всей матрицы. Сочетание групп сохраняется, если в нём больше minr строк, больше minc столбцов и разброс ниже порога. Логарифмическая нормализация сдвигает матрицу так, чтобы минимум стал равен 1, затем логарифмирует и вычитает средние по строкам и столбцам. IRRC и бистохастизация сдвигают матрицу в неотрицательную область, и при обоих методах первый сингулярный вектор отражает фон и отбрасывается, поэтому запрос одного сингулярного вектора с предупреждением повышается до двух.
Шаг 2: выбор числа кластеров
Тенденция к кластеризации проверяется только в режиме наблюдений, только на числовой матрице и в евклидовой геометрии при любой метрике. Статистика Хопкинса опирается на ограничивающий параллелепипед кластеризуемых объектов, а эталонное распределение Дуды — Харта зависит от их размерности, поэтому на транспонированной матрице режима переменных обе величины описывали бы число наблюдений, а не переменные. У таблицы смешанных типов нет ни параллелепипеда, ни Je, поэтому при метрике Гауэра не выполняется ни один из двух критериев. Для статистики Хопкинса берётся по одной пробной точке на десять объектов (не меньше двух и не больше n − 1), равномерно по размаху каждой переменной в матрице анализа, и столько же реальных объектов, всё при зерне воспроизводимости. Статистика равна по евклидовым расстояниям до ближайшего соседа — в той ориентации, где 0.5 соответствует случайному разбросу, а 1 — кластерам. Критерий Дуды — Харта делит данные выбранным алгоритмом при k = 2, вычисляет Je(2)/Je(1) как отношение евклидовой внутрикластерной суммы квадратов этого разбиения к общей и проверяет его односторонним критерием по нормальной аппроксимации Дуды и Харта (среднее , дисперсия для n объектов и p переменных) на уровне значимости из раздела Настройки.
Проверка. На трёх стандартизованных переменных кафе статистика Хопкинса по 19 пробным точкам равна 0.811, а лучшее разбиение k-средних на два кластера оставляет 410.51 из 573 единиц вариации: отношение Дуды — Харта 0.716, p = 0.078.
Статистика разрыва вычисляется по собственной мере различия расчёта, а k выбирается по правилу самих Тибширани, Уолтера и Хасти. cluster::clusGap измеряет W(k) в евклидовом расстоянии, поэтому модуль вычисляет разрыв сам. W(k) складывает по кластерам суммы квадратов внутрикластерных различий, делённые на размер кластера, в той мере различия, по которой идёт кластеризация. Каждая из B = 100 эталонных выборок генерируется равномерно в параллелепипеде главных компонент данных (второй эталон из статьи, повторяющий разброс и ориентацию самих данных) и кластеризуется тем же алгоритмом, всё при зерне воспроизводимости. Стандартная ошибка равна . k выбирается функцией maxSE(method = "Tibs2001SEmax"), то есть как наименьшее k, для которого Gap(k) ≥ Gap(k + 1) − s(k + 1), — это правило статьи, в отличие от принятого в cluster по умолчанию firstSEmax. Правило проходит по k от 1 до максимума диапазона, каким бы ни был минимум, поэтому может вернуть 1 или k ниже диапазона. Отметка в таблице и выделение на графике следуют правилу, а не максимуму столбца. Для корреляции между профилями равномерный параллелепипед не служит бесструктурной нулевой моделью, поэтому при 1 − r и 1 − |r| выводится примечание. При метрике Гауэра параллелепипеда нет, и статистика разрыва не вычисляется.
Проверка. На трёх стандартизованных переменных кафе при К-средних разрыв при k от 2 до 10 держится между 0.976 и 1.032, и правило возвращает k = 1. Карточка сообщает, что статистика разрыва указывает на отсутствие кластерной структуры, тогда как силуэт рекомендует 3.
Центроидные индексы евклидовы при любой метрике, а индексы разделения работают с мерой различия самого расчёта. Внутрикластерная СК, объяснённая дисперсия, индексы Калинского — Харабаша и Хартигана — суммы квадратов относительно средних кластеров, и общепринятой неевклидовой формы у них нет. Поэтому при другой метрике они вычисляются по центроидам в матрице анализа, а примечание под таблицей предупреждает, что в таблице смешаны две геометрии. При метрике Гауэра среднее взять не из чего, и эти величины не выводятся (индекс Калинского — Харабаша — вместе с флажком). Силуэт и индекс Данна вычисляются по той же матрице различий, на которой шла кластеризация. Индекс Дэвиса — Болдина его авторы определяют через расстояние, поэтому в этой метрике измеряются и разбросы, и расстояния между центроидами, а при метрике Гауэра — по медоидам.
Для внутрикластерной СК и объяснённой дисперсии лучшее значение не отмечается, а индекс Хартигана сравнивается с 10. Обе величины монотонно улучшаются с ростом k, и отметка указывала бы на наибольшее проверенное k, а именно этого таблица для выбора k советовать не должна. Индекс Хартигана равен . По эмпирическому правилу Хартигана значение больше 10 — повод добавить кластер. Индексу нужно следующее k, поэтому в последней строке диапазона стоит Н/Д.
Проверка. На переменных кафе индекс Хартигана при k = 2 равен 50.7 = (410.51 ⁄ 323.67 − 1) × 189 и остаётся больше 10 вплоть до k = 9.
Рекомендация по силуэту — простой максимум, но с оговоркой о диапазоне. Рекомендуемое k по силуэту — k с наибольшим средним силуэтом среди тех, где силуэт определён (2 ≤ k < n). Если таких k нет, берётся минимум диапазона. Лучшее значение в диапазоне ничего не говорит о том, есть ли в нём структура, поэтому, если ни одно k не достигает 0.25 (порог Кауфмана и Руссеу, ниже которого существенной структуры нет), под рекомендациями выводится оговорка с лучшим значением.
Рекомендуемое число бикластеров — самое устойчивое k, а если устойчивость не измерялась, — локоть объяснённой дисперсии. У плоской кривой локтя нет. Если устойчивость измерена, рекомендуется k с наибольшим средним бутстрэп-коэффициентом Жаккара. Иначе рекомендуется последнее k перед тем, как прирост объяснённой дисперсии падает ниже половины предыдущего прироста. Но одно такое отношение срабатывает и на шуме: BiMax на тестовом наборе, которым проверяется модуль, набирает по k от 2 до 8 всего 0.45 процентного пункта, и первая же смена знака принимается за локоть, тогда как алгоритм Ченга — Чёрча на том же диапазоне набирает 37 пунктов. Поэтому правило применяется, только если кривая поднялась хотя бы на 5 пунктов выше начального значения. В противном случае рекомендуется наименьшее проверенное k, а карточка называет применённое правило. Если устойчивость была запрошена, но ни для одного k не вычислена, предупреждение называет запасное правило.
Проверка. BiMax на шести пунктах анкеты при k от 2 до 6 даёт приросты 6.28, 4.72 и 0.93 пункта, то есть после k = 4 прирост обрывается, и рекомендуется 4.
Устойчивость бикластеров оценивается по бутстрэп-выборкам наблюдений при неизменных переменных. Решение при каждом k (то, что показано в таблице) оценивается покластерным бутстрэп-коэффициентом Жаккара Хеннига (Hennig, 2007). Строки извлекаются с возвращением, алгоритм заново подгоняется на бутстрэп-выборке, каждый исходный бикластер, ограниченный попавшими в бутстрэп-выборку строками, сопоставляется с самым похожим новым по коэффициенту Жаккара для множеств ячеек, и результаты усредняются по бикластерам и повторам. Переменные выбирает сам пользователь, поэтому они в бутстрэп-выборки не попадают. Бутстрэп-выборки зависят от зерна бутстрэпа из раздела Настройки, по умолчанию не заданного. Столбец Устойчивость в таблице бикластеров на шаге 3 — то же вычисление для одного решения. Повтор, в котором подгонка не удалась, даёт 0.
Диагностика F-статистики носит только описательный характер. Она берётся из функции diagnoseColRow пакета biclust, нулевое распределение которой строится по ячейкам, извлекаемым независимо из всей матрицы. Такая нулевая модель разрушает любую структуру строк и столбцов, вместо того чтобы сохранить блок: при ней значимо любое реальное различие между переменными, а бикластеры как раз и подбирались так, чтобы максимизировать эти эффекты. F-значения для строк и столбцов усредняются по бикластерам, в которых не меньше двух строк и двух столбцов, доли считают p-значения ниже уровня значимости из раздела Настройки, а примечание в таблице и на графике предупреждает, что обе величины завышены по построению. Линия 80% — ориентир, а не критерий. Параметрические p-значения computeObservedFstat проверяют выбранный блок относительно него самого и по той же причине ≈ 0, поэтому используются только наблюдаемые F-значения.
Консенсусная оценка сравнивает повторные запуски в обе стороны и доступна только для стохастических алгоритмов. При каждом k алгоритм заново подгоняется столько раз, сколько указано в поле повторов, каждый раз со своим зерном, выведенным из зерна бутстрэпа (или из энтропии, если оно не задано). Каждая пара запусков оценивается по коэффициенту Жаккара для ячеек лучших соответствий в обе стороны, и две оценки усредняются. Поэтому запуск, нашедший два бикластера, не может получить высокую оценку в паре с запуском, нашедшим десять, а неудачный запуск не выпадает, а даёт пары с нулевой оценкой. При 100 запусках по умолчанию это 4950 пар на каждое k. BiMax не использует случайных чисел, его запуски одинаковы, поэтому при нём этот параметр скрыт.
Шаг 3: полный анализ
Силуэт оценивается по шкале Кауфмана и Руссеу, а кофенетическая корреляция для дерева выводится всегда. Средняя ширина силуэта получает оценку «сильная структура» от 0.70, «достаточная» от 0.50, «слабая» от 0.25 и «нет существенной структуры» ниже этого — по Кауфману и Руссеу. Средние силуэты кластеров в таблице размеров берутся из того же вызова cluster::silhouette() и потому выводятся всякий раз, когда силуэт вычисляется, для строки метрик качества или для силуэтного графика. Кофенетическая корреляция показывает, насколько точно высоты слияний дерева воспроизводят исходные различия. Она относится к дереву, а не к разрезу, поэтому флажка у неё нет, для всех k она одинакова и сравнивается с общепринятым порогом 0.75.
Устойчивость кластеров вычисляется по clusterboot Хеннига: уникальные объекты бутстрэп-выборки, с отнесением обратно каждого объекта. Силуэт, индекс Калинского — Харабаша и другие индексы оценивают разбиение на той же выборке, по которой оно построено, и не отличают настоящую группу от удачного разреза. Бутстрэп по Хеннигу (Hennig, 2007) проверяет, появился бы кластер снова в другой выборке. В каждом повторе объекты извлекаются с возвращением и оставляются только уникальные, потому что повторённый объект тянул бы центроид к себе. Затем они заново кластеризуются тем же алгоритмом с теми же настройками, и каждый объект полных данных относится к ближайшему бутстрэп-представителю в мере различия расчёта: к центроиду k-средних, медоиду PAM, среднему кластера иерархического разреза или его медоиду при метрике Гауэра. После этого каждый исходный кластер сопоставляется с самым похожим новым по коэффициенту Жаккара для целых множеств, так что каждый повтор с успешной подгонкой оценивает все кластеры. Среднее оценивается по шкале Хеннига: от 0.85 — очень устойчивый, от 0.75 — устойчивый, от 0.60 — закономерность с сомнительным составом, ниже — доверия не заслуживает. Повтор с коэффициентом 0.5 и ниже считается распадом кластера, с 0.75 и выше — восстановлением. В режиме переменных перевыбираются наблюдения, по которым строится мера различия, поэтому все переменные присутствуют в каждом повторе и бутстрэп-решение само является распределением по кластерам. Повторы зависят от зерна бутстрэпа из раздела Настройки.
Консенсусное разбиение строится методом накопления свидетельств, режется методом средней связи и перенумеровывается по приведённому решению. Используются те же повторы, что и для устойчивости, поэтому обе диагностики вместе обходятся одним циклом. Как в методе Фреда и Джейна (Fred and Jain, 2005), распределение по кластерам в каждом повторе добавляется к счётчику коассоциаций каждой пары объектов, а в конце счётчик делится на число повторов с успешной подгонкой. Консенсусное разбиение — это hclust(as.dist(1 − C), "average"), разрезанное при k расчёта. Разрез дерева нумерует кластеры по положению в дереве, поэтому каждый консенсусный кластер получает номер приведённого кластера, с которым больше всего пересекается (начиная с наибольшего пересечения), и совпадения оказываются на диагонали таблицы сопряжённости. Согласованность — скорректированный индекс Рэнда по этой таблице, оцениваемый по шкале Стейнли (Steinley, 2004): 0.90 — отличная, 0.80 — хорошая, 0.65 — умеренная. Ниже 0.65 выводится предупреждение. PAC Шенбабаоглу и др. (Şenbabaoğlu et al., 2014) — доля пар с коассоциацией строго между 0.1 и 0.9. Средняя коассоциация кластера — это консенсус кластера по Монти и др. (Monti et al., 2003), среднее по всем различным парам внутри кластера. Для кластера из одного объекта средняя коассоциация не определена.
Суммы квадратов вычисляются в евклидовой геометрии по матрице анализа, а F в таблице вклада переменных — описательная величина. Внутри- и межкластерная суммы квадратов — центроидные величины, поэтому при любой метрике они считаются относительно средних кластеров той матрицы, по которой шла кластеризация (стандартизованной, если стандартизация включена, и транспонированной в режиме переменных). При метрике Гауэра среднего нет, и суммы квадратов не выводятся. Таблица вклада — однофакторный дисперсионный анализ каждой переменной в исходной шкале по кластерам, с df k − 1 и n − k и . Разбиение выбиралось так, чтобы разделить именно эти переменные, поэтому F велика по построению, и p-значение не выводится.
Многомерное шкалирование сообщает ту меру соответствия, которую оптимизирует метод, и ограничено 1000 объектами. Классическое шкалирование — это cmdscale(d, k = 2, eig = TRUE), а мера соответствия — GOF[1], вычисленная по модулям собственных чисел, чтобы отрицательные собственные числа неевклидовой меры различия уменьшали долю, приходящуюся на две оси, а не отбрасывались. Стресса у классического шкалирования нет. Неметрическое шкалирование — это MASS::isoMDS() с классическим решением в качестве начального. Стресс-1 этого метода возвращается в процентах, делится на 100 и оценивается по шкале Краскала (Kruskal, 1964): 0.025 — отлично, 0.05 — хорошо, 0.1 — удовлетворительно. Оба метода выполняются в скомпилированном коде, который не проверяет запросы на прерывание, поэтому размер задачи ограничивается заранее, а не отменой: при 1000 объектах cmdscale на движке работал около 1.6 с, isoMDS — около 14 с, и время обоих растёт как куб числа объектов. Подписей выводится не больше 40, начиная с самых удалённых от центра: при большем числе они уже ничего не проясняют.
Биплот строится с сохраняющим расстояния масштабированием Габриэля по матрице кластеризации. prcomp() выполняется с центрированием, но без масштабирования, по той же матрице, по которой шла кластеризация, так что картинка передаёт геометрию расчёта. Точки — значения главных компонент, стрелки — исходные нагрузки. Из семейства масштабирований Габриэля (Gabriel, 1971) только такое сохраняет расстояния: расстояния между точками на графике приближают те, с которыми работала кластеризация. Альтернатива, масштабирование с сохранением ковариаций, позволяет читать длины стрелок как стандартные отклонения, но искажает именно то разделение, ради которого строится график. Нагрузки — единичные векторы, а значения компонент выражены в единицах данных, поэтому общий множитель растягивает самую длинную стрелку до 85% размаха облака: длины стрелок сравнимы между собой, но не с осями. При метрике Гауэра биплот не строится, так как у таблицы смешанных типов нет главных компонент, и в режиме переменных, где точками уже служат сами переменные. Если наблюдений больше 5000, облако прореживается с постоянным шагом внутри каждого кластера, так что доли кластеров сохраняются.
Оптимальный порядок листьев находится точно. Это алгоритм Бар-Джозефа и др. (Bar-Joseph et al., 2001), и модуль вычисляет его в R: повороты всех ветвей подбираются так, чтобы суммарное различие между соседними листьями было минимальным. Алгоритм точный, а не эвристический, поэтому на 200 объектах занимает около секунды, к 800 объектам — десятки секунд, а при методе одиночной связи — больше.
Таблицы бикластеров, кроме профилей, строятся по матрице анализа, а MSR сравнивается с MSR самой матрицы. Среднее значение, MSR и столбцы когерентности вычисляются по матрице, с которой работал алгоритм (в z-оценках, если была стандартизация), тогда как профили бикластеров приводят средние в исходной шкале. Поэтому в этом случае столбец сводки помечен (z). MSR зависит от шкалы матрицы, и само по себе его значение ничего не говорит. Его отношение к MSR всей матрицы меньше 1 ровно тогда, когда блок согласованнее матрицы. Ср. r не учитывает строк, постоянных на столбцах блока, поскольку их корреляция не определена.
Дендрограммы тепловой карты — деревья Уорда по евклидовым расстояниям, а размер сетки ограничен. Строки и столбцы упорядочиваются функцией hclust(dist(x), "ward.D2") по матрице анализа и по транспонированной матрице, какой бы ни была бикластеризация. Сетка нужна для обзора, поэтому при числе строк больше 320 (сетка в 640 пикселей по два пикселя на строку) строки равномерно прореживаются в порядке дендрограммы, причём каждая нарисованная ячейка остаётся реальным наблюдением. Дендрограмма строк в этом случае скрывается, так как её листья не совпали бы со строками. При числе нарисованных ячеек больше 8000 подсказки для отдельных ячеек обходятся дороже, чем приносят пользы, и отключаются.
Пропущенные данные
Наблюдения исключаются построчно по переменным, которые сохраняет расчёт, что бы ни было выбрано в настройке пропущенных данных. Все результаты расчёта опираются на один набор объектов: разбиение, его силуэты и суммы квадратов, критерии тенденции, счётчики совместных попаданий в бутстрэпе. К тому же kmeans(), prcomp() и cmdscale() пропусков не допускают вовсе. Матрица, собранная попарно (каждое различие по переменным, общим для двух наблюдений, или, в режиме переменных, каждая корреляция по своим полным наблюдениям), сводила бы вместе расстояния, измеренные на разных данных, и индексы качества описывали бы выборку, которой не существует. Поэтому модуль отбирает полные наблюдения ещё до обращения к R, и попарный вариант, выбранный в настройке по умолчанию, здесь не действует. Отбор идёт по переменным, которые сохраняет расчёт: при обычной метрике — после отделения нечисловых, при метрике Гауэра — по всем выбранным. Столбец, который расчёт не использует, наблюдений не отнимает, а постоянные переменные определяются уже на полных наблюдениях. Построчное удаление и импутация учитываются, потому что действуют на набор данных до того, как его получит модуль. Полные наблюдения сохраняют исходные номера строк на протяжении всего расчёта, поэтому добавленная принадлежность к кластеру или бикластеру попадает на нужное наблюдение, а исключённое наблюдение остаётся пустым.