Базовые графические возможности R: категоризованные графики

Функция coplot(), входящая в базовую "комплектацию" R, предназначена для построения т.н. "conditioning plots". Это название можно перевести на русский язык как "категоризованные графики". Суть здесь сводится к тому, что анализируемые данные разбиваются на отдельные категории (например, в соответствии с уровнями какого-то фактора), для каждой из которых строится свой график (= панель) определенного типа. Все эти графики затем объединяются на одном рисунке, что существенно облегчает выявление статистических закономерностей и структур в данных (подробности можно найти в книге Cleveland (1993) Visualizing Data. AT&T Bell Laboratories).

Функция coplot() обладает большим количеством управляющих аргументов (см. ?coplot), основными из которых являются следующие:

• formula - формула, описывающая взаимодействие между анализируемыми переменными. Формула вида y ~ x | a  показывает, что графики зависимости у от x должны быть построены для каждого уровня пемеренной а. В свою очередь формула вида y ~ x| a * b показывает, что графики зависимости у от x должны быть одновременно построены для каждого уровня как переменной а, так и переменной b. Все три или четыре переменные могут быть как количественными, так и качественными (факторами). Если x или y являются факторами, то их уровни будут автоматически преобразованы в численные значения (при помощи функции as.numeric()).
• data - таблица данных, содержащая значения переменных, указанных в formula. 
• panel - функция, позволяющая задать тип и настроить внешний вид отдельных панелей категоризованного графика. По умолчанию эти панели представляют собой диаграммы рассеяния (см. примеры ниже).
• rows и columns - Панели категоризованного графика размещаются в виде матрицы с количеством строк и столбцов, соответствующим значениям аргументов rows и columns. Изменяя значения этих аргументов, можно упорядочить панели на графике необходимым образом.
• show.given - логическое значение (или вектор из двух логических значений в случае двух категориальных переменных, по которым разбиваются данные), позволяющее включать (TRUE) или отключать (FALSE) изображение "вывески" графика (см. ниже).
• number - количество интервалов, на которые разбиваются переменные а и b в случаях, когда они (эти переменные) не являются факторами.
• overlap - число (< 1), определяющее область перекрытия между данными, которые группируются в соответствии с уровнями количественных переменных a и b. Если overlap < 0, соответсвующая доля наблюдений на "стыках" групп не будет изображаться.

Кроме перечисленных аргументов, функция coplot() принимает также такие стандартные графические параметры, как col, pch, xlim, ylim, и др. (см. здесь и здесь).

