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

Pages:     | 1 || 3 | 4 |   ...   | 5 |

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

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

Если ожидается, что тот или иной метод построения предсказательной модели чувствителен к наличию многомерных выбросов, исходные данные можно преобразовать при помощи метода пространственных признаков (spatial sign – Serneels et al., 2006). С математической точки зрения, этот метод проецирует значения предикторов на поверхность многомерной сферы, в результате чего отдельные наблюдения становятся нормированными m относительно центра этой сферы: X ij = X ij / X ij, где m – это число ' 2

i =1

преобразуемых предикторов.

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

preProcess(х, method = c("center", "scale", "spatialSign")) чем будет достигнут примерно одинаковый их вклад в величину квадрата расстояния. Кроме того, важно уточнить, что поскольку это преобразование применяется одновременно к группе предикторов, то последующее изменение их состава может исказить смысл проведенной трансформации.

Как было показано в разделе 2.4, одним из наиболее популярных методов снижения размерности исходного набора данных является анализ главных компонент PCA (principle components analysis). Этот метод позволяет сформировать ортогональную систему координат, обеспечивающую оптимальное расположение точек объектов относительно осей главных компонент нового редуцированного пространства.

Пересчет из исходной системы координат в p-мерное (m p) пространство главных компонент осуществляется с использованием m линейных ортогональных преобразований переменных Xi: Tk = Pik ( X i µ i ), i =1 где Pik – нагрузки (компоненты собственного вектора, соответствующего k-му собственному значению), k = 1, 2, …, p. Основной критерий качества такого преобразования – достаточно высокая доля общей вариации исходных данных, объясняемая p первыми собственными значениями (например, 80% или 95%).

Для реализации PCA-процедуры с помощью preProcess() необходимо определить параметр thresh – кумулятивную долю дисперсии исходных данных, которая должна содержаться в главных компонентах (по умолчанию thresh = 0.95) или pcaComp – максимальное число главных компонент (по умолчанию этот параметр имеет значение NULL, т.е. он выключен, и оптимальное число главных компонент выбирается на основе thresh).

Построим две модели логистической регрессии на основе "очищенных" y - factor(GermanCredit$Class) gcred - GermanCredit[, -10] nz - nearZeroVar(gcred) gcred.clean - gcred[, -nz] highCor - findCorrelation(cor(gcred.clean), cutoff = 0.75) gcred.clean - gcred.clean[, -highCor] linCombo - findLinearCombos(gcred.clean) gcred.clean - gcred.clean[, -linCombo$remove] dim(gcred.clean) (см.

раздел 3.2) данных таблицы GermanCredit и после РСА-преобразования:

[1] 1000 43 Определим оптимальную размерность пространства главных компонент по критерию Кайзера-Гуттмана, который рекомендует оставить только те главные компоненты, собственные значения которых превышают их среднее library(vegan) mod.pca - rda(gcred.clean, scale=TRUE) ev - mod.pca$CA$eig # Иллюстрация Критерия Кайзера-Гуттмана barplot(ev, col="bisque", las=2) abline(h=mean(ev), col="red") legend("topright", "Средние собственных значений", lwd=1, col=2, bty="n") c(length(ev[ev mean(ev)]), sum(ev[ev mean(ev)])/sum(ev)) (рис. 3.2):

[1] 19.0000000 0.7267951

Рис. 3.2. Оценка числа главных компонент по критерию Кайзера-Гуттмана





Сформируем новую таблицу предикторов из 19 главных компонент, которые объясняют 72.7% общей дисперсии:

prePCA - preProcess(gcred.clean, method= c("center", "scale", "pca"), pcaComp=19) gcred.pca - predict(prePCA,gcred.clean) Попробуем сравнить точность прогноза двух моделей логистической регрессии с использованием функции train(), которая будет предметом нашего подробного рассмотрения ниже. Для тестирования моделей будем многократно (times = 100) случайным образом делить всю выборку на обучающую (800 объектов или 80%) и контрольную (200 объектов или 20%), для чего с помощью функции createDataPartition() создадим соответствующую Метод тестирования "заготовку" train.index.

method="LGOCV" (многократное разбиение на обучающую и контрольную выборки) и train.index определим в специальном объекте trControl.

На функцию train() подадим данные для построения модели, тип модели и условия тестирования:

train.index - createDataPartition(y, p =.8, times = 100) trControl = trainControl(method="LGOCV", index = train.index) print("Модель на основе исходного набора из 43 предикторов") modSource - train(gcred.clean, y, "glm", family=binomial, trControl = trControl) print("Модель на основе 19 главных компонент") (modPCA - train(gcred.pca, y, "glm", family=binomial, trControl = trControl)) [1] Модель на основе исходного набора из 43 предикторов’ Resampling: Repeated Train/Test Splits (100 reps, 0.75%) Resampling results Accuracy Kappa Accuracy SD Kappa SD 0.753 0.37 0.0251 0.0603 [1] Модель на основе 19 главных компонент Resampling results Accuracy Kappa Accuracy SD Kappa SD 0.749 0.362 0.0224 0.0554 Нельзя сказать, что преобразовав переменные, мы получили более точную модель, но зато обошлись значительно меньшим числом переменных.

При значении параметра method="ica" функции preProcess() используется другой метод снижения размерности пространства переменных

– анализ независимых компонент ICA (Independent Component Analysis), который выполняет ту же задачу, что и РСА, однако основан на несколько иных математических концепциях и процедурах. В частности, цель ICA состоит в нахождении таких новых компонент (нового пространства латентных переменных), которые взаимно независимы в полном статистическом смысле.

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

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

3.3):

plot(varImp(modSource, scale = FALSE))

–  –  –

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

–  –  –

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

Проблема "борьбы с пропусками" столь же сложна, как и сама статистика, поскольку в этой области существует впечатляющее множество подходов. В русскоязычных книгах по использованию R (Кабаков, 2014;

Мастицкий, Шитиков, 2015) бегло представлены только некоторые функции пакета mice, который, несмотря на свою "продвинутость", мало удобен для практической работы с данными умеренного и большого объема. Хорошей альтернативой являются методы и "knnImpute", "bagImpute" "medianImpute" функции preProcess() из пакета caret, которую мы рассмотрели в разделе 3.3 как инструмент для трансформации данных.

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

Каждое из 200 наблюдений содержит информацию о 18 переменных, в том числе:

три номинальных переменных, описывающих размеры size = c("large", "medium", "small") и скорость течения реки speed = c("high", "low", "medium"), а также время года season = c("autumn", "spring", "summer", "winter"), сопряженное с моментом взятия проб;

8 переменных, составляющих комплекс наблюдаемых гидрохимических показателей: максимальное значение рН mxPH (1), минимальное содержание кислорода mnO2 (2), хлориды Cl (10), нитраты NO3 (2), ионы аммония NH4 (2), орто-фосфаты oPO4 (2), общий минеральный фосфор PO4 (2) и количество хлорофилла а Chla (12) (в скобках приведено число пропущенных значений);

средняя численность каждой из 7 групп водорослей a1-a7 (видовой состав не идентифицировался).

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

3.4):

library(DMwR) library(ggplot2) summary(algae) # вывод не приводится ggplot(data=algae[!is.na(algae$mnO2),], aes(speed, mnO2)) + geom_violin(aes(fill = speed), trim = FALSE, alpha = 0.3) + geom_boxplot(aes(fill = speed), width = 0.2, outlier.colour = NA) + theme(legend.position = "NA")

–  –  –

На рис. 3.4 мы получили так называемую "скрипичную диаграмму" (violin plot), которая объединяет в себе идеи диаграмм размахов и кривых распределения вероятности. Суть достаточно проста: продольные края "ящиков с усами" (для сравнения приведены тоже) замещаются кривыми плотности вероятности. В итоге, например, легко можно выяснить не только тот факт, что в потоках с быстрым течением (high) содержание кислорода выше, но и ознакомиться с характером распределения соответствующих значений.

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

Такие диаграммы, или "панели" (от англ. panels, facets или multiples), определенным образом упорядочиваются и размещаются на одной странице.

Из графиков, представленных на рис. 3.5, легко увидеть, что численность водорослей группы а1 падает с увеличением концентрации фосфатов.

qplot(PO4, a1, data = algae[!is.na(algae$PO4), ]) + facet_grid(facets = ~ season)+ geom_smooth(color="red", se = FALSE) Рис. 3.5.

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

algae[!is.na(algae$PO4), ].

Если в обрабатываемой таблице обнаружены недостающие данные, то в общих чертах можно избрать одну из следующих возможных стратегий:

удалить строки с неопределенностями;

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

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

постараться обходить эту неприятную ситуацию, используя, например, формальный параметр na.rm некоторых функций.

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

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

# Число строк с пропущенными значениями nrow(algae[!complete.cases(algae),]) # Их удаление algae - na.omit(algae) [1] 16

–  –  –

61 6.50 10.4 5.970 10.0000 2.000 14.0000 5.475 32.730 62 6.40 2.675 103.1665 40.150 14.0000 5.475 9.8 32.730 63 7.83 11.7 4.083 1.328 18.0000 3.333 6.6670 5.475 116 9.70 10.8 0.222 0.406 10.0000 22.444 10.1110 5.475 161 9.00 5.8 0.900 142.0000 102.000 186.0000 68.050 32.730 184 8.00 10.9 9.055 0.825 40.0000 21.083 56.0910 5.475 199 8.00 7.6 32.730 2.675 103.1665 40.150 103.2855 5.475 Альтернативой "наивному" подходу является учет структуры связей между переменными. Например, можно воспользоваться тем, что между двумя формами фосфора oPO4 и PO4 существует тесная корреляционная связь. Тогда, например, можно избавиться от некоторых неопределенностей в показателе PO4, вычислив его пропущенные значения по уравнению простой регрессии:

data(algae) lm(PO4 ~ oPO4, data = algae)

Coefficients:

(Intercept) oPO4 42.897 1.293 # Функция вывода значений PO4 в зависимости от оPO4 fillPO4 - function(oP) {if (is.na(oP)) return(NA) else return(42.897 + 1.293 * oP) } # Восстановление пропущенных значений PO4 algae[is.na(algae$PO4), 'PO4'] sapply(algae[is.na(algae$PO4), 'oPO4'], fillPO4) algae[ind, 10] [1] 48.069 3.000 6.000 6.500 1.000 4.000 6.000 11.000 6.000 [10] 14.000 14.000 6.667 10.111 186.000 56.091 NA Одно из пропущенных значений удалось восстановить.

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

Использование метода "bagImpute" осуществляет для каждой из имеющихся переменных построение множественной бутстреп-агрегированной модели, или бэггинг-модели (bagging), на основе деревьев регрессии, используя все остальные переменные в качестве предикторов.

Этот метод мудр и точен, но требует значительных затрат времени на вычисление, особенно при работе с данными большого объема:

data(algae) pPbI - preProcess(algae[, 4:11], method = 'bagImpute') algae[, 4:11] - predict(pPbI, algae[, 4:11]) (Imp.Bag - algae[ind, 4:11]) Наконец, третий метод функции preProcess() для заполнения пропусков – "knnImpute" – основан на простейшем, но чрезвычайно эффективном алгоритме k ближайших соседей (k-nearest neighbours) или kNN.

В основе метода kNN лежит гипотеза о том, что тестируемый объект d будет иметь примерно тот же набор признаков, что и обучающие объекты в локальной области его ближайшего окружения (рис. 3.6).

Рис. 3.6. Интерпретация метода k ближайших соседей Если речь идет о классификации, то неизвестный класс объекта определяется голосованием k его ближайших соседей (на рис. 3.6 k = 5).

kNN-регрессия оценивает значение неизвестной координаты Y, усредняя известные ее величины для тех же k соседних точек.

Одна из важных проблем kNN – выбор метрики, на основе которой оценивается близость объектов. Наиболее общей формулой для подсчета расстояния в m-мерном пространстве между объектами X1 и X2 является мера

Минковского:

DS(X1, X2) = [ |x1i – x2i|p ]1/r, где i изменяется от 1 до m, а r и p – задаваемые исследователем параметры, с помощью которых можно осуществить нелинейное масштабирование расстояний между объектами. Мера расстояния по Евклиду получается, если принять в метрике Минковского r = p = 2, и является, по-видимому, наиболее общей мерой расстояния, знакомой всем со школы по теореме Пифагора. При r = p = 1 имеем "манхэттенское расстояние" (или "расстояние городских кварталов"), не столь контрастно оценивающее большие разности координат x.

Вторая проблема метода kNN заключается в решении вопроса о том, на мнение какого числа соседей k нам целесообразно положиться? В свое время мы обсудим этот вопрос детально, а сейчас будем ориентироваться на значение k = 5, используемое по умолчанию:

data(algae) pPkI - preProcess(algae[, 4:11], method = 'knnImpute') alg.stand - predict(pPkI, algae[, 4:11]) Получив в результате применения predict() матрицу переменных с пропущенными значениями, заполненными этим методом, мы с удивлением обнаруживаем, что данные оказались стандартизованными (т.е.

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

m - pPkI$mean sd - pPkI$std algae[, 4:11] - t(apply(alg.stand, 1, function (r) m + r * sd)) (Imp.Knn - algae[ind, 4:11]) Наконец, зададимся следующим закономерным вопросом: а какой метод лучше? Обычно эта проблема не имеет теоретического решения, и исследователь полагается на собственную интуицию и опыт. Но мы можем оценить, насколько расходятся между собой результаты, полученные каждым способом заполнения. Для этого сформируем блок данных из 3 * 16 = 48 строк исходной таблицы с пропусками, заполненными тремя методами ("Med", "Bag", "Knn"), и выполним снижение размерности пространства переменных методом главных компонент в двумерное (см. раздел 2.4).

Посмотрим, как "лягут карты" на плоскости:

ImpVal - rbind(Imp.Med, Imp.Knn) ImpVal - rbind(ImpVal, Imp.Bag) Imp.Metod - as.factor(c(rep("Med", 16), rep("Knn", 16), rep("Bag", 16))) library(RANN) library(vegan) Imp.M - rda(ImpVal ~ Imp.Metod, ImpVal) plot(Imp.M, display = "sites", type = "p") ordihull(Imp.M, Imp.Metod, draw = "polygon", alpha = 67, lty = 2, col = c(1, 2, 3), label = TRUE) На рис. 3.7 мы выделили контуром (hull), проведенным через крайние точки, области каждого из трех блоков данных и поместили метку метода в центры тяжести полученных многоугольников. Понятно, что медианное заполнение характеризуется меньшей вариацией результатов, поскольку игнорирует специфичность свойств каждого объекта. Оба других метода, учитывающих внутреннюю структуру данных, дали приблизительно похожие результаты.

Рис. 3.7. Ординационная диаграмма блоков данных таблицы algae с пропущенными значениями, заполненными разными способами

–  –  –

Как подчеркивалось ранее в разделе 3.1, пакет caret был разработан как эффективная надстройка, позволяющая унифицировать и интегрировать использование множества различных функций и методов построения моделей классификации и регрессии, представленных в других пакетах R. При этом происходит всестороннее тестирование и оптимизация настраиваемых коэффициентов и гиперпараметров моделей. Эта разработанная единая технология реализована в функции использующей train(), полуавтоматические интеллектуальные подходы и универсальные критерии качества с применением алгоритмов ресэмплинга.

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

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

trainControl(method = "boot", number = ifelse(grepl("cv", method), 10, 25), p = 0.75, repeats = ifelse(grepl("cv", method), 1, number), search = "grid", initialWindow = NULL, horizon = 1, fixedWindow = TRUE, verboseIter = FALSE, returnData = TRUE, returnResamp = "final", savePredictions = FALSE, classProbs = FALSE, summaryFunction = defaultSummary, selectionFunction = "best", seeds = NA, preProcOptions = list(thresh = 0.95, ICAcomp = 3, k = 5), sampling = NULL, index = NULL, indexOut = NULL, timingSamps = 0, predictionBounds = rep(FALSE, 2), adaptive = list(min = 5, alpha = 0.05, method = "gls", complete = TRUE), trim = FALSE, allowParallel = TRUE)

Остановимся на описании наиболее важных аргументов этой функции:

метод ресэмплинга: "boot", "boot632", "cv", method – "repeatedcv", "LOOCV", "LGOCV" (для повторяющихся разбиений на обучающую и проверочную выборки), "none" (исследуется только модель на обучающей выборке), "oob" (для таких алгоритмов, как случайный лес, бэггинг деревьев и др.), "adaptive_cv", "adaptive_boot" или "adaptive_LGOCV";

number – задает число итераций ресэмплинга, в частности, число блоков (folds) при перекрестной проверке;

repeats – число повторностей (только для k-кратной перекрестной проверки);

p – доля обучающей выборки от общего объема при выполнении перекрестной проверки;

verboseIter – TRUE означает, что в ходе вычислений caret будет показывать, на каком этапе они находятся (это удобно для оценки оставшегося времени вычислений, которое часто бывает очень большим);

search – способ перебора параметров модели (по сетке – "grid", или случайным назначением – "random");

returnResamp и savePredictions – определяют условия сохранения результатов выполнения ресэмплинга и прогнозируемых значений, т.е.

"none", только для итоговой модели "final" или все "all";

classProbs – TRUE означает, что в процессе вычислений алгоритм будет сохранять данные о вероятностях попадания объекта в каждый класс, а не только конечные метки класса;

summaryFunction – определяет функцию, которая вычисляет метрику качества модели при ресэмплинге;

selectionFunction – определяет функцию выбора оптимального значения настраиваемого параметра;

preProcOptions – список опций, который передается функции предобработки данных preProcess().

Например, при создании объекта ctrl ctrl - trainControl(method = "repeatedcv", number = 10, repeats = 10) параметры перекрестной проверки будут иметь следующий смысл:

method = ‘repeatedcv’ означает, что будет использоваться повторная перекрестная проверка (также возможна перекрестная проверка, leaveone-out проверка, и т.д.);

number = 10 означает, что в процессе перекрестной проверки выборку надо разбивать на 10 равных частей;

repeats = 10 означает, что повторная перекрестная проверка должна быть запущена 10 раз.

Теперь перейдем непосредственно к описанию функции train(), которая имеет следующий формат вызова:

train(x, y, method = "rf", preProcess = NULL, weights = NULL, metric = ifelse(is.factor(y), "Accuracy", "RMSE"), maximize = ifelse(metric %in% c("RMSE", "logLoss"), FALSE, TRUE), trControl = trainControl(), tuneGrid = NULL, tuneLength = 3) Исходные данные задаются, как всегда, либо матрицей предикторов x и вектором отклика y, либо объектом formula с указанием таблицы данных data. Следующий аргумент – method – это, в сущности, название модели классификации или регрессии, которую необходимо построить и протестировать. Если выполнить команду names(getModelInfo()), то можно увидеть список из 233 доступных методов (это количество постоянно расширяется).

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

http://topepo.github.io/caret/modelList.html или http://topepo.github.io/caret/bytag.html Можно просмотреть названия всех моделей, которые имеют, например, отношение к линейной регрессии:

library(caret) ls(getModelInfo(model = "lm")) [1] "bayesglm" "elm" "glm" "glmboost" "glmnet" [6]"glmStepAIC" "lm" "lmStepAIC" "plsRglm" "rlm" С каждым методом связан набор гиперпараметров, подлежащих оптимизации.

Можно убедиться в том, что простая линейная регрессия (method = "lm") не имеет параметров для настройки:

modelLookup("lm") model parameter label forReg forClass probModel 1 lm parameter parameter TRUE FALSE FALSE В свою очередь, деревья решений rpart, которые мы рассмотрим позднее, имеют один параметр для настройки – Complexity Parameter (аббревиатура – cp):

modelLookup("rpart") model parameter label forReg forClass probModel 1 rpart cp Complexity Parameter TRUE TRUE TRUE Из сообщений функции modelLookup() можно также увидеть, что линейная регрессия не используется для классификации (forClass= FALSE), тогда как rpart (от "Recursive Partitioning and Regression Trees") можно применять как для построения деревьев регрессии, так и классификации. В последнем случае модель осуществляет не только предсказание класса, но и оценивает апостериорные вероятности (probModel = TRUE).

Метод перекрестной проверки, заданный объектом trControl = trainControl(), обеспечивает сканирование настраиваемых параметров и оценку эффективности модели на каждой итерации по определенным критериям качества. При построении каждой частной модели предварительно осуществляется предобработка данных с использованием методов, перечисленных в preProcess (и с учетом опций preProcOptions объекта trControl).

По умолчанию аргумент metric использует в качестве критерия качества точность предсказания ("Accuracy") в случае классификации (см.

пример, рассмотренный в разделе 3.3) и корень из среднеквадратичного отклонения прогнозируемых значений от наблюдаемых ("RMSE") для регрессии. Логический аргумент maximize уточняет, должен ли этот критерий быть максимизирован или минимизирован. Другие значения metric вместе с различными значениями аргументов summaryFunction и объекта обеспечивают широкие selectionFunction trControl возможности пользовательского назначения метрик для оценки эффективности моделей.

Количество перебираемых значений настраиваемого параметра задается аргументом tuneLength. Например, чтобы задать 30 повторов оценки параметра Ср модели rpart вместо исходных трех, необходимо указать tuneLength=30. Другой вариант – определить последовательность этих значений в списке, задающем сетку: например, tuneGrid = expand.grid(.cp = 0.5^(1:10)), если заранее известен диапазон, к которому принадлежит оптимизируемое значение.

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

ls(mytrain) [1]"bestTune" "call" "coefnames" "control" [5]"dots" "finalModel" "maximize" "method" [9]"metric" "modelInfo" "modelType" "perfNames" [13]"pred" "preProcess" "resample" "resampledCM" [17]"results" "terms" "times" "trainingData" [21]"xlevels" "yLimits" Приведем краткий пример на тему подбора полиномиальной регрессии для описания зависимости электрического сопротивления (Ом) мякоти фруктов киви от процентного содержания в ней сока, который подробно разбирался нами в разделе 2.1. Позже (раздел 2.2) мы показали, как найти оптимальную степень полинома d = 4 с использованием написанной нами функции скользящего контроля.

К сожалению, функция train() селекцию предикторов для метода "lm" не выполняет и оптимальную степень полинома не настраивает. Но мы можем выполнить любую перекрестную проверку для оценки динамики изменения критериев качества модели – среднеквадратичной ошибки RMSE и среднего коэффициента детерминации RSquared (рис.

3.8):

library(DAAG) data("fruitohms") set.seed(123) max.poly - 7 degree - 1:max.poly RSquared - rep(0,max.poly) RMSE - rep(0,max.poly) # Выполним 10-кратную кросс-проверку с 10 повторностями fitControl - trainControl(method = "repeatedcv", number = 10,repeats = 10) # Тестируем модель для различных степеней for ( d in degree) { f - bquote(juice ~ poly(ohms,.(d))) LinearRegressor - train(as.formula(f), data = fruitohms, method = "lm",trControl =fitControl) RSquared[d] - LinearRegressor$results$Rsquared RMSE[d]- LinearRegressor$results$RMSE } library(ggplot2) Degree.RegParams = data.frame(degree,RSquared,RMSE) ggplot(aes(x = degree,y = RSquared), data = Degree.RegParams) + geom_line() ggplot(aes(x = degree,y = RMSE), data = Degree.RegParams) + geom_line() Рис. 3.8. Поиск степени функции полиномиальной регрессии с использованием функции train() В отличие от ранее проведенных расчетов, минимум ошибки и максимум коэффициента детерминации имеют место при d = 5.

Если не задавать непосредственно объект trControl, то по умолчанию вместо кросс-проверки функция train() осуществляет бутстреппинг критериев качества подгонки RMSE и RSquared:

Poly5 - train(ohms ~ poly(juice,5), data = fruitohms, method = "lm") summary(Poly5$finalModel) summary(lm(ohms ~ poly(juice, 5), data = fruitohms))

Coefficients:

Estimate Std. Error t value Pr(|t|) (Intercept) 4360.0 83.2 52.406 2e-16 *** poly(juice, 5)1 -16750.2 941.3 -17.796 2e-16 *** poly(juice, 5)2 4750.6 941.3 5.047 1.59e-06 *** poly(juice, 5)3 3945.6 941.3 4.192 5.27e-05 *** poly(juice, 5)4 -3371.2 941.3 -3.582 0.000492 *** poly(juice, 5)5 1057.3 941.3 1.123 0.263509

--Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 941.3 on 122 degrees of freedom Multiple R-squared: 0.7539, Adjusted R-squared: 0.7439 F-statistic: 74.76 on 5 and 122 DF, p-value: 2.2e-16 Poly5 Resampling: Bootstrap (25 reps) RMSE Rsquared RMSE SD Rsquared SD 966 0.734 141 0.0569 Мы получили в точности те же коэффициенты модели, что и при использовании базовой функции lm(), но для критериев RMSE и RSquared были найдены стандартные ошибки, позволяющие рассчитать доверительные интервалы этих статистик. Оценка проводилась по результатам 25 бутстрепитераций, выполняемых функцией train() по умолчанию. Обратите внимание, что несмещенное бутстреп-значение коэффициента детерминации (0.734) несколько меньше, чем рассчитанное функцией summary() для финальной модели (0.754).

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

–  –  –

В последующих разделах мы приведем описание большинства перечисленных методов и продемонстрируем процесс оптимизации их параметров с помощью функции train().

4. ПОСТРОЕНИЕ РЕГРЕССИОННЫХ МОДЕЛЕЙ РАЗЛИЧНОГО ТИПА

4.1. Селекция оптимального набора предикторов линейной модели Модель множественной линейной регрессии является, безусловно, весьма полезной и широко применяемой для прогнозирования количественного отклика. Считается, что наиболее эффективный путь улучшения качества регрессии – исключение незначимых коэффициентов, или, выражаясь точнее, отбор информативного комплекса Sq из q переменных (q m). Причины, по которым стоит проводить селекцию "оптимального подмножества" предикторов, Дж. Фаравэй (Faraway, 2006) видит в следующем:

1. Принцип бритвы Оккама утверждает, что из нескольких вероятных объяснений явления лучшим является самое простое. Компактная модель, из которой удалены избыточные предикторы, лучше объясняет имеющиеся данные.

2. Ненужные предикторы добавляют шум к оценке влияния других интересующих нас факторов. Иначе степени свободы (часто ограниченные) будут тратиться впустую.

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

4. Если модель используется для прогнозирования, то можно сэкономить время и/или деньги, не измеряя избыточные переменные.

Алгоритмы выбора оптимального подмножества Sq обычно основаны на последовательном "переборном" процессе, при котором многократно создаются модели с различными наборами предикторов, лучшая из которых определяется по некоторому критерию эффективности (ошибка или точность модели, R2, АIC и проч.). Рассмотрим технику применения и оценим полученные результаты с использованием трех таких методов: пошаговый (forward/backward) метод селекции, рекурсивное исключение и генетический алгоритм, представленные в пакете caret.

В качестве примера рассмотрим построение регрессионных моделей, прогнозирующих обилие водорослей группы a1 в зависимости от гидрохимических показателей воды и условий отбора проб в различных водотоках (см.

подробное описание таблицы переменных в разделе 3.4):

library(DMwR) data(algae) library(caret) # Заполним пропуски в данных на основе алгоритма бэггинга pPbI - preProcess(algae[, 4:11], method='bagImpute') algae[,4:11] - predict(pPbI, algae[,4:11]) # Сохраним таблицу для использования в дальнейшем save(algae, file="algae.RData") Важные предварительные выводы можно сделать, сформировав корреляционную матрицу предикторов (рис.

4.1):

Mcor -cor(algae[,4:11]) library(corrplot) corrplot(Mcor, method="color", addCoef.col=" darkgreen", addgrid.col = "gray33", tl.col = "black") Рис. 4.1. Корреляционная матрица показателей качества воды в реках Очевидно, что между предикторами существуют корреляционные связи умеренной силы, а коллинеарность, в целом, выражена слабо.

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

