Статистическая среда R

Самым универсальным программным средством для анализа биоразнообразия является среда R^[1]. R — это язык программирования для статистической обработки данных и работы с графикой, свободная программная среда. В R реализованы самые современные методы анализа данных. Эта бесплатная программа работает на различных операционных системах (Windows, MacOS, Linux).

Основная трудность в использовании R — это строковый командный интерфейс. Пользователь не выбирает команды из меню, а вводит их с клавиатуры. Профессионалы считают, что ввод

команд с клавиатуры способствует осмысленному применению статистических методов, вдумчивому анализу данных. Для этого необходимо запомнить некоторое количество команд — функций языка R. Названия команд созданы на основе английского языка и в основном интуитивно понятны. Примеры: plot(), read.table(), summary(). Синтаксис данного языка достаточно логичен. Так, в скобках () заключен объект, над которым выполняется действие. Знак < — означает присваивание. Последовательность команд, необходимая для того или иного анализа, записывается в текстовый файл-скрипт. Фактически скрипт — это программа, которая загружается в среду R и автоматически выполняется в ней.

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

Самый простой способ ввода данных — это импорт из *.csv-файла, который предварительно нетрудно создать в MS Excel.

Помимо статистической обработки, R может выполнять самые разные операции с таблицами, текстами, картами и другими объектами. В целом по затратам времени и компьютерных ресурсов командный интерфейс выгоднее кнопочного, ведь записанный скрипт можно воспроизводить многократно, меняя только исходные данные. Однако для создания красивых рисунков в R надо затратить большие усилия, чем, например, в STATISTIC А.

Существуют специальные программы, делающие работу с R очень удобной. Например, RStudio позволяет разделить окно на несколько частей и одновременно просматривать исходные данные, выполненные команды (консоль), полученные в результате их выполнения рисунки, а также обеспечивает подсветку синтаксиса (подсказывает правильное написание команды при вводе первых букв) и т.д. (рис. 18).

File Edit

View Workspace Plots Tools Help

Рис. 18. Анализ данных в R с помощью RStudio

Workspace History

t^S 3 To Console To Source О % Q

гаатіг<аагагдtzyJ) n

radfit(datat4[3,])

radf it(datat4 [4 ,])

radfit(datat4[5,])

radfit(datat4[6,])

plot(rad.preempt(datat4[6, ]))

plot (rad. zipfbrot(datat4[4, ])) oat f д Г c *1

Для оценки биоразнообразия предназначены пакеты Vegan¹¹⁰ и BiodiversityR. Приведем пример вычислений с использованием первого из них (построчное описание скрипта).

• 110 ^{Oksanen R.} Multivariate Analysis of Ecological Communities in R: vegan tutorial.
• 2011. URL: http://cc.oulu.fi/~jarioksa/opetus/metodi/vegantutor.pdf

Пример анализа сообществ в пакете Vegan

Подготовка данных.

Предварительно в программе MS Excel создаем файл, в котором строки соответствуют сообществам, а столбцы — видам, сохраняем файл в формате csv. В программе:

Определяем рабочую директорию (в какой папке находятся наши данные).

Создаем объект data. Для этого импортируем данные из файла data 12.csv, который мы предварительно создали в MS Excel и поместили в рабочую папку. Названия строк row.names берем из первого столбца.

Просматриваем данные, чтобы убедиться в правильности импорта.

Привязываем данные, чтобы можно было выполнять все операции.

Загружаем пакет Vegan из библиотеки (если он был раньше установлен, в противном случае нужно воспользоваться командой установки пакета).

> setwd("D:/flaHHbie/R analyse/R mezoph")

> data<-read.csv2("datal2.csv", row.names=l)

> view(data)

> attach(data)

> UbraryC'vegan")

Вычисление показателей альфа-разнообразия.

Для каждого столбца подсчитываем число видов:

> specnumber(data)

Вычисляем индекс Шеннона (он установлен по умолчанию):

> diversity(data) и полидоминантности Симпсона: > diversity(data, "simpson")

Для других показателей можно написать более или менее простые функции. Например, показатель выравненное™ по Шеннону находится как отношение индекса Шеннона к числу видов: >3<-H/log(specnumber(data))

Статистическое оценивание индексов можно проводить методом бутстреп^[2].

Подсчитаем общее количество видов во всех сообществах: > speepool(data)

Ранговое распределение видовых обилий.

Строим и тестируем модели ранговых распределений для первого сообщества — первого столбца из нашей таблицы данных:

> radfit(data[l,])

В отличие от предыдущей программы, соответствие модели фактическим данным оценивается по информационному критерию Акаике (AIC). Наименьшее значение критерия свойственно строке Log normal , следовательно, наше сообщество ближе к логнормальной модели, чем к другим, вычисленным здесь.

Строим график рангового распределения видовых обилий:

> plot(radfit(data [1,]))

ні

Оказалось, что большое количество кривых мешает восприятию рисунка. Наверное, лучше оставить на рисунке только логнормальную модель:

> plot(rad.lognormal(data[1, ]))

Отклонение точек от линии слишком заметное, поэтому возникает потребность в графике без линий модели. Это можно задать, указав как тип графика «точки»: > plotfrad.lognormal(data[l, ]), type="p")

Наконец, можно не писать выражение для каждого столбца отдельно, а вычислить сразу для всей таблицы (аргумент 1 означает, что функцию radfit нужно применить к строкам таблицы): > apply(data, 1, radfit)

Оценка сходства сообществ.

При помощи функции vegdist вычисляем матрицу дистанций на основе коэффициента Жаккара с учетом обилия (если бы обилие учитывать было не нужно, мы бы написали bi nary=TRUE): > dis<-vegdist(data, method="jaccard", binary=FALSE)

Отображаем полученную матрицу: > dis

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

> summary(dis)

Выясним, зависит ли степень сходства сообществ от типа местообитания, при помощи рандомизации. Предварительно добавим переменную, обозначающую тип биотопа, импортировав для этого таблицу с такими же строками, в которой есть столбец «Тип»:

> env<-read.csv2("dataenv.csv", row.names=l)

> attach(env)

> mod<-betadisper(dis, Тип)

> permutest(mod)

Permutation test for homogeneity of multivariate dispersions

No. of permutations: 999 **** STRATA ****

Permutations are unstratified **** SAMPLES ****

Permutation type: free

Mirrored permutations for Samples?: No

Response: Distances

Df Sum Sq Mean Sq F N.Perm Pr(>F) Groups 3 0.047385 0.015795 4.9251 999 0.006 **

Residuals 34 0.109039 0.003207

Signif. codes: 0 '***’ 0.001 ***’ 0.01 **’ 0.05 0.1 ‘ ’ 1

Значит, различия между группами не случайны. Чтобы оценить, между какими именно типами местообитаний существуют различия в степени сходства сообществ, используем процедуру множественных сравнений Тьюки:

> TukeyHSD(mod)

Tukey multiple comparisons of means 95% family-wise confidence level

Fit: aov(formula = distances ~ group, data = df) $group

diff Iwr upr p adj лес-двор -0.008158986 -0.086521563 0.0702035922 0.9920967 луг-двор 0.083333948 -0.003859879 0.1705277755 0.0653247 сад-двор 0.008573326 -0.076737043 0.0938836952 0.9928777 луг-лес 0.091492934 0.025264468 0.1577214000 0.0037040 сад-лес 0.016732312 -0.046996059 0.0804606829 0.8927870 сад-луг -0.074760622 -0.149080035 -0.0004412089 0.0482071

Итак, для двух пар из трех с участием луга получены значимые результаты, для пары «луг и двор» лишь ненамного превышает принятое в статистике критическое значение 0,05. Уровень значимости несколько превышает критический. Но можно считать, что луг отличается более дифференцированным животным населением, чем иные типы биотопов.

Проведем кластерный анализ сообществ:

Создаем объект «кластер» hell на основе полученной матрицы с объединением по методу полной связи (дальнего соседа).

Выводим дендрограмму.

Статистическую оценку полученной классификации дадим при помощи однофакторного дисперсионного анализа.

Создаем новую переменную «номер кластера» псі при помощи функции cutree (буквально «обрезаем дендрограмму на уровне четырех кластеров»).

Запускаем (aov): изучаем влияние кластерной принадлежности на каждую переменную, которую мы получаем как столбец таблицы данных. Можно видеть, что есть статистически значимые различия между кластерами по обилию таксонов Lumbricidae и Isopoda (вывод результатов анализа для других таксонов мы опускаем):

> hcll<-hclust (dis, method="complete")

> plot(hcll)

нажмите <Ввод>, чтобы увидеть следующий график:rect.hclust(hcl1, 4)

> ncl<-cutree(hcll, 4)

> summary(aov(as.matrix(data)~ncl, data))

Проведем ординацию изученных сообществ методом многомерного неметрического шкалирования. Предварительно подключаем необходимый для этого пакет:

> library("MASS")

> іnvMDS<-isoMDS(dis)

Строим ординационную диаграмму при помощи функции:

> figl<-ordiplot(invMDS)

и затем подписываем названия сообществ:

> text(figl, "sites").

• [1] См.: Шипунов А. Б., Балдин Е. М. Анализ данных с R // Linux Format" январь-апрель 2008. URL: http://www.inp.nsk.su/~baldin/DataAnalysis/; Габидуллин Р. Статистика в R. Начало. URL: http://voliadis.ru/R-intro; Мастицкий С. R: Анализ и визуализация данных. URL: http://r-analytics.blogspot.com/; Шипунов А. Б. R — объектно-ориентированная статистическая среда. URL: http://herba.msu.rU/shipunov/software/r/r-ru.htm. 2 The R Project for Statistical Computing. URL: http://www.r-project.org/.
• [2] Скрипт см.: Шитиков В. К., Розенберг Г. С. Указ. соч. С. 36-38.