WWW.BOOK.LIB-I.RU
БЕСПЛАТНАЯ  ИНТЕРНЕТ  БИБЛИОТЕКА - Электронные ресурсы
 


Pages:     | 1 |   ...   | 2 | 3 || 5 |

«В.К. Шитиков, С.Э. Мастицкий Классификация, регрессия и другие алгоритмы Data Mining с использованием R Тольятти, Лондон - 2017 Шитиков В.К., Мастицкий С.Э. (2017) Классификация, ...»

-- [ Страница 4 ] --

Настройка сети представляет достаточно сложную задачу: необходимо подобрать количество промежуточных слоев и число нейронов в каждом из них, установить типы активационной функции и оценить коэффициенты wi, ci каждого нейрона, которые минимизировали бы ошибку прогноза, выдаваемого сетью. Нахождение глобального оптимума аналитическими методами тут невозможно, но разработаны вполне эффективные методы исследования поверхностей ошибок и движения по градиенту (Горбань и др., 1998). Выбор нейросети "правильной" сложности, как и в рассмотренных ранее случаях, сводится к двум рецептам: использование контрольных выборок и экспериментирование.

Для обучения ИНС в среде R используются два пакета – neuralnet и nnet, – обеспечивающие гибкие функциональные возможности построения моделей классификации и регрессии на основе многослойного персептрона.

Вернемся снова к данным по ирисам Фишера и разделим ее на обучающую и тестовую выборки в соотношении 7:3, а также добавим в исходную таблицу три столбца, содержащих значения TRUE/FALSE для каждого вида растений:

data(iris) ind = sample(2, nrow(iris), replace = TRUE, prob=c(0.7, 0.3)) trainset = iris[ind == 1,] testset = iris[ind == 2,] trainset$setosa = trainset$Species == "setosa" trainset$virginica = trainset$Species == "virginica" trainset$versicolor = trainset$Species == "versicolor" Построим сеть с тремя слоями, содержащими три нейрона в скрытом слое:

library(neuralnet) net.iris = neuralnet(versicolor + virginica + setosa ~ Sepal.

Length + Sepal.Width + Petal.Length + Petal.Width, trainset, hidden=3) net.iris$result.matrix error 2.753538364739 reached.threshold 0.009945475689 steps 59246.000000000000 Из сводки итоговых результатов видно, что процесс обучения нуждался в 59246 шагах, пока максимумы всех частных производных функции ошибок reached.threshold по абсолютной величине не оказались ниже чем 0.01 (порог по умолчанию). Ошибка error оценивается как уровень правдоподобия по критерию AIC. Остальные величины, представленные в net.iris$result.matrix и net.iris$generalized.weights, содержат весовые коэффициенты w и значения сигмоид с для каждого нейрона, которые можно должным образом проанализировать и интерпретировать. Это легко сделать, получив изображение ИНС (рис.

7.12):

plot(net.iris)

–  –  –

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

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

net.prob = compute(net.iris, testset[-5])$net.result [,1] [,2] [,3] 3 -0.003883942574 0.009369070501 0.99453152167152 7 -0.017620571917 0.022468391113 0.99518743983979 11 0.019655437423 -0.010584995810 0.99095013953826 pred = c("versicolor", "virginica", "setosa") [apply(net.prob, 1, which.max)] table(Факт=testset$Species, Прогноз= pred) Acc = mean(pred == testset$Species) paste("Точность=", round(100*Acc, 2), "%", sep="")

–  –  –

8.1. Модель логита для порядковой переменной Если категории отклика являются упорядоченными, то можно воспользоваться этой дополнительной информацией и построить потенциально более качественную модель логистической регрессии с более простой интерпретацией результатов, чем модель, основанную на использовании категориальных классов отклика. Пусть для произвольной порядковой случайной величины Y, изменяющейся на интервале от 1 до J, справедливо неравенство P(Y 1) P(Y 2) … P(Y J) = 1, определяющее процесс накопления вероятности P(Y 1) = 1+...+ j, j = 1, …, J.





Тогда модель для оценки логарифма отношения накопленных шансов, или кумулятивного логита (cumulative logit), будет иметь вид:

1 +... + j logit[ P(Y j )] = log = j + x, j = 1, …, J – 1, j +1 +... + J где параметр j определяет величину, на которую увеличивается логарифм отношения шансов при включении вероятности j, а коэффициенты линейной модели не зависят от номера группы j и отражают эффекты воздействия независимых переменных. На рис. 8.1 показан пример такой модели с одинаковым эффектом независимой переменной х на каждую из трех функций накопленных вероятностей отклика с четырьмя категориями: видно, что происходит пропорциональный сдвиг кривых вправо, определяемый величиной коэффициента j.

–  –  –

Рис. 8.2. Корреляционная матрица морфометрических показателей Рис. 8.3. Зависимость числа колец от длины раковины Разделим интервал варьирования rings на диапазоны и выделим четыре возрастные группы. С точки зрения теории информации назначение границ предпочтительнее всего осуществить из условия равной численности групп. Рассмотрим предварительно функцию плотности распределения и границы квартилей для возможного разделения (рис.

8.4):

ggplot(abalone) + aes(rings, fill = пол) + geom_density(position = "stack")+ geom_vline(xintercept = quantile(abalone$rings, p = c(0.25, 0.5, 0.75)),colour = "blue", linetype = 5, size = 1.5) table (cut(abalone$rings,breaks=quantile(abalone$rings, c(0,.25,.50,.75,1), include.lowest=TRUE))) [1,8] (8,9] (9,11] (11,29] Рис. 8.4. Плотность распределения числа колец раковины и границы квартилей Однако из практических соображений мы несколько скорректируем диапазоны, обозначенные квартилями, переместив объекты с 8 возрастными кольцами в группу более старых особей.

Добавим новый столбец Возраст с четырьмя возрастными категориями в исходную таблицу abalone и сохраним ее в преобразованном виде в файле abalone.RData:

abalone$Возраст - cut(abalone$rings, breaks=c(0,7,9,11,29), labels=c("Q1","Q2","Q3","Q4"), include.lowest=TRUE) save(abalone, file="abalone.RData") table (abalone$Возраст ) Q1 Q2 Q3 Q4 Сама по себе идея перейти от счетной переменной к категориальной переменной-фактору может показаться дискуссионной, но она дает нам возможность приступить к построению моделей с порядковым откликом.

Модель кумулятивного логита получим с использованием функции polr() из пакета MASS.

library(MASS) CL.aq - polr(Возраст ~., data = abalone[,-9]) summary(CL.aq, digits = 3)

Coefficients:

Value Std. Error t value полI -0.9424 0.0934 -10.090 полM 0.0542 0.0748 0.724 длина 1.5643 1.6722 0.935 диаметр 9.2270 2.0472 4.507 высота 13.4736 2.1206 6.354 вес.общ 6.3932 0.7897 8.095 вес.тела -15.3736 0.9126 -16.846 вес.внут -7.1420 1.2760 -5.597 вес.рак 7.3099 1.2157 6.013

Intercepts:

Value Std. Error t value Q1|Q2 4.056 0.326 12.447 Q2|Q3 6.540 0.341 19.206 Q3|Q4 8.487 0.344 24.644 Residual Deviance: 8224.074 AIC: 8248.074 # Оценка доверительных интервалов коэффициентов confint.default(CL.aq) 2.5 % 97.5 % полI -1.12547618 -0.7593602 полM -0.09236741 0.2006869 длина -1.71315347 4.8417381 диаметр 5.21461784 13.2394216 высота 9.31738606 17.6298645 вес.общ 4.84535656 7.9410052 вес.тела -17.16223162 -13.5850167 вес.внут -9.64297569 -4.6411202 вес.рак 4.92703181 9.6926746

Заметим, что блок коэффициентов состоит из двух разделов:

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

Выше отмечалась высокая мультиколлинеарность данных, т.е. высокая степень взаимной зависимости морфометрических показателей. Считается, что наиболее эффективный путь устранения мультиколлинеарности – исключение из регрессионной модели малозначимых коэффициентов, или, выражаясь точнее, отбор информативного комплекса из q переменных (q m). Пошаговый регрессионный анализ, выполняемый функцией stepAIC(), представляет собой последовательную процедуру включения и исключения отдельных предикторов в модель, пока не будет достигнута оптимальная регрессия по критерию Акаике.

В результате проведенной пошаговой процедуры был исключен показатель длина, а AIC-критерий уменьшился с 8248 до 8247. Однако последующие итерации не привели к его дальнейшей оптимизации (мы опускаем итоги расчетов, поскольку они незначительно отличаются от приведенных выше).

Тест на адекватность полученной модели в целом проведем, сравнивая отношения правдоподобия нуль-модели без параметров (1) и модели CLs.aq (2):

CLs.aq - stepAIC (CL.aq) summary(CLs.aq, digits = 3) confint.default(CLs.aq) CL0.aq - polr(Возраст ~ 1, data = abalone[,-9]) anova(CL0.aq, CLs.aq) Likelihood ratio tests of ordinal regression models Model Resid. df Resid. Dev Test Df LR stat. Pr(Chi) 1 4174 11484.659 2 4166 8224.946 1 vs 2 8 3259.712 0 Модель имеет высокую статистическую значимость, что практически всегда бывает, если число наблюдений превышает несколько тысяч.

Рассмотрим теперь эффективность использования построенной модели для прогнозирования возраста моллюсков. С помощью функции fitted() найдем значения вероятностей для каждой из 4 рассматриваемых градаций, которые вычисляются на основе оценок коэффициентов j и для модели CLs.aq. Отнесение каждого объекта выборки к конкретной возрастной категории будем осуществлять по максимальной вероятности. Как всегда, функция table() поможет нам выделить несовпадение действительных и предсказанных значений возрастных категорий.

Probs - fitted(CLs.aq) head(Probs, 3) Q1 Q2 Q3 Q4 1 0.16982105 0.5401767 0.23502603 0.054976171 2 0.42821924 0.4714134 0.08472626 0.015641138 3 0.03303907 0.2572052 0.45143960 0.258316097 # Построение матрицы неточностей "Факт-Прогноз" pred=apply(Probs,1,function(x) colnames(Probs)[which(x==max(x))]) table(Факт=abalone$Возраст,Прогноз=pred) Acc = mean(pred == abalone$Возраст) paste("Точность=", round(100*Acc, 2), "%", sep="") Прогноз Факт Q1 Q2 Q3 Q4 Q1 589 233 14 3 Q2 113 773 306 65 Q3 22 379 446 274 Q4 6 144 292 518 [1] "Точность=55.69%" Можно отметить, что модель недостаточно хорошо справилась с предсказанием смежных возрастных категорий.

Многомерность параметров модели не позволяет показать графически изменение вероятности классов в зависимости от значений всех предикторов одновременно. Распространено мнение, что это можно сделать для отдельных предикторов, если зафиксировать значения остальных независимых переменных на уровне их средних значений. Выполним такой частный анализ, например, для диаметра раковины (см. рис.

8.5):

# Подготовка данных для графика VM - apply(abalone[,2:8],2, mean) d.plot - as.data.frame(matrix(VM, ncol=7, nrow=50, byrow=TRUE, dimnames = list(1:50,names(VM)))) d.plot - cbind(пол=rep("F",50), d.plot) d.plot$диаметр - seq(min(abalone$диаметр), max(abalone$диаметр), len = 50) d.pplot - cbind(d.plot, predict(CL.aq, newdata = d.plot, type = "probs", se = TRUE)) # Прорисовка компонент графика plot (1,1, xlim=c(0,max(abalone$диаметр)),ylim=c(0,0.8), type='n', xlab="Диаметр раковины", ylab="Вероятность Р") lines(d.pplot[,c(3,9)],lwd=2, col="green") lines(d.pplot[,c(3,10)],lwd=2, col="blue") lines(d.pplot[,c(3,11)],lwd=2, col="black") lines(d.pplot[,c(3,12)],lwd=2, col="red") legend("topright", c("Q1","Q2","Q3","Q4"), lwd=2, col=c(3,4,1,2)) Рис. 8.5. Вероятности отнесения моллюсков к возрастным категориям в зависимости от диаметра раковины На рис. 8.5 видно, что экстремальные значения диаметра раковины соответствуют максимальным вероятностям отнесения к младшим и старшим возрастным категориям соответственно слева и справа. Промежуточные значения диаметра, вероятнее всего, приведут к средним возрастам Q4 и Q3.

8.2. Настройка параметров нейронных сетей средствами пакета caret В разделе 7.5 мы рассмотрели основные принципы построения искусственных нейронных сетей на примере многослойного персептрона.

Основная проблема обучения ИНС заключается в необходимости предварительно исследовать поверхность ошибок и задать архитектуру сети и параметры оптимизации, по возможности, приводящие в окрестность глобального минимума. С использованием функции train() из пакета caret путем перекрестной проверки можно оценить оптимальные значения числа скрытых нейронов size и параметр "ослабления весов" decay, который осуществляет регуляризацию точности подстройки коэффициентов (при decay = 0 стремление к точности может перерасти в эффект переусложнения модели).

Покажем, как это можно сделать, на примере оценки возрастной категории морских ушек:

library(nnet); library(caret) set.seed(123); load (file="abalone.RData") train.aba - train(Возраст ~., data = abalone[,c(3:8,10)], method = "nnet", trace = F, linout = 1, tuneGrid = expand.grid(.decay = c(0,0.05,0.2),.size = 4:9), trControl = trainControl(method = "cv")) 4177 samples 6 predictors 4 classes: 'Q1', 'Q2', 'Q3', 'Q4' Resampling: Cross-Validation (10 fold)

Resampling results across tuning parameters:

size decay Accuracy Kappa Accuracy SD Kappa SD 5 0.2 0.589 0.446 0.0211 0.0291 6 0 0.593 0.452 0.0255 0.0348 6 0.05 0.592 0.45 0.0251 0.0346 6 0.2 0.59 0.447 0.0231 0.0318 7 0 0.595 0.455 0.0287 0.0391 7 0.05 0.591 0.449 0.0245 0.034 7 0.2 0.59 0.447 0.0226 0.0312 8 0 0.59 0.448 0.0257 0.035 8 0.05 0.588 0.445 0.0256 0.0353 Accuracy was used to select optimal model using largest value.

The final values used for model were size = 7 and decay = 0.