lm.a1 - lm(a1 ~., data = algae[, 1:12]) summary(lm.a1)

Coefficients:

Estimate Std. Error t value Pr(|t|) (Intercept) 41.2439725 23.9452133 1.722 0.08667.

seasonspring 3.5927401 4.1371119 0.868 0.38630 seasonsummer 0.3093462 3.9942858 0.077 0.93835 seasonwinter 3.5401487 3.8521670 0.919 0.35930 sizemedium 3.6080130 3.7240594 0.969 0.33390 sizesmall 9.8601114 4.1192992 2.394 0.01769 * speedlow 3.2806448 4.6757407 0.702 0.48380 speedmedium -0.2221952 3.2139617 -0.069 0.94496 mxPH -3.3723566 2.6959345 -1.251 0.21256 mnO2 1.0421490 0.7045909 1.479 0.14083 Cl -0.0377924 0.0337294 -1.120 0.26398 NO3 -1.5215484 0.5501278 -2.766 0.00626 ** NH4 0.0016244 0.0010029 1.620 0.10702 oPO4 -0.0007233 0.0391875 -0.018 0.98529 PO4 -0.0559043 0.0300688 -1.859 0.06459.

Chla -0.0629710 0.0781788 -0.805 0.42159

--Signif. codes: 0 С***Т 0.001 С**Т 0.01 С*Т 0.05 С.Т 0.1 С Т 1 Residual standard error: 17.64 on 184 degrees of freedom Multiple R-squared: 0.3686, Adjusted R-squared: 0.3172 F-statistic: 7.162 on 15 and 184 DF, p-value: 2.94e-12 Отметим, что при всей внушительной адекватности модели в целом (по F-критерию), почти все коэффициенты оцениваются как статистически незначимые. Выполним стандартную пошаговую процедуру включений с исключениями "слабых" предикторов, используя функцию step():

lm_step.a1 - step(lm.a1, trace = 0) summary(lm_step.a1)

Coefficients:

Estimate Std. Error t value Pr(|t|) (Intercept) 55.7204080 20.5365454 2.713 0.00727 ** sizemedium 3.3401764 3.3721698 0.991 0.32316 sizesmall 10.9864235 3.7636536 2.919 0.00393 ** mxPH -3.7829230 2.4260402 -1.559 0.12056 NO3 -1.5247600 0.4931425 -3.092 0.00228 ** NH4 0.0014816 0.0009336 1.587 0.11416 PO4 -0.0691868 0.0102496 -6.750 1.68e-10 ***

--Residual standard error: 17.49 on 193 degrees of freedom Multiple R-squared: 0.3494, Adjusted R-squared: 0.3292 F-statistic: 17.27 on 6 and 193 DF, p-value: 5.882e-16 Доля значимых предикторов стала существенно выше. Проверим также, можно ли считать статистически значимой некоторое увеличение ошибки модели:

anova(lm.a1, lm_step.a1) Analysis of Variance Table Res.Df RSS Df Sum of Sq F Pr(F) 2 193 59006 -9 -1743.4 0.6224 0.777 Выполним тестирование обоих моделей функцией train() из пакета

caret (см. раздел 3.5) с использованием 10-кратной перекрестной проверки:

lm.a1.cv - train(a1 ~., data = algae[, 1:12], method = 'lm', trControl = trainControl(method = "cv")) Resampling: Cross-Validation (10 fold) RMSE Rsquared RMSE SD Rsquared SD 19.0968 0.2955792 4.273678 0.1369611 lm_step.a1.cv - train(a1 ~ size + mxPH + mnO2 + NO3 + NH4 + PO4, data = algae[, 1:12], method = 'lm', trControl = trainControl(method = "cv")) Resampling: Cross-Validation (10 fold) RMSE Rsquared RMSE SD Rsquared SD 17.82252 0.3376223 3.372695 0.1312212 Преимущества модели, полученной с использованием пошаговой процедуры отбора предикторов, вполне очевидны.

Рекурсивное исключение переменных Алгоритм рекурсивного исключения (RFE – Recursive Feature Elimination) выполняется следующим образом. Вначале строится модель по всем предикторам, которые ранжируются по их важности. Далее рассматривается последовательность подмножеств S (S1 S2, и т.д.) из переменных наивысшего ранга. На каждой итерации подмножества Si ранги предикторов пересматриваются, а модели пересчитываются. Итоговая модель основывается на подмножестве Si, обеспечивающем оптимум заданному критерию качества.

Функция rfe() из пакета caret включает алгоритм RFE в процедуру ресэмплинга, и тогда цикл рекурсивного исключения приобретает следующий вид:

1 Цикл для каждой итерации ресэмплинга 2 | Данные разделяются на обучающие и проверочные в зависимости, 2| например, от количества блоков, участвующих в перекрестной проверке 3 | Построение модели на обучающем множестве с использованием всех предикторов 4 | Предсказание на проверочной выборке 5 | Вычисление важности предикторов и их ранжирование 6 | Цикл для каждого подмножества Si, i = 1…S 7 | | Заполнение Si наиболее важными переменными 8 | | Построение модели на обучающем множестве с использованием всех | | предикторов Si 9 | | Предсказание на проверочной выборке 10 | | Повторное вычисление рангов важности для каждого предиктора 11 | Конец цикла 12 Конец цикла 13 Вычисление профилей эффективности для каждого Si с использованием проверочной выборки 14 Определение оптимального числа предикторов 15 Уточнение итогового списка наилучших предикторов 16 Построение итоговой модели на основе оптимального Si с использованием исходной обучающей выборки Синтаксис параметров функций rfe() и rfeControl() в целом похож на таковой у функций train() и trainControl() (см. раздел 3.5). Отметим, что функция rfe() не осуществляет оптимизацию модели, включающей номинальные предикторы, поэтому их необходимо преобразовать в индикаторные переменные с использованием функции model.matrix().

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

x - model.matrix(a1 ~., data = algae[, 1:12])[, -1] set.seed(10) ctrl - rfeControl(functions = lmFuncs, method = "cv", verbose = FALSE, returnResamp = "final") lmProfileF - rfe(as.data.frame(x), algae$a1, sizes = 1:10, rfeControl = ctrl)

–  –  –

Рис. 4.2. Изменение корня из среднеквадратичной ошибки и коэффициента детерминации в зависимости от числа предикторов (по результатам перекрестной проверки) [1] "NO3" "sizesmall" "PO4" "NH4" "mnO2" "mxPH" [7] "Cl" "seasonwinter" "sizemedium"

–  –  –

sizemedium 3.2637287 3.3692541 0.969 0.33394 Cl -0.0354895 0.0324232 -1.095 0.27509

--Signif. codes: 0 С***Т 0.001 С**Т 0.01 С*Т 0.05 С.Т 0.1 С Т 1 Residual standard error: 17.44 on 190 degrees of freedom Multiple R-squared: 0.3627, Adjusted R-squared: 0.3325 F-statistic: 12.02 on 9 and 190 DF, p-value: 6.177e-15 Нельзя утверждать, что нам удалось достигнуть серьезных успехов, применив метод RFE: пошаговая модель lm_step.a1 была компактнее и чуть лучше (по результатам перекрестной проверки), хотя и уступала модели RFE по величине ошибки на обучающей выборке.

Генетический алгоритм Генетический алгоритм, позаимствованный у природных аналогов и разработанный Дж. Холландом (Holland, 1975), отличается от большинства иных процедур селекции тем, что поиск оптимального решения развивается не сам по себе, а с учетом предыдущего опыта.

Смысл его заключается в следующем:

формируемое решение кодируется как вектор х0, который называется хромосомой и соответствует битовой маске, т.е. двоичному представлению набора исходных переменных;

инициализируется исходная "популяция" 0 ( x10... x ) потенциальных решений, состоящая из некоторого количества хромосом ;

каждой хромосоме в популяции присваиваются две оценки: значение эффективности µ( xi0 ) в соответствии с заданной функцией оптимальности и вероятность воспроизведения P( xi0 ), которая зависит от "перспективности" этой хромосомы;

в соответствии с вероятностями воспроизведения хромосомы производят потомков, используя операции кроссинговера (обмена фрагментами между хромосомами) и мутации (замена значения бита 0/1 на противоположный) – см. рис. 4.3;

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

формирование новой популяции многократно повторяется и осуществляется поиск субоптимальных моделей;

процесс останавливается, если получено удовлетворительное решение, либо исчерпано все отведенное на эволюцию время.

Синтаксис пары функций gafs() и gafsControl(), реализующих генетический алгоритм в пакете caret, в целом аналогичен таковому у функций rfe() и rfeControl(), но требует задания дополнительных аргументов, регулирующих скорость эволюции.

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

