В предыдущих уроках мы говорили об иерархической кластеризации и о дендрограмме как наиболее естественном способе представления таких кластеров. Однако важно помнить, что результат кластеризации зависит от нескольких параметров:
число предикторов (например, наиболее частотных слов для разных произведений);
расстояние, которое используется для попарных сравнений (евклидово, косинусное, др.)
метод связи (метод полной связи, метод средней связи, метод Уорда и др.);
Иногда добавление или удаление буквально одного слова меняет структуру дерева.
17.1 Танглграммы
Поэтому бывает необходимо сравнить разные деревья. Один из способов это сделать визуально – построить tanglegram, например, при помощи пакета dendextend.
library(stylo)library(dendextend)data("galbraith")data <-as.data.frame.matrix(galbraith)distmx <-dist(scale(data))d1 <-as.dendrogram(hclust(distmx, method ="average")) |>set("labels_col", value =c("skyblue", "orange", "grey40"), k=3) |>set("branches_k_color", value =c("skyblue", "orange", "grey40"), k =3)d2 <-as.dendrogram(hclust(distmx, method ="ward.D2")) |>set("labels_col", value =c("skyblue", "orange", "grey40"), k=3) |>set("branches_k_color", value =c("skyblue", "orange", "grey40"), k =3)dlist <-dendlist(d1, d2)par(family ="Arial Bold")tanglegram(dlist, common_subtrees_color_lines =FALSE, highlight_distinct_edges =TRUE, highlight_branches_lwd=FALSE, margin_inner=10, lwd=2, axes=FALSE, main_left ="Cредняя", main_right ="Уорд", lab.cex =1.3)
На картинке видно, что книги группируются немного по-разному в зависимости от метода связи, хотя для кластеризации использовалась одна и та же матрица расстояний.
17.2 Консенсусные деревья в stylo
Консенсусное дерево позволяет “обобщить” произвольное число дендрограмм. В stylo за консенсусные деревья отвечает метод BCT (Bootstrap Consensus Tree), к которому можно обратиться через GUI (но здесь мы показываем решение без него).
Работать через GUI удобно, но есть нюансы. Во-первых, не получится кастомизировать внешний вид дерева, а, во-вторых, в Stylo реализована достаточно специфическая процедура бутстрепа (повторных выборок).
Under the FEATURES tab, users can define the minutes of the MFW division and sampling procedure, using the increment, the minimum and maximum parameters. For minimum = 100, maximum = 3000, and increment = 50, stylo will run subsequent analyses for the following frequency bands: 100 MFW, 50–150 MFW, 100–200 MFW, …, 2900–2950 MFW, 2950–3000 MFW.
Для консенсуса нужно много деревьев, и Stylo будет строить эти деревья в заданном интервале. Это значит, что последние деревья будут построены уже не на основе самой частотной лексики, т.е. скорее всего на них отразится тематика текстов, входящих в корпус.
В некоторых случаях это работает неплохо, но, возможно, у вас есть другие идеи для консенсуса. Разные расстояния. Разные методы кластеризации. Случайные выборки из первых двух сотен слов или еще что-то. Тогда придется самим строить сразу множество деревьев.
17.3 Целый лес c purrr
Если изучить изнанку функции stylo(), которая вызывает GUI в одноименном пакете, то можно заметить, что за консенсусное дерево там отвечает пакет для работы с филогенетическими данными под названием Ape 🦍
Что делает машина, когда вы заказываете у нее консенсусное дерево? Принимает на входе матрицу с 1 … n столбцами, в которых хранится частотность для слов. Потом отбирает первые сколько-то слов (скажем, сотню или сколько скажете), считает расстояние, строит на основе матрицы расстояний дерево, складывает его в корзинку. Потом берет следующую сотню слов, считает расстояние, строит дерево, складывает в корзинку… Ну вы поняли. Получается целый лес.
Звучит как итерация, а такие задачи в R решаются при помощи цикла for или пакета purrr. Функции map() из пакета purrr надо вручить другую функцию – у нас это будет пользовательская get_tree(). Она берет случайные 100 столбцов в таблице с частотностями galbraith из пакета Stylо, считает евклидово расстояние между документами и строит дерево.
Запускаем ее 13 раз при помощи map(), получаем список из тринадцати деревьев. Если хочется на них посмотреть по отдельности, то функцией walk() печатаем сразу несколько деревьев одной строчкой кода 🧙🪄
library(purrr)set.seed(123)trees_result <-map(1:13, ~get_tree(data))# отдельные деревьяpar(mfrow =c(2, 2), mar =c(1,1,1,1))walk(trees_result[1:4], plot)
Так можно построить и 100, и 1000 деревьев. Но сравнивать их вручную мы не будем, а вместо этого посчитаем консенсус. Но сначала разберемся, что это такое.
17.4plot.phylo(): простой пример
Допустим, у нас есть три дерева. Создадим их с использованием формата Ньюика, т.е. просто-напросто комбинации скобок и запятых.
Кластеры 1-2, 3-4 встречаются в двух деревьях, остальные лишь в одном. Задача — найти наиболее устойчивые кластеры методом простого большинства. Для этого считаем консенсус, причем аргумент p указывает, что кластер должен быть представлен не менее, чем в половине деревьев. Также уточняем, что наши деревья укоренены.
cons <-consensus(list(tr1, tr2, tr3), p =0.5, rooted =TRUE)
Значение p не может быть меньше 0.5, потому что конфликтующие сплиты не могут быть представлены вместе в одном дереве.
Теперь изобразим консенсус в виде дерева; дополнительно для узлов укажем силу консенсуса (2/3 = 0.67):
В таких случаях на помощь приходит консенсусная сеть. Построим сеть с использованием пакета phangorn. На входе отдаем объект класса multiPhylo, это по сути просто три дерева в одном “букете”.
Объект cons.nw относится к классу networx. Его можно изобразить как в двух, так и в трех измерениях. Вот 2D.
set.seed(16092024)par(mar =c(0,0,0,0))plot(cons.nw, type ="2D")
А вот 3D.
library(rgl) plot(cons.nw, "3D")# create animated gif file movie3d(spin3d(axis=c(0,1,0), rpm=3), duration=10, dir =".", type ="gif")
Теперь попробуем понять, что это значит (иллюстрация и объяснение отсюда).
Рассмотрим неукорененные деревья в середине: их внутренние ветви определяют расщепления (splits), а именно 12|34, 13|24 и 14|23, которые явно не могут наблюдаться в одном дереве и, следовательно, все они несовместимы. Сеть в левом верхнем углу представляет одновременно два первых дерева с прямоугольником, символизирующим две внутренние ветви. Чтобы представить все три расщепления, нам нужен куб, как показано справа.
17.7consensusNet(): galbraith
Итак, у нас есть 13 деревьев для данных galbraith.
# рисуемset.seed(16092024)par(mar =c(0,0,0,0))plot(cons.nw, type ="2D", direction ="axial",tip.color = cons.nw$col,edge.color ="grey30",edge.width =1,)
Таким образом, consensusNet() строит консенсусную сеть на основе набора деревьев: это позволяет визуализировать степень поддержки различных связей, найденных в наборе деревьев. Подход полезен для выявления областей неопределенности в филогенетических отношениях, когда несколько разных деревьев одинаково хорошо соответствуют данным.
17.8neighborNet(): простой пример
Еще один алгоритм для визуализации неоднозначных филогенетических отношений в R назвается neighborNet. Он подходит для тех случаев, когда мы подозреваем нарушения в древовидной структуре (в генетике это может быть горизонтальный перенос генов, а в литературе – например, отношения подражания или т.п.).
Ключевое различие по сравнению с consensusNet() заключается в том, что neighborNet() строит сеть непосредственно из данных, а не на основе набора деревьев.
NeighborNet работает в два шага:
Сначала строит круговую раскладку для таксонов таким образом, чтобы минимизировать расстояния между парами кластеров, каждый из которых включает в себя 1 или 2 таксона.
Потом считает веса для сплитов. На этом этапе некоторые ребра удаляются, а другие вытягиваются сообразно весам. Чем длиннее ребро, тем больше вес сплита.
Рассмотрим это на простом примере. Представьте, что у нас есть следующая матрица расстояний.
Аналогичным образом считается длина вертикальных ребер. Формула сработает максимум для четырех таксонов, для более сложных структур понадобится метод наименьших квадратов. Все вычисления делает функция neighborNet из пакета phangorn.
Если аргументу edge.label оставить значение по умолчанию, то на картинке увидите номер сплита.
У каждого сплита есть свой вес (рассчитанный методом наименьших квадратов). Его можно достать из объекта nnet.
w =attr(nnet$splits, "weights")w
[1] 0.01 0.03 0.03 0.02 0.06 0.02
Это можно понять так: чтобы попасть из пунка b в пункт d, нужно сложить веса для сплитов 4, 5 и 1:
w[4] + w[5] + w[1]
[1] 0.09
Это вернет нам 0.09. Сверяемся с матрицей расстояний — все верно!
Сплит — это разбиение совокупности таксонов на два непустых множества. Посмотрим, какие сплиты возможны для 4 таксонов из нашего примера:
as.matrix(nnet$splits)
d c a b
[1,] 1 0 0 0
[2,] 0 1 0 0
[3,] 0 0 1 0
[4,] 0 0 0 1
[5,] 1 1 0 0
[6,] 1 0 0 1
Первые четыре сплита довольно заурядны: мы просто откусываем по одному углу от нашего прямоугольника. Пятый сплит делит прямоугольник поперек, а шестой — вдоль. Дальше алгоритм для каждого сплита считает, какие пары таксонов оказались с разных сторон сплита. На основе матрицы сплитов А и исходной матрицы расстояний D рассчитывается длина ребра таким образом, чтобы кратчайшие пути между таксонами были максимально приближены к исходной матрице расстояний.
Как уже говорилось, для 4-х таксонов соответствие может быть полным. Это легко проверить, достав атрибут RSS (Residual Sum of Squares, остаточная сумма квадратов) из объекта nnet, который мы создали.
В статье “Untangling Our Past: Languages, Trees, Splits and Networks” создатели алгоритма NeighborNet объясняют, как правильно интерпретировать подобный граф на примере дерева германских языков.
Конфликтующие сигналы передаются за счет “ретикулярной структуры” (квадратиков, проще говоря). Там, где конфликта нет, мы видим дерево.
Каждый сплит представлен несколькими параллельными линиями, и если эти параллели удалить, то граф распадется на два. Чем длиннее ребро, тем “весомее” сплит.
На графе видно, что креольский язык сранан-тонго обладает сходством и с английским, и с нидерландским (граф можно разрезать по зеленым линиям двояко).
Более слабый конфликтующий сигнал прослеживается между немецким, нидерландским и фламандским, с одной стороны, и пенсильванским немецким, с другой (синий разрез).
Рассмотренные филогенетические методы (консенсусные сети, консенсусные деревья и neighborNet) ничего не говорят о происхождении одного текста от другого. Филограмма, полученная дистанционными методами, не отражает эволюционный процесс, а показывает степень дивергенции таксонов.
Это значит, что модель NeighborNet не делает никаких допущений о происхождении, однако в каком-то смысле она вполне способна показывать то, что называют «конфликтующими сигналами». В биологии это рекомбинация, гибридизация и т.п., а в гуманитарных науках — жанровые и диалектные особенности, отношения подражания, заимствования и все то, что способно влиять на результат классификации текстов, помимо авторства.
17.10 Сетевой анализ в stylo
Пакет styloне создает сетей как таковых, однако он генерирует таблицы ребер/узлов (или только ребер), используя два алгоритма Эдера для установления связей между узлами. Таблицу можно загрузить в Gephi (https://gephi.org) или прочитать в R (что мы сделаем дальше). Чтобы получить такую таблицу, вызовите функцию stylo() с аргументом network=TRUE и, по желанию, с некоторыми другими аргументами.
library(ggraph)# нормализация весовweights <- (E(galbraith_graph)$Weight -min(E(galbraith_graph)$Weight)) / (max(E(galbraith_graph)$Weight) -min(E(galbraith_graph)$Weight))E(galbraith_graph)$Weight <- weights# атрибут с именем автораlabels =str_remove(V(galbraith_graph)$name, "_.+$")V(galbraith_graph)$label <- labels# графggraph(galbraith_graph, layout ="kk") +geom_edge_link(aes(alpha = Weight), linewidth =1.1,show.legend =FALSE, color ="grey70") +geom_node_point(aes(color = label),size =3, shape =21, fill ="white", show.legend =FALSE) +geom_node_label(aes(label = name, color = label), vjust =-1, cex =2,show.legend =FALSE) +labs(x =NULL, y =NULL) +theme_void()
У нас получились три не связанные между собой подсети. О том, как работать с такого рода объектами в R, и как их интерпретировать, мы поговорим уже в следующий раз 🧙♂️.