Мы выполнили 10-кратную перекрестную проверку 18 нейросетевых моделей с числом нейронов в скрытом слое от 4 до 9 и разных значениях "ослабления". При найденных значениях size = 7 и decay = 0, приводящих к максимальной точности Accuracy, построим далее модель с помощью функции nnet(). Для визуализации сети (рис. 8.6) применим функцию из скрипта nnet_plot_update.r, которая имеет ряд полезных опций (можно скачать с https://www.r-bloggers.com).

source("nnet_plot_update.r") # plot.nnet(train.aba) - можно все извлечь из train-объекта # nn.aba - train.aba$finalModel nn.aba - nnet(Возраст ~., data = abalone[, c(3:8, 10)], decay = 0, size = 7, niter=200) plot.nnet(nn.aba) summary(nn.aba) Рис. 8.6. Персептрон для оценки возрастной категории морских ушек Мы получили сеть из 17 нейронов (6 7 4), для которых рассчитано 81 нижеприведенных весовых коэффициентов и порогов b.

a 6-7-4 network with 81 weights b-h1 i1-h1 i2-h1 i3-h1 i4-h1 i5-h1 i6-h1

-8.67 10.81 0.22 -4.89 7.14 -1.40 16.56 b-h2 i1-h2 i2-h2 i3-h2 i4-h2 i5-h2 i6-h2 3.03 -3.80 -4.03 4.59 -4.55 7.01 -14.24 b-h3 i1-h3 i2-h3 i3-h3 i4-h3 i5-h3 i6-h3

-5.70 34.35 6.02 -2.66 -1.47 1.00 -5.93 b-h4 i1-h4 i2-h4 i3-h4 i4-h4 i5-h4 i6-h4 7.87 -18.03 -0.85 18.30 -24.87 -6.31 -21.53 b-h5 i1-h5 i2-h5 i3-h5 i4-h5 i5-h5 i6-h5

-2.20 9.24 -16.03 -23.37 41.34 10.69 16.89 b-h6 i1-h6 i2-h6 i3-h6 i4-h6 i5-h6 i6-h6

-7.12 -8.78 -5.86 1.73 -5.04 18.64 5.11 b-h7 i1-h7 i2-h7 i3-h7 i4-h7 i5-h7 i6-h7

-2.62 10.64 1.02 -7.28 11.45 1.54 -13.55 b-o1 h1-o1 h2-o1 h3-o1 h4-o1 h5-o1 h6-o1 h7-o1 7.05 0.24 -0.42 -15.02 5.22 4.68 -1.05 11.34 b-o2 h1-o2 h2-o2 h3-o2 h4-o2 h5-o2 h6-o2 h7-o2 4.49 -2.53 0.54 -3.46 -1.97 1.39 -6.17 2.48 b-o3 h1-o3 h2-o3 h3-o3 h4-o3 h5-o3 h6-o3 h7-o3

-3.14 0.81 1.71 5.18 -2.48 -1.39 -1.04 -3.68 b-o4 h1-o4 h2-o4 h3-o4 h4-o4 h5-o4 h6-o4 h7-o4

-7.55 1.40 -1.86 12.61 -1.04 -3.73 8.48 -10.25 Рассмотрим теперь, насколько хороша полученная нейросетевая модель при выполнении предсказаний:

pred = predict(nn.aba, abalone[,3:8], type="class") nn.table = table(abalone[,10], pred) confusionMatrix(nn.table) Confusion Matrix and Statistics pred Q1 Q2 Q3 Q4 Q1 660 166 7 6 Q2 177 764 231 85 Q3 46 367 492 216 Q4 10 129 239 582 Overall Statistics Accuracy : 0.598 95% CI : (0.583, 0.613) No Information Rate : 0.3414 P-Value [Acc NIR] : 2.2e-16 Kappa : 0.4591 Mcnemar's Test P-Value : 2.33e-13

Statistics by Class:

Class: Q1 Class: Q2 Class: Q3 Class: Q4 Sensitivity 0.7391 0.5358 0.5077 0.6547 Specificity 0.9455 0.8208 0.8039 0.8850 Pos Pred Value 0.7867 0.6078 0.4389 0.6063 Neg Pred Value 0.9302 0.7733 0.8439 0.9046 Prevalence 0.2138 0.3414 0.2320 0.2128 Detection Rate 0.1580 0.1829 0.1178 0.1393 Detection Prevalence 0.2009 0.3009 0.2684 0.2298 Точность предсказания возрастных категорий несколько возросла по сравнению с кумулятивным логитом (0.598 против 0.557), даже несмотря на то, что ковариату пол мы не использовали.

Рассмотрим кратко еще две функции из пакета caret, предназначенные для работы с искусственными нейронными сетями. Серьезной проблемой при обучении ИНС может стать высокая мультиколлинеарность предикторов.

Одним из вариантов ее решения является проецирование исходной выборки в пространство главных компонент (см. метод PCR в разделе 4.2). Если задаться некоторым порогом (аргумент thresh), то значения p главных компонент, для которых накопленная доля объясненной дисперсии превышает этот порог, можно использовать как входы в нейронную сеть и настроить ее обычным способом как для регрессии, так и для классификации.

Функция pcaNNet() является своеобразной оберткой для совместного выполнения предобработки данных, анализа главных компонент, и запуска функции nnet() для обучения сети на полученных главных компонентах.

Выполним обучение сети для предсказания возрастной категории морских ушек с использованием параметра thresh = 0.975 и числа нейронов в скрытом слое size = 7:

pcaNNet.Fit - pcaNNet(abalone[,3:8], abalone[,10], size = 7, thresh = 0.975, linout = TRUE, trace = FALSE) Neural Network Model with PCA Pre-Processing PCA needed 4 components to capture 97.5 % of the variance a 4-7-4 network with 67 weights При значении thresh = 0.975 на вход сети подается 4 главных компоненты, отвечающих этому условию. Для тестируемых данных такое же преобразование (основанное на факторных нагрузках для обучающего множества) применяется к новым значениям предикторов.

pred - predict(pcaNNet.Fit,abalone[,3:8], type="class") table(Факт=abalone$Возраст,Прогноз=pred) Acc = mean(pred == abalone$Возраст) paste("Точность=", round(100*Acc, 2), "%", sep="") Прогноз Факт Q1 Q2 Q3 Q4 Q1 654 170 8 7 Q2 187 694 266 110 Q3 54 302 511 254 Q4 15 137 224 584 [1] "Точность=58.49%" В разделе 4.4 была описана процедура бэггинга, как общего метода агрегирования моделей произвольной структуры, построенных на бутстрепвыборках из исходного набора данных. Функция avNNet() осуществляет обучение заданного множества моделей нейронной сети на одном и том же наборе данных. Модели могут различаться как из-за случайного дрейфа стартовых параметров калибровки сети (Ripley, 1996), так и вследствие использования бутстреп-выборок, извлеченных из исходного обучающего множества. Для моделей классификации функция avNNet() оценивает среднее значение вероятностей классов на основе частных прогнозов каждой из моделей созданного ансамбля и далее производит заключительное предсказание класса.

Для рассматриваемого примера сформируем 10 экземпляров моделей

ИНС с использованием бэггинга:

avNNet.Fit - avNNet(Возраст ~., data = abalone[,c(3:8,10)], size = 7, repeats = 10, linout = TRUE, trace = FALSE, bag = TRUE) pred - predict(avNNet.Fit,abalone[,3:8], type="class") Acc = mean(pred == abalone$Возраст) paste("Точность=", round(100*Acc, 2), "%", sep="") Model Averaged Neural Network with 10 Repeats a 6-7-4 network with 81 weights [1] "Точность=59.48%"

8.3. Методы комплексации модельных прогнозов Существует два возможных подхода к получению некоторого обобщенного предсказания на основе формального или неформального объединения частных прогнозов. В разделе 4.4 нами подробно рассматривалось агрегирование результатов моделирования, когда каждая модель имеет фактически одну и ту же структуру, но многократно обучается на слегка модифицированных исходных данных или при небольшой вариации управляющих параметров (bagging). Другой путь – построить на одних и тех же исходных данных некоторый "коллектив" (ensemble) из принципиально разнотипных моделей и выполнить комбинирование прогнозов (forecast combinations).

По аналогии с методами коллективного решения, столь эффективно используемым в человеческом обществе, считается, что суммарная эффективность системы коллективного прогнозирования теоретически будет в среднем значительно выше любого из его членов. Этот оптимизм мало чем подтверждается на практике, но нельзя не признать, что комбинированный прогноз значительно более устойчив к не вполне объяснимым "провалам", которые свойственны отдельным индивидуальным методам, и поэтому может быть эффективным инструментом страхования на случай возможных стохастических рисков (в мире трейдинга это называется хеджированием, англ. hedging).

Пусть имеется m прогнозов y1, y2,..., ym для переменной Y, полученных с помощью m различных моделей. Под коллективным прогнозом g будем понимать прогноз той же переменной Y как некоторую функции от индивидуальных прогнозов yk: g = F(y1, y2,..., ym; X).

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

Коллективный прогноз g, построенный без адаптации, представляют в виде линейной комбинации из базового или суженного (наиболее m g = y k k, информативного) множества частных прогнозов k =1 где – вектор неизвестных весовых коэффициентов, компоненты которого неизменны на всем диапазоне варьирования исходных данных X. Тогда задача комплексации эквивалентна определению совокупности значений вектора k, удовлетворяющих заданным ограничениям и минимизирующих некоторый критерий качества. Очевидно, что качество (надежность) прогноза g тем выше, чем "ближе" расчетные значения Y к фактическим.

Обратимся вновь к примеру из предыдущих разделов, но будем теперь строить модели прогнозирования не для возрастных категорий, а непосредственно для числа колец в раковинах rings тасманских морских ушек (см. рис. 8.4). Заменим фактор пол на две метрические переменные полM и полI, после чего разделим исходную выборку на обучающую и проверочную в соотношении 2:1.

load(file="abalone.RData") # Формирование матрицы предикторов X - model.matrix(rings ~.,data=abalone[,1:9])[,-1] Y - abalone$rings # Разделение на обучающую и проверочную выборки train - runif(nrow(X)) =.66 На объектах обучающей выборки (2748 наблюдений) построим шесть моделей оценки rings с использованием функции train() из пакета

caret:

library(caret); library(party) myControl - trainControl(method='cv', number=10, savePredictions=TRUE, returnData=FALSE, verboseIter=TRUE) PP - c('center', 'scale') # Обучение избранных моделей # Регрессия на k-ближайших соседей m.knn - train(X[train,], Y[train], method='knn', trControl=myControl, preProcess=PP) # Линейная модель m.lm - train(X[train,], Y[train], method='lm', trControl=myControl, preProcess=PP) # Гребневая регрессия с регуляризацией m.rlm - train(X[train,], Y[train], method='glmnet', trControl=myControl, preProcess=PP) # Модель опорных векторов m.svm - train(X[train,], Y[train], method='svmRadial', trControl=myControl, preProcess=PP) # Метод случайного леса m.rf - train(X[train,], Y[train], method='rf', trControl=myControl) # Бэггинг деревьев условного вывода m.ctr - train(X[train,], Y[train], "bag", B = 10, bagControl = bagControl(fit = ctreeBag$fit, predict = ctreeBag$pred, aggregate = ctreeBag$aggregate)) Обратим внимание на модель m.ctr, построенную методом "bag". При построении этой модели используется специальная функция bag() из пакета caret, позволяющая выполнить бэггинг широкого класса различных моделей (daBag, plsBag, nbBag, ctreeBag, svmBag, nnetBag). Мы осуществили бэггинг деревьев условного вывода (см. раздел 4.3), агрегируя каждый раз по 10 деревьев.

Упакуем результаты моделирования в объект класса list и извлечем значения среднеквадратичной ошибки RMSE каждой индивидуальной модели на обучающей выборке:

# Создание списка всех моделей all.models - list(m.knn, m.lm, m.rlm, m.svm, m.rf, m.ctr) names(all.models) - sapply(all.models, function(x) x$method) sort(sapply(all.models, function(x) min(x$results$RMSE))) rf svmRadial glmnet lm bag knn 2.136772 2.153634 2.172855 2.178365 2.222380 2.233542 В целом все модели показали близкие результаты, хотя можно отметить традиционное лидерство методов случайного леса (rf) и опорных векторов (svmRadial).

Сформируем таблицу частных прогнозов для объектов тестовой выборки (1429 наблюдений) и оценим ошибки индивидуальных моделей

RMSE:

preds.all - data.frame(sapply(all.models, function(x){predict(x, X)})) head(preds.all) knn lm glmnet svmRadial rf bag 1 7.111111 7.732803 7.781386 7.389098 8.417767 8.031008 2 9.000000 9.614801 9.622614 9.537425 9.006100 9.451514 3 11.888889 13.043430 12.944758 12.409405 14.279467 12.646068 4 11.888889 10.927220 10.908910 11.563814 12.673133 10.957280 5 9.555556 8.813379 8.783278 9.038425 8.866833 8.803584 6 8.111111 8.360550 8.393360 8.198748 8.211833 8.467949 rmse - function(x,y){sqrt(mean((x - y)^2))} print("Среднеквадратичная ошибка на тестовой выборке") print(sort(apply(preds.all[!train,], 2, rmse, y = Y[!train])), digits= 3) [1] “Среднеквадратичная ошибка на тестовой выборке" rf svmRadial bag lm glmnet knn 2.17 2.20 2.21 2.28 2.28 2.36 В целом незначительное (2-4%) увеличение ошибки RMSE на свежих наблюдениях свидетельствует о хорошем качестве моделей.

Рассмотрим теперь возможные алгоритмы комплексации частных прогнозов, реализованные в пакете ForecastCombinations.

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

простое среднее (simple) при k = const = 1/m;

по лучшему частному прогнозу (best), т.е. k = 1 для max(yk) и k = 0 для остальных yk.

Комплексация, основанная на дисперсиях (variance based), требует оценки на некоторой обучающей выборке значений среднеквадратичных отклонений MSEк для каждой индивидуальной модели и тогда k = (1 / MSEk ) / k =1 (1 / MSEk ).

m Другим способом комплексации прогнозов разных моделей является использование линейной модели, в которой коэффициенты играют роль весов, определяющих вклад каждой из объединяемых моделей. С легкой руки букмекеров эту процедуру стали часто называть стэкингом (stacking).

Разумеется, обычный метод наименьших квадратов (ols) является одним из претендентов на создание такой "модели моделей". Однако, поскольку прогнозы частных моделей сильно коррелируют между собой, то для получения комплексного результата предлагают использовать иные методы регрессии, которые, по мнению авторов, могут лучше справиться с этой напастью. В пакет ForecastCombinations включены метод наименьших абсолютных отклонений (robust) и наименьших квадратов с ограничениями (cls - Constrained Least Squares). Во втором случае на значения коэффициентов регрессии накладываются ограничения: все они неотрицательны и в сумме равны 1, т.е. по сути отражают вклад каждой частной модели в общий прогноз. Коэффициенты CLS-модели рассчитываются с использованием метода квадратичного программирования, представленного функцией solve.QP() из пакета quadprog.

На вход функции Forecast_comb() из пакета ForecastCombinations подаются объекты для обучения (матрица fhat со значениями частных прогнозов y1, y2,..., ym и вектор obs с соответствующими наблюдаемыми величинами y) и прогноза (матрица fhat_new со значениями частных прогнозов, которые следует объединить с использованием метода Averaging_scheme). На выходе формируются векторы прогнозов pred и весов weights.

С учетом выполненного нами разбиения исходной таблицы abalone рассчитаем сравнительную оценку RMSE, полученную каждым методом комплексации на проверочной выборке:

library(ForecastCombinations) scheme= c("simple", "variance based", "ols", "robust", "cls", "best") combine_f - list() for (i in scheme) { combine_f[[i]] - Forecast_comb(obs = Y[train], fhat = as.matrix(preds.all[train, ]), fhat_new = as.matrix(preds.all[!train, ]), Averaging_scheme = i) tmp - round(sqrt(mean((combine_f[[i]]$pred Y[!train])^2)), 3) cat(i, ":", tmp, "\n") } simple : 2.16 variance based : 2.143 ols : 2.266 robust : 2.266 cls : 2.166 best : 2.166 Рассмотрим веса, которые получили каждые модели при реализации каждой схемы усреднения (приведем их предварительно к единой шкале):

w.list - sapply(combine_f, function(sp.list) sp.list$weights) w.list$ols - w.list$ols[-1] w.list$robust - w.list$robust[-1] weights - as.data.frame.list(w.list) rownames(weights) - scheme wstan - function (x) abs(x)/sum(abs(x)) weights[3,] - wstan(weights[3,]) weights[4,] - wstan(weights[4,]) print(round(weights,3)) knn lm glmnet svmRadial rf bag simple 0.167 0.167 0.167 0.167 0.167 0.167 variance based 0.111 0.094 0.094 0.105 0.489 0.108 ols 0.056 0.077 0.063 0.088 0.660 0.056 robust 0.043 0.137 0.128 0.086 0.562 0.043 cls 0.000 0.000 0.000 0.000 1.000 0.000 best 0.000 0.000 0.000 0.000 1.000 0.000 Комплексные прогнозы, созданные по схемам cls и best, эквивалентны мнению "непогрешимого" авторитета (а таковым оказался метод rf). Соответственно они дали ошибку, в точности равную погрешности наиболее качественной индивидуальной модели. Схемы ols и robust потерпели относительную неудачу, тогда как наиболее точным оказался прогноз, основанный на дисперсии (variance based).

П.

Полищук (http://www.qsar4u.com) предлагает применять метод неотрицательных наименьших квадратов (non-negative least squares), используя пакет nnls:

require(nnls) m.nnls - nnls(as.matrix(preds.all[train, ]), Y[train]) coef(m.nnls) [1] 0.000000 0.000000 0.000000 0.000000 1.008773 0.000000 Этот метод дал те же значения весов, что и CLS-регрессия (впрочем, эти методы идентичны, по-видимому, и по своей сути).

Кроме линейного стэкинга для комплексации можно применять любые изощренные алгоритмы статистического моделирования. Например, набор частных прогнозов можно выразить через небольшое число главных компонент и на последующих шагах применить любой другой алгоритм взвешивания (Горелик, Френкель, 1983)3. П. Брусиловский и Г. Розенберг (1983) использовали для комплексации многорядный алгоритм МГУА, строили самообучающиеся иерархические модели и назвали этот метод "модельным штурмом". Дело устойчиво идет к проблематике создания "коллектива коллективов".

Пакет также позволяет реализовать ForecastCombinations эксперименты подобного рода. Из 6 векторов частных прогнозов, полученных нами по индивидуальным моделям, можно построить 53 модели регрессии с различными наборами переменных от одного до 5. Функция Forecast_comb_all() выполняет настройку всех этих моделей и сводит результаты в новую таблицу (т.е. вместо 6 прогнозов предлагается иметь дело с 53-мя прогнозами на основе прогнозов).

Далее их можно просто усреднить:

combine_f_all - Forecast_comb_all(Y[train], fhat = as.matrix(preds.all[train, ]), fhat_new = as.matrix(preds.all[!train, ])) # Усредняем комбинированные прогнозы по всем регрессиям Combined_f_all_simple - apply(combine_f_all$pred, 1, mean) print(sqrt(mean((Combined_f_all_simple - Y[!train])^2)), digits= 3 ) [1] 2.16 Для моделей регрессии по каждому подмножеству отдельных моделей функция Forecast_comb_all() рассчитывает также набор информационных критериев (AIC, BIC и Mallow’s Cp), которые могут быть использованы для расчета взвешенного среднего:

# Комбинирование прогнозов с использованием Mallow's Cp:

Combined_f_all_mal t( combine_f_all$mal %*% t(combine_f_all$pred) ) print(sqrt(mean((Combined_f_all_mal - Y[!train])^2)), digits= 3) Информация о всех литературных источниках этого раздела приведена в (Розенберг и др., 1994) [1] 2.24 Итогом всех этих достаточно трудоемких вычислительных процедур оказался весьма скромный результат.

Альтернативой механическому усреднению прогнозов является учет структурных связей в ансамбле таким образом, чтобы положительные свойства той или иной модели (метода) доминировали при формировании коллективного отклика, а отрицательные – отфильтровывались (Растригин, Эренштейн, 1981). В частности, существенным влиянием должны обладать самые "непохожие" между собой модели-индивидуумы, являющиеся носителями нетривиальных тенденций моделируемого процесса. Основным направлением развития в этой области является разработка методов адаптации, в основе которых лежат алгоритмы Бейтса-Гренджера (Bates, Grander, 1965), Ньюболда-Гренджера (Newbald, Grander, 1974), Ершова (1975) и др. К сожалению, нам не удалось обнаружить пакеты R, в которых были бы реализованы вменяемые методы комплексации прогнозов с адаптацией.

8.4. Обобщенные линейные модели для счетных данных В предыдущих главах нами рассматривались примеры использования методов построения моделей для различных типов распределения отклика.

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

Этот выбор, в первом приближении, делается априори на основании имеющихся знаний относительно отклика. Например, если моделируется присутствие (1) и отсутствие (0) животных на некоторой изучаемой площади как функция нескольких факторов, то выбор прост: необходимо использовать биномиальное распределение. Однако это, вероятно, единственный сценарий, где ситуация столь очевидна.

Обычно, если это позволяет характер вариации переменной Y, прибегают к аппроксимации счетных данных тем или иным непрерывным распределением – нормальным, гамма или Вейбулла, – что мы в конце концов и сделали с числом колец в раковинах морских ушек (рис. 8.4). Хорошему результату подгонки параметров модели и гарантированной неотрицательности модельных значений может способствовать преобразование исходных данных с помощью различных функций, таких как, логарифмирование или трансформация Бокса-Кокса (см. раздел 3.3).

Однако зачастую дискретность отклика оказывается столь очевидной (yi = 0, 1, 2, …), что аппроксимация непрерывным распределением оказывается неудовлетворительной и это сильно влияет на адекватность регрессии с использованием МНК. Если счетные данные имеют тенденцию к гетерогенности (т.е. в большинстве случаев можно не встретить ни одного экземпляра, реже попадутся единичные особи, и лишь иногда их большие скопления) и фиксированный верхний предел у них отсутствует, то хорошим вариантом является аппроксимация наблюдений распределением Пуассона (Zuur et al., 2009).

Распределение Пуассона P() имеет случайная величина Y, отражающая количество событий, произошедших за некоторый промежуток времени, когда эти события независимы и происходят с постоянной интенсивностью, (0, ). Предполагается, что среднее µ = E(Y) = и дисперсия Var(Y) =, а функция плотности вероятности p(k) = e- k/k!, k = 1, 2, …,, где k! – факториал от числа событий.

Пуассоновскую регрессию применяют, когда отклик является Y счетной переменной, имеющей такое распределение. Соответствующая модель имеет вид ln( ) = 0 + j x j p где x1, …, xp – набор из p независимых переменных, 0 – математическое ожидание Y при равенстве нулю всех предикторов xj, j – коэффициенты независимых переменных.

Нередкими проблемами подгонки распределением Пуассона являются:

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

эксцесс – избыток некоторых наблюдаемых значений по сравнению с ожидаемыми частотами;

наличие большой гетерогенности данных, вследствие чего оценка дисперсии начинает значимо превышать среднее, Var(Y) µ (overdispersion).

Если наблюдается отчетливая избыточная дисперсия ("перерассеяние"), то альтернативой распределению Пуассона для счетных данных является отрицательное биномиальное распределение.

Отрицательное биномиальное распределение NB(r, p), также называемое распределением Паскаля, – это распределение дискретной случайной величины, которая отражает количество произошедших неудач в последовательности испытаний Бернулли с вероятностью успеха p, проводимой до r-го успеха (https://ru.wikipedia.org). Его обычно интерпретируют как смесь распределений гамма Г(.) и Пуассона, т.е. функция ( у + ) µ y y f ( y; µ, ) = плотности вероятности имеет вид:, ( у ) ! µ + ( y + ) где µ = E(Y) – среднее, а – параметр формы распределения.

Включение в функцию распределения дополнительного параметра позволяет учесть превышение дисперсии над средним:

Var(Y) = µ + µ2 /.

При отрицательное биномиальное распределение сводится к распределению Пуассона (Agresti, 2007; Zeileis et al., 2008).

Обобщенная линейная модель для отрицательного биномиального распределения оценивает зависимость µ от совокупности объясняющих переменных с использованием функции связи (link function), как и в случае с лог-линейной моделью Пуассона. Принято считать, что параметр формы постоянен на всей области определения модели, что аналогично предположению о гомоскедастичности для моделей с нормальным распределением остатков.

Если аналитик колеблется в выборе между двумя конкурирующими статистическими моделями, то основанием для выбора может служить оценка согласия между наблюдаемыми данными и GLM с соответствующим теоретическим распределением. Рассмотрим эту задачу на примере двухлетнего подсчета числа амфибий, раздавленных автомобилями на 52 участках дорожной сети в южной части Португалии. Файл с данными RoadKills.txt можно скачать с полным комплектом данных к книге Zuur et al. (2009) на сайте http://www.highstat.com/book2.htm.

Само по себе распределение отклика TOT.N (числа погибших животных на каждом участке, экз.) трудно принять как пуассоновское (рис.

8.7):

RK - read.table(file = "RoadKills.txt", header = TRUE, dec = ".") y - RK$TOT.N library(vcd) gf-goodfit(table(y),type= "poisson",method= "ML") # Визуализация подогнанных данных гистограммой plot(gf,ylab="Частоты", xlab="Число классов") Goodness-of-fit test for poisson distribution X^2 df P( X^2) Likelihood Ratio 970.7709 29 1.290368e-185 Рис. 8.7. Гистограмма отклонений эмпирических частот от теоретических частот распределения Пуассона (красные кружки) Модель регрессии Пуассона дает нам возможность объяснить эти отклонения за счет влияния внешних факторов.

На первом этапе оценим общий вид моделируемых зависимостей, для чего применим одну независимую переменную – D.PARK (расстояние до заповедника), которая нам представляется наиболее экологически значимой:

M1 - glm(TOT.N ~ D.PARK, family = poisson, data = RK) summary(M1)

Coefficients:

Estimate Std. Error z value Pr(|z|) (Intercept) 4.316e+00 4.322e-02 99.87 2e-16 D.PARK -1.059e-04 4.387e-06 -24.13 2e-16 Null deviance: 1071.4 on 51 degrees of freedom Residual deviance: 390.9 on 50 degrees of freedom AIC: 634.29 Была построена модель log(µi) = 4.13 – 0.0000106 D.PARKi, или в другой форме, µi = exp(4.13 – 0.0000106 D.PARKi), где µi - прогнозируемая величина пуассоновской частоты с учетом значения расстояния D.PARKi.

Построим график полученной зависимости:

require(ggplot2) MyData=data.frame(D.PARK=seq(from=0,to=25000,by=1000)) G - predict(M1, newdata=MyData, type="link", se=T) MyData$FIT - exp(G$fit) MyData$FSEUP - exp(G$fit+1.96*G$se.fit) MyData$FSELOW - exp(G$fit-1.96*G$se.fit) ggplot(MyData, aes(D.PARK, FIT)) + geom_ribbon(aes(ymin = FSEUP, ymax = FSELOW), fill = 3, alpha =.25) + geom_line(colour = 3) + labs(x = "Расстояние до парка", y = "Убийства на дорогах")+ geom_point(data=RK, aes(D.PARK, TOT.N) ) Обращает на себя внимание очень узкая полоса 95%-го доверительного интервала относительно линии регрессии на рис. 8.8. Это обусловлено, возможно, как неверной спецификацией построенной модели, так и заниженной оценкой остаточной дисперсии.

–  –  –

Включим в модель дополнительные предикторы с целью повышения ее адекватности. В исходной таблице данных для этого имеется 17 полезных переменных-кандидатов. Выполним, однако, предварительно две операции – преобразование данных и удаление высоко коррелированных признаков.

Поскольку 8 переменных имеют размерность площади (га), а другие заданы длиной (расстоянием в м), то разумно привести их к единой функциональной форме и извлечь квадратный корень из значений первой группы признаков. Далее оценим степень мультиколлинеарности 17-мерного комплекса переменных и вычислим для каждого предиктора фактор инфляции дисперсии (Variance Inflation Factor, VIF).

RK[, 7:14] - sqrt(RK[, 7:14]) library(car) res - sort(vif(glm(TOT.N ~., data=RK[,c(5,7:23)], family = poisson))) ; print(res,4) RK.11 - cbind(RK$D.WAT.RES, RK$D.WAT.COUR, RK$WAT.RES, RK$D.PARK, RK$L.D.ROAD, RK$L.P.ROAD, RK$MONT.S, RK$SHRUB, RK$L.WAT.C, RK$POLIC, RK$URBAN) D.WAT.RES D.WAT.COUR D.PARK WAT.RES L.P.ROAD L.D.ROAD 1.876 2.315 2.630 2.667 4.323 4.538 MONT.S L.WAT.C SHRUB POLIC URBAN OPEN.L 4.801 5.196 5.469 5.786 6.822 18.955 OLIVE L.SDI N.PATCH P.EDGE MONT 20.382 24.535 25.856 26.869 36.951

–  –  –

Напомним, что распределение разности девиансов асимптотически приближается к 2-распределению с р1-р2 степенями свободы D1 – D2 ~ 2р1-р2, чем мы воспользовались, чтобы оценить статистическую значимость однородности ошибок двух регрессионных моделей с помощью команды anova(M2, M3, test = "Chi").

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

anova(M3, test = "Chi") Terms added sequentially (first to last) Df Deviance Resid. Df Resid. Dev Pr(Chi) NULL 51 1071.44 D.WAT.RES 1 109.55 50 961.89 2.2e-16 *** WAT.RES 1 0.28 49 961.61 0.5947268 D.PARK 1 595.60 48 366.01 2.2e-16 *** L.P.ROAD 1 7.17 47 358.84 0.0074114 ** MONT.S 1 16.35 46 342.49 5.272e-05 *** SHRUB 1 20.79 45 321.71 5.138e-06 *** L.WAT.C 1 51.96 44 269.75 5.666e-13 *** POLIC 1 12.80 43 256.95 0.0003475 *** URBAN 1 17.97 42 238.98 2.240e-05 *** Очевидно, что наибольший вклад в снижение ошибки модели вносят расстояния до лесного массива, реки или другого водного резервуара (D.PARK, D.WAT.RES, L.WAT.C), в меньшей мере это относится к наличию кустарников поликультурных насаждений или (SHRUB), (POLIC) урбанизированных территорий (URBAN), а доля девианса, связанного с площадью водного резервуара (WAT.RES), не отличается от нуля.

Однако следует задаться вопросом, насколько точно была оценена нами дисперсия аппроксимирующего распределения. Поскольку остаточный девианс имеет 2-распределение с (n – p) степенями свободы, то одним из признаков избыточной дисперсии является значение = D/(n – p). Если это отношение приблизительно равно 1, то можно благополучно предположить, что избыточной дисперсии нет. В рассматриваемом примере = 239/42 = 5.69 и опасность недооценки дисперсии существует.

Мы можем компенсировать возможную избыточную дисперсию, приняв квази-распределение Пуассона и задав в функции glm() параметр family = quasipoisson.

Для этого типа распределения среднее и дисперсия для Y будет равно E(Y) = µ и var(Y) = µ:

M4 - glm(Y ~ D.WAT.RES + WAT.RES + D.PARK + L.P.ROAD + MONT.S + SHRUB + L.WAT.C + POLIC + URBAN, family = quasipoisson, data = RK.11) summary(M4) (Dispersion parameter for quasipoisson family taken 5.181428) Null deviance: 1071.44 on 51 degrees of freedom Residual deviance: 238.98 on 42 degrees of freedom AIC: NA Другой формой команды anova() является функция drop1(), которая поочередно удаляет из модели по одной объясняющей переменной и каждый раз пересчитывает значение остаточного девианса. Значимость различий отдельных моделей оценивается по 2-критерию.

drop1(M4, test = "Chi") Single term deletions Df Deviance scaled dev. Pr(Chi) none 238.98 D.WAT.RES 1 252.26 2.564 0.1092996 WAT.RES 1 258.82 3.829 0.0503756.

D.PARK 1 888.32 125.321 2.2e-16 *** L.P.ROAD 1 284.12 8.712 0.0031606 ** MONT.S 1 268.17 5.633 0.0176229 * SHRUB 1 285.09 8.899 0.0028526 ** L.WAT.C 1 301.24 12.017 0.0005271 *** POLIC 1 247.95 1.731 0.1882444 URBAN 1 256.95 3.469 0.0625355.

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

значения коэффициентов модели и оценки ее девианса не меняются;

р-величины коэффициентов увеличиваются и некоторые параметры модели оцениваются уже как статистически незначимые;

величина AIC-критерия рассчитана быть не может.

Заметим также, что ширина доверительного интервала модели на рис. 8.8 должна быть увеличена в = 7.63 раз.

Мы получили регрессионную модель, достаточно убедительную с экологической точки зрения и позволяющую ранжировать и интерпретировать степень влияния отдельных предикторов. Однако насколько она хороша для прогнозирования на новых данных? Загрузим пакет caret и проведем отбор информативного комплекса переменных методом RFE (Recursive Feature Elimination – см.

раздел 4.1):

library(caret) glmFuncs - lmFuncs glmFuncs$fit - function (x, y, first, last,...) { tmp - as.data.frame(x) tmp$y - y glm(y ~., data = tmp, family=quasipoisson(link='log')) } set.seed(13) ctrl - rfeControl(functions = glmFuncs, method = "cv", verbose = FALSE, returnResamp = "final") lmProfileF - rfe(x = RK.11, y = RK$TOT.N, sizes = 1:10, rfeControl = ctrl) Recursive feature selection Outer resampling method: Cross-Validated (10 fold) Variables RMSE Rsquared RMSESD RsquaredSD Selected 1 30.82 0.6908 10.85 0.2590 * 2 30.88 0.6144 10.83 0.2984 3 30.92 0.5726 10.81 0.3393 4 30.98 0.5394 10.80 0.3345 5 30.98 0.4872 10.80 0.3166 6 30.96 0.5077 10.80 0.3137 7 30.94 0.5203 10.81 0.2751 8 30.91 0.5650 10.79 0.3059 … … … 11 30.91 0.5619 10.78 0.2890 The top 1 variables (out of 1): D.PARK

Проверим справедливость сделанного выбора с помощью функции train():

train(TOT.N ~ D.PARK, data=RK, method='glm', family = quasipoisson, trControl = trainControl(method = "cv")) Resampling: Cross-Validated (10 fold) 1 predictor RMSE Rsquared 15.3978 0.7070051 15.3979 train(TOT.N ~ D.WAT.RES + WAT.RES + D.PARK + L.P.ROAD + MONT.S + SHRUB + L.WAT.C + POLIC + URBAN, data=RK, method='glm', family = quasipoisson, trControl = trainControl(method = "cv")) Resampling: Cross-Validated (10 fold) 9 predictor RMSE Rsquared 17.29405 0.6148829 Тестирование показало, что модель с 1 предиктором при перекрестной проверке дает меньшую ошибку, чем множественная регрессия на основе 9 параметров.

Выполним теперь построение модели на основе отрицательного биномиального распределения (NB), для чего воспользуемся функцией

glm.nb() из пакета MASS:

library(MASS) M5 - glm.nb(Y ~., data = RK.11) summary(M5) (Dispersion parameter for Negative Binomial(5.5397) family taken to be 1) Null deviance: 214.316 on 51 degrees of freedom Residual deviance: 54.764 on 40 degrees of freedom AIC: 396.92 M6 - step(M5, trace=0) summary(M6)

Coefficients:

Estimate Std. Error z value Pr(|z|) (Intercept) 3.943e+00 2.435e-01 16.194 2e-16 *** D.WAT.RES 4.478e-04 2.217e-04 2.020 0.0434 * WAT.RES 2.546e-01 1.719e-01 1.482 0.1385 D.PARK -1.447e-04 1.385e-05 -10.449 2e-16 *** L.P.ROAD 2.820e-01 1.419e-01 1.987 0.0469 * MONT.S 2.024e-01 8.882e-02 2.279 0.0227 * SHRUB -6.690e-01 3.074e-01 -2.176 0.0296 * L.WAT.C 1.644e-01 8.351e-02 1.968 0.0491 * Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

--Dispersion parameter for Negative Binomial(5.2047) Null deviance: 204.391 on 51 degrees of freedom Residual deviance: 54.054 on 44 degrees of freedom AIC: 390.57 Theta: 5.20 Std. Err.: 1.37 2 x log-likelihood:

-372.571 В дополнение к обычному набору показателей в протокол включены оценка дисперсии и параметр распределения Theta, обозначенный выше как. В отношении построенных моделей NB можно также выполнить анализ девиансов, как это мы делали выше на примере пуассоновской регрессии.

При сравнении моделей, построенных для разных распределений, необходимо использовать величины логарифмов функции правдоподобия или значения AIC-критерия: значение последнего для NB-модели AIC = 390.57 существенно меньше, чем для оптимальной модели регрессии Пуассона AIC = 498.37.

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

train(TOT.N ~ D.PARK, data=RK, method='glm.nb', trControl = trainControl(method = "cv")) Negative Binomial Generalized Linear Model 1 predictor Resampling: Cross-Validated (10 fold) link RMSE Rsquared identity NaN NaN log 31.23333 0.7155445 sqrt 29.53605 0.7155445 RMSE was used to select optimal model using smallest value.

The final value used for the model was link = sqrt.

train(TOT.N ~ D.WAT.RES + WAT.RES + D.PARK + L.P.ROAD + MONT.S + SHRUB + L.WAT.C, data = RK, method='glm.nb', trControl = trainControl(method = "cv")) 7 predictor link RMSE Rsquared identity NaN NaN log 30.99842 0.7018376 sqrt 29.16430 0.7072774 The final value used for the model was link = sqrt.

В этом случае обе модели показали весьма близкие результаты и окончательный вывод зависит от решения аналитика. Заметим, что функция train() попутно сравнила два вида функции связи (link) и более точным оказалось использование квадратного корня.

8.5. ZIP- и барьерные модели счетных данных Как обсуждалось выше, реальные данные часто характеризуются существенными отклонениями от теоретического однопараметрического распределения Пуассона: избытком нулевых значений и/или нарушениями предположения о параметре дисперсии Var(Y) E(Y). Учесть повышенную разреженность данных можно с использованием модели Пуассона со смешанными параметрами, учитывающей избыток нулей (Zero Inflated Poisson, ZIP), или специальной двухкомпонентной модели с измененными нулями (Zero-altered Poisson, ZAP), которую часто называют "барьерной" моделью (hurdle model). Если наряду с избытком нулей имеет место отчетливая избыточная дисперсия, то можно воспользоваться соответствующими моделями на основе отрицательного биномиального распределения: ZINB (Zero-inated Negative Binomial) и ZANB (Zero-altered Negative Binomial). Эти четыре модели будут предметом нашего дальнейшего обсуждения (подробности см. в Zeileis et al., 2008; Zuur et al., 2009; Cameron, Trivedi, 2013; Кшнясев, 2010).

Рассмотрим в качестве примера таблицу NMES1988 из пакета AER, содержащий данные о визитах 4406 респондентов старше 66 лет в учреждения бесплатного медицинского обслуживания США в 1988 г.

Число посещений visits примем в качестве отклика, а три метрических переменных и три фактора используем в качестве предикторов:

library(AER) data("NMES1988") # Загружаем данные и отбираем столбцы nmes - NMES1988[, c(1, 6:8, 13, 15, 18)] plot(table(nmes$visits)) sum(nmes$visits 1) # наблюдаемое число нулей [1] 683

Рис. 8.9. Наблюдаемые частоты посещений врача

Вычислим теоретическую частоту нулевых значений:

mod1 - glm(visits ~., data = nmes, family = "poisson") mu - predict(mod1, type = "response") # среднее Пуассона exp - sum(dpois(x = 0, lambda = mu)) # теоретическая частота round(exp) # нулевых значений [1] 47 Модели ZIP и ZINB называют моделями со смешанными параметрами (mixture models), поскольку присутствие нулей считается итогом двух различных процессов: биномиального процесса и процесса формирования счетной случайной величины. Здесь биномиальный процесс моделирует вероятность измерения ложно-положительного результата против всех других типов данных (включая истинные нули и конкретные счетные значения).

Иными словами, мы делим данные на три воображаемые группы:

первая группа содержит только ложные нули – наблюдения с вероятностью появления i, вторая группа – истинные нули с вероятностью (1 - i), и третья

–  –  –

M.ZINB - zeroinfl(visits ~., data = nmes, dist = "negbin", link="logit") summary(M.ZINB) Theta = 1.484 Log-likelihood:

-1.209e+04 on 17 Df Визуально сравнить степень согласия эмпирических и оцененных моделью частот можно с использованием гистограмм (поскольку для лучшего восприятия графика из частот извлекают квадратный корень, их называют еще "рутограммами" – от англ. "root", корень):

library(countreg) rootogram(M.ZIP) rootogram(M.ZINB) Рис. 8.10. Рутограмма частот для модели ZIP В результате введения в модель на основе отрицательного биномиального распределения ZINB дополнительного параметра Theta, компенсирующего избыточную дисперсию, оценки коэффициентов и их уровней значимостей были уточнены, максимум логарифма функции правдоподобия увеличился, а гистограмма на рис. 8.11 приобретает более плавный характер.

Рис. 8.11. Рутограмма частот для модели ZINB

Заметим, что мы построили простейший вариант моделей, исходя из предположения о том, что параметр i равен постоянному значению (т.е.

= const). Изменяя структуру правой части объекта formula, мы можем задать любую зависимость для от имеющихся ковариат (Zuur et al., 2009).

Барьерная модель по своему смыслу является уже не смешанной, а двухкомпонентной моделью, т.е.

задается как комбинация двух самостоятельных уравнений регрессии: логит-регрессии для моделирования вероятности Pr(yi = 0) на основе биномиального распределения, и множественной регрессии для ненулевой средней численности в предположении одно- или двухпараметрического распределения для стохастической компоненты, например:

–  –  –

Count model coefficients (truncated poisson with log link):

Estimate Std. Error z value Pr(|z|) (Intercept) 1.406459 0.024180 58.167 2e-16 *** hospital 0.158967 0.006061 26.228 2e-16 *** healthpoor 0.253521 0.017708 14.317 2e-16 *** healthexcellent -0.303677 0.031150 -9.749 2e-16 *** chronic 0.101720 0.004719 21.557 2e-16 *** gendermale -0.062247 0.013055 -4.768 1.86e-06 *** school 0.019078 0.001872 10.194 2e-16 *** insuranceyes 0.080879 0.017139 4.719 2.37e-06 ***

Zero hurdle model coefficients (binomial with logit link):

Estimate Std. Error z value Pr(|z|) (Intercept) 0.043147 0.139852 0.309 0.757688 hospital 0.312449 0.091437 3.417 0.000633 *** healthpoor -0.008716 0.161024 -0.054 0.956833 healthexcellent -0.289570 0.142682 -2.029 0.042409 * chronic 0.535213 0.045378 11.794 2e-16 *** gendermale -0.415658 0.087608 -4.745 2.09e-06 *** school 0.058541 0.011989 4.883 1.05e-06 *** insuranceyes 0.747120 0.100880 7.406 1.30e-13 ***

--Log-likelihood:

-1.613e+04 on 16 Df M.ZANB - hurdle (visits ~., data = nmes, dist = "negbin", link="logit") summary(M.ZANB) rootogram(M.ZANB) Theta = 1.484 Log-likelihood:

-1.209e+04 on 17 Df

–  –  –

На вопрос о том, какая модель эффективнее, можно дать очевидный ответ (M.ZANB), основанный на оптимальных значениях логарифма максимального правдоподобия и информационного критерия:

–  –  –

9.1. Преобразование данных и вычисление матрицы расстояний В большинстве случаев статистика оперирует с результатами предметно-ориентированной системы наблюдений, включающей два X Y, множества переменных:

где X – фиксированные или случайно варьируемые факторы (в экологии, например, условия среды или наличие ресурсов), которые потенциально определяют свойства изучаемых объектов Y (обилие видов, показатели здоровья и т.д.). С формальной точки зрения результаты мониторинга еще не представляют собой строго определенного информативного пространства, для создания которого необходимо в рамках поставленной задачи задать структуру данных и количественную меру отношений между объектами.

Как упоминалось в главе 1 при обсуждении природы многомерного отклика, современные сложные системы (прежде всего, экономические и экологические) часто включают связанные ансамбли данных, состоящие из некоторого количества однородных компонент S, что определяет ряд особенностей построения моделей на их основе. Если отсортировать список S={x1, x2,..., xs} по частотам встречаемости элементов каждого класса N(хi ), то получим функцию рангового распределения Ф(r) = N(x), оценивающую вероятности N(r)/N в зависимости от ранга r(х). Все устойчивые плотности таких распределений стартуют с больших величин N(r) и круто убывают при увеличении r приблизительно как гиперболы, имеющие длинные правосторонние "хвосты". Математики по этому поводу обнаружили шутливую закономерность: «20% жителей выпивают 80% пива».

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

Поясним смысл сказанного на примере. Пусть мы планируем оценить сходство двух водоемов по обилию донных беспозвоночных организмов, взяв соответствующие гидробиологические пробы. Животный мир бентоса состоит из самых различных видов: от крупных моллюсков до мелких нематод, индивидуальный вес которых отличается в тысячи раз. Сравним водоемы, рассчитав евклидово расстояние двух векторов биомасс {x1, x2, …, xm} для обнаруженных m видов. Мы увидим, что это расстояние целиком определяется разностью масс 1-2 видов крупнейших моллюсков, а остальные (m – 1) видов даже не имело никакого смысла отлавливать. Если мы выполним, например, стандартизацию биомассы каждого вида на диапазоне [0, 1] и снова рассчитаем расстояние Евклида, то оно также не будет способствовать выяснению сути дела, поскольку доминирующие и функционально важные виды будут иметь ту же относительную значимость, что и случайные, маргинальные или редкие виды, часто играющие ничтожную роль в экосистеме. Читатель, далекий от проблем беспозвоночных, может мысленно трансформировать этот пример на сравнение уровней материальной обеспеченности жителей двух регионов, включая олигархов и дворников, или любой иной вариант.

П. Лежандр с соавторами (Legendre, Gallagher 2001; Legendre, Legendre,

2012) разработали общие правила многомерного анализа данных и построения хорошо интерпретируемых ординационных диаграмм. Основные способы трансформации и стандартизации данных, рекомендуемые ими (применительно к подсчету численности видов в экологических исследованиях), представлены функцией decostand() в пакете vegan:

decostand(x, method, MARGIN), где MARGIN = 1, если операция применяется к строкам таблицы x, и 2 – если к столбцам (обычно именно это значение принимают по умолчанию).

Параметр method может принимать следующие значения:

normalize Сумма квадратов значений по строкам делается равной 1 total Деление на суммы по строкам hellinger Корень квадратный из значений по методу total max Деление на максимумы по столбцам freq Деление на максимумы по столбцам и умножение на число ненулевых компонент chi.square См. формулу (9.1) ниже log Логарифмическая трансформация (не требует добавления 1) range Данные по столбцам стандартизуются на диапазоне [0, 1] standardize Обычная стандартизация x к нулевому среднему и единичной дисперсии; см. scale(x, center = TRUE, scale = TRUE) pa Приведение x к бинарной шкале (0/1) В представленной выше ситуации с бентосом, рекомендуется либо простейшая логарифмическая трансформация, либо преобразование, приводящее к 2-дистанции, которая является, видимо, наиболее разумным компромиссом при учете как роли ведущих компонент, так и вклада длинного правого "хвоста". Такое преобразование имеет вид xij xij = x + + ', (9.1) xi + x+ j где xi+ – сумма по строкам, x+j – сумма по столбцам, x++ – общая сумма элементов таблицы x.

Вторым важным этапом является спецификация метрики отношений между объектами.

В частном случае при применении статистических методов информативное пространство может интерпретироваться как вероятностное:

тогда пара векторов действительных чисел {x1, x2, …, xm} и {y1, y2, …, ym}, описывающих произвольные объекты x и y, будут трактоваться как выборочные реализации m-мерной случайной величины. В этом случае в качестве мер сходства между объектами могут выступать оценки ковариации cov( x, y ) = i =1 ( xi mx )( yi m y ), коэффициент корреляции m rxy = cov( x, y ) / x y или произвольное ковариационное отношение K = [cov(x, y) - covmin]/(covmax covmin), где m – математическое ожидание, – стандартное отклонение, covmin и covmax – экстремальные значения ковариации для теоретической ("эталонной") выборки (Воробейчик, 1993).

В общем случае использование вероятностных представлений совершенно не обязательно.

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

а) тождества (x, y) = 0 при x = y, б) симметрии (x, y) = (y, x) и в) правила треугольника (x, y) + (y, z) (x, z). Конкретной дефиницией функции может быть, например, обобщенная мера Минковского (см. раздел 3.4), наиболее популярными реализациями которой являются манхэттенская ("manhattan") или евклидова ("euclidean") дистанции, а также расстояние Хемминга ("binary"), равное числу совпавших единиц для двух бинарных кодов. Эти опции являются основными для базовой функции R dist(x, method = "euclidean"), вычисляющей матрицу дистанций между всеми парами объектов, которые представлены строками таблицы x.