Рис. 4.3. Схема кроссинговера и мутации во время применения генетического алгоритма set.seed(10) ctrl - gafsControl(functions = rfGA, method = "cv", verbose = FALSE, returnResamp = "final") lmProfGafs - gafs(algae[, 1:11], algae$a1, iters = 10, # 10 генераций gafsControl = ctrl) lmProfGafs Genetic Algorithm Feature Selection Maximum generations: 10 Population per generation: 50 Crossover probability: 0.8 Mutation probability: 0.1 Elitism: 0 Internal performance values: RMSE, Rsquared Subset selection driven to minimize internal RMSE External performance values: RMSE, Rsquared Best iteration chose by minimizing external RMSE External resampling method: Cross-Validated (10 fold)

During resampling:

* the top 5 selected variables (out of a possible 11):

Cl (100%), PO4 (100%), size (100%), mxPH (50%), oPO4 (40%) * on average, 5.2 variables were selected (min = 3, max = 7)

In the final search using the entire training set:

* 4 features selected at iteration 1 including:

mxPH, Cl, NO3, PO4 * external performance at this iteration is RMSE Rsquared 15.381 0.519 lm_gafs.a1 - lm(a1 ~ size + mxPH + Cl + NO3 + PO4, data = algae[, 1:12]) train(a1 ~ size + mxPH + Cl + NO3 + PO4, data = algae[, 1:12], method = 'lm', trControl = trainControl(method = "cv")) RMSE Rsquared RMSE SD Rsquared SD 0.3430528 2.838075 0.09572028 17.72904 Несмотря на большие затраты вычислительных ресурсов (поиск решения продолжался более 30 мин.), была получена комбинация предикторов, превосходящая по критериям RMSE и Rsquared все три предыдущие модели.

Тестирование моделей с использованием дополнительного набора данных Возникает естественный вопрос: а какой из четырех полученных моделей следует все же отдать предпочтение при прогнозировании?

Хорошую возможность ответить на него предоставляет нам Л. Торго (Torgo, 2011), подготовивший на сайте своей книги (http://www.dcc.fc.up.pt) специально предназначенный для этого набор данных из 140 наблюдений (см.

файл Eval.txt с предикторами и Sols.txt со значениями отклика).

Пропущенные значения заполним с использованием алгоритма бэггинга:

Eval - read.table('Eval.txt', header=F, dec='.', col.names=c('season','size','speed','mxPH','mnO2','Cl', 'NO3','NH4','oPO4','PO4','Chla'), na.strings=c('XXXXXXX')) Sols - read.table('Sols.txt', header=F, dec='.', col.names=c('a1','a2','a3','a4','a5','a6','a7'), na.strings=c('XXXXXXX')) ImpEval - preProcess(Eval[, 4:11], method='bagImpute') Eval[,4:11] - predict(ImpEval, Eval[,4:11]) # Cохраним данные для дальнейшего использования save(Eval,Sols, file="algae_test.RData") y - Sols$a1 EvalF - as.data.frame(model.matrix(y ~.,Eval)[,-1]) Выполним прогноз для набора предикторов тестовой выборки и оценим точность каждой модели по трем показателям: среднему абсолютному отклонению (MAE), корню из среднеквадратичного отклонения (RSME) и квадрату коэффициента детерминации Rsq = 1 - NSME, где NSME – относительная ошибка, равная отношению средних квадратов отклонений от регрессии и от общего среднего:

# Функция, выводящая вектор критериев ModCrit - function (pred, fact) { mae - mean(abs(pred-fact)) rmse - sqrt(mean((pred-fact)^2)) Rsq - 1-sum((fact-pred)^2)/sum((mean(fact)-fact)^2) c(MAE=mae, RSME=rmse, Rsq = Rsq ) } y - Sols$a1 EvalF - as.data.frame(model.matrix(y ~.,Eval)[,-1]) Result - rbind( lm_full = ModCrit(predict(lm.a1,Eval),Sols[,1]), lm_step = ModCrit(predict(lm_step.a1,Eval),Sols[,1]), lm_rfe = ModCrit(predict(lmProfileF$fit,EvalF),Sols[,1]), lm_gafs = ModCrit(predict(lm_gafs.a1,Eval),Sols[,1])) Result MAE RSME Rsq lm_full 12.64484 16.80890 0.3268512 lm_step 12.47868 16.62016 0.3418829 lm_rfe 12.38730 16.52863 0.3491115 lm_gafs 12.74591 17.19192 0.2958234 Вероятно, нет смысла делать серьезные выводы из того, что рейтинг лучших моделей при перекрестной проверке и при независимом тестировании на объектах, не участвовавших в построении моделей, оказался столь несовпадающим. Во-первых, разброс критериев точности моделей находится в пределах доверительных интервалов, поэтому их ранжирование, строго говоря, можно трактовать как обусловленное случайными причинами.

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

4.2. Регуляризация, частные наименьшие квадраты и kNN-регрессия Регрессия по методу "лассо" Регрессия по методу наименьших квадратов (МНК) часто может стать неустойчивой, то есть сильно зависящей от обучающих данных, что обычно является проявлением тенденции к переобучению. Избежать такого переобучения помогает регуляризация – общий метод, заключающийся в наложении дополнительных ограничений на искомые параметры, которые могут предотвратить излишнюю сложность модели. Смысл процедуры заключается в "стягивании" вектора коэффициентов в ходе их настройки таким образом, чтобы они в среднем оказались несколько меньше по абсолютной величине, чем это было бы при оптимизации по МНК.

Метод регрессии "лассо" (LASSO, Least Absolute Shrinkage and Selection Operator) заключается во введении дополнительного слагаемого регуляризации в функционал оптимизации модели, что часто позволяет получать более устойчивое решение.

Условие минимизации квадратов ошибки при оценке параметров выражается следующей формулой:

n m = arg min[ ( yi j xij ) 2 + | |], i =1 j =1 где – параметр регуляризации, имеющий смысл штрафа за сложность.

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

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

Обобщением регрессии с регуляризацией можно считать модель "эластичных сетей" (elastic net – Zou, Hastie, 2005).

Эта модель устанавливает сразу два типа штрафных параметров 1 и 2, объединяет гребневую регрессию (при 2 =0) и регрессию "лассо" (при 1 =0):

n m = arg min[ i =1 ( yi j =1 j xij ) 2 + 1 () 2 + 2 | |].

Продолжим рассмотрение примеров построения регрессионных моделей, прогнозирующих обилие водорослей группы a1 в зависимости от гидрохимических показателей воды и условий отбора проб в различных водотоках (см. подробное описание таблицы переменных в разделе 3.4).

Поскольку даже ориентировочная величина параметра регуляризации нам неизвестна, то на первом этапе с использованием функции glmnet() из одноименного пакета проведем анализ изменения значений коэффициентов в зависимости от в ее широком диапазоне от 0.1 до 1000:

load(file="algae.RData") # Загрузка таблицы algae - раздел 4.1 grid=10^seq(10,-2,length=100) library( glmnet) lasso.a1 - glmnet(x,algae$a1, alpha=1, lambda=grid) plot(lasso.a1, xvar = "lambda", label = TRUE, lwd=2)

–  –  –

Функция train() для метода glmnet выполняет оптимизацию сразу двух моделей регуляризации: при alpha = 1 подгоняется модель по методу лассо, а при alpha = 0 – гребневая регрессия.

Примем лассо-модель, и выполним тонкую настройку значения lambda в интервале от 0.5 до 4.5:

x - model.matrix(a1 ~.,data=algae[,1:12])[,-1] grid.train = seq(0.5,4.5,length=15) lasso.a1.train - train(as.data.frame(x), algae$a1, method='glmnet', tuneGrid = expand.grid(.lambda = grid.train,.alpha = 1), trControl = trainControl(method = "cv"))

Resampling results across tuning parameters:

lambda RMSE Rsquared RMSE SD Rsquared SD 0.5000000 17.86789 0.3122892 3.148867 0.11004143 0.7857143 17.70163 0.3226125 3.119873 0.10317362 1.0714286 17.60777 0.3300379 3.064511 0.09865072 1.3571429 17.58081 0.3339868 3.010830 0.09807188 1.6428571 17.58690 0.3368441 2.974649 0.09840106 1.9285714 17.61382 0.3394520 2.945534 0.09889217 2.2142857 17.66452 0.3408752 2.914110 0.09916679 2.5000000 17.73367 0.3414334 2.892287 0.09909826 2.7857143 17.81998 0.3409785 2.889505 0.09957978 … 4.5000000 18.49634 0.3249896 2.953191 0.10310683 Tuning parameter 'alpha' was held constant at a value of 1 RMSE was used to select optimal model using smallest value.

The final values used for model were alpha = 1 and lambda = 1.357143.

Выведем протокол со значениями коэффициентов модели:

coef(lasso.a1.train$finalModel, lasso.a1.train$bestTune$lambda) 16 x 1 sparse Matrix of class "dgCMatrix" (Intercept) 38.37192105 seasonspring.

seasonsummer.

seasonwinter.

sizemedium.

sizesmall 7.42569282 speedlow.

speedmedium.

mxPH -1.78557347 mnO2 0.11517464 Cl -0.03960609 NO3 -0.51584616 NH4.

oPO4.

PO4 -0.05193156 Chla -0.02286621 Коэффициенты регрессии, отличные от 0, составляют информативный набор предикторов.

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

Напомним основные шаги построения модели регрессии на главные компоненты (PCR):

стандартизация матрицы исходных предикторов X и отклика Y;

выбор числа собственных значений p и формирование матрицы главных компонент (часто называемой таблицей счетов – scores) Tnr = Xnm Pтmp, где P – матрица нагрузок (loadings);

оценка вектора коэффициентов А множественной регрессии на главные Ar = (TтT)-1ТтY;

компоненты: Y = TA, пересчет коэффициентов регрессии на главные компоненты обратно в коэффициенты регрессии для исходных предикторов: Bm = PA.

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

Метод частных наименьших квадратов PLS (Partial Least Squares, или Projection into Latent Structure) также использует разложение исходных предикторов по осям главных компонент, но дополнительно выделяет подмножество латентных переменных, в пространстве которых связь между зависимой переменной и предикторами достигает максимального значения.

В случае одномерного отклика сначала оценивается корреляционная связь между предикторами X и Y, осуществляется сингулярное разложение матрицы XтY, и формируется вектор W первого направления PLS (direction), при вычислении которого особое внимание уделяется переменным, которые наиболее тесно связаны с откликом. Исходные переменные X ортогонально проецируются на ось, задаваемую вектором W, а затем последовательно рассчитываются значения T счетов и нагрузок P.

Для нахождения второго направления PLS вычисляются остатки, которые остались необъясненными первым PLS-направлением. Оценивается направление новой оси наибольшей корреляции и оцениваются значения счетов с использованием ортогонализированных остатков. Такой итеративный подход можно применить p раз пока модель не достигнет оптимальной сложности.

Реализация PLS-регрессии в среде R представлена в пакете pls, который включает в себя большое количество функций для построения регрессионных моделей, создания диагностических графиков и извлечения информации из результатов вычислений:

library(pls) M.pls - plsr(algae$a1 ~ x, scale = TRUE, Validation = "CV", method = "oscorespls") summary(M.pls) Data: X dimension: 200 15 Y dimension: 200 1 Fit method: oscorespls Number of components considered: 15

–  –  –

Оптимальная модель в нашем случае основывается только на 1-м направлении PLS, поскольку при включении новых осей проецирования ошибка перекрестной проверки монотонно возрастает – см. рис. 4.5. Заметим, что накопленная доля объясненной дисперсии для матрицы предикторов X равномерно возрастает при увеличении числа компонент, тогда как почти вся вариация отклика algae$a1 сконцентрирована вдоль главного направления.

Выполним подбор размерности пространства латентных переменных с использованием функции train(). По умолчанию диапазон значений оптимизируемого параметра ncomp принимает значения от 1 до

tuneLength:

set.seed(100) ctrl - trainControl(method = "cv", number = 10) (plsTune.a1 - train(x, algae$a1, method = "pls", tuneLength = 14, trControl = ctrl, preProc = c("center", "scale"))) Resampling: Cross-Validation (10 fold) ncomp RMSE Rsquared RMSE SD Rsquared SD 1 17.8 0.349 3.97 0.167 2 18.6 0.299 4.01 0.151 3 18.9 0.287 4.27 0.158...

13 18.8 0.3 4.29 0.163 14 18.8 0.3 4.29 0.163 RMSE was used to select optimal model using smallest value.

The final value used for the model was ncomp = 1.

Выполним для сравнения подбор числа компонент для модели регрессии на главные компоненты PCR:

set.seed(100) ctrl - trainControl(method = "cv", number = 10) (pcrTune.a1 - train(x, algae$a1, method = "pcr", tuneLength = 14, trControl = ctrl, preProc = c("center", "scale"))) ncomp RMSE Rsquared RMSE SD Rsquared SD 1 17.5 0.358 3.7 0.173 2 17.6 0.352 3.59 0.169 3 17.6 0.355 3.68 0.165...

13 18.4 0.321 4.34 0.162 14 18.5 0.324 4.63 0.157 The final value used for the model was ncomp = 1.

Никаких преимуществ метода PLS по сравнению с PCR на нашем примере обнаружено не было; более того, ошибка в целом модели на главные компоненты оказалась несколько ниже.

Регрессия по методу k ближайших соседей Ключевая идея, лежащая в основе метода k ближайших соседей, как уже обсуждалось в разделе 3.4, состоит в формулировке модели в терминах евклидовых расстояний в исходном многомерном пространстве признаков.

Задача заключается в том, чтобы для каждой тестируемой точки x0 найти такую -окрестность (многомерный эллипсоид), чтобы в ней поместилось ровно k точек с известными значениями y. Тогда прогноз f(x0) можно получить, усредняя значения отклика всех обучающих наблюдений из.

Функция train() оптимизирует число соседей k, используя другую функцию – knnreg() из пакета caret (рис.

4.6):

knnTune.a1 - train(x, algae$a1, method = "knn", preProc = c("center", "scale"), trControl = ctrl, tuneGrid = data.frame(.k = 4:25)) plot( knnTune.a1)

Resampling results across tuning parameters:

k RMSE Rsquared RMSE SD Rsquared SD 4 18.5 0.307 3.25 0.157 5 18.2 0.322 3.64 0.179 6 18.2 0.313 3.57 0.176 7 17.6 0.343 3.74 0.166...

18 17.1 0.395 3.59 0.115 19 16.9 0.41 3.69 0.116 20 16.9 0.413 3.74 0.121 21 16.8 0.42 3.52 0.119 22 16.9 0.418 3.48 0.113 23 17 0.409 3.62 0.117 24 17.1 0.407 3.54 0.118 25 17.2 0.397 3.57 0.104 RMSE was used to select optimal model using smallest value.

The final value used for the model was k = 21.

Рис. 4.6. Зависимость ошибки перекрестной проверки от числа ближайших соседей Тестирование моделей с использованием дополнительного набора данных Используем для прогноза набор данных (Torgo, 2011) из 140 наблюдений, который мы уже применяли в предыдущем разделе. Данные, подготовленные для тестирования с восстановленными пропущенными значениями – таблицы Eval с предикторами и Sols со значениями отклика – мы сохранили в файле algae_test.RData (см раздел 4.1).

Выполним прогноз для набора предикторов проверочной выборки и оценим точность каждой модели по трем показателям: среднему абсолютному отклонению (MAE), корню из среднеквадратичного отклонения (RSME) и квадрату коэффициента детерминации Rsq = 1 - NSME, где NSME – относительная ошибка, равная отношению среднего квадрата отклонений от модельных значений и от общего среднего:

load(file="algae_test.RData") # Загрузка таблиц Eval, Sols y - Sols$a1 EvalF - as.data.frame(model.matrix(y ~.,Eval)[,-1]) # Функция, выводящая вектор критериев ModCrit - function (pred, fact) { mae - mean(abs(pred-fact)) rmse - sqrt(mean((pred-fact)^2)) Rsq - 1-sum((fact-pred)^2)/sum((mean(fact)-fact)^2) c(MAE=mae, RSME=rmse, Rsq = Rsq ) } Result - rbind( lasso = ModCrit(predict(lasso.a1.train,EvalF),Sols[,1]), pls_1c = ModCrit(predict(plsTune.a1, EvalF),Sols[,1]), pcr_1c = ModCrit(predict(pcrTune.a1, EvalF),Sols[,1]), kNN_21 = ModCrit(predict(knnTune.a1, EvalF),Sols[,1])) Result MAE RSME Rsq lasso 12.74513 17.13765 0.3002625 pls_1c 12.76695 17.27928 0.2886493 pcr_1c 12.85117 17.29990 0.2869504 kNN_21 12.21037 16.30965 0.3662442 Отметим, что регрессия на главные компоненты и PLS в условиях слабой мультиколлинеарности данных не приносит никаких ощутимых преимуществ по сравнению с классическими линейными моделями раздела

4.1. Модель по методу лассо, отлично зарекомендовавшая себя при перекрестной проверке, сработала хуже на свежих данных, что связано, видимо, с дрейфом параметра регуляризации от своего оптимального значения. И, наконец, пока лучшей из всех исследованных оказалась методически простейшая модель регрессии по 21 ближайшему соседу. В то же время эта непараметрическая модель имеет важнейший недостаток – она не предоставляет никакой информации для графической или содержательной интерпретации.

4.3. Построение деревьев регрессии Деревья решений (Breiman at al., 1984; Quinlan, 1986) осуществляют разбиение пространства объектов в соответствии с некоторым набором правил разбиения (splitting rule). Эти правила являются логическими утверждениями в отношении той или иной переменной и могут быть истинными или ложными. Ключевыми здесь являются три обстоятельства:

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

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

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

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

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

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

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

5. Деревья решений можно эффективно применять к данным с пропущенными значениями, что очень полезно при решении практических задач, где наличие пропущенных значений – это, скорее, правило, чем исключение.

6. Деревья решений одинаково хорошо применимы как к количественным, так и к качественным зависимым переменным.

К недостаткам этого класса моделей иногда относят нестабильность и невысокую точность предсказаний, что, как будет показано ниже, не всегда подтверждается. По своей сути, деревья используют "наивный подход" (naive approach) в том смысле, что они исходят из предположения о взаимной независимости признаков. Поэтому модели регрессионных деревьев статистически наиболее работоспособны, когда комплекс анализируемых переменных является не слишком мультиколлинеарным или имеется регулярная внутренняя множественная альтернатива в исходной комбинации признаков.

Алгоритм CART (Classification and Regression Tree) рекурсивно делит исходный набор данных на подмножества, которые становятся все более и более гомогенными относительно определенных признаков, в результате чего формируется древовидная иерархическая структура. Деление осуществляется на основе традиционных логических правил в виде ЕСЛИ (А) ТО (В), где А – некоторое логическое условие, а В – процедура деления подмножества на две части, для одной из которых условие А истинно, а для другой – ложно.

Примеры условий: Xi==F, Хi = V; Хi = V и др., где Хi – один из предикторов исходной таблицы, F – выбранное значение категориальной переменной, V – специально подобранное опорное значение (порог).

На первой итерации корневой узел дерева связывается с наиболее оптимальным условным суждением, и все множество объектов делится на две группы. От каждого последующего узла-родителя к узлам-потомкам также может отходить по две ветви, в свою очередь связанные c граничными значениями других наиболее подходящих переменных и определяющие правила дальнейшего разделения (splitting criterion). Конечными узлами дерева являются "листья", соответствующие найденным решениям и объединяющие все разделенные на группы объекты обучающей выборки.

Общее правило выбора опорного значения для каждого узла построенного дерева можно сформулировать следующим образом: «выбранный признак должен разбить множество Х* так, чтобы получаемые в итоге подмножества Х*k, k = 1, 2, …, p, состояли из объектов, принадлежащих к одному классу, или были максимально приближены к этому».

Описанный процесс относится к так называемым "жадным" алгоритмам, стремящимся, не считаясь ни с чем, построить максимально "кустистое" дерево (также "глубокое дерево", deep tree). Естественно, чем обширнее и кустистее дерево, тем лучше будут результаты его тестирования на обучающей выборке, но не столь успешными – на проверочной выборке.

Поэтому построенная модель должна быть еще и оптимальной по размерам, т.е. содержать информацию, улучшающую качество распознавания, и игнорировать ту информацию, которая его не улучшает. Для этого обычно проводят "обрезание" дерева (tree pruning) – отсечение ветвей там, где эта процедура не приводит к серьезному возрастанию ошибки.

Невозможно подобрать объективный внутренний критерий, приводящий к хорошему компромиссу между безошибочностью и компактностью, поэтому стандартный механизм оптимизации деревьев основан на перекрестной проверке (Loh, Shih, 1997). Для этого обучающая выборка разделяется, например, на 10 равных частей: 9 частей используется для построения дерева, а оставшаяся часть играет роль проверочной совокупности. После многократного повторения этой процедуры из некоторого набора деревьев-претендентов, у которых имеется практически допустимый разброс критериев качества модели, выбирается дерево, показавшее наилучший результат при перекрестной проверке.

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

Функция rpart() из одноименного пакета выполняет рекурсивный выбор для каждого следующего узла таких разделяющих значений, которые приводят к минимальной сумме квадратов внутригрупповых отклонений Dt для всех t узлов дерева.

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

штраф за сложность модели (cost complexity), включающий штрафной множитель за каждую неотсечённую ветвь СС(T) = Dt + t ;

t

–  –  –

ошибка перекрестной проверки (CVer) с разбиением на k блоков (обычно k = 10), также отнесенная к девиансу нуль-дерева D0; CVer, как правило, больше, чем RELer;

стандартное отклонение (SE) ошибки перекрестной проверки.

Лучшим считается дерево, состоящее из такого количества ветвей t, для которого сумма (CVer + SE) является минимальной.

В качестве примера рассмотрим построение дерева CART, прогнозирующего обилие водорослей группы a1 в зависимости от гидрохимических показателей воды и условий отбора проб в различных водотоках (см. разделы 3.4 и 4.1-4.2). Используем сначала пакет rpart, для работы с которым обычно применяется двухшаговая процедура: функция устанавливает связи между зависимой и независимыми rpart() переменными и формирует бинарное дерево, а функция prun() выполняет обрезание лишних ветвей.

load(file="algae.RData") # Загрузка таблицы algae - раздел 4.1 (rt.a1 - rpart::rpart(a1 ~., data = algae[, 1:12])) n= 200 node), split, n, deviance, yval * denotes terminal node

1) root 200 90694.880 16.923500