Работу функции coplot() я проиллюстрирую на примере собственных данных о плотности популяции (экз./м²) моллюска Dreissena polymorpha в озере Нарочь, Беларусь (подробнее см. в оригинальной статье Mastitsky and Veres 2010). Учеты плотности популяции дрейссены были выполнены в мае и октябре 2005 г. на пяти глубинах озера вдоль восьми трансект. В приведенном ниже примере для простоты использована лишь часть данных (для трех из восьми трансект). Загрузить текстовый файл с данными (из моего Dropbox'а) можно при помощи простой команды:

density <- read.delim(file =
    "http://dl.dropbox.com/u/7521662/Dreissena_in_Naroch_Lake.txt",
    header = TRUE)

Созданная указанным образом таблица density содержит такие переменные, как Transect (Трансекта), Month (Месяц), Depth (Глубина) и собственно Density (Плотность популяции):

str(density)
'data.frame': 79 obs. of  4 variables:
 $ Transect: Factor w/ 3 levels "A","B","C": 1 1 1 1 1 1 1 1 1 1 ...
 $ Month   : Factor w/ 2 levels "May","October": 1 1 1 1 1 1 1 1 1 1 ...
 $ Depth   : num  0.8 0.8 0.8 0.8 0.8 0.8 0.8 0.8 2 2 ...
 $ Density : int  36 340 40 40 24 56 28 0 1560 40 ...

Для начала изобразим зависимость плотности популяции дрейсены от глубины для каждого месяца:

coplot(Density ~ Depth | Month,
       xlab = c("Depth", "Month"), ylab = "Density", data = density)

Воспользовавшись аргументами pch и col, мы можем изменить внешний вид и цвет точек в каждой панели графика:

coplot(Density ~ Depth | Month,
       pch = 19, col = "blue",
       xlab = c("Depth", "Month"), ylab = "Density", data = density)

Прямоугольную область графика с названиями уровней фактора, по которым разбиваются данные, по-английски иногда называют "shingle" ("вывеска"). При необходимости можно изменить заливку прямоугольников, в которые вписаны названия уровней фактора (см. аргумент bar.bg):

coplot(Density ~ Depth | Month,
       bar.bg = c(fac = "coral"), # fac - сокращение от "factor"
       pch = 19, col = "blue",
       xlab = c("Depth", "Month"), ylab = "Density",
       data = density)

Как это часто бывает с данными, характеризующими численность того или иного организма, данные по плотности популяции дрейссены в озере Нарочь распределены резко асимметрично - большинство наблюдений сосредоточено в области малых значений плотности, и лишь некоторые наблюдения имеют очень высокие значения. Как результат, плотное расположение точек в области низких значений затрудняет вывод о характере анализируемой связи. Для лучшей визуализации этой зависимости мы можем добавить сглаживающие кривые в каждую панель графика. Для этого воспользуемся аргументом panel, который, как отмечалось выше, позволяет позволяет выполнять тонкую настройку внешнего вида панелей путем прописывания соответствующей графической функции (здесь panel.smooth()):

coplot(Density ~ Depth | Month,
       panel = function(x, y, ...) {
          panel.smooth(x, y, lty = 2, pch = 19, col = "blue", span = 0.6)},
          # span регулирует кривизну сглаживающей кривой 
       bar.bg = c(fac = "coral"),
       xlab = c("Depth", "Month"), ylab = "Density",
       data = density)

Пока закономерность вырисовывается нечетко - и все из-за уже упомянутой асимметричности распределения значений плотности популяции дрейссены. Мы можем несколько снизить эту ассиметричность, преобразовав данные путем логарифмирования. Поскольку некоторые значения плотности равны нулю, перед логарифмированием добавим по единице к каждому наблюдению:

coplot(log(Density + 1) ~ Depth | Month,
       panel = function(x, y, ...) {
          panel.smooth(x, y, lty = 2, pch = 19, col = "blue", span = 0.6)},
       bar.bg = c(fac = "coral"),
       xlab = c("Depth", "Month"), ylab = "Density",
       data = density)

На полученном графике гораздо лучше прослеживается тенденция к возрастанию плотности популяции дрейссены от глубины 0.8 к 2 м и к последующей стабилизации на более глубоких участках озера. При этом данная тенденция имела место как в мае, так и в октябре.

Вспомним, что учёты плотности популяции дрейссены проводились вдоль нескольких трансект. При помощи функции coplot() мы можем легко изобразить полученные данные в отношении не только времени отбора проб, но и трансект, вдоль которых эти пробы отбирались (обратите внимание на формулу log(Density + 1) ~ Depth | Month * Transect):

coplot(log(Density + 1) ~ Depth | Month * Transect,
       panel = function(x, y, ...) {
           panel.smooth(x, y, lty = 2, pch = 19, col = "blue", span = 0.6)},
       bar.bg = c(fac = "coral"),
       xlab = c("Depth", "Month"), ylab = c("Density", "Transect"),
       data = density)

Полученный категоризованный график дает весьма полную пространственно-временную характеристику изучаемой зависимости.

Как отмечалось выше, в формуле вида y ~ x | a * b  любая из переменных может быть количественной, в том числе и переменные, стоящие по правую сторону от |. Например, мы можем разбить наши данные не по месяцам, а по интервалам глубин:

coplot(log(Density+1) ~ Month | Depth,
       bar.bg = c(num = "coral"), # num - сокращение от numeric
       pch = 19, col = "blue",
       xlab = c("Month", "Depth"), ylab = c("Density", "Transect"),
       data = density)

Важно понимать, как на полученном графике располагаются отдельные панели: они выстроены в порядке возрастания значений плотности популяции дрейссены. Первой является левая нижняя панель, за которой следуют средняя нижняя и правая нижняя панели. Во втором верхнем ряду порядок аналогичный (т.е. "читать" график нужно слева направо). "Вывеска" в верхней части графика является своего рода легендой, подсказывающей, какие именно данные изображены в каждой панели. Так, в левом нижнем углу представлены данные для диапазона глубин от 0 до примерно 1.5 м; средняя нижняя панель содержит данные для диапазона глубин от 0 до примерно 2.5 м, и т.д. Практически все автоматически выделенные программой диапазоны глубин перекрываются. Количество диапазонов и уровень их перерекрытия регулируются при помощи аргументов number и overlap соответственно (см. выше).

coplot(log(Density+1) ~ Month | Depth,
       number = 2, overlap = 0.1,
       bar.bg = c(num = "coral"),
       pch = 19, col = "blue",
       xlab = c("Month", "Depth"), ylab = c("Density", "Transect"),
       data = density)