Более широкими возможностями обладает функция из пакета vegan vegdist(x, method = "bray", binary = FALSE), где параметр method охватывает значительную часть "изобретений" экологов в области мер сходства/расстояния: "manhattan", "euclidean", "canberra","bray", "kulczynski", "jaccard", "gower", "altGower", "morisita", "horn", "mountford", "raup", "binomial", "chao" и "cao".

По умолчанию используется мера Брея–Кёртиса :

m ab m m M xy = 1 2 min[ xi, yi ] /( i =1 xi + i =1 yi ), i =1

–  –  –

Другие его названия: индекс Ренконена, процентное подобие, коэффициент общности, индекс Штейнгауза, количественная мера сходства Чекановского и т.д.

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

Отметим также, что некоторые из вышеперечисленных мер (например, коэффициенты Жаккара и Съеренсена) являются коэквивалентными, т.е. они порождают одну и ту же предупорядоченность анализируемых объектов.

Более сильное влияние на характер результатов оказывают различия в функциональной форме разностей x и y, такие как абсолютная, квадратичная или логарифмическая.

9.2. Непараметрический дисперсионный анализ матриц дистанций Результаты наблюдений, представленные в виде многомерной таблицы объемом mp, обычно можно разделить по строкам на r групп, (например, в соответствии со списком регионов, где оценивается покупательский спрос, или шириной водотока, где выполняются гидробиологические пробы). Тогда, в соответствии с моделью дисперсионного анализа, общую изменчивость каждого k-го фиксированного показателя xk (k = 1, 2, …, p) можно разложить Var xk = Var + Var, на компоненты:

где Var – вариация, обусловленная влиянием группирующего фактора, Var – изменчивость, связанная с воздействием неконтролируемых (в том числе, случайных) факторов.

Однако, если количество p одновременно варьируемых показателей велико (например, номенклатура из 300-400 продаваемых товаров), то дисперсионный анализ не дает возможности оценить совокупную изменчивость всего изучаемого многомерного комплекса объектов под воздействием группирующего фактора. Один из приемов избежать "проклятия размерности" – выполнить дисперсионный анализ симметричной матрицы D размерностью mm, элементами dij которой являются коэффициенты расстояния между каждой парой изучаемых объектов i и j.

Метод, известный как непараметрический многофакторный дисперсионный анализ (npMANOVA), осуществляет разложение многомерной изменчивости, заключенной в матрице расстояний, в соответствии с уровнями влияния изучаемых факторов. В первом приближении, как это показано на рис.

9.1, можно рассчитать общую SST и внутригрупповую SSW суммы квадратов dij, после чего оценить их дисперсионное отношение:

SS (m r ) 1 m 1 m 1 m 1 m SST = i =1 j =i +1 d ij ; SSW = i =1 j =i +1 d ij ij ; F = W, (9.2) SST (r 1) r m где ij равны 1, если объекты i и j принадлежат одной группе, и 0 в противном случае (Anderson, 2001).