2) PO4=43.818 148 31359.210 8.918919

4) Cl=7.8065 141 21678.580 7.439716

8) oPO4=51.118 85 3455.770 3.801176 *

9) oPO4 51.118 56 15389.430 12.962500

18) mnO2=10.05 24 1248.673 6.716667 *

19) mnO2 10.05 32 12502.320 17.646880

38) NO3=3.1875 9 257.080 7.866667 *

39) NO3 3.1875 23 11047.500 21.473910

78) mnO2 8 13 2919.549 13.807690 *

79) mnO2=8 10 6370.704 31.440000 *

5) Cl 7.8065 7 3157.769 38.714290 *

3) PO4 43.818 52 22863.170 39.705770

6) mxPH 7.87 28 11636.710 32.875000

12) oPO4=3.1665 14 1408.304 23.978570 *

13) oPO4 3.1665 14 8012.309 41.771430 *

7) mxPH=7.87 24 8395.785 47.675000

14) PO4=15.177 12 3047.517 38.183330 *

15) PO4 15.177 12 3186.067 57.166670 * Приведенной командой мы построили полное дерево без обрезания ветвей, состоящее из 9 узлов и 10 листьев, обозначенных в приведенном протоколе разбиения символом *. В каждой строке представлены по порядку:

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

Например, перед первой итерацией общее множество из 200 наблюдений имеет среднее значение m =16.92 при девиансе D = 90694. При PO4=43.8 это множество делится на две части: 2) 148 наблюдений (m =8.92, D = 31359) и 3) 52 наблюдения с высоким уровнем обилия водорослей (m =39.7, D = 22863). Дальнейшие разбиения каждой из этих двух частей аналогичны.

Разумеется, лучший вариант – представить дерево графически.

Популярны три варианта визуализации с использованием различных функций: plot(), prettyTree() из пакета DMwR и prp() из чрезвычайно продвинутого пакета rpart.plot:

DMwR::prettyTree(rt.a1)

Рис. 4.7. Дерево rpart без обрезания ветвей

Полезно также проследить изменение перечисленных выше статистических критериев по мере выращивания дерева:

printcp(rt.a1)

Variables actually used in tree construction:

[1] Cl mnO2 mxPH NO3 oPO4 PO4 Root node error: 90695/200 = 453.47 n= 200 CP nsplit rel error xerror xstd 1 0.402145 0 1.00000 1.01569 0.13062 2 0.071921 1 0.59785 0.71689 0.11912 3 0.031241 2 0.52593 0.71641 0.11979 4 0.031211 3 0.49469 0.74315 0.12153 5 0.024435 4 0.46348 0.73140 0.12122 6 0.023840 5 0.43905 0.74064 0.11900 7 0.018065 6 0.41521 0.73495 0.11532 8 0.016291 7 0.39714 0.74622 0.11461 9 0.010000 9 0.36456 0.75839 0.11467 Функция rpart() и другие функции из пакета rpart имеют собственные возможности выполнить перекрестную проверку и оценить ее ошибку при различных значениях штрафа за сложность модели cp:

# Снижаем порог штрафа за сложность с шагом.005 rtp.a1 - rpart(a1 ~., data=algae[, 1:12], control=rpart.control(cp=.005)) # График зависимости относительных ошибок от числа узлов plotcp(rtp.a1) with(rtp.a1, {lines(cptable[,2]+1, cptable[,3], type="b", col="red") legend("topright", c("Ошибка обучения", "Ошибка крос-проверки (CV)", "min(CV ошибка)+SE"), lty=c(1,1,2), col=c("red","black","black"), bty="n") })

–  –  –

На рис. 4.8 видно, что минимум относительной ошибки при перекрестной проверке приходится на значение cp=0.029.

Выполним обрезку дерева при этом значении:

rtp.a1 - prune(rtp.a1, cp=0.029) prettyTree(rtp.a1) Рис. 4.9. Дерево rpart c обрезанием ветвей при cp=0.029

–  –  –

RMSE was used to select optimal model using smallest value.

The final value used for the model was cp = 0.0418.

Рис. 4.10. Оценка параметра cp с использованием функции train() rtt.a1 - rt.a1.train$finalModel prettyTree(rtt.a1) Рис. 4.11. Дерево rpart c обрезанием ветвей при cp=0.0418 При cp = 0.0418 было получено существенно урезанное дерево, которое, правда, значительно потеряло в своей объясняющей ценности.

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

Стандартный механизм проверки статистического гипотез, который предотвращает переусложнение модели, реализован в функции ctree(), использующей метод построения деревьев на основе "условного вывода" (conditional inference). Алгоритм принимает во внимание характер распределения отдельных переменных и осуществляет на каждом шаге рекурсивного разделения данных несмещенный выбор влияющих ковариат, используя формальный тест на основе статистического критерия c(tj, µ j, j), j = 1,..., m, где µ, – соответственно среднее и ковариация (Hothorn et al., 2006). Оценка статистической значимости с-критерия выполняется на основе перестановочного теста, в результате чего формируются компактные деревья, не требующие процедуры обрезания.

library(party) # Построение дерева методом "условного вывода" (ctree.a1 - ctree(a1 ~., data = algae[, 1:12])) plot(ctree.a1) Conditional inference tree with 4 terminal nodes Number of observations: 200

1) PO4 = 43.5; criterion = 1, statistic = 47.106 2)* weights = 52

1) PO4 43.5

3) oPO4 = 51.111; criterion = 0.985, statistic = 10.234

4) size == {small}; criterion = 0.993, statistic = 14.65 5)* weights = 14

4) size == {large, medium} 6)* weights = 49

3) oPO4 51.111 7)* weights = 85 Оптимизацию параметра mincriterion выполним с использованием функции train() при тех же условиях перекрестной проверки:

ctree.a1.train - train(a1 ~., data = algae[, 1:12], method = "ctree", tuneLength = 10, trControl = cvCtrl) ctreet.a1 - ctree.a1.train$finalModel plot(ctreet.a1) Resampling: Cross-Validation (10 fold, repeated 3 times) mincriterion RMSE Rsquared RMSE SD Rsquared SD 0.01 16.9 0.423 4.58 0.193 0.119 16.8 0.432 4.42 0.19 0.228 16.5 0.436 4.58 0.194 0.337 16.4 0.45 4.63 0.197 0.446 16.2 0.462 4.21 0.18 0.554 16.1 0.464 4.25 0.18 0.663 15.6 0.487 4.16 0.186 0.772 15.6 0.485 4.17 0.174 0.881 15.6 0.487 4.36 0.182 0.99 16 0.443 4.65 0.157 RMSE was used to select optimal model using smallest value.

The final value used for the model was mincriterion = 0.881.

Рис. 4.12. Дерево cpart без оптимизации параметра mincriterion Рис. 4.13. Дерево cpart после оптимизации параметра mincriterion Здесь имел место обратный процесс: число узлов дерева было предложено увеличить с 7 до 11. Обратим также внимание на то, что в дереве появились категориальные переменные (размер и скорость течения реки), которые были проигнорированы rpart-деревьями.

Тестирование моделей с использованием дополнительного набора данных Используем для прогноза набор данных (Torgo, 2011) из 140 наблюдений, который мы уже применяли в предыдущем разделе. Данные с восстановленными пропущенными значениями мы сохранили ранее в файле algae_test.RData (см раздел 4.1).

Оценим точность каждой модели на этом наборе данных по трем показателям: среднему абсолютному отклонению (MAE), корню из среднеквадратичного отклонения (RSME) и квадрату коэффициента детерминации Rsq = 1 - NSME, где NSME – относительная ошибка, равная отношению средних квадратов отклонений от модельных значений и от общего среднего:

load(file="algae_test.RData") # Загрузка таблиц Eval,Sols # Функция, выводящая вектор критериев ModCrit - function (pred, fact) { mae - mean(abs(pred-fact)) rmse - sqrt(mean((pred-fact)^2)) Rsq - 1-sum((fact-pred)^2)/sum((mean(fact)-fact)^2) c(MAE=mae, RSME=rmse, Rsq=Rsq) } Result - rbind( rpart_prune = ModCrit(predict(rtp.a1,Eval),Sols[,1]), rpart_train = ModCrit(predict(rt.a1.train,Eval),Sols[,1]), ctree_party = ModCrit(predict(ctree.a1,Eval),Sols[,1]), ctree_train = ModCrit(predict(ctree.a1.train,Eval),Sols[,1]) ) Result MAE RSME Rsq rpart_prune 11.16546 16.09485 0.3828278 rpart_train 10.72834 15.36578 0.4374751 ctree_party 11.32286 16.53470 0.3486336 ctree_train 11.25551 16.40532 0.3587876 Можно с разумной осторожностью сделать вывод о том, что деревья, построенные rpart(), немного точнее, чем деревья условного вывода ctree(). При этом все деревья решений оказались существенно эффективней для прогнозирования свежих данных, чем все ранее построенные модели.

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

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

Рис. 4.14. Ансамбль из пяти линейных классификаторов: каждый сегмент пространства объектов отличается средними вероятностями предсказания классов (подробности см. Флах, 2015, с. 344) В разделе 2.2 был описан бутстреп – процедура генерации повторных случайных выборок из исходного набора данных. Бутстреп-выборки производятся равномерно и с возвращением, поэтому некоторые исходные примеры будут отсутствовать, а другие – дублироваться: в среднем одна такая выборка содержит около 2/3 уникальных исходных наблюдений.

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

Подобно тому, как усреднение нескольких наблюдений снижает оценку дисперсии данных, так и разумным способом снижения дисперсии прогноза является извлечение большого количества порций данных из генеральной совокупности, построение предсказательной модели по каждой обучающей выборке и усреднение полученных предсказаний. Если вместо отдельных обучающих выборок (которых нам, как правило, всегда не хватает) выполнить бутстреп и на основе сгенерированных псевдо-выборок построить В деревьев регрессии, то средний коллективный прогноз B f bag ( x) ={f (x) + f (x) +...+ f (x)}/B будет обладать более низкой дисперсией. Эта процедура и называется бэггингом (bagging, сокр. от bootstrap aggregating). Как нами будет показано в последующих главах, бэггинг можно проводить не только в отношении деревьев регрессии, но и иных моделей: опорных векторов, нейронных сетей, линейных дискриминантов, байесовских вероятностей и др.

Метод случайного леса (Random Forest) представляет собой дальнейшее улучшение бэггинга деревьев решений, которое заключается в устранении корреляции между деревьями. Как и в случае с бэггингом, мы строим несколько сотен деревьев решений по обучающим бутстреп-выборкам.

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

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

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

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

Естественно, если случайный лес строится с использованием m = р, то вся процедура сводится к простому бэггингу.

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

mtry = ncol(x):

load(file="algae.RData") # Загрузка таблицы algae - раздел 4.1 x - as.data.frame(model.matrix(a1~.,data=algae[,1:12])[,-1]) library(randomForest) randomForest(x, algae$a1, mtry = ncol(x)) Type of random forest: regression Number of trees: 500 No. of variables tried at each split: 15 Mean of squared residuals: 271.6539 % Var explained: 40.09 Как видно из полученных результатов, прогнозирование выполнялось по 500 деревьям, в которых было использовано только 40% исходных переменных.

Оценить эффективность этой модели при перекрестной проверке можно с использованием функции train() из пакета caret:

bag.a1 - train(x, algae$a1, preProc = c('center', 'scale'), method = 'rf', trControl = trainControl(method = "cv"), tuneGrid = expand.grid(.mtry = ncol(x))) Resampling: Cross-Validated (10 fold) RMSE Rsquared RMSE SD Rsquared SD 16.27656 0.4503684 4.495015 0.1934045 Модель случайного леса можно построить этой же процедурой, задав последовательность значений mtry для оптимизации:

ranfor.a1 - train(x, algae$a1, preProc = c('center', 'scale'), method = 'rf', trControl = trainControl(method = "cv"), tuneGrid = expand.grid(.mtry = 2:10), importance = TRUE)

Resampling results across tuning parameters:

mtry RMSE Rsquared RMSE SD Rsquared SD 2 15.58681 0.4973347 3.339887 0.2115185 3 15.64973 0.4893305 3.613366 0.2112102 4 15.87789 0.4819333 3.467073 0.2047955 5 15.92607 0.4760140 3.492468 0.2048518 6 15.94487 0.4722841 3.618424 0.2049526 7 16.09182 0.4666949 3.623971 0.2080536 8 16.12491 0.4651339 3.684320 0.2122767 9 16.41238 0.4516508 3.713966 0.2038540 10 16.30965 0.4522984 3.696520 0.2039949 RMSE was used to select optimal model using smallest value.

The final value used for the model was mtry = 2.

Заметим, что бутстреп дает хорошую возможность провести специальную процедуру перекрестной проверки, называемую "тестом по наблюдениям, не попавшим в сумку" (out–of–bag observations). Поскольку ключевая идея бэггинга состоит в многократном построении моделей по наблюдениям из бутстреп-выборок, то каждое конкретное дерево строится на основе примерно двух третей всех наблюдений. Остальная треть наблюдений не используется в обучении, но вполне может быть использована для независимого тестирования: ошибка на таких оставшихся данных (out–of–bag error) является состоятельной оценкой ошибки на контрольной выборке (Джеймс и др., 2016).

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

Их графики легко построить при помощи функции varImpPlot():

varImpPlot(ranfor.a1$finalModel)

–  –  –

На рис. 4.15 приведены два показателя важности: %IncMSE основан на среднем снижении точности предсказания на оставшихся данных, а IncNodePurity – мера среднего увеличения "чистоты узла" дерева (node purity) в результате разбиения данных по соответствующей переменной. В случае деревьев регрессии чистота узла выражается через ошибку RSS.

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

plot(ranfor.a1$finalModel, col="blue", lwd=2) plot(bag.a1$finalModel, col="green", lwd=2, add=TRUE) legend("topright",c("Bagging", "RandomForrest"), col=c("green","blue"), lwd=2) Рис. 4.16. Зависимость ошибки на обучающей выборке от числа агрегируемых деревьев при бэггинге и использовании алгоритма "случайный лес" Бустинг Другим методом улучшения предсказаний является бустинг (boosting), идея которого заключается в итеративном процессе последовательного построения частных моделей. Каждая новая модель обучается с использованием информации об ошибках, сделанных на предыдущем этапе, а результирующая функция представляет собой линейную комбинацию всего ансамбля моделей с учетом минимизации любой штрафной функции.

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

Здесь мы ограничимся обсуждением градиентного бустинга в контексте деревьев регрессии.

Бутстреп-выборки в ходе реализации бустинга не создаются, но вместо этого каждое дерево строится по набору данных {X, r}, который на каждом шаге модифицируется определенным образом. На первой итерации по значениям исходных предикторов строится дерево f1(x) и находится вектор остатков r1. На последующем этапе новое регрессионное дерево f2(x) строится уже не по обучающим данным X, а по остаткам r1 предыдущей модели.

Линейная комбинация прогноза по построенным деревьям дает нам новые остатки r2 r1 + f2(x), и этот итерационный процесс повторяется B раз.

Благодаря построению неглубоких деревьев по остаткам, прогноз отклика медленно улучшается в областях, где одиночное дерево работает не очень хорошо. Такие деревья могут быть довольно небольшими, лишь с несколькими конечными узлами. Параметр сжатия регулирует скорость этого процесса, позволяя создавать комбинации деревьев более сложной формы для "атаки" остатков. Итоговая модель бустинга представляет собой B ансамбль f ( x) = f b ( x).

b =1 В среде R для построения бустинг-моделей на основе деревьев решений можно использовать функцию gbm() из пакета gbm (Generalized Boosted Models).

Процесс моделирования проходит под управлением трех гиперпараметров:

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

2. Параметр сжатия (shrinkage), который корректирует величину вклада каждого дополнительного дерева и контролирует скорость, с которой происходит обучение модели при реализации бустинга. Типичные значения варьируют от 0.01 до 0.001, и их оптимальный выбор зависит от решаемой проблемы. Для достижения хорошего качества предсказаний очень низкие значения требуют очень большого значения B.

3. Число внутренних узлов d (interaction.depth) в каждом дереве, которое контролирует сложность получаемого в результате бустинга ансамбля моделей. По своей сути, параметр d отражает глубину взаимодействий между предикторами в итоговой модели. Если эти взаимодействия не слишком выражены, то хорошо работает d = 1, и тогда дополнительные деревья представляют собой просто "пни" (stump), т.е.

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

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

для решения задач регрессии задается значение "gaussian" – квадратичный штраф, или "laplace" – штраф по абсолютной величине отклонения;

для задач бинарной классификации используют значение "bernoulli" – функция кросс-энтропии, или "adaboost" – экспоненциальный штраф.

Используем значение shrinkage = 0.001, установленное функцией gbm() по умолчанию.

Функция summary() в отношение этого метода выводит список предикторов и соответствующие им значения показателя важности:

library(gbm) set.seed(1) xd - cbind(a1 = algae$a1, x) boost.a1=gbm(a1 ~., data= xd, distribution="gaussian", n.trees=1000,interaction.depth=3) summary(boost.a1)

–  –  –

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

bst:

modelLookup("bstTree") 1 bstTree mstop # Boosting Iterations TRUE TRUE FALSE 2 bstTree maxdepth Max Tree Depth TRUE TRUE FALSE 3 bstTree nu Shrinkage TRUE TRUE FALSE Параметры, оптимизируемые методом bstTree, имеют несколько отличающиеся названия, но фактически эквивалентный содержательный смысл.

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

library(bst) boostFit.a1 - train(a1 ~., data= xd, method='bstTree', trControl=trainControl(method = "cv"), preProc=c('center','scale')) plot(boostFit.a1) Рис. 4.17. Зависимость ошибки от числа агрегируемых деревьев при бустинге (по результатам перекрестной проверки)

–  –  –

3 50 16.31444 0.4480353 3.605141 0.1762307 3 100 16.67521 0.4249761 3.483872 0.1652263 3 150 17.06778 0.4033277 3.421522 0.1589667 Tuning parameter 'nu' was held constant at a value of 0.1 RMSE was used to select optimal model using smallest value.

The final values used for the model were mstop = 50, maxdepth = 2 and nu = 0.1.

Тестирование моделей с использованием дополнительного набора данных Используем для прогноза набор данных (Torgo, 2011) из 140 наблюдений, который мы уже применяли в предыдущем разделе. Данные, с восстановленными пропущенными значениями мы сохранили ранее в файле algae_test.RData (см раздел 4.1).

Выполним прогноз для набора предикторов тестовой выборки и оценим точность каждой модели по трем показателям: среднему абсолютному отклонению (MAE), корню из среднеквадратичного отклонения (RSME) и квадрату коэффициента детерминации Rsq = 1 - NSME, где NSME – относительная ошибка, равная отношению среднему квадрату отклонений от модельных значений и от общего среднего:

load(file="algae_test.RData") # Загрузка таблиц Eval, Sols y - Sols$a1 EvalF - as.data.frame(model.matrix(y ~.,Eval)[,-1]) # Функция, выводящая вектор критериев ModCrit - function (pred, fact) { mae - mean(abs(pred-fact)) rmse - sqrt(mean((pred-fact)^2)) Rsq - 1-sum((fact-pred)^2)/sum((mean(fact)-fact)^2) c(MAE=mae, RSME=rmse, Rsq=Rsq) } Result - rbind( bagging = ModCrit(predict(bag.a1, EvalF), Sols[,1]), ranfor = ModCrit(predict(ranfor.a1, EvalF), Sols[,1]), bst.gbm = ModCrit(predict(gbmFit.a1, EvalF), Sols[,1]), bst.bst = ModCrit(predict(boostFit.a1, EvalF), Sols[,1])) Result

MAE RSME Rsq bagging 10.101308 14.53114 0.4969259 ranfor 9.968411 14.22920 0.5176151 bst.gbm 10.034274 14.34638 0.5096377 bst.bst 9.971036 14.38706 0.5068524

4.5. Сравнение построенных моделей и оценка информативности предикторов Разделы 4.1-4.4 содержат подробную информацию о результатах тестирования различных типов моделей регрессии в идентичных условиях и на одном и том же примере, обобщенную в файле Models.txt.

Сравнительная точность прогноза, оцениваемая по квадрату коэффициента детерминации Rsquared при 10-кратной перекрестной проверке (ось Y) и на контрольной выборке из 140 наблюдений (ось X), представлена на рис.

4.18:

Models - read.delim('Models.txt',header=T) plot(Models$Rsq,Models$Rsquared,pch=CIRCLE-16, col=8-Models$col, cex=2.5, xlab="Rsquared на дополнительной выборке", ylab="Rsquared при кросс-проверке") text(Models$Rsq,Models$Rsquared,rownames(Models), pos=4, font=4,cex=0.8) legend('bottomright',c('Бэггинг/бустинг','Деревья', 'Регрессия kNN','PLS/PCR','Лассо','Линейные модели'), col=2:7, pch=CIRCLE-16, cex=1) Рис. 4.18. Результаты тестирования точности моделей регрессии различного класса при кросс-проверке и на внешнем дополнении Хотя использованный тестовый пример представляется вполне типичным, отсутствие повторностей нашего вычислительного эксперимента не дает нам права делать далеко идущие выводы. Однако некоторые достоинства и недостатки отдельных типов прогнозирующих моделей проявились достаточно четко.

Бесспорными лидерами по точности прогноза явились модели случайного леса (12), бустинга (13-14) и бэггинга (11), основанные на ансамблях деревьев решений. Неплохо себя проявили также одиночные деревья и регрессия k ближайших соседей. Модели, основанные на обобщенных характеристиках обучающей выборки или ее преобразованиях, такие как регрессия на главные компоненты (7), PLS, лассо, дерево условного вывода (10), показали сравнительно неплохие результаты при перекрестной проверке, но оказались не столь хорошими "предсказателями" на "свежих" данных, возможно, не столь похожих на обучающие.

Метод случайного леса не только позволяет построить превосходные модели прогнозирования, но и выполнить такую работу, как селекция набора всех информативных признаков (finding all relevant variables).

В общем случае эта проблема решается с использованием трех возможных подходов (https://habrahabr.ru/post/264915/):

методы фильтрации (filter methods), которые рассматривают каждую переменную независимо и, в некоторой степени, изолированно, оценивая ее по тому или иному показателю (информационные или статистические критерии, минимальная избыточность при максимальной релевантности mRmR и др.);

адаптационные методы осуществляющие (wrapper methods), направленный перебор разных подмножеств признаков и оценка их по заданному критерию (в разделе 4.1 рассматривались три таких алгоритма:

пошаговый, генетический и RFE);

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

В адаптационных методах требуется регрессионная модель (или классификатор), которая используется как черный ящик, возвращая признаки, ранжированные по какому-нибудь удобному критерию – см. рис. 4.15. По практическим соображениям эта модель должна быть в вычислительном отношении быстрой, эффективной и простой, а также мало зависящей от параметров пользователя. Пакет Boruta (бог леса в славянской мифологии), включающий одноименную функцию, реализует адаптационный алгоритмдля модели случайного леса (Kursa, Rudnicki, 2010).

Как и функция varImp() из пакета caret, функция Boruta() оценивает меру информативности каждой переменной в виде дополнительной ошибки регрессии, вызванной исключением этой переменной из модели. Среднее µ этой дополнительной ошибки и его стандартное отклонение рассчитываются по всем деревьям в лесу, которые используют оцениваемый признак для прогнозирования. Оценка Z = µ/ может непосредственно использоваться для ранжирования признаков, однако она не является мерой статистической значимости, поскольку не распределена нормально.

Для того, чтобы оценить, является ли ценность признака существенной, а не обусловленной случайными флуктуациями (т.е. проверить гипотезу Н0: Z = 0), алгоритм Boruta использует внешнее дополнение, полученное в ходе рандомизации. Исходная таблица переменных расширяется таким образом, что в пару каждому предиктору создается соответствующий "теневой" (shadow) признак, вектор которого получен случайным перемешиванием значений основного признака между строками. Для таких признаков корреляция с откликом отсутствует. Далее запускается алгоритм множественного построения моделей с использованием этой вдвое расширенной таблицы и вычисляется ценность всех признаков Z.

Информативная ценность теневого признака может отличаться от нуля только из-за случайных флуктуаций, поэтому множество значений Z теневых признаков служит эталоном того, чтобы решить, какие признаки действительно информативны. Для этого вычисляется "теневой порог" MZSA (maximum Z score among shadow attributes) и признаки, для которых Z MZSA объявляются значимо важными (important), в то время как остальные – незначимыми (unimportant). Поскольку Boruta – высокозатратный с вычислительной точки зрения алгоритм и не всегда удается в условиях сокращенного числа итераций Random Forrest достичь полной ясности, некоторые признаки могут быть обозначены как неопределенные (tentative).

Используем метод Boruta для оценки важности предикторов, определяющих обилие водорослей в реках разного типа.

load(file="algae.RData") # Загрузка таблицы algae - раздел 4.1 library(Boruta) algae.mod - as.data.frame(model.matrix(a1 ~., data=algae[,1:12])[,-1]) algae.mod - cbind(algae.mod, a1=algae$a1) algae.Boruta - Boruta(a1 ~., data = algae.mod, doTrace = 2, ntree = 500) getConfirmedFormula(algae.Boruta) Boruta performed 99 iterations in 1.376562 mins.

8 attributes confirmed important: Chla, Cl, mxPH, NH4, NO3, Sizesmall, oPO4, PO4, ;

5 attributes confirmed unimportant: seasonspring, seasonsummer, seasonwinter, sizemedium, speedlow;

2 tentative attributes left: mnO2, speedmedium;

a1 ~ sizesmall + mxPH + Cl + NO3 + NH4 + oPO4 + PO4 + Chla Таким образом, в результате 99 итераций создания моделей Random Forrest по 500 деревьев в каждой статистически значимыми были признаны 8 переменных, которые были включены в объект formula.

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

attStats(algae.Boruta) meanImp medianImp minImp maxImp normHits decision seasonspring 0.3855680 0.09735517 -0.8022415 2.1263346 0.01010101 Rejected seasonsummer -0.3086411 0.05693780 -2.1205250 0.9013821 0.00000000 Rejected seasonwinter -0.4294279 -0.45213668 -1.8695717 1.7123137 0.01010101 Rejected sizemedium 1.0810475 0.92457355 -0.5714529 3.0446745 0.04040404 Rejected sizesmall 3.4921206 3.52485556 0.5738518 5.9942249 0.73737374 Confirmed speedlow 1.5411720 1.51556358 -0.8321455 3.6901298 0.16161616 Rejected speedmedium 1.9939529 2.00780207 -0.4444101 3.6953856 0.42424242 Tentative mxPH 4.6936750 4.65861782 2.6808697 6.7698339 0.90909091 Confirmed mnO2 2.2344718 2.34487402 -0.3452873 4.7319744 0.39393939 Tentative Cl 12.2215939 12.31395995 10.0974482 15.2918406 1.00000000 Confirmed NO3 4.5193051 4.55940735 1.8457780 7.5044228 0.90909091 Confirmed NH4 11.4733743 11.46492878 9.3601708 13.2503691 1.00000000 Confirmed oPO4 14.0490915 13.97411465 12.2834284 16.0710432 1.00000000 Confirmed PO4 16.0725740 15.97896519 14.2029194 18.3010180 1.00000000 Confirmed Chla 7.5156806 7.57219524 5.4705756 9.8583923 1.00000000 Confirmed Разумеется, более наглядно результаты выглядят на диаграмме (рис.

4.19). К сожалению, разработчики пакета сделали трудночитаемым перечень предикторов по оси Х, и нам пришлось немного повозиться, чтобы исправить этот недостаток:

plot(algae.Boruta, xlab = "", xaxt = "n") lz-lapply(1:ncol(algae.Boruta$ImpHistory),function(i) algae.Boruta$ImpHistory[is.finite(algae.Boruta$ImpHistory[,i]),i]) names(lz) - colnames(algae.Boruta$ImpHistory) Labels - sort(sapply(lz,median)) axis(side = 1,las=2,labels = names(Labels), at = 1:ncol(algae.Boruta$ImpHistory), cex.axis = 0.7) Рис. 4.19. Ранжирование предикторов с использованием алгоритма Boruta;

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

4.6. Деревья регрессии с многомерным откликом В разделе 4.3 мы рассмотрели деревья CART, прогнозирующие конкретное значение одной зависимой переменной. Развитием этих идей являются деревья многомерной классификации и регрессии (MRT, Multivariate Regression Trees – De’Ath, 2002). Если в случае обычной регрессии строится модель, прогнозирующая значения одномерного вектора, то многомерный отклик задается в виде двумерной таблицы, содержащей несколько столбцов наблюдаемых признаков. Ставится задача оценить, какие предикторы и в какой степени влияют на совокупную изменчивость количественных соотношений между отдельными компонентами отклика.

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

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

Искомым критерием, минимизирующим внутригрупповые различия, может быть, например, сумма квадратов отклонений SS D = ( yij y j ) 2, где yij ij

– значение показателя отклика j для наблюдения i; y j – средние значения этого показателя для формируемого кластера, куда включается i-е наблюдение. Геометрически SSD можно представить как сумму евклидовых расстояний от объединяемых объектов до центра их группировки.

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

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

Для переменных a1-a7, описывающих численности 7 групп водорослей, выполним преобразование Бокса-Кокса и шкалирование, что даст нам возможность соизмерить между собой значения различных показателей и корректно вычислить статистики SSD:

load(file="algae.RData") # Загрузка таблицы algae - раздел 4.1 library(caret) Transal = preProcess(algae[,12:18], method = c("BoxCox", "scale")) Species = predict(Transal, algae[,12:18]) Построение дерева MRT будем осуществлять с использованием функции mvpart() из пакета mvpart.

В левой части формулы, задающей структуру модели, представим матрицу из трансформированных численностей семи групп водорослей Species, а в правой части – независимые гидрохимические показатели рек algae[, 1:11]:

library(mvpart) spe.mvpart - mvpart(data.matrix(Species) ~., algae[,1:11], xv="pick", xval=nrow(Species), xvmult = 1, margin=0.08, which=4, bars = TRUE) Здесь параметры xval и xvmult соответствуют числу блоков и повторностей перекрестной проверки, т.е. при xval = nrow(Species) осуществляется скользящий контроль. Аргументы margin (ширина полей), which (расположение текста) и bars (вывод столбиковых диаграмм) задают атрибуты визуализации дерева.

Важный параметр xv определяет условие выбора числа листьев дерева: при значении "1se" оно определяется автоматически по результатам перекрестной проверки, а при "pick" размеры дерева можно выбрать интерактивно, выполнив щелчок мышью по следующему графику:

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

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

Поэтому, чтобы выполнить содержательный анализ зависимости свойств найденных групп водорослей от внешних факторов, разделим все множество наблюдений на 4 кластера и изобразим полученное дерево графически:

plot(spe.mvpart) ; text(spe.mvpart)

Рис. 4.21. Построенное дерево с многомерным откликом

Рассмотрим более внимательно некоторые нюансы итеративной процедуры построения модели MRT:

summary(spe.mvpart) # Протокол приведен с сокращениями CP nsplit rel error xerror xstd 1 0.09564230 0 1.0000000 1.0102813 0.09977143 2 0.04933011 1 0.9043577 0.9492214 0.10004843 3 0.03391020 2 0.8550276 0.9647968 0.09961182 4 0.03286729 3 0.8211174 0.9886724 0.10069152 Node number 1: 200 observations, complexity param=0.0956423 left son=2 (47 obs) right son=3 (153 obs)

Primary splits:

Cl 9.0275 to the left, improve=0.09564230, (0 missing) PO4 31.97933 to the left, improve=0.09327022, (0 missing) oPO4 27.375 to the left, improve=0.08805637, (0 missing) Node number 2: 47 observations Means=1.89,0.1381,0.6381,0.3959,0.07981,0.06642,0.1394 Node number 3: 153 observations, complexity param=0.04933011 left son=6 (102 obs) right son=7 (51 obs)

Primary splits:

mnO2 7.85 to the right, improve=0.06307793, (0 missing) NH4 7588.8 to the left, improve=0.05783644, (0 missing) oPO4 448.125 to the left, improve=0.05719325, (0 missing) Node number 6: 102 observations, complexity param=0.0339102 left son=12 (86 obs) right son=13 (16 obs)

Primary splits:

NO3 6.194 to the left, improve=0.07251793, (0 missing) mxPH 7.55 to the right, improve=0.05984036, (0 missing) Cl 44.0875 to the left, improve=0.05770594, (0 missing) Node number 7: 51 observations Means=0.307,0.943,1.1,0.8536,0.3646,0.2085,0.7994, Summed MSE=7.241111 Node number 12: 86 observations Means=0.5884,0.8402,0.3861,0.3035,0.9321,0.6551,0.5444 Node number 13: 16 observations Means=0.2169,0.5259,0.2977,0.1231,2.044,2.013,0.1636 На первом шаге Node number 1 рассматриваются все варианты разбиения исходной выборки на две части при разных опорных значениях независимых факторов и выбирается такой из них (в нашем случае хлориды Cl 9.0275), который в наибольшей мере обеспечивает статистическую гомогенность формируемых кластеров (complexity param=0.0956423).

Кластер 2 из 47 наблюдений на последующих шагах не разбивается, а 153 наблюдения слева вновь делятся на два подмножества по условию mnO2

7.85 в соотношении 51 + 102. Последнее подмножество по условию NO3

6.194 в свою очередь делится на два кластера, на чем итерации завершаются.

Анализ многомерного отклика сообщества водорослей с помощью деревьев MRT предоставляет исследователю много дополнительных возможностей интерпретации результатов.

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

# Относительная доля групп водорослей в кластерах groups.mrt - levels(as.factor(spe.mvpart$where)) leaf.sum - matrix(0, length(groups.mrt), ncol(Species)) colnames(leaf.sum) - colnames(Species) rownames(leaf.sum) - groups.mrt for(i in 1:length(groups.mrt)){ leaf.sum[i,] - apply(Species[which (spe.mvpart$where==groups.mrt[i]),], 2, sum) } leaf.sum a1 a2 a3 a4 a5 a6 a7 2 88.81706 6.492446 29.991924 18.608213 3.750967 3.121594 6.552211 5 50.60338 72.260194 33.201232 26.101301 80.158567 56.334476 46.815353 6 3.47099 8.414790 4.763592 1.969482 32.704161 32.202155 2.617007 7 15.65459 48.094876 56.083746 43.532352 18.594652 10.634001 40.767158 Разумеется, эти суммы не соответствуют значениям абсолютных численностей водорослей, поскольку набор исходных данных подвергался преобразованию Бокса-Кокса и шкалированию.

opar - par() # Вывод диаграммы типа "разрезанный пирог" par(mfrow = c(2,2)) ; for(i in 1:length(groups.mrt)){ pie(leaf.sum[i, which(leaf.sum[i,]0)], radius = 1, main = paste("Кл. №", groups.mrt[i])) } par(opar)

Рис. 4.22. Диаграммы долей средних численностей групп водорослей

Мы получили диаграммы долей усредненных численностей групп водорослей для четырех кластеров, составивших "листья" МRТ. Их визуальный анализ показывает, что подмножество 2 составлено с явным доминированием группы a1, подмножество 6 – с доминированием групп a5 и a6, а два остальных кластера не имеют столь характерных признаков. Видимо, перекрестная проверка все же имела все основания рекомендовать нам разбиение только на два кластера.

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

# Вывод диаграммы РСА rpart.pca(spe.mvpart) Рис. 4.23. Диаграмма четырех групп водорослей в пространстве двух главных компонент Из центра координат диаграммы проводятся дополнительные оси ординации, косинусы углов между которыми соответствуют коэффициентам корреляции между каждой парой из 7 групп водорослей. Проекции точек на каждую ось ординации определяют характер распределения показателя по кластерам с разными внешними факторами. Например, на ось а1 проецируются в основном красные точки кластера 2, объединяющего наблюдения с низким содержанием хлоридов Cl 9.0275.

Настоящий раздел можно рассматривать в качестве "реквиема" по этому весьма оригинальному и практически полезному пакету: mvpart был удален из репозитория CRAN и для версий R выше 3.1 отсутствует. Это уже не первый случай, когда авторы перестают поддерживать важные и популярные функциональные компоненты: например, та же участь постигла удачный пакет lmRepm, в котором функции построения линейных моделей lm() и anova() при оценке статистических критериев используют перестановочные алгоритмы и не столь жестко связаны с предположениями о характере распределения данных.

Выхода из подобной ситуации может быть два. Во-первых, ничто не мешает установить на компьютер несколько версий R (так, один авторов этой книги при любой возможности отдает предпочтение старенькой, но симпатичной версии 2.12). Другая возможность – осуществить компиляцию последней версии пакета mvpart_1.6-2 из архива CRAN в вашей текущей версии R (см. рекомендации в сообщении на http://stackoverflow.com).

5. БИНАРНЫЕ МАТРИЦЫ И АССОЦИАТИВНЫЕ ПРАВИЛА

5.1. Классификация в бинарных пространствах с использованием классических моделей Большое количество задач связано с нахождением закономерностей в пространствах бинарных переменных: расшифровка геномного кода, анализ пиков на спектрофотометрах, обработка социологических анкет с ответами "Да/Нет" и т.д. Методы анализа таких выборок рассмотрим на не слишком серьезном примере, фигурирующем в книге Дюка и Самойленко (2001).

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

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

Прежде всего, выделим признаки, характеризующие изображенные лица.

Это следующие характеристики:

х1 (голова) – круглая (1) или овальная (0);

х2 (уши) – оттопыренные (1) или прижатые (0);

х3 (нос) – круглый (1) или длинный (0);

х4 (глаза) – круглые (1) или узкие (0);

х5 (лоб) – с морщинами (1) или без морщин (0);

x6 (носогубная складка) – есть (1) или нет (0);

х7 (губы) – толстые (1) или тонкие (0);

х8 (волосы) – есть (1) или нет (0);

х9 (усы) – есть (1) или нет (0);

х10 (борода) – есть (1) или нет (0);

х11 (очки) – есть (1) или нет (0);

х12 (родинка на щеке) – есть (1) или нет (0);

х13 (бабочка) – есть (1) или нет (0);

х14 (брови) – подняты кверху (1) или опущены книзу (0);

х15 (серьга) – есть (1) или нет (0);

х16 (курительная трубка) – есть (1) или нет (0).

Загрузим исходную матрицу данных, соответствующую изображенным лицам, cтроки которой соответствуют объектам (n = 16), а столбцы – выделенным бинарным признакам (m = 16). Объекты с номерами 1-8 относятся к классу 1, а с номерами 9-16 – к классу 2.

DFace - read.delim(file="Faces.txt", header=TRUE, row.names=1) head(DFace,8) x1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 x12 x13 x14 x15 x16 Class Попытаемся решить поставленную задачу с помощью классических статистических методов. Множественная логистическая регрессия – один из распространенных подходов к моделированию значений отклика, заданного альтернативной шкалой 0/1 (Мастицкий, Шитиков, 2015).

Если x – вектор значений m независимых переменных, то вероятность P того, что отклик y примет значение 1, может быть описана моделью P(y = 1|x) p(x) = 0 + 1 x1 +... + m xm,

–  –  –

Прогноз Факт 1 2 [1] "Точность=100.00%" Рис. 5.2. Диаграмма вероятностей, предсказанных логит-регрессией Все члены обучающей выборки с использованием полученной регрессионной модели были верно классифицированы и с максимальной вероятностью (см. рис. 5.2). Перекрестную проверку построенной модели можно выполнить, например, с использованием функции cv.glm() из пакета boot.

Если параметр K, определяющий число блоков разбиения, опущен, то выполняется скользящий контроль:

library(boot) cv.glm(DFace, logit.step)$delta [1] 2.564444e-17 2.404152e-17 В результате получены среднеквадратичное отклонение по координате логита от проверочных наблюдений до гиперплоскости, найденной на обучающих выборках, а также его скорректированное значение (adjusted cross-validation estimate of prediction error). Поскольку ошибки перекрестной проверки ничтожно малы и все объекты в ее ходе были правильно опознаны, то в целом полученную модель можно считать высоко эффективной.

Мы убедились в том, что логистическая регрессия и функция glm() отлично работают с индикаторными переменными. Однако строгие аналитики (Дюк, Самойленко, 2001) считают, что «вряд ли такое правило способно удовлетворить разработчика интеллектуальной системы. Оно формально, не дает нового знания, а попытка дать интерпретацию коэффициентам регрессии приводит к сомнительным выводам. Непонятно, например, почему вес формы носа (х3) в два раза превышает вес оттопыренности ушей (х2) и т. д.». Несмотря на высказанные сомнения, мы достигли желаемой цели – достоверное и эффективное правило классификации построено, а из анализа коэффициентов модели легко сделать вывод, что признаки x3, x4, x5 и x14 склонны иметь демократы, а показатели x2, x7 и x13 характеризуют патриотические устремления.

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

library("MASS") Fac.lda - lda(DFace[,1:16], grouping = DFace$Class) pred - predict(Fac.lda, DFace[,1:16]) table(Факт=DFace$Class, Прогноз=pred$class) Acc = mean(DFace$Class == pred$class) paste("Точность=", round(100*Acc, 2), "%", sep="")

Warning message:

In lda.default(x, grouping,...) : variables are collinear Прогноз Факт 1 2 [1] "Точность=75.00%" Обратим внимание на неприятное сообщение о коллинеарности переменных и мало впечатляющие результаты прогноза.

Стохастические модели многомерных двоичных случайных величин основываются на распределении Бернулли X = {x1, …, xm} ~ Bem(µ, ), где µ

– вектор математических ожиданий, а включает (2m – m – 1) параметров, оценивающих взаимодействие между переменными xj. Как и в одномерном случае, дисперсии полностью определяются через средние Var(xj) = µj (1 - µj).

Во многих случаях, особенно для предсказания на основе выборок малого объема, но с большим числом предикторов, представляется разумным следовать "наивному" байесовскому подходу, т.е. игнорировать зависимости между индивидуальными переменными xj. Тогда количество параметров многомерной модели Бернулли существенно сокращается, и функция объединенной плотности вероятности определяется произведением µj и диагональной ковариационной матрицей Var(X) = diag{µj (1 - µj)}.

С учетом этих положений разработан метод бинарного дискриминантного анализа (binDA – binary discriminant analysis; Gibb, Strimmer, 2015). Дискриминантная функция, полученная на основе теоремы

Байеса, для каждого класса с меткой y имеет следующий вид:

dy(x) = log y + log( Pr(x|y)) + C, где y - априорная вероятность класса y, C – константа, а совокупная условная вероятность предикторов в зависимости от значения отклика вычисляется по формуле

–  –  –

Тестируемый объект относится к тому классу y, дискриминантная функция которого dy(x) максимальна. Параметром алгоритма является – степень ослабления дисбаланса в частотах классов (lambda.freqs – shrinkage intensity for the frequencies). При lambda.freqs=0 используются эмпирические частоты классов y = ny/n, а при lambda.freqs=1 все вероятности принимаются равными.

Важным разделом алгоритма является устранение подмножества предикторов, влияние которых на значение метки y статистически незначимо.

Разность между значениями логарифма функции правдоподобия полной модели и нуль-модели без параметров равна относительной энтропии, или дивергенции D Кульбака-Лейблера (Kullback-Leibler). Относительные вклады Sj каждого индивидуального предиктора в величину D определяют меру значимости j-й переменной в результате классификации. Отсюда легко найти матрицу t-значений, соответствующих разностям между общим средним и средними для каждой переменной относительно каждого класса.

Алгоритм binDA реализован в пакете binda для R.

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

library(binda) Xtrain - as.matrix(DFace[,1:16]) is.binaryMatrix(Xtrain) # Проверяем бинарность матрицы Class = as.factor(DFace$Class) binda.fit = binda(Xtrain,Class,lambda.freq=1) pred - predict(binda.fit, Xtrain) table(Факт=DFace$Class, Прогноз=pred$class) Acc = mean(DFace$Class!= pred$class) paste("Точность=", round(100*Acc, 2), "%", sep="") Прогноз Факт 1 2 [1] "Точность=93.75%" При равных апостериорных вероятностях классов 0.5/0.5 объект №4 был ошибочно отнесен ко второму классу.

Выполним теперь оценку информативности предикторов по t-критерию (рис.

5.3):

# ранжируем предикторы binda.x - binda.ranking(Xtrain, Class) plot(binda.x, top=40, arrow.col="blue", zeroaxis.col="red", ylab="Предикторы", main="")

Рис. 5.3. Оценка значимости предикторов модели binda



Pages:     | 1 || 3 | 4 |   ...   | 5 |



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

«Инструкция по установке модуля интеграции "1С-Битрикс" и службы доставки ShopLogistics Оглавление Описание модуля Требования Установка модуля 1 шаг. Установка модуля штатными средствами "1С-Битрикс" 2 шаг. Заполнение ключей доступа, сопоставление Особеннос...»

«ІДЕОЛОГІЯ І ПОЛІТИКА ИДЕОЛОГИЯ И ПОЛИТИКА IDEOLOGY AND POLITICS Наши Авторы Доктор наук, старший консультант, dutchTest, Нидерланды; Стивен Беккер консультант Международного альянса USETI. Сфера научных интересов:...»

«МИНИСТЕРСТВО ОБРАЗОВАНИЯ И НАУКИ РОССИЙСКОЙ ФЕДЕРАЦИИ Федеральное государственное образовательное учреждение высшего образования "Сибирский государственный университет геосистем и технологий" (СГУГиТ) Рассмотрено Утверждаю на заседании Ученого совета СГУГиТ Ректор _ А.П. Карпик 26 января 2...»

«Касьянова Людмила Олеговна (автор) Книга № 1 "Диалог Маленькой души" ВВЕДЕНИЕ Дорогой читатель! Книга – откровение попала к Вам в руки и это неслучайно, ведь случайностей не бывает. Её автор искренне молил...»

«УДК 82-93 ББК 84(2Рос-Рус) Л 49 Оформление серии О. Горбовской Лермонтов М. Ю. Герой нашего времени / Михаил Лермонтов. — Л 49 М. : Эксмо, 2014. — 320 с. — (Классика в школе). ISBN 978-5-699-72936-4 Перед вами книга из серии "Классика в школе", в которой собраны все произведения, изучающи...»

«RT-G32 беспроводный роутер Руководство пользователя RT-_UM_rus.indd 1 1//0 :1:1 AM R4264 первая редакция Ноябрь 2008 Copyright © 2008 ASUSTeK Computer Inc. Все права защищен...»

«1 ПОЯСНИТЕЛЬНАЯ ЗАПИСКА Рабочая программа разработана в соответствии со ст. 2, 12, 13, 47, 48 Закона Российской Федерации "Об образовании в Российской Федерации" от 29.12.2012г 273-ФЗ; ФГОС НОО, утвержденный приказом Минобрнауки РФ от 06.10.20...»

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

«МИНИСТЕРСТВО ОРАЗОВАНИЯ СТАВРОПОЛЬСКОГО КРАЯ ГОСУДАРСТВЕННОЕ БЮДЖЕТНОЕ ОБРАЗОВАТЕЛЬНОЕ УЧРЕЖДЕНИЕ СРЕДНЕГО ПРОФЕССИОНАЛЬНОГО ОБРАЗОВАНИЯ "ГЕОРГИЕВСКИЙ РЕГИОНАЛЬНЫЙ КОЛЛЕДЖ "ИНТЕГРАЛ" ОСНОВНАЯ ПРОФЕССИОНАЛЬНАЯ ОБРАЗОВАТЕЛЬНАЯ ПРОГРАММА ГБОУ СПО "Георгиевский региональный колледж "Интеграл" по профессии среднего пр...»

«ЗАКРЫТОЕ АКЦИОНЕРНОЕ ОБЩЕСТВО ТРАНСЛЯЦИОННЫЙ УСИЛИТЕЛЬ МОЩНОСТИ ЗВУКОВОЙ ЧАСТОТЫ РА450Тм ТУ 6587-059-02963487-07 (Руководство по эксплуатации) РОССИЯ 630003, г. Новосибирск, ул. Владимировская, 1а ЗАО "НОЭМА" те...»

«области Выпусн I Январь М ар т 1985 Свердловск 1985 t '•) о /, ЦЛЦу листок К ОН ТРОЛЬН Ы Й С РО К О В В О ЗВ РА ТА КНИ ГА Д О Л Ж Н А БЫ ТЬ В О ЗВ РА Щ ЕН А НЕ П О ЗЖ Е У К А ЗА Н Н О ГО З Д Е С Ь С РО К А К олич. пред. в ы...»

«СОДЕРЖАНИЕ ПРЕДИСЛОВИЕ ВВЕДЕНИЕ Раздел I. МОДЕЛИРОВАНИЕ. ИНСТРУМЕНТЫ ОБЩЕГО НАЗНАЧЕНИЯ Тема 1. СИСТЕМНАЯ ГАРАНТИЯ РЕЗУЛЬТАТА МОДЕЛИРОВАНИЯ 1.1 Неопределенность построения моделей 1.2 О методе гарантированного результата Тема 2. ОСНОВЫ СИСТЕМНОГО ПОДХОДА 2.1 Основные категории теории систем 2.2 Полиструктурность и свойства объектов 2.3 К...»

«Подлежит публикации,Л СОГЛАСОВАНО в открытой печати „ Зам. генерального директора Регистраторы электрических Внесены в Государственный процессов цифровые “Парма РП4.08” реестр средств измерений 2Ji S § {~0 Регистрационный № Взамен № Выпускаются по ТУ 4222-008-31920409-01. НАЗНАЧЕНИЕ И ОБЛАСТЬ ПРИМЕНЕНИЯ Регистра...»

«2 Программа подготовки специалистов среднего звена (далее ППССЗ) государственного бюджетного образовательного учреждения среднего профессионального образования Ростовской области "Ростовский технологический техникум сервиса" составлена на основе Федерального государст...»

«Постановление администрации городского округа "Город Йошкар-Ола" от 30.11.2015 № 2223 Об утверждении административного регламента предоставления муниципальной услуги "Перевод жилого помещения в нежилое помещение или нежилого помещения в жилое помещение" В соответ...»

«МИНИСТЕРСТВО ОБРАЗОВАНИЯ РОССИЙСКОЙ ФЕДЕРАЦИИ УТВЕРЖДАЮ Заместитель Министра образования Российской Федерации В.Д.Шадриков 3 марта 2000г. Регистрационный номер 3 – тех/дс ГОСУДАРСТВЕННЫЙ ОБРАЗОВАТЕЛЬНЫЙ СТАНДАРТ ВЫСШЕГО ПРОФЕССИОНАЛЬНОГО ОБРАЗОВАНИЯ Направление подготовки дипломированного специал...»

«Исполнительское творчество Марии Юдиной. Подготовила Соколова Ю.С. Мое выступление посвящено такому явлению в музыкальном мире, как Мария Вениаминовна Юдина. В наш прогрессивный электронный 21 век не трудно найти краткую, лаконичную информацию о любом и...»

«Р. Н. Абалуев, Н. Е. Астафьева, Н. И. Баскакова, Е. Ю. Бойко, О. В. Вязовова, Н. А. Кулешова, Л. Н. Уметский, Г. А. Шешерина ИНТЕРНЕТ-ТЕХНОЛОГИИ В ОБРАЗОВАНИИ ИЗДАТЕЛЬСТВО ТГТУ Министерство образования Российской Федерации Тамбовский региональный центр Федерации Интернет-Образования Тамбовский областной институт повышения квалифика...»

«ПЛАЗМЕННЫЕ ТЕХНОЛОГИИ УПРОЧНЕНИЯ ПРИ РЕМОНТЕ ПОДВИЖНОГО СОСТАВА С.В Петров., д-р техн. наук Эксплуатация нагруженных сопряженных и трущихся узлов подвижного состава сопровождается контактно-усталостными повреждениями поверхностей, проявляющимися в виде мелких или крупных участков уноса мате...»

«Муниципальное бюджетное образовательное учреждение дополнительного образования детей "Детская школа искусств № 1 имени Александра Семеновича Розанова" П РИ Н Я ТО П едагогическим советом М БО У Д О Д "ДШ И № 1 им. А.С. Розан...»








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

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