Рис. 9.1. Схема группировки матрицы расстояний Более корректный способ вычисления тестовой статистики F использует ряд конструкций линейной модели ANOVA (McArdle, Anderson, 2001). При применении этого подхода на основе исходной матрицы расстояний D с использованием преобразований Гувера рассчитывается матрица G, описывающая "центроиды", т. е. центры распределения значений dij в пределах каждой группы. Другая проекционная матрица H (или "шляпа", hat) вычисляется с использованием вектора ортогональных контрастов, связанных с уровнями группировочных факторов. Тогда тестовая статистика для проверки нулевой гипотезы об отсутствии различий между r группами, или псевдо-F-критерий, находят как tr (HGH )(m r ) F=, (9.3) tr[(I H )G (I H )](r 1) где tr – сумма диагональных элементов матрицы, I – единичная матрица. При выполнении требований к конфигурации данных и использовании евклидова расстояния обе формулы (9.2) и (9.3) становятся эквивалентными и воспроизводят традиционное дисперсионное отношение Фишера.

Распределение многомерной псевдо-F-статистики (9.2-9.3) при справедливости нулевой гипотезы теоретически неизвестно, и поэтому для его аппроксимации используется рандомизационная процедура. Ее практическая стратегия состоит в том, чтобы получить большое подмножество значений F*, извлеченных случайно, независимо и с равной вероятностью из нуль-распределения F*, полученного в результате многократного хаотичного перемешивания наблюдений относительно выделенных групп. Р-значение, соответствующее верной нулевой гипотезе, вычисляется как доля значений F*, которые больше или равны значению F, найденному на реальных данных.

Непараметрический дисперсионный анализ Андерсона реализован в пакете TraMineR (см.

раздел 5.5) или пакете vegan следующей функцией:

adonis(formula, data, permutations = 999, method = "bray"), где объект formula чаще всего имеет форму Y ~ A + B*C, в которой Y – таблица исходных данных, которая пересчитывается в матрицу дистанций по выражению method, а A, B и C – независимые ("объясняющие") факторы или непрерывные переменные. Если факторные пространства являются непересекающимися по какому-то показателю, то, указав его в качестве параметра strata, можно выполнить гнездовой дисперсионный анализ (nested ANOVA).

Рассмотрим проведение анализа на примере таблицы dune, которая содержит данные по 20 пробным участкам, заложенным в дюнах одного датского острова, на которых было обнаружено 30 видов луговой травянистой растительности. Обилие каждого вида оценивалось в баллах проективного покрытия от 0 до 7.

Параллельно для каждого участка измерялась толщина гумусового слоя почвы A1 и оценивалась одна из 4-х категорий Management агротехнических мероприятий по улучшению травостоя (таблица dune.env):

library(vegan) data(dune) data(dune.env) results - adonis(dune ~ Management*A1, data=dune.env, permutations=499) Terms added sequentially (first to last) Df SumsOfSqs MeanSqs F.Model R2 Pr(F) Management 3 1.4686 0.48953 3.2629 0.34161 0.002 ** A1 1 0.4409 0.44089 2.9387 0.10256 0.018 * Management:A1 3 0.5892 0.19639 1.3090 0.13705 0.212 Residuals 12 1.8004 0.15003 0.41878 Total 19 4.2990 1.00000

--Signif. codes: 0 С***Т 0.001 С**Т 0.01 С*Т 0.05 С.Т 0.1 С Т 1 Традиционная таблица дисперсионного анализа полученных результатов показывает статистическую значимость обоих основных факторов и незначимость эффекта их парного взаимодействия. Например, рзначение, равное 0.018, означает, что только в 9 случаях из 500 имитируемое значение F*, полученное при справедливости нулевой гипотезы, превысило значение F, вычисленное по (9.3) для эмпирических данных.

Для однофакторного дисперсионного анализа и оценки однородности многомерной вариации в группах наблюдений удобно использовать также другую функцию vegan – betadisper(), на вход которой подается произвольная матрица дистанций d.

Выполним расчеты с применением этой функции на основе группировочного фактора group = Management:

# Предварительно рассчитываем матрицу расстояний Брея-Кёртиса D - vegdist(dune) mod - betadisper(D, dune.env$Management, type = "centroid") plot(mod, hull= FALSE, main="") ; legend("topright", levels(dune.env$Management), pch=1:4, col=1:4) # Выполнение рандомизационного теста и парных сравнений permutest(mod, pairwise = TRUE) Permutation test for homogeneity of multivariate dispersions Response: Distances Df Sum Sq Mean Sq F N.Perm Pr(F) Groups 3 0.12457 0.041525 3.8687 999 0.025 *

–  –  –

Рис. 9.2. Центроиды групп в пространстве двух главных координат В целом перестановочный тест показал значимое (р = 0.025) влияние агротехнических мероприятий на вариацию видов растительности. Однако выполненные парные сравнения выявили значимые отличия только для метода NM (оставить луга в естественном состоянии) от биологического BF или рекреационного HF вмешательства. Относительные расстояния между центроидами групп можно приблизительно оценить по ординационной диаграмме (рис. 9.2).

Для выполнения множественных сравнений наиболее распространённым и рекомендуемым в литературе является тест Тьюки, использующий критерий "подлинной" значимости (Honestly Significant Difference, HSD). HSD оценивает наименьшую величину разности математических ожиданий в группах, которую можно считать значимой. Тест Тьюки-HSD реализуется в R для модели betadisper с использованием базовой функции TukeyHSD(). На рис. 9.3 представлены рассчитанные этой функцией одновременные 95%-ные доверительные интервалы разности между каждой парой групп. Поскольку этот интервал не включает 0 только для одной сравниваемой пары (NM–BF), то значимые различия в видовом составе растительности имеются лишь между этими группами.

# Тест Тьюки HSD и график доверительных интервалов (mod.HSD - TukeyHSD(mod)) ; plot(mod.HSD) Рис. 9.3. Множественные сравнения методов улучшения травостоя с использованием критерия HSD Тьюки

9.3. Методы ординации объектов и переменных:

построение и сравнение диаграмм Ординация в узком смысле представляет собой нахождение таких координатных осей на плоскости, относительно которых можно выполнить оптимальное проецирование многомерных анализируемых объектов.

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

К настоящему времени в среде R разработано и реализовано значительное количество алгоритмов сжатия информационного пространства, в частности:

метод главных компонент PCA (см. раздел 2.4), оперирующий с ковариационной (корреляционной) матрицей;

метод главных координат PCoA и неметрическое многомерное шкалирование NMDS (nonmetric multidimensional scaling), выполняющие последовательную процедуру преобразования любой матрицы дистанций;

анализ соответствий или корреспондентный анализ CA (correspondence analysis), основанный на итерационной процедуре встречного усреднения взвешивающих коэффициентов для объектов и переменных (Джогман и др., 1999);

функции многомерного факторного анализа, включая версии для иерархических и смешанных данных, представленные в пакете FactoMineR.

Метод главных координат PCoA, или многомерное шкалирование (MDS, multidimensional scaling), во многом похож на РСА, но вместо корреляционной матрицы выполняет вычисление собственных значений и собственных векторов произвольной квадратной симметричной матрицы расстояний D. Так можно компенсировать некоторые отклонения от предпосылок в отношении статистического распределения данных, принятых для корреляционного анализа, но одновременно возникает проблема выбора подходящей метрики дистанции. На странице gastonsanchez.com описана работа с 7 функциями различных пакетов R, в которых реализован метод PCoA, но обычно используют базовую функцию cmdscale().

В качестве примера рассмотрим достаточно представительный набор геоботанических данных bryceveg из пакета labdsv, в котором 165 видов травянистой растительности было обнаружено на 160 пробных площадках, заложенных на территории Брайс-Каньо;н (Bryce Canyon, штат Юта, США):

library(labdsv) data(bryceveg) ; data(brycesite) # Рассчитываем три матрицы расстояний по разным формулам dis.mht - dist(bryceveg,"manhattan") dis.euc - dist(bryceveg,'euclidean') dis.bin - dist(bryceveg,"binary") # Создаем объекты cmdscale mht.pco - cmdscale(dis.mht, k=10,eig=TRUE) euc.pco - cmdscale(dis.euc, k=10,eig=TRUE) bin.pco - cmdscale(dis.bin, k=10,eig=TRUE) # График изменения относительных собственных чисел plot(1:10, (mht.pco$eig/sum(mht.pco$eig))[1:10], type="b",xlab="Число координат", ylab="Относительные собственные значения") lines((euc.pco$eig/sum(euc.pco$eig))[1:10],type="b",col=2) lines((bin.pco$eig/sum(bin.pco$eig))[1:10],type="b",col=3) legend ("topright",c("Манхэттен","Евклид", "Хемминг"), lwd=1, col=1:3) Рис. 9.4. Уменьшение собственных значений матриц с увеличением числа координат Возникает естественный вопрос, какая из трех наиболее характерных метрик расстояния приводит к наилучшей ординации? Если ориентироваться на величину первых двух собственных значений матрицы D (рис. 9.4), которые соответствуют доле объясненной дисперсии видовой структуры участков, то наименее адекватной оказывается ординация, полученная при использовании расстояния Хемминга, т.е. когда учитывалось только наличие или отсутствие видов растений. Второй способ – оценить коэффициент корреляции двух матриц расстояния: исходной и рассчитанной на основе первых двух главных координат, что соответствует доле информации, сохраненной после сжатия 165-мерного пространства до двумерного.

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

plot(euc.pco$points[,1:2], col = 4, pch=17, xlab = "PCO1", ylab = "PCO2", main="Расстояние Евклида") plot(0-bin.pco$points[,1:2], col = 3,pch=17, xlab = "PCO1", ylab = "PCO2", main="Расстояние Хемминга") CCor -c(cor(dis.euc, dist(euc.pco$points[,1:2],'euclidean')), cor(dis.mht, dist(mht.pco$points[,1:2],'manhattan'))) names(CCor) - c("Евклидово","Манхэттенское") Евклидово Манхэттенское 0.8171305 0.8256334

Рис. 9.5. Распределение точек ординации PCoA для разных метрик дистанции

Важной косвенной оценкой, насколько полученная ординация адекватна задаче исследования, является ее тесная взаимосвязь с ведущими внешними факторами. Функция ordtest() из пакета labdsv выполняет оценку р значимости такой связи с использованием рандомизационного теста, рассчитывая сумму квадратов разностей двух матриц расстояния (исходной и построенной на основе фактора). Другой подход основан на аппроксимации двумерной сглаживающей поверхностью значений предиктора в точках, определенных главными координатами. Подробности изложены в курсе лекций проф. Д. Робертса студентам университета штата Монтана (http://ecology.msu.montana.edu). Ниже используется функция surf.pco(), которая выполняет сглаживание сплайнами и визуализацию на диаграмме линий изоном фактора.

Ее скрипт приведен в лекции № 8 и скопирован нами в файл surf_pco.r:

data(brycesite) ; H - brycesite$elev mht.pco - pco(dis.mht, k = 2) # аналог функции cmdscale() source("surf_pco.r") plot(mht.pco, pch=16, col="red") surf.pco(mht.pco, H, col="blue") ordtest(mht.pco,H)$p var ~ s(x, y)

Estimated degrees of freedom:

16.077 total = 17.07685 GCV score: 189195.3 D^2 = 0.5999 [1] 0.001 Рис. 9.6. Ординация геоботанических данных с использованием манхэттенского расстояния и ее связь с высотой пробных площадок Для построения сглаживающей поверхности высоты участков над уровнем моря по осям двух главных координат использовалась функция gam() из пакета mgcv, выполняющая аппроксимацию сплайном s(x, y).

Качество подгонки может быть оценено с помощью псевдо-коэффициента детерминации, который для произвольной GAM-модели m рассчитывается как D2 = (m$null.deviance - m$deviance)/m$null.deviance.

Если построить различные варианты ординации с использованием разных методов или метрик, то наилучшей будет считаться та, которая доставляет максимум D2. При этом, разумеется, необходимо выбрать для тестирования ведущий внешний фактор. Например, если вместо высоты над уровнем моря использовать угол наклона площадки к горизонту brycesite$av, то его статистическая связь с изменчивостью видового состава растительности будет незначимой р=0.544, а величина D2 будет мала D^2 = 0.06873.

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

Для этого проводится ряд последовательных итераций подгонки масштаба осей:

сравниваемые точки помещаются внутри единичного круга || = 1;

вся структура вращается относительно центра координат;

находится минимум суммы квадратов расстояний между двумя ординациями m 2 = [x1 T (x 2 )]2 min или максимум прокрустовой корреляции r = 1 m.

Используя функции из пакета vegan, выполним сравнение расположения площадок на двух диаграммах, построенных методами главных координат (mht.pco на рис.

9.6) и анализом соответствий CA:

library(vegan) ca - cca(bryceveg) pro.comp - procrustes(mht.pco, ca) plot(pro.comp) Procrustes sum of squares: 1.164e+04 Рис. 9.7. Сравнение двух диаграмм после прокрустового преобразования

–  –  –

где Fs(i) и Gs(k) – координаты i-го объекта и k-й переменной на s-й оси компоненты, а s – собственное значение, связанное с s-й осью.

В случае анализа соответствий такая диаграмма для 165 видов на территории Брайс-Каньо;н выглядит как на рис.

9.8:

plot(ca, type="n") # Отрисовка «пустой» диаграммы ordilabel(ca, dis = "sp", cex=0.7, font=3,add = TRUE) Рис. 9.8. Ординационная диаграмма видов растительности, полученная методом анализа соответствий Если вид встречается только в одном месте, то его положение на диаграмме совпадает с точкой, обозначающей соответствующую пробную площадку. В противном случае местоположение вида в системе осей СА1-СА2 определяется средневзвешенными координатами нескольких его возможных местообитаний. По-видимому, эта точка определяет экологический оптимум вида, где его появление наиболее вероятно, либо обилие максимально.

Функции ordipointlabel() и ordilabel() из пакета vegan дают возможность расположить на графике метки данных наилучшим образом и задать последовательность их визуализации: в верхнем слое диаграммы на рис. 9.8 показаны экологически важные виды с наибольшей суммой обилия по столбцам.

Важное место в современном статистическом анализе занимает построение внешне элегантных и информационно насыщенных графиков.

Ведущая роль при этом отводится системе визуализации ggplot2 (Мастицкий, 2016). В предыдущих главах мы приводили ординационные диаграммы, полученные с использованием функций ggplot (см. рис. 2.10, 2.11, 2.13, 7.3, 7.7), когда обсуждали метод главных компонент и дискриминантный анализ. Но читатель, вероятно, смог заметить, что построение графиков ggplot2 – непростой, хотя и увлекательный процесс.

Возможность построения элегантных графиков "одной строкой" предоставляется специализированными пакетами, использующими "в темную" инструментарий графической системы ggplot2. К таким относится пакет factoextra, предназначенный для визуализации результатов многомерного и кластерного анализа (см. также главу 10). Подробный рассказ об этом пакете см. на сайте http://www.sthda.com (STHDA, Statistical tools for high-throughput data analysis).

Функция viz_pca() из пакета factoextra осуществляет построение ординационных диаграмм для объектов и/или переменных (отдельно или совмещенно) на основе PCA-объектов, полученных функциями prcomp(), princomp(), PCA() из пакета FactoMineR и некоторых других.

Набор данных decathlon из пакета FactoMineR (Le et al., 2008) содержит результаты 41 ведущего десятиборца, показанные в 2004 г. на олимпийских играх и турнире "ДекаСтар".

Выполним его PCA-ординацию:

library(FactoMineR) library("factoextra") data("decathlon") colnames(decathlon) - c("100м","Длина.прыжок","Ядро", "Высота.прыжок","400м","110м.барьер","Диск","Шест.прыжок", "Копье","1500м", "Место","Очки", "Соревнование" ) res.pca - PCA(decathlon[, 1:10], scale.unit = TRUE, ncp = 5, graph = FALSE) head(res.pca$eig) eigenvalue percent of var cumulative percent of var comp 1 3.2719055 32.719055 32.71906 comp 2 1.7371310 17.371310 50.09037 comp 3 1.4049167 14.049167 64.13953 comp 4 1.0568504 10.568504 74.70804 comp 5 0.6847735 6.847735 81.55577 comp 6 0.5992687 5.992687 87.54846 Первые две компоненты объясняют только 50% статистического разброса результатов. Вклад (contribution) каждой переменной при формировании каждой главной компоненты может быть оценен как отношение квадрата факторной нагрузки (score) к соответствующему собственному значению (Abdi, Williams, 2010). Построим ординационную диаграмму признаков (т.е.

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

res.pca$var$contrib fviz_pca_var(res.pca, col.var="contrib", gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"), repel = TRUE) # Раскрасим также текст Dim.1 Dim.2 Dim.3 100м 18.34376957 2.016090 2.42049891 Длина.прыжок 16.82246707 6.868559 2.36319121 Ядро 11.84353954 20.606785 0.03890276 Высота.прыжок 9.99788710 7.063694 4.79362526 400м 14.11622887 18.666374 1.23027094 110м.барьер 17.02011495 3.013382 0.61083225 Диск 9.32848615 21.162245 0.13131711 Шест.прыжок 0.07745541 1.872547 34.06090024 Копье 2.34696326 5.784369 10.80714169 1500м 0.10308808 12.945954 43.54331962

Рис. 9.9. Ординационная диаграмма видов соревнований в десятиборье

По сути примера видно, что первая компонента связана с темповыми видами спорта (бег на короткие дистанции), а вторая – с силовыми (метание диска и ядра).

Оценим теперь, имеются ли различия между результатами, показанными на Олимпиаде и на турнире "ДекаСтар".

Для этого в структуру объекта PCA включим в качестве вспомогательных переменных признаки, не участвующие в построении ординации:

res.pca - PCA(decathlon, quanti.sup = 11:12, quali.sup = 13, graph=FALSE) fviz_pca_ind(res.pca, geom = "point", habillage = 13, addEllipses =TRUE, ellipse.level = 0.68) + scale_color_brewer(palette="Dark2") + theme_minimal() Рис. 9.10. Ординационная диаграмма участников соревнований в десятиборье Нетрудно убедиться в том, что по оси первой главной компоненты обозначился "зазор" между групповыми центроидами и результаты, показанные на Олимпиаде, можно считать существенно более высокими. В то же время, разброс показателей на турнире "ДекаСтар" оказался меньше, т.е.

состав участников был более однородным.

9.4. Оценка связи ординации с внешними факторами Ранее мы рассматривали ординационный анализ только одной матрицы Y, чтобы получить отображение в ортогональной системе координат некоторых частных структурных особенностей изучаемой системы в форме графических проекций объектов и их признаков. Однако это – пассивный путь, поскольку он ограничивается фиксацией некоторого распределения точек на плоскости и не дает никаких объяснений этому феномену. Такую ординацию называют также непрямой, или неограниченной (unconstrained), в отличие от прямой, или ограниченной (constrained) ординации, которая ставит задачу связать внутреннюю изменчивость матрицы Y с теми или иными внешними воздействиями X. Нам кажется, что непосредственный перевод приведенных терминов не вполне отражает суть дела и было бы точнее назвать эти методы ординации "необъясняющим" и "объясняющим" соответственно.

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

Y = f(X), где Y матрица размером nm содержит значения отклика уij, измеренные по m признакам (столбцы) для n объектов (строки); X матрица размером nq, в которой i-я строка содержит значения объясняющих переменных xik.

Наиболее часто используются две математические процедуры канонической ординации, известные под аббревиатурами RDA и CCA. В обоих случаях формой канонизации является вычисление собственных чисел и собственных векторов.

Анализ избыточности (RDA, redundancy analysis Rao, 1964; Legendre, Legendre, 2012), общий формализм и алгебра которого были описаны в разделе 2.5, является расширенной канонической формой, объединяющей линейную модель множественной регрессии и метод главных компонент РСА. Канонический анализ соответствий (CCA, canonical correspondence analysis ter Braak, 1986) является модификацией RDA, выполняющей "взвешивание" ординационных компонент пропорционально их вкладу в статистику 2, которая оценивает расстояния между объектами в многомерном пространстве. В связи с этим ССА в экологических приложениях придает существенно бо;льшее значение редко встречающимся объектам, чем это делает RDA, основанный на евклидовых дистанциях.

В целом, оба метода RDA и CCA осуществляют аппроксимацию данных многомерной гауссовой регрессией и выполняют подгонку таких коэффициентов a = (ai) для переменных матрицы Y, i = 1,..., m, и c = (ck) для факторов среды X, k = 1,..., q, которые делают максимальной корреляцию между z* = Y'a и z = X'c. Расчеты в среде R обычно проводят с помощью функций rda() или cca() из пакета vegan. С этими функциями связано еще не менее двух десятков сервисных функций, выполняющих вычисление различных статистик, проверку гипотез о значимости как внутренних взаимодействий, так и о влиянии внешних факторов, и т.д.

Рассмотрим еще один геоботанический пример – набор varespec, описывающий относительное обилие 44 видов травянистых растений на 24 пробных площадках. Параллельный набор varechem содержит данные о 14 показателях химического состава почвы на тех же площадках.

Оценим предварительно характер распределения частот встречаемости объектов по рангам r списка видов S. Осуществим аппроксимацию данных моделью Ципфа Nr = Np1 r с двумя оцениваемыми параметрами: p1 – долей самого многочисленного класса в общей численности N и – коэффициентом затухания гиперболической зависимости.

Правда, мы имеем дело не с числом экземпляров растений, а с проективным покрытием в %, но нет никакой сложности превратить формально непрерывные значения обилия в ряд положительных целочисленных величин:

library(vegan) data(varechem) ; data(varespec) B - rev(sort(colSums(varespec))) BF - as.matrix(t(round(B/min(B)))) rad - rad.zipf(BF) barplot(BF) lines(1:ncol(BF), rad$fitted, col="green", lwd=2) RAD model: Zipf Family: poisson No. of species: 44 Total abundance: 24175 p1 gamma Deviance AIC BIC 0.34724 -1.29554 9036.36623 9308.19192 9311.76030

Рис. 9.11. Аппроксимация рангового распределения видов моделью Ципфа

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

Оценим влияние отдельных химических элементов в почве на встречаемость видов с использованием канонического анализа соответствий CCA (рис.

9.12):

varespec.chi - decostand(varespec, method="chi.square") mod1.cca - cca(varespec.chi ~., varechem) spec - as.data.frame(mod1.cca$CCA$v[,1:2]) # Виды sites - as.data.frame(mod1.cca$CCA$u[,1:2]) # Площадки fact - as.data.frame(mod1.cca$CCA$biplot[,1:2]) # Факторы library(ggplot2) ggplot() + geom_point(data=sites,aes(x=CCA1,y=CCA2),size=3,shape=17, color="red") + # добавляем точки площадок geom_text(data=spec,aes(x=CCA1,y=CCA2, label = abbreviate(rownames(spec),4)),colour="darkgreen", fontface = 3) + # добавляем метки видов geom_segment(data=fact,aes(x=0,xend=CCA1,y=0,yend=CCA2), arrow = arrow(length = unit(0.5, "cm")), colour="darkblue") + # добавляем стрелки факторов geom_text(data = fact, aes(x = CCA1, y = CCA2, label = rownames(fact), size = 6), hjust = 1.2, fontface = 2) + theme_bw() Inertia Proportion Rank Total 3.2761 1.0000 Constrained 2.2603 0.6899 14 Unconstrained 1.0158 0.3101 9

Eigenvalues for constrained axes:

CCA1 CCA2 CCA3 CCA4 CCA5 CCA6 CCA7 … CCA14 0.4624 0.3017 0.2815 0.2651 0.2333 0.2028 0.1476 … 0.0190

Eigenvalues for unconstrained axes:

CA1 CA2 CA3 CA4 CA5 CA6 … CA9 0.29087 0.21535 0.17088 0.09812 0.07670 0.06552 … 0.02342

Рис. 9.12. Ординационный триплот, полученный методом CCA

Из результатов следует, что всего было сформировано 23 (т.е. n 1) линейно независимых оси ординации, в том числе 14 осей объясняющих переменных (Constrained). 69% совокупной многомерной дисперсии (Inertia) матрицы связаны с химическим составом почвы, а 31% со случайными или неучтенными факторами.

Результирующая диаграмма ("триплот" см. рис. 9.12) отражает не только изменчивость видового состава относительно двух осей проецирования CCA1-CCA2, но и статистические связи между видами и каждой независимой переменной xji. Для этого из центра координат стрелками синего цвета проведены дополнительные оси (constrained axes), ориентация которых зависит от значений канонических собственных векторов. Относительная длина стрелки пропорциональна силе связи внешнего фактора с осями ординации CCA1-CCA2, а косинус угла между каждой парой стрелок приблизительно равен коэффициенту корреляции между факторами (положительной или отрицательной).

Как и в случае биплота, координаты каждого вида находятся вблизи точек тех площадок, где его обилие наиболее высоко. Если виды представить стрелками, исходящими из начала координат, то косинус угла между стрелкой соответствующего вида и стрелкой фактора среды приблизительно равен коэффициенту корреляции между ними. Например, обилие вида Callvulg (Cllv) в значительной степени коррелирует с содержанием в почве молибдена, в отличие от Cladphyl (Cldp). Проекция точки вида на каждую стрелку показывает его экологический оптимум (точнее, центр тяжести распределения обилия) относительно этого химического компонента.

Нетрудно убедиться из графика на рис. 9.12, что показатели состава почвы очень высоко коррелирут друг с другом, о чем свидетельствует фактор инфляции дисперсии VIF.

library(car) sort(vif(mod1.cca)) paste("AIC=", extractAIC(mod1.cca)[2]) N Baresoil Mn Humdepth Mo P pH 1.637648 1.966942 5.163652 5.258233 5.760964 6.398181 6.593768 Zn Fe Ca Mg K S Al 8.446984 9.100191 9.295130 10.531890 11.705305 20.424393 21.546563 [1] "AIC= 65.58352" Дисперсионный анализ anova() в отношении объектов ССА или RDA имеет важный параметр by: при by = "axis" оценивается значимость объясняющих осей, значение by дает возможность = "term" последовательно проверить по F- и -критериям значимость отдельных коэффициентов регрессионной модели, а в случае by = "margin" проверяется наличие маргинального эффекта. Если параметр by не задан, то оценивается значимость ординационной модели в целом. Естественно, что все проверки осуществляются с использованием перестановочного теста, т.е.

проводится многократное случайное перемешивание строк исходных таблиц

–  –  –

При выполнении пошаговой процедуры уже после 4-го шага началось увеличение AIC-критерия до значения 65.58 при полном составе переменных. В итоге была получена модель с высокой степенью статистической значимости в целом, однако использование уже третьей факторной оси возможно оказалось излишним. Судя по всему, алгоритм не смог четко распорядиться высокой корреляционной зависимостью между содержанием фосфора и калия (см. рис. 9.11).

9.5. Неметрическое многомерное шкалирование и построение распределения чувствительности видов Перспективным методом ординации, находящим все большее применение в различных предметных областях, является алгоритм неметрического многомерного шкалирования (NMDS, nonmetric multidimensional scaling – Краскел, 1986; McCune, Grace, 2002). Его главным преимуществом является то, что он не требует от исходных данных никаких априорных предположений о характере статистического распределения.

Считается, что этот метод дает наиболее адекватные результаты, особенно для больших матриц с сильными "шумами".

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

Для этого обычно используют величину "стресса":

n d ij d ij, где dij и dij – расстояния между объектами i и j в исходном и = d ij i, j =1 редуцированном пространствах, и – задаваемые коэффициенты.

Пакет vegan содержит ряд функций для проведения расчетов методом NMDS и построения объясняющих ординаций. Рассмотрим пример, включающий наборы данных varespec и varechem, которые использовались нами в предыдущем разделе для ординации каноническим методом CCA, основанным на линейных преобразованиях. На первом этапе определим наилучшую метрику расстояния, рассчитав для каждого из ее вариантов коэффициент корреляции Спирмена между двумя матрицами дистанций (в пространствах видов растительности и химического состава почвы).

Эту операцию выполняет функция rankindex():

library(vegan) data(varechem) ; data(varespec) varespec.chi - decostand(varespec, method="chi.square") # Коэффициенты корреляции Спирмена для различных метрик rankindex(varechem, varespec.chi, c("euc","man","bray","jac","kul")) euc man bray jac kul 0.2228233 0.1894385 0.3022988 0.3022988 0.3071600 Используем для расчета матрицы дистанций формулу Кульчицкого.

Далее функция metaMDS() осуществляет неметрическое шкалирование и минимизацию стресса, а функция envfit() расчет коэффициентов корреляции каждого из показателей химического состава почвы с осями ординации NMDS1-NMDS2 и оценку статистической значимости p этих коэффициентов на основе перестановочного теста.

Функция MDSrotate() вращает систему координат таким образом, чтобы совместить ось NMDS1, например, с осью концентрации алюминия:

vare.mds - metaMDS(varespec,distance = "kul", trace = FALSE) ord.mds - MDSrotate(vare.mds, varechem$Al) ef - envfit(vare.mds, varechem, permu = 999) global Multidimensional Scaling using monoMDS Data: wisconsin(sqrt(varespec)) Distance: kulczynski Stress: 0.1825658

---NMDS1 NMDS2 r2 Pr(r) N -0.43946 0.89826 0.2536 0.043 * P -0.15541 -0.98785 0.1939 0.106 K -0.35324 -0.93553 0.1809 0.115 Ca -0.24020 -0.97072 0.4119 0.004 ** Mg -0.17161 -0.98517 0.4270 0.004 ** S 0.31435 -0.94931 0.1752 0.123 Al 1.00000 0.00000 0.5269 0.001 *** Fe 0.98839 0.15193 0.4450 0.002 ** Mn -0.99114 0.13279 0.5231 0.002 ** Zn -0.15262 -0.98828 0.1879 0.117 Mo 0.99764 0.06863 0.0610 0.521 Baresoil -0.99256 -0.12179 0.2508 0.052.

Humdepth -0.98969 -0.14321 0.5201 0.002 ** pH 0.93818 -0.34615 0.2308 0.053.

--Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Number of permutations: 999 Построим ординационные диаграммы с использованием стандартных графических средств R (как это можно сделать на основе ggplot2 см.

или github.com/gavinsimpson/ggvegan, chrischizinski.github.io stackoverflow.com). Первая диаграмма (рис. 9.13) включает ординацию площадок и оси девяти внешних факторов, для которых p 0.1.

plot(vare.mds, type = "t", disp="sites", font=2) points(vare.mds, disp="sites",cex=3, col="red") abline(h = 0, lty = 3) ; abline(v = 0, lty = 3) plot(ef, p.max = 0.1, col="blue", font=2, lwd=2, cex=1)

Рис. 9.13. Ординация площадок, полученная методом NMDS

Прежде чем построить вторую диаграмму, для концентрации в почве того или иного химического компонента yхэ по его эмпирическим значениям для каждой площадки подгонялись обобщенные аддитивные модели (Wood, yхэ = s (v1, v2) +, 2006): (9.4) где s – функция двумерного сплайна с k степенями свободы от NMDSкоординат v1 и v2. Построение этой модели осуществляется функцией ordisurf(), которая параллельно отрисовывает изолинии трехмерной сглаживающей поверхности на ординационной диаграмме.

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

На рис. 9.14 показано, как на ординацию видов растительности могут быть наложены изолинии по моделям сглаживания содержания Al и Ca.

plot(vare.mds, type="n") ordipointlabel(vare.mds, display = "sp", font=4, col="darkgreen", add=TRUE) abline(h = 0, lty = 3) ; abline(v = 0, lty = 3) with(varechem, ordisurf(vare.mds, Al, add = TRUE, col =2)) with(varechem, ordisurf(vare.mds, Ca, add = TRUE, col =4)) Рис. 9.14. Ординация видов методом NMDS и изолинии сглаживающих поверхностей для Al и Ca Вероятно, можно предположить (по крайней мере, в рамках имеющейся выборки), что каждая точка ординации, соответствующая тому или иному виду, определяет координаты наилучших природных условий существования этого вида, а прогнозное значение модели (9.4) оценивает его экологический оптимум по сглаживаемому фактору. Функция calibrate(mod, newdata), аналогичная функции predict(), осуществляет с использованием модели mod вычисление для произвольной таблицы newdata прогнозных значений и их стандартных ошибок se.fit.

Оценим экологические оптимумы val содержания алюминия Al в почве для каждого обнаруженного вида растительности. Вычислим также верхние 95%-ные доверительные границы предсказания valC по формуле valC = val + tc*se.fit, где tc – квантиль распределения Стьюдента при = 95%.

Чтобы блокировать появление отрицательных значений, значение внешнего фактора предварительно прологарифмируем.

Al.log -log(varechem$Al) fit.Al - ordisurf(vare.mds, Al.log) a - as.data.frame(vare.mds$species) Sp.Al - cbind(a, calibrate(fit.Al, newdata = a, se.fit = TRUE)) Sp.Al$val - exp(Sp.Al$fit) tc - qt(0.975, nrow(Sp.Al)-1) Sp.Al$valC - exp(Sp.Al$fit)+ exp(Sp.Al$se.fit)*tc head(Sp.Al)

–  –  –

Одним из методов нормирования техногенных загрязнений на исследуемой территории является построение статистических моделей распределения чувствительности видов SSD (Species Sensitivity Distribution – Posthuma et al., 2001). Кривая SSD рассчитывается как интегральная функция некоторого теоретического распределения плотности вероятности встречаемости видов, параметры которого оцениваются по выборке того или иного показателя жизнедеятельности таксономических групп (чаще всего для этого используются показатели токсикометрии NOEC, DС50, LС50 и т.д.).

Поскольку для травянистых растений параметры токсичности отсутствуют, то в качестве предположительно опасных значений содержания алюминия в почве для j-го вида примем настолько большие его величины, которые маловероятны в рамках сглаживающей модели GAM, т.е. столбец valC сформированной таблицы Sp.Al. Если отсортировать этот показатель по его величине, то функция ppoints() сгенерирует последовательность эмпирических вероятностей frac, т.е. долю значений valC, которые не превышают величину valC для каждого j-го вида.

# Сортировка по возрастанию val df - Sp.Al[,5:6] df - df[order(df$val), ] df$frac - ppoints(df$val, 0.5) head(df) val valC frac Hylosple 5.272568 8.87682 0.01136364 Nepharct 6.750738 10.92165 0.03409091 Descflex 6.925953 10.27115 0.05681818 Rhodtome 8.281667 11.68154 0.07954545 Polycomm 10.531262 13.60922 0.10227273 Betupube 12.113196 15.33292 0.12500000 Для аппроксимации данных теоретическим распределением воспользуемся функцией fitdist() из пакета fitdistplus. Наилучшая аппроксимация, процедуру подбора которой мы опускаем, соответствует логнормальному распределению. Найдем параметры этого распределения, т.е.

оценки среднего meanlog и дисперсии sdlog, а также некоторые квантильные значения Al.

library(fitdistrplus) fit - fitdistr(df$valC, "lognormal") ep1s - fit$estimate[1]; ep2s - fit$estimate[2] hcs - qlnorm(c(0.05, 0.1, 0.2), meanlog=ep1s,sdlog=ep2s) meanlog sdlog 4.290294 1.007003 13.928272 20.080861 31.273734 Для построения кривой распределения чувствительности видов (SSD) будем использовать скрипт, представленный в блоге «Data in Environmental Science and Ecotoxicology» (http://edild.github.io/ssd/, автор E. Szcs). Оценка доверительных интервалов теоретической кривой может быть осуществлена с использованием параметрического бутстрепа, который использует предположение о том, что исходные выборочные данные представляют собой случайные реализации вероятностного процесса, определяемого параметрами заданного распределения (Davison, Hinkley, 2006). Подробные комментарии к тексту скрипта см. в нашей книге (Шитиков, 2016).

Определим предварительно две функции, необходимые нам в основных расчетах и выполняющие генерацию псевдо-выборок:

# 1. Функция для нахождения p-квантили случайной выборки из # логнормального распределения с параметрами fit myboot - function(fit, p){ xr - rlnorm(fit$n, meanlog = fit$estimate[1], sdlog = fit$estimate[2]) fitr - fitdist(xr, 'lnorm') hc5r - qlnorm(p, meanlog = fitr$estimate[1], sdlog = fitr$estimate[2]) return(hc5r) } # 2. Функция, возвращающая значения вероятностей для случайной # выборки из логнормального распределения с параметрами fit myboot2 - function(fit, newxs){ xr - rlnorm(fit$n, meanlog = fit$estimate[1], sdlog = fit$estimate[2]) fitr - fitdist(xr, 'lnorm') pyr - plnorm(newxs, meanlog = fitr$estimate[1], sdlog = fitr$estimate[2]) return(pyr) } Выполним построение 1000 бутстреп-кривых из заданного распределения и определим координаты графиков основной функции распределения и ее 95%-ных доверительных огибающих:

set.seed(1234) # Установка генератора случайных чисел require(reshape2) # новые данные для построения плавной кривой newxs - 10^(seq(log10(0.01), log10(max(df$valC)), length.out = 1000)) # получение матрицы для построения 1000 кривых boots - replicate(1000, myboot2(fit, newxs)) bootdat - data.frame(boots) ; bootdat$newxs - newxs bootdat - melt(bootdat, id = 'newxs') # извлечение доверительных интервалов cis - apply(boots, 1, quantile, c(0.025, 0.975)) rownames(cis) - c('lwr', 'upr') # Формирование итоговой таблицы подогнанных значений и cis pdat - data.frame(newxs, py = plnorm(newxs, meanlog = fit$estimate[1], sdlog = fit$estimate[2])) pdat - cbind(pdat, t(cis)) # координаты x для названия видов df$fit - 10^(log10(qlnorm(df$frac, meanlog = fit$estimate[1], sdlog = fit$estimate[2])) -0.4) Осуществим вывод полноценного графика SSD с использованием пакета ggplot2 (рис.

9.15):

library(ggplot2) ggplot() + geom_line(data = bootdat, aes(x = newxs, y = value, group = variable), col = 'steelblue', alpha = 0.05) + geom_point(data = df, aes(x = valC, y = frac)) + geom_line(data = pdat, aes(x = newxs, y = py), col = 'red') + geom_line(data = pdat, aes(x = newxs, y = lwr), linetype = 'dashed') + geom_line(data = pdat, aes(x = newxs, y = upr), linetype = 'dashed') + geom_text(data = df, aes(x = fit, y = frac, label = rownames(df)), hjust = 1, size = 4) + scale_x_log10(breaks = c(5, 10, 50, 100, 300), limits = c(1, max(df$valC))) + labs(x = 'Содержание алюминия в почве', y = 'Доля видов с неблагоприятными условиями') + theme_bw() Рис. 9.15. Распределение чувствительности видов в зависимости от содержания алюминия в почве Если задаться произвольной критической вероятностью на оси ординат SSD (см. рис. 9.15), то с использованием кумулятивной кривой распределения можно оценить величину потенциально опасного уровня воздействия внешнего фактора (Шитиков, 2016). Из представленных результатов следует, что уже при содержании алюминия в почве Al = {13.9,

20.1 и 31.3} из травянистого сообщества могут исчезнуть соответственно доли р = {5, 10, и 20%} видов от их общего числа.

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

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

10. КЛАСТЕРНЫЙ АНАЛИЗ

10.1. Алгоритмы кластеризации, основанные на разделении Алгоритмы неиерахического разделения (partitioning algorithms) осуществляют декомпозицию набора данных, состоящего из n наблюдений, на k групп (кластеров) с заранее неизвестными параметрами. При этом выполняется поиск центроидов – максимально удаленных друг от друга центров сгущений точек Ck с минимальным разбросом внутри каждого кластера. К разделяющим алгоритмам относятся:

метод k средних Мак-Кина (k-means clustering, MacQueen, 1967), в котором каждый из k кластеров представлен центроидом;

разделение вокруг k медоидов или PAM (Partitioning Around Medoids – Kaufman, Rousseeuw, 1990), где медоид – это центроид, координаты которого смещены к ближайшему из исходных объектов данных;

алгоритм CLARA (Clustering Large Applications) – метод, весьма похожий на PAM и используемый для анализа больших наборов данных.

Метод k средних выполняет кластеризацию следующим образом:

1. Назначается число групп (k), на которые должны быть разбиты данные.

Случайно выбирается k объектов исходного набора как первоначальные центры кластеров.

2. Каждому наблюдению присваивается номер группы по самому близкому центроиду, т.е. на основании наименьшего евклидова расстояния между объектом и точкой Ck.

3. Пересчитываются координаты центроидов µk всех k кластеров и вычисляются внутригрупповые разбросы (within-cluster variation) W (Ck ) = ( xi µ k ). Если набор данных включает р переменных, то µk

–  –  –

min, для чего шаги 2 и 3 повторяются многократно, пока назначения групп не прекращают изменяться или не достигнуто заданное число итераций iter.max. Предельное число итераций для минимизации Wtotal, установленное функцией kmeans() по умолчанию, составляет iter.max = 10.

В разделе 2.6 мы уже подробно разбирали пример кластеризации 50 штатов США по криминогенной обстановке на 5 групп с использованием функции kmeans() из пакета cluster. Пример сопровождается скриптами, графической интерпретацией разбиения и приводятся конкретные вычисленные значения вышеприведенных статистик. В связи с этим в дальнейшем мы будем обсуждать только нюансы практического использования метода.

Объединение в кластеры методом k средних – очень простой и эффективный алгоритм, имеющий, однако, две существенные проблемы. Вопервых, итоговые результаты чувствительны к начальному случайному выбору центров групп. Возможное решение этой проблемы состоит в многократном выполнении алгоритма с различным случайным назначением начальных центроидов. Итерация с минимальным значением Wtotal отбирается как конечный вариант кластеризации.

Число таких итераций можно задать параметром nstart функции kmeans():

library(cluster) data("USArrests") df.stand - as.data.frame(scale(USArrests)) c(kmeans(df.stand,centers=5, nstart=1)$tot.withinss, kmeans(df.stand,centers=5, nstart=25)$tot.withinss) [1] 50.72761 48.94420 Нам удалось за 25 повторов алгоритма несколько уменьшить значение критерия оптимальности.

Вторая проблема – необходимость априори задавать фиксированное число кластеров для разбиения, которое, безусловно, далеко не всегда выбирается оптимальным. Поэтому одной из задач кластерного анализа является подбор оптимального значения k, для которой существует несколько версий решения.

Метод "локтя" (elbow method) рассматривает характер изменения разброса Wtotal с увеличением числа групп k. Объединив все n наблюдений в одну группу, мы имеем наибольшую внутрикластерную дисперсию, которая будет снижаться до 0 при k n. На каком-то этапе можно усмотреть, что снижение этой дисперсии замедляется – на графике это происходит в точке, называемой "локтем" (родственник "каменистой осыпи" для анализа главных компонент).

Построить такой график можно в результате прямого перебора, либо с использованием функции fviz_nbclust() из прекрасного пакета factoextra, предназначенного для визуализации результатов кластерного анализа на основе графической системы ggplot2:

k.max - 15 # максимальное число кластеров wss - sapply(1:k.max, function(k) {kmeans(df.stand, k, nstart=10 )$tot.withinss}) plot(1:k.max, wss, type="b", pch = 19, frame = FALSE, xlab="Число кластеров K", ylab="Общая внутригрупповая сумма квадратов") # Формируем график с помощью fviz_nbclust() library(factoextra) fviz_nbclust(df.stand, kmeans, method = "wss") + geom_vline(xintercept = 4, linetype = 2) Ограничимся представлением только одного графика, полученного функцией fviz_nbclust() (рис. 10.2). Два вдумчивых аналитика, безусловно, могут поспорить о месте расположения "максимально согнутого локтевого сустава": для числа групп k = 2 или k = 4.

Рис. 10.1. Выбор оптимального числа кластеров по методу "локтя" Альтернативой субъективно понимаемым графикам "локтя" является использование статистики разрыва (gap statistic – Tibshirani et al., 2001), которые генерируются на основе ресэмплинга и имитационных процедур Монте-Карло. Пусть En {log( Wk* )} обозначает оценку средней дисперсии Wk*, *

–  –  –

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

При сравнительном анализе последовательности значений Gapn(k), k = 2, …, Kmax наибольшее значение статистики соответствует наиболее полезной группировке, дисперсия которой максимально меньше внутригрупповой дисперсии кластеров, собранных из случайных объектов исходной выборки:

set.seed(123) gap_stat - clusGap(df.stand, FUN = kmeans, nstart = 10, K.max = 10, B = 50) # Печать и визуализация результатов print(gap_stat, method = "firstmax") fviz_gap_stat(gap_stat) logW E.logW gap SE.sim [1,] 3.458369 3.633279 0.1749100 0.04956376 [2,] 3.135112 3.370352 0.2352405 0.04192860 [3,] 2.977727 3.234819 0.2570927 0.04085179 [4,] 2.826221 3.120886 0.2946656 0.04153549 [5,] 2.738868 3.022296 0.2834282 0.04377128 [6,] 2.669860 2.935234 0.2653739 0.04201917 [7,] 2.612957 2.855663 0.2427060 0.04029606 [8,] 2.545027 2.782767 0.2377398 0.04031744 [9,] 2.471474 2.715964 0.2444898 0.03998257 [10,] 2.397200 2.654907 0.2577068 0.04074831

Рис. 10.2. График GAP-статистики для выбора оптимального числа кластеров

Параметр FUNcluster использованной выше функции clusGap() указывает на необходимый метод кластеризации и при FUN = pam осуществляется разделение вокруг k медоидов:

set.seed(123) gap_stat - clusGap(df.stand, FUN = pam, K.max = 7, B = 100) print(gap_stat, method = "firstmax") (k.pam - pam(df.stand, k=4)) logW E.logW gap SE.sim [1,] 3.458369 3.632850 0.1744804 0.03737864 [2,] 3.135112 3.380017 0.2449056 0.04476136 [3,] 2.981802 3.250681 0.2688793 0.04601583 [4,] 2.834744 3.141994 0.3072506 0.04709690 [5,] 2.745099 3.049445 0.3043454 0.04729300 [6,] 2.685908 2.962081 0.2761727 0.04687566 [7,] 2.634287 2.886058 0.2517713 0.04560393

Medoids:

ID Murder Assault UrbanPop Rape Alabama 1 1.2425641 0.7828393 -0.5209066 -0.003416473 Michigan 22 0.9900104 1.0108275 0.5844655 1.480613993 Oklahoma 36 -0.2727580 -0.2371077 0.1699510 -0.131534211 New Hampshire 29 -1.3059321 -1.3650491 -0.6590781 -1.252564419 Метод PAM во многом идентичен алгоритму k средних за исключением того, что вместо вычисления центроидов осуществляется поиск k наиболее представительных объектов (или медоидов) среди анализируемых наблюдений. Внутрикластерный разброс оценивается при этом по манхэттенскому, а не евклидовому расстоянию, как в kmeans(). Мы прибегли к простейшему варианту вычислений – использованию функции pam() из пакета cluster, которая выделила в качестве центральных представителей классов Алабаму, Мичиган, Оклахому и Нью-Гэмпшир.

Поскольку параметром функции pam() также является число кластеров k, проверим оптимальность его значения третьим из наиболее популярных методов – по средней ширине силуэта (average silhouette width). Для каждого найденного кластера может быть вычислена "ширина силуэта" b(i ) a (i ), где a(i) – среднее расстояние между объектами i-го кластера, si = max[b(i ), a (i )] b(i) – среднее расстояние от объектов i-го кластера до другого кластера, самого близкого к i-му.

Диаграмму силуэтов можно построить с использованием функции fviz_silhouette() из пакета factoextra:

fviz_silhouette(silhouette(k.pam))

Рис. 10.3. Диаграмма силуэтов при разбиении на 4 кластера

На рис. 10.3 для каждого из 50 штатов показана разность s средних расстояний до объектов своего кластера и чужого кластера, ближайшего к анализируемому объекту. Общее среднее из этих значений Smean определяет качество выполненной кластеризации. Функция fviz_nbclust() выполняет сканирование разбиений с различными значениями k, и строит график зависимости Smean от k, подобный представленным на рис. 10.1-10.2.

fviz_nbclust(df.stand, pam, method = "silhouette")

–  –  –

Метод силуэтов оценил как наилучшее разбиение на два кластера k = 2.

Отметим, что в условиях неопределенности различные алгоритмы могут порождать конкурирующие решения. Колеблющийся читатель может в таких случаях пользоваться функцией pamk() из пакета fpc, которая не требует задавать число классов, а оценивает его самостоятельно.

Алгоритм CLARA мало чем отличается от метода PAM, но приспособлен для обработки больших массивов данных (миллионы наблюдений и более). Соответствующая функция clara() из пакета cluster специально оптимизирована для сокращения времени расчетов и рационального использования оперативной памяти.

Значительный интерес представляет построение двумерных диаграмм (ординационных "биплотов") распределения наблюдений по кластерам, которые формируются с предварительным приведением исходного пространства признаков к двум главным компонентам (см. рис. 2.13 из раздела 2.6). Обычно для этого используется функция clustplot(), но мы предлагаем вниманию читателей функцию fviz_cluster() из пакета factoextra, которая использует для создания диаграмм графическую систему ggplot2.

Эта функция может быть использована для визуализации результатов по методам k средних, PAM, CLARA и Fanny.

Ее простейший формат имеет вид:

fviz_cluster(object, data = NULL, stand = TRUE, geom = c("point", "text"), frame = TRUE, frame.type = "convex") где:

object: объект класса "partition", созданный функциями pam(), clara() или fanny() из пакета cluster, или объект, возвращаемый функцией kmeans();

data: таблица с исходными данными для кластеризации (этот аргумент необходим только, если используется объект kmeans);

stand = TRUE: данные стандартизуются перед выполнением анализа главных компонент;



Pages:     | 1 |   ...   | 2 | 3 || 5 |



Похожие работы:

«1. ЦЕЛИ И ЗАДАЧИ ПРЕДДИПЛОМНОЙ ПРАКТИКИ Целью преддипломной практики является закрепление и углубление теоретической подготовки студента и приобретение им практических навыков и компетенций в сфере профессиональной деятельности, при...»

«ISSN 2411–6602 (Online) ISSN 1607–2855 (Print) С. 27 – 37 Том 12 • № 1 • 2016 УДК 523.4 Периодические изменения активности процессов в атмосфере Юпитера А.П. Видьмаченко Главная астрономическая обсерватория НАН Украины Вариации йовимагнитной широты Земл...»

«Том 7, №1 (январь февраль 2015) Интернет-журнал "НАУКОВЕДЕНИЕ" publishing@naukovedenie.ru http://naukovedenie.ru Интернет-журнал "Науковедение" ISSN 2223-5167 http://naukovedenie.ru/ Том 7, №1 (2015) http://naukovedenie.ru/index.php?p=vol7-1 URL статьи: http://naukovedenie.ru/PDF/106TVN115.pdf DOI: 10....»

«4 Аннотация В данной диссертационной работе предложена методика обессоливания водных растворов низкомолекулярных солей, основанная на использовании полимерных гидрогелей в режиме сжатия набухания. Показана разработка системы, способной действовать...»

«СЛОВАРЬ КВИЖНОЙ МАЛОРУССКОЙ Р'М И ПО РУКОПИСИ ХЛ Б. 17 Выяснятся— озаряюся. Высокій розумъ— мудрованіе. В даніе— свденіе. Высокодгішкурую— великому дрВдаю— вмъ, спмъ, разумю. ствую, велерчую. В едом ст во— вжство. Высокой мысли— высокоумный. Вдомость даю— возвщаго. Высоко о себ розумю— высоВдомща— свдитель. комудрствую...»

«2015, № 2 ВОПРОСЫ РУССКОЙ ЛИТЕРАТУРЫ, УДК: 82-14 (477) ПЕЙЗАЖНЫЙ ДИСКУРС РУССКОЙ ЛИРИКИ КОНЦА 1960-Х – СЕРЕДИНЫ 1980-Х ГОДОВ: ТЕОРЕТИКО-МЕТОДОЛОГИЧЕСКИЙ ПОДХОД Остапенко Ирина Владимировна, д. филол. н., профессор кафедры русской и...»

«Карапец Мария Валерьевна, Черноситова Татьяна Леонидовна КОМПАРАТИВНЫЙ АНАЛИЗ ПЕРЕВОДОВ ПОЭМЫ РЕКВИЕМ А. АХМАТОВОЙ НА АНГЛИЙСКИЙ И ФРАНЦУЗСКИЙ ЯЗЫКИ В статье на материале сравнительного анализа трех переводов на английский и французский языки отрывка из Реквиема А. Ахматовой рассматривается проблема коммуникативной эквивалентности. Показывается,...»

«ДЕВЯТНАДЦАТЫЙ АРБИТРАЖНЫЙ АПЕЛЛЯЦИОННЫЙ СУД ПОСТАНОВЛЕНИЕ от 10 июня 2009 г. по делу N А48-182/2009 Резолютивная часть постановления объявлена 05 июня 2009 г. Постановление в полном объеме изготовлено 10 июня 2009 г.Девятнадцатый арб...»

«МИНИСТЕРСТВО ОБРАЗОВАНИЯ РОССИЙСКОЙ ФЕДЕРАЦИИ Государственное образовательное учреждение высшего профессионального образования "Оренбургский государственный университет" И.Т. КАЗАРМЩИКОВ ПРОИЗВОДСТВО МЕТАЛЛИЧЕСКИХ КОНСТРУКЦИОННЫХ МАТЕРИАЛОВ Рекомендовано Учебным Советом Оренбургского Го...»

«Сводка изменений Руководства кандидата Отражены изменения с версии от 19 сентября 2011 года до версии от 11 января 2012 года Раздел Тема Изменение текста Обоснование и комментарии Модуль 1 1.1.2.4 Заблаговременное Одновременно с проведением 60-дневного периода Добавлено уточнение для указания на то, что датой предупреждение со комме...»

«Министерство образования Республики Беларусь Учебно-методическое объединение по образованию в области горнодобывающей промышленности УТВЕРЖДАЮ Первый заместитель Министра образования Ред]^лики Беларусь РегйстраіЬгонный № ТД-і. •/ОбУтип. МЕТАЛЛИЧЕСКИЕ ПОЛЕЗНЫЕ ИСКОПАЕМЫ...»

«UNW/2013/3 Организация Объединенных Наций Исполнительный совет Distr.: General Структуры Организации Russian Объединенных Наций Original: English по вопросам гендерного равенства и расширения прав и возможностей же...»

«Viasat History Russia – Неделя 51, 2012 Понедельник 17 Декабрь 10:00 Команда времени (4) (Time Team YR X) Великобритания, документальный сериал, 2003 Эта серия: "Команда времени" берет курс на Шетландские острова, где ей предстои...»

«БОКЛЕВСКИЙ массовая библиотека „ИСКУССТВО Н. НИКИФОРАКИ ПЕТР МИХАЙЛОВИЧ БОКЛЕВСКИЙ 1816-1897 ГОСУДАРСТВЕННОЕ ИЗДАТЕЛЬСТВО П. И. Чичикоп "ИСКУССТВО" Иллюстрация к поэме Н. В. Гоголя "Мертвые души" Москеа 79,52 Петр Михайлович...»

«Опубликовано: “25 ”09. 2013 ИНФОРМАЦИОННАЯ БРОШЮРА С р о ч н ы й в к л ад ф и з ич еск и х л и ц К о н в ер т и р уе мы й Описание вклада, основные положения 1.1.1. Вклад вносится физическими лицами – резидентами или нерезидентами.1.2. Вклад с неограниченным количеством пополнений...»

«ДОГОВОР НОМИНАЛЬНОГО СЧЕТА (заключается с физическими лицами опекунами/попечителями/родителями) г. Магнитогорск "Кредит Урал Банк" (Акционерное общество), именуемый в дальнейшем "Банк", с одной стороны, и Владелец счета, реквизиты которого указаны в пункте 10.1. настоящего договора, с другой стороны, заключили насто...»

«ОБОСТРЕНИЕ ГАЗОПРОВОДНОЙ КОНКУРЕНЦИИ НА КАСПИЙСКО-ЧЕРНОМОРСКОЙ ЭНЕРГЕТИЧЕСКОЙ СЦЕНЕ: "ЮЖНЫЙ ПОТОК" ПРОТИВ "НАБУККО" А.К. Магомедов доктор политических наук, профессор заведующий кафедрой "Связи с общест...»

«Рабочая программа дисциплины Б1.В.ДВ.5.2 "Профессиональное консультирование в профориентации" Содержание стр.1. Цели и задачи дисциплины 3 2. Место дисциплины в структуре ОПОП. 3 3. Требования к результатам освоения дисциплины 3 4. Объем дисциплины и виды учебной работы 4 5. Содерж...»

«МИНОБРНАУКИ РОССИИ ФЕДЕРАЛЬНОЕ ГОСУДАРСТВЕННОЕ БЮДЖЕТНОЕ ОБРАЗОВАТЕЛЬНОЕ УЧРЕЖДЕНИЕ ВЫСШЕГО ПРОФЕССИОНАЛЬНОГО ОБРАЗОВАНИЯ "САМАРСКИЙ ГОСУДАРСТВЕННЫЙ АЭРОКОСМИЧЕСКИЙ УНИВЕРСИТЕТ ИМЕНИ АКАДЕМИКА С.П.КОРОЛЕВА (НАЦИОНАЛЬНЫЙ ИССЛЕДОВАТЕЛЬСКИЙ УНИВЕРСИТ...»

«CLASIC VI Relax раскладной диван Диваны раскладные Мягкая мебель, диваны Мебель PREMIUM класса. Высокое качество. Европейские ткани. Деревянная основа. Европейский сертификат. mebeles.buv.l...»

«СОВЕТ ПО ЖЕЛЕЗНОДОРОЖНОМУ ТРАНСПОРТУ ГОСУДАРСТВ-УЧАСТНИКОВ СОДРУЖЕСТВА Утверждены Советом по железнодорожному транспорту государствучастников Содружества, Протокол от 21-22 мая 2009 г. № 50 ПРАВИЛА ПЕРЕВОЗОК ЖИДКИХ ГРУЗОВ НАЛИВОМ В ВАГОНАХ-ЦИСТЕРНАХ И ВАГОНАХ БУНКЕРНОГО ТИПА ДЛЯ ПЕРЕВОЗКИ НЕ...»

«Русский баталист Василий Васильевич Верещагин (1842 — 1904) Библиографический список материалов к выставке I. Книги: 1) Библиография русского художника В. В. Верещагина на Кубе (1889 — 1992).Череповец, 1992. 48 с. КХ, Шифр Л 277/649 2) Булгаков Ф. И. В. В. Верещагин и его произведения. СПб., 1905. с. КХ, Шифр А II 6/67, РК...»

«РАСШИРЕНИЕ: новости Ключ расширения : tt_news Copyright 2000-2005, Rupert Germann, rupi@gmx.li Этот документ публикуется в соответствии с Open Content License доступной на http://www.opencontent.org/opl.shtml Содержимое этого документа относится к TYPO3 GNU/GPL CMS/Framework доступной с www.typo3.com Оглавление Ключевые...»

«Приложение к Решению Совета депутатов Муниципального образования Переволоцкий район Оренбургской области от 17.06.2015г. № 387 ВНЕСЕНИЕ ИЗМЕНЕНИЙ В ГЕНЕРАЛЬНЫЙ ПЛАН МУНИЦИПАЛЬНОГО ОБРАЗОВАНИЯ "АБРАМОВСКИЙ СЕЛЬСОВЕТ" ПЕРЕВОЛОЦКОГО РАЙОНА ОРЕНБУРГСКОЙ ОБЛАСТИ ТОМ 2 МАТЕРИАЛЫ ПО ОБОСНОВАНИЮ Заказчик: ЗАО "Центр наукоёмких те...»

«АДМИНИСТРАЦИЯ МАГАДАНСКОЙ ОБЛАСТИ ПОСТАНОВЛЕНИЕ от 26 июля 2012 года N 531-па Об областной целевой программе Развитие малого и среднего предпринимательства в Магаданской области на 2013-2016 годы (Текст изменен) {Изменено: Постановлением о...»

«УТВЕРЖДЕНО Общим собранием акционеров ЗАО "Монокристалл" "30" июня 2010 г. (Протокол от "05" июля 2010 г.) ПОЛОЖЕНИЕ О комитете Совета директоров открытого акционерного общества "Монокристалл" ПО АУДИТУ (редакция 1) г. Ставрополь 2010 год Страница 1 из 7 Настоящ...»

«Илья РЕЙДЕРМАН У моря жить, у моря умереть. Мой город.на выученный наизусть бульвар. Иосиф БРОДСКИЙ Да, выученный наизусть бульвар. Гуляй с женою, тещею, детьми. Морской вокзал, гостиницы кошмар и море, что уморено людьми. Вдруг понимаешь, как тоскует Дюк, что повернулся к городу спиной. Должно быть,...»

«Вступительное слово председателей: на пути к окончательному варианту предложения 22 апреля 2015 г. Сегодня рабочая группа ICANN, состоящая из представителей нескольких сообществ и о...»








 
2017 www.book.lib-i.ru - «Бесплатная электронная библиотека - электронные ресурсы»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.