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

Pages:   || 2 | 3 | 4 | 5 |

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

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

В.К. Шитиков, С.Э. Мастицкий

Классификация, регрессия

и другие алгоритмы Data Mining

с использованием R

Тольятти, Лондон - 2017

Шитиков В.К., Мастицкий С.Э. (2017) Классификация, регрессия и

другие алгоритмы Data Mining с использованием R. 351 с. Электронная

книга, адрес доступа: https://github.com/ranalytics/data-mining

Описана широкая совокупность методов построения статистических моделей

классификации и регрессии для откликов, измеренных в альтернативной,

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

Особое внимание уделяется сравнительной оценке эффективности и поиску оптимальных областей значений гиперпараметров моделей с использованием пакета caret для статистической среды R. Рассматриваются такие алгоритмы Data Mining, как генерация ассоциативных правил и анализ последовательностей.

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

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

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

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

Данная работа распространяется в рамках лицензии Creative Commons «Атрибуция – Некоммерческое использование – На тех же условиях

4.0 Всемирная». Согласно этой лицензии, Вы можете свободно копировать, распространять и видоизменять данное произведение при условии точного указания его авторов и источника. При изменении этого произведения или использовании его в своих работах, Вы можете распространять результат только по такой же или подобной лицензии. Запрещается использовать эту работу в коммерческих целях без согласования с авторами. Более подробная информация о лицензии представлена на сайте www.creativecommons.com © 2017, Владимир Кириллович Шитиков, Сергей Эдуардович Мастицкий

Контактная информация:

В.К. Шитиков Институт экологии Волжского бассейна РАН г. Тольятти, ул. Комзина, 10 Самарская обл., 445003, Россия Сайты авторов: http://www.ievbras.ru/ecostat/Kiril E-mail: stok1946@gmail.com rtutorialsbook@gmail.com http://r-analytics.blogspot.com Информация о книге: http://www.ievbras.ru/ecostat/Kiril/Article/A44/DM.htm Иллюстрация на обложке худ. А. Банных (г. Владимир) СОДЕРЖАНИЕ

–  –  –

От статистического анализа разового эксперимента к Data Mining Экспериментальные данные, представленные в компьютерном формате в виде взаимосвязанных таблиц, нуждаются в таких процедурах обработки анализа и интерпретации, которые, во-первых, делают очевидными потенциально возможные закономерности и связи между отдельными компонентами и, во-вторых, дают возможность предсказать новые факты. В узком плане речь может идти об оценке значения целевого признака y (отклика) для любого объекта а по его описанию х – набору независимых переменных (предикторов).





Однако в более широком смысле затрагиваются традиционно ключевые вопросы многомерного анализа систем:

Можно ли считать идентичными анализируемые объекты и за счет каких признаков можно объяснить их возможные отличия?

Как можно объединить отобранные объекты в группы?

Существует ли пространственная или временная изменчивость описанных объектов и каковы ее структурные особенности?

Изменение каких признаков приводит к систематическим причинноследственным изменениям других?

Как можно осуществить прогноз состояния или поведения анализируемого объекта?

–  –  –

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

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

Попытки ранжирования статистических методов по их эффективности предпринимались неоднократно.

Например, одна из недавно опубликованных статей (Fernandez-Delgado et al., 2014) так и называлась «Do We Need Hundreds of Classifiers to Solve Real World Leskovec et al., Classification Problems?» («Нужны ли нам сотни классификаторов для решения практических проблем классификации?»).

В этом обстоятельном исследовании была изучена эффективность работы 179 методов классификации из 17 "семейств" на 121 наборе данных.

Читатель может рассматривать эту интересную статью даже просто как справочник по методам распознавания.

В проведенном исследовании для каждого метода классификации оценивалась общая верность предсказаний (overall accuracy) и другие показатели эффективности моделей.

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

1. "Случайный лес" (Random Forest);

2. Машины опорных векторов (Support Vector Machines);

3. Искусственные нейронные сети (Artificial Neural Networks);

4. Бустинговые ансамбли моделей (Boosting Ensembles).

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

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

Интересное обсуждение того, что простые методы (на примере классификаторов) при решении практических задач часто превосходят более сложные алгоритмы, можно найти также в широко цитируемой работе Д. Xэнда (Hand, 2006).

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

Типы и характеристики групп моделей Data Mining Предварительно, не вторгаясь в терминологические детали, отметим, что все три дисциплины – Data Mining, машинное обучение (Machine Learning) и статистический анализ – работают на одном и том же предметном поле и используют фактически одни и те же алгоритмы. Общие для них черты и темы гораздо обширней, чем различия, которые носят характер того или иного "уклона": Data Mining включает в сферу своей компетенции практически необходимые приемы работы с большими массивами данных, сетями и проч., машинное обучение склонно к красивой фразеологии по поводу искусственного интеллекта, а статистический анализ ищет обоснование своих процедур в теоретико-вероятностном подходе. В рамках нашего изложения эти дисциплины рассматриваются как синонимы.

Представленное выше многообразие методов построения моделей порождает разнообразие подходов к их группировке.

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

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

Алгоритмы индуктивного распознавания с обучением (Ripley, 1996;

Vapnik, 1995; Agresti, 2007; Hastie et al., 2007; Zaki M., Meira, 2014) предполагают наличие априори заданной выборки прецедентов, позволяющей построить модели статистической связи x y, где y Y, Y – наблюдаемые реализации отклика (responses) или моделируемой случайной величины;

x X, X – множество переменных (predictors, independent variables, features), с помощью которых предполагается объяснить изменчивость переменной y.

Большинство моделей с учителем устроено таким образом, что их можно записать в виде y = f(x, ) +, где f – математическая функция, выбранная из некоторого произвольного семейства, – вектор параметров этой функции, а – ошибки модели, которые обычно сгенерированы несмещенным, некоррелированным случайным процессом. В ходе построения модели по фиксированным выборочным значениям y минимизируют некоторую функцию остатков Q(y, ) и находят – вектор с оптимальными оценками параметров модели.

Варьируя вид функций f() и Q(), можно получать разные модели, из которых предпочтение отдается наиболее эффективным – т.е. моделям, которые дают несмещенные, точные и надежные прогнозы отклика y. Под смещением (bias) понимается различие между рассчитанным по модели прогнозом (например, среднее из оценок всех возможных выборок, которые могут быть взяты из совокупности) и истинным неизвестным значением моделируемой переменной. Точность (accuracy) – различие между оценками отклика, основанного на выборочных данных и истинным значением реализации y. Надежность (precision) – различие между оценкой y по разовой выборке и средним из прогнозов по всему множеству выборок, которые могут быть взяты из той же совокупности.

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

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

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

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

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

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

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

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

Обеспечивая точность предсказания, модели прогнозирования в большинстве случаев существенно теряют в интерпретируемости. Часто бывает затруднительным сделать предметный анализ нескольких десятков коэффициентов моделей МГУА или нейронных сетей. Очевидна стремительно набирающая темп современная тенденция на создание ансамблей из сотен моделей, осуществляющих коллективное предсказание (bagging, random forest, boosting), когда традиционные регрессионные методы, такие как анализ ковариаций, стандартных отклонений, доверительных интервалов коэффициентов и проч. вообще не имеют смысла.

Природа многомерного отклика и его моделирование При изучении современных информационных, биологических и социально-экономических систем характерна оценка взаимных зависимостей между комплексами многомерных переменных, и тогда задачей построения статистической модели является объяснение изменчивости многомерного отклика Y. Часто такой отклик может быть представлен в форме некоторой относительно замкнутой системы элементов N, относящихся к множеству S = {y} различных типов этих элементов (Арапов и др., 1975; Шитиков и др., 2012, глава 5). В разных предметных областях это может быть каталог продаваемых автомобилей, категории обеспеченности населения, словарь встретившихся словоформ, список экологических видов, обнаруженных в водоеме, перечень кандидатов, за которых голосовали на выборах и т.д. С помощью классификатора S объекты разбиваются на подпопуляции (классы), т. е. каждой рубрике y Є S соответствует подмножество N(y), определяющее частоту всех вхождений объектов этого типа в N. Закон Ципфа–Парето, описывающий характер таких распределений, играет практически ту же универсальную роль, что и закон Гаусса в обычных стохастических процессах с конечной дисперсией (Кудрин, 2002) – подробности см. в главе 9.

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

Если модель обычной регрессии оценивает условную вероятность математического ожидания y, то в случае многомерного отклика Y анализируется влияние предикторов на структурную изменчивость корреляционной матрицы (или иной матрицы дистанций). К настоящему времени с использованием такого подхода разработаны методы непараметрического дисперсионного анализа (McArdle, Anderson, 2001), анализа избыточности и канонического корреспондентного анализа (Legendre, Legendre, 2012), построения деревьев с многомерным откликом (De’Ath, 2002) и др. В экологических исследованиях процедуры оптимального проецирования многомерных данных по осям главных компонент и одновременной оценки коэффициентов линейных канонических моделей объединяются под названием прямая ординация.

–  –  –

Stanton, 2013 Список только признанных фундаментальных монографий типа «Используем R!» (Use R!) насчитывает уже несколько сотен томов. В многочисленных фундаментальных руководствах детально описаны функции и пакеты R, использующие широкий набор статистических методов, таких как классические линейные и нелинейные регрессионные модели, различные формы дисперсионного анализа, алгоритмы иерархической кластеризации, анализ временных рядов, деревья классификации и регрессии (Faraway, 2006;

Fox, 2008; Zuur et al., 2009; Myers, Patil, 2012; Mount, Zumel, 2014, 2016; Kuhn, Johnson, 2013; Stanton, 2013; Кабаков, 2014; Сhiu, 2015, Kassambara, 2017).

Отдельно можно отметить библиографические источники, в названии которых имеется термин "Data Mining": Torgo (2011), Zhao (2012), Zhao, Cen (2014). Большинство из этих книг с разной степенью легальности доступно для скачивания в Интернете (список 9 свободно распространяемых пособий приводится, в частности, на www.kdnuggets.com).

–  –  –

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

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

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

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

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

с GitHub-репозитория книги https://github.com/ranalytics/data-mining или сайта http://www.ievbras.ru/ecostat/Kiril/R/DM/Scripts_data.zip Мы хотели сделать свое изложение, по возможности, разноплановым и иллюстрированным. Например, убедительным доказательством наличия закономерности в данных часто может служить, их качественная графическая иллюстрация (Зиновьев, 2000; Maindonald, Braun, 2010). Поэтому, где было возможно, мы старались включить в наши скрипты примеры визуализации данных с использованием пакета qqplot2 (Маcтицкий, 2016) или других пакетов, поддерживающих эту графическую систему.

В одной из дискуссий в Интернете промелькнула мысль «Data Mining – это обработка данных без применения статистики». Видимо тут имелось в виду скупое употребление псевдо-теоретической фразеологии вроде «истинных параметров некой мифической случайной величины или теоретической функции регрессии», которые необходимо оценить. Избегнуть подобных оборотов постарались и мы, чем рискуем навлечь недовольство ортодоксальных статистиков.

К счастью, большинство описанных ниже методов не апеллируют к вероятностной природе обрабатываемых данных и их использование не всегда связано с (безосновательными) надеждами на закон нормального распределения. Однако это нисколько не мешает легко сгенерировать с их помощью продуктивные научные гипотезы, которые могут оказаться нетривиальными, практически полезными и удобными для интерпретации. В тех же случаях, когда проверка статистической значимости необходима, мы (следуя разработчикам пакетов R) активно использовали методы ресэмплинга: бутстреп, рандомизацию и перекрестную проверку (Эфрон, 1988; Edgington, 1995; Davison, Hinkley, 2006; Шитиков, Розенберг, 2014).

Эффективными средствами получения оценок параметров и их апостериорных распределений является байесовский подход и итерационные методы Монте-Карло, использующие цепи Маркова (MCMC – Monte Carlo Markov chain). Их разработка иногда трактуется как наиболее существенный прорыв в статистике за последние несколько десятков лет. Однако мы не стали касаться этой важной темы, требующей подробного и обоснованного обсуждения.

В настоящее пособие включена совокупность методических сообщений, опубликованных авторами за последнее время в блоге «R: Анализ и визуализация данных» (http://r-analytics.blogspot.com). Нам показалась целесообразной идея представить для удобства читателей весь этот несколько разобщенный материал в концентрированной форме, а также расширить некоторые разделы для полноты изложения.

Прокомментируем кратко содержание глав книги:

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

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

Глава 3 содержит подробное описание функций пакета caret, а разделы главы 4 показывают, как с его помощью выполнить построение и тестирование различных моделей регрессии.

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

Главы 6 и 7 объединяют описание и особенности реализации различных алгоритмов классификации на два и более класса. Как и в главе 4, особое внимание уделено процедурам тестирования и подстройки гиперпараметров моделей.

Глава 8 посвящена сразу трем проблемам: методам классификации порядковых данных, созданию комплексных прогнозов на основе "коллектива" моделей и построению различных версий моделей для счетных данных по материалам сайта http://www.highstat.com и книги А. Зуура с коллегами (Zuur et al., 2009).

Глава 9 описывает различные алгоритмы непрямой и прямой ординации многомерных наблюдений. Предлагается процедура построения кумулятивной кривой распределения чувствительности видов (SSD).

Глава 10 (за исключением раздела о самоорганизующихся картах Кохонена) представляет собой изложение материалов сайта http://www.sthda.com/ и книги А. Кассамбары (Kassambara, 2017) по кластерному анализу в среде R.

Демонстрации возможностей пакета rattle посвящена целиком глава 11, основанная на публикациях Г. Вильямса (Williams, 2009, 2011).

Как видно из рубрикатора глав, в книге довольно много заимствований из различных источников и, вероятно, она имела бы трудности при проверке в системе «Антиплагиат». Однако все в науке в той или иной степени плагиат, поскольку подавляющее большинство исследований лишь опровергает или подтверждает уже высказанные гипотезы. Вспомним, например, историю об Амбруазе Паре, который многое переписал у Везалия, а потом придумал блестящую формулу: «Я просто зажег свою свечку от соседней горящей свечи».

По вполне понятным причинам мы не смогли охватить в книге весь перечень методов построения статистических моделей и соответствующих пакетов. Например, в спектральном анализе и хемометрике (https://ru.wikipedia.org) чрезвычайно популярны методы, основанные на технике условного математического ожидания, и, наряду с методом частных наименьших квадратов (PLS) и регрессией на главные компоненты (PCR), рекомендуется использовать параллельный факторный анализ (PARAFAC), анализ Тюкера (Tucker) и др.

Мы также совсем не затронули такие обширные разделы Data Mining, как анализ временных рядов и моделирование пространственнораспределенных данных. Мы сожалеем, что уделили столь важным темам недостаточное внимание, но это вряд ли можно считать существенным недостатком книги, поскольку детальный разговор планируется впереди, в том числе на страницах блога http://r-analytics.blogspot.com.

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

2. СТАТИСТИЧЕСКИЕ МОДЕЛИ: КРИТЕРИИ И МЕТОДЫ ИХ

ОЦЕНИВАНИЯ

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

1. Разведочный анализ данных (exploratory data analysis), главная цель которого – изучение статистических свойств имеющихся в наличии выборок переменных, наличие выбросов, необходимость (распределение трансформации и др.) и выявление характера взаимосвязей между откликом и предикторами (Mount, Zumel, 2014).

2. Выбор методов построения моделей и спецификация систематической части последних. Например, модели типа "доза-эффект" в биологии на одном и том же исходном материале могут быть построены с использованием самых различных функций: логистической, экспоненциальной, Вейбулла, Гомперца, Михаелиса-Ментен, Брейна-Кузенса и т.д. (Шитиков, 2016).

3. Оценка параметров моделей и их диагностика (Мастицкий, Шитиков, 2014). Диагностика и оценка валидности (model validation) включает в себя ряд стандартных процедур. Так, в случае с классическими регрессионными моделями, подгоняемыми по методу наименьших квадратов, выполняются: а) проверка статистической значимости модели в целом и анализ неопределенности оцененных коэффициентов; б) проверка допущений в отношении остатков модели; в) обнаружение необычных наблюдений и выбросов; г) построение графиков, позволяющих оценить соответствие модели структуре анализируемых данных.

4. Анализ вклада отдельных предикторов и селекция оптимальной комбинации из них (model selection). Оценка качества каждой моделипретендента (model evaluation) по совокупности объективных критериев эффективности, включая тестирование на порции "свежих" данных, не участвовавших в процессе оценивания коэффициентов.

5. Ранжирование нескольких альтернативных моделей и, при необходимости, подстройка их важнейших параметров (model tuning).

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

E (Y x1, x2,…, xm ) = f(, x1, x2,…, xm) +, где остатки отражают ошибку модели, т.е. необъяснимую случайную вариацию наблюдаемых значений зависимой переменной относительно ожидаемого среднего значения.

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

Смещение (bias), или систематическая ошибка модели (сдвиг предсказываемых значений на некоторую трудно объяснимую величину);

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

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

Переусложнение (overfitting) модели, которое «так же вредно, как и ее недоусложнение» (Ивахненко, 1982).

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

James et al., 2013):

E[y0f(x0)]2 = Var[f(x0)] + [Bias(f(x0))]2+Var(), где Bias означает смещение, а Var – дисперсию. Здесь мы предположили, что неизвестная истинная функция f(…) была получена на большом числе обучающих выборок, а отклонения y0 вычислялись по каждой из множества моделей с последующим усреднением результатов.

Из приведенного уравнения следует, что для минимизации ожидаемой ошибки прогноза мы должны подобрать такую модель, для которой одновременно достигаются низкое смещение и низкая дисперсия. Обратите внимание, что дисперсия никогда не может быть ниже некоторого уровня неустранимой ошибки Var() (irreducible error).

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

В качестве простого примера используем данные по электрическому сопротивлению (Ом) мякоти фруктов киви в зависимости от процентного содержания в ней сока. Таблица с этими данными (fruitohms) входит в состав пакета DAAG, который является приложением к книге Maindonald (2010). В рассматриваемом примере у нас есть лишь один предиктор – содержание сока в мякоти фруктов juice, и было бы вполне логичным на первоначальном этапе рассмотреть простую линейную регрессию.

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

2.1):

# Загрузка и визуализация данных в виде диаграммы рассеяния:

library(DAAG) data("fruitohms") library(ggplot2) ggplot(fruitohms, aes(x = juice, y = ohms)) + geom_point() + stat_smooth(method = "lm") + xlab("Содержание сока, %") + ylab("Сопротивление, Ом") Рис. 2.1. График линейной зависимости электрического сопротивления мякоти плодов киви от содержания сока Все критерии оценки качества классических моделей регрессии f так или иначе основаны на анализе остатков (residuals), т.е. разностей между прогнозируемыми f(x[i,]) и наблюдаемыми y[i] значениями, где x – матрица независимых переменных. Базовыми критериями являются сумма квадратов остатков RSS (residual sum of squares), корень из среднеквадратичной ошибки RMSE (root mean square error) и стандартное отклонение остатков RSE (residual standard error). Для унификации обозначений в листингах кода перейдем от предметных обозначений переменных к их формальным эквивалентам (y и x соответственно).

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

x - fruitohms$juice y - fruitohms$ohms/1000 # Строим модель с одним предиктором n - dim(as.matrix(x))[1]; m - dim(as.matrix(x))[2] M_reg - lm(y ~ x); pred - predict(M_reg) RSS - sum((y - pred) * (y - pred)) RMSE - sqrt(RSS/n) RSE - sqrt(RSS/(n - m - 1)) c(RSS, RMSE, RSE) [1] 158.706891 1.113507 1.122309 Мы не можем сказать с определенностью, малы или велики эти ошибки, поскольку все зависит от шкалы измерения y, и поэтому необходимо определиться с набором "эталонов" для сравнения. Тогда, выбрав некоторый подходящий критерий качества, мы сможем оценить, насколько эффективность тестируемой модели по этому критерию отличается от эталона.

Такими эталонными моделями являются (Mount, Zumel, 2014):

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

2. Модель с минимальным уровнем байесовской ошибки (Bayes rate model)

– самая лучшая модель для данных, имеющихся под рукой. Она, как правило, основывается на всех имеющихся переменных (т.е. является максимально "насыщенной" – saturated model) и ее ошибка определяется только набором наблюдений с разными значениями отклика y при одних и тех же x1,..., хn.

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

3. Модель с одной наиболее информативной переменной. Если в процессе селекции более сложные модели по своей эффективности не превосходят модель с одной переменной, то включение дополнительных переменных вряд ли имеет смысл.

Традиционными оценками качества аппроксимации данных являются коэффициент детерминации Rsquared, дисперсионное отношение Фишера F и соответствующее ему р-значение.

Отметим, что F-критерий интерпретируется как мера превышения точности предсказания отклика у построенной модели над нуль-моделью:

Rsquared - 1 - RSS/sum((mean(y) - y)^2) F - (sum((mean(y) - pred)^2)/m)/(RSS/(n - m - 1)) p - pf(q = F, df1 = m, df2 = (n - m - 1), lower.tail = FALSE) c(Rsquared, F, p) [1] 6.387089e-01 2.227492e+02 1.234110e-29 Разумеется, все эти величины можно получить и с помощью стандартной функции summary(), но нам показалось интересным показать всю "кухню" их расчетов.

summary(M_reg)

Coefficients:

Estimate Std. Error t value Pr(|t|) (Intercept) 7.519404 0.233779 32.16 2e-16 *** x -0.089877 0.006022 -14.93 2e-16 ***

--Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.122 on 126 degrees of freedom Multiple R-squared: 0.6387, Adjusted R-squared: 0.6358 F-statistic: 222.7 on 1 and 126 DF, p-value: 2.2e-16 anova(M_reg) Analysis of Variance Table Response: y Df Sum Sq Mean Sq F value Pr(F) x 1 280.57 280.57 222.75 2.2e-16 *** Residuals 126 158.71 1.26 Мы убедились в том, что линейная модель со статистически значимыми коэффициентами по всем формальным критериям вполне адекватна полученным данным. Для проверки условия однородности дисперсии остатков в пакете car имеется функция ncvTest() (от "non-constant variance test" – "тест на непостоянную дисперсию"), которая позволяет проверить нулевую гипотезу о том, что дисперсия остатков никак не связана с предсказанными моделью значениями.

ncvTest(M_reg) Non-constant Variance Score Test Variance formula: ~ fitted.values Chisquare = 4.647123 Df = 1 p = 0.03110564 Те же данные мы можем использовать для построения обобщенной линейной модели (general linear model, GLM), задав в качестве аргумента family функции glm() гауссовый характер распределения данных "gaussian". Вместо минимизации суммы квадратов отклонений эта модель ищет экстремум логарифма функции наибольшего правдоподобия (log maximum likelihood), вид которой зависит от заданного распределения.

В общем случае, значение функции правдоподобия численно равно вероятности того, что модель правильно предсказывает любое предъявленное ей наблюдение из заданной выборки. Логарифм функции правдоподобия LL будет всегда отрицательным, и, поскольку log(1) = 0, то для оптимальной модели, прогнозирующей наименее ошибочное значение отклика, значение LL будет стремиться к 0. Часто для оценки расхождений между наблюдаемыми и прогнозируемыми данными вместо LL используют остаточный девианс1 (deviance), который определен как D = -2(LL - S), где S – правдоподобие "насыщенной модели" с минимально возможным уровнем неустранимой ошибки. Обычно принимают S = 0, поскольку наилучшая модель выполняет, как правило, верное предсказание.

Чем меньше выборочный остаточный девианс D, тем лучше построенная модель. Аналогично можно рассчитать логарифм правдоподобия и девианс для нулевой модели D.null. Тогда адекватность модели определяется соотношением девианса остатков D и нуль-девианса D.null путем вычисления псевдо-коэффициента детерминации PRsquare.

Девиансанализ является обобщением техники дисперсионного анализа и статистическую значимость разности двух значений девианса (D.null – D) можно оценить по критерию 2:

M_glm - glm(y ~ x) lgLik - logLik(M_glm) D.null - M_glm$null.deviance D - M_glm$deviance df - with(M_glm, df.null - df.residual) p - pchisq(D.null - D, df, lower.tail = FALSE) PRsquare = 1 - D/D.null c(lgLik, D, D.null, p, PRsquare) [1] -1.9538e+02 1.5870e+02 4.3927e+02 5.6410e-63 6.3870e-01

Те же величины можно получить с использованием функции summary():

summary(M_glm)

Coefficients:

Estimate Std. Error t value Pr(|t|) (Intercept) 7.519404 0.233779 32.16 2e-16 *** x -0.089877 0.006022 -14.93 2e-16 ***

--Null deviance: 439.28 on 127 degrees of freedom Residual deviance: 158.71 on 126 degrees of freedom AIC: 396.77 with(M_glm, null.deviance - deviance) [1] 280.57 anova(M_glm, glm(x ~ 1), test = "Chisq") Analysis of Deviance Table Df Deviance Resid. Df Resid. Dev P(|Chi|) NULL 127 439.28 x 1 280.57 126 158.71 2.2e-16 *** Нетрудно заметить, что классическая и обобщенная линейные модели в случае нормального распределения дают идентичные результаты и отличаются лишь по характеру используемой терминологии.

Этот термин в русскоязычной литературе пока не устоялся: используются также девиация (Ллойд, Ледерман, 1990) и аномалия (International Statistical Institute) Рассмотрим теперь обоснованность выбора систематической части модели. Поверим утверждениям, что наиболее совершенным инструментом экспертизы нелинейности пока остается глаз человека, и отметим на рис.

2.1 некоторую асимметрию распределения точек относительно линии регрессии:

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

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

2.2):

library(hexbin) ggplot(fruitohms, aes(x=juice, y=ohms)) + geom_hex(binwidth=c(3, 500)) + geom_smooth(color="red", se=F) + xlab("Содержание сока, %") + ylab("Сопротивление, Ом") Рис. 2.2. Кривая сглаживания, описывающая зависимость электрического сопротивления мякоти плодов киви от содержания сока Поскольку функциональная форма истинной модели нам неизвестна, сделаем предположение, что удовлетворительная аппроксимация данных может быть выполнена полиномиальной зависимостью. Отметим, что при увеличении степени m полинома любые критерии эффективности, основанные на остатках (RSE, Rsquared, F) будут монотонно улучшаться, поскольку каждая новая модель будет полнее отражать характер зависимости, имеющий место в обучающей выборке (рис. 2.3).

Рис. 2.3. Аппроксимация данных полиномами разных степеней

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

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

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

Если D – выборочный остаточный девианс, k – число свободных параметров модели, а n – объем обучающей выборки, то можно определить семейство следующих весьма популярных критериев, которые обобщены под названием информационных (Burnham, Anderson, 2002):

классический критерий Акаике AIC = D + 2*k;

байесовский критерий Шварца BIC = D + k*ln(n);

скорректированный критерий Акаике AICс= AIC+2k(k+1)/(n-k-1).

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

k - extractAIC(M_glm)[1]; AIC - extractAIC(M_glm)[2] AIC - AIC(M_glm); AICc - AIC(M_glm) + 2*k*(k+1)/(n-k-1) BIC - BIC(M_glm); BIC - AIC(M_glm, k=log(n)) c(AIC, AICc, BIC) [1] 396.7719 396.8679 405.3280 Оптимальной считается такая модель, которой соответствуют субэкстремальные значения критериев качества (например, минимум AIC).

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

2.4):

# Построение моделей со степенью полинома от 1 до 7 max.poly - 7 # Создание пустой таблицы для хранения значений AIC и BIC, # рассчитанных для всех моделей, и ее заполнение AIC.BIC - data.frame(criterion = c(rep("AIC", max.poly), rep("BIC", max.poly)), value = numeric(max.poly*2), degree = rep(1:max.poly, times = 2)) for(i in 1:max.poly) { AIC.BIC[i, 2] - AIC(lm(y ~ poly(x, i))) AIC.BIC[I + max.poly, 2] - BIC(lm(y ~ poly(x, i))) } # График AIC и BIC для разных степеней полинома qplot(degree, value, data = AIC.BIC, geom = "line", linetype = criterion) + xlab("Степень полинома") + ylab("Значение критерия")

–  –  –

Минимальные значения информационных критериев соответствуют выводу, что наилучшая модель – полином 4-й степени.

2.2. Использование алгоритмов ресэмплинга для тестирования моделей и оптимизации их параметров Фундаментальным для статистики является рассуждение о соотношении свойств случайных выборок и генеральной совокупности, из которых они были извлечены. Любой эксперимент в принципе ограничен некоторым (чаще всего, небольшим) множеством наблюдений, причем никакие сверхинтеллектуальные методы обработки не являются панацеей от влияния неучтенных факторов или систематических погрешностей. Один из возможных путей надежного оценивания свойств данных заключается в использовании компьютерно-интенсивных (computer-intensive) технологий, объединенных общим термином "численный ресэмплинг" (англ. resampling) или, как их иногда называют в русскоязычной литературе, "методов генерации повторных выборок" (Эфрон, 1988; Edgington, 1995; Davison, Hinkley, 2006). Ключевая особенность этой технологии заключается в том, что повторные выборки при классическом сценарии извлекаются из генеральной совокупности, а псевдоповторные выборки при выполнении ресэмплинга – из самой эмпирической выборки (хорошо известна фраза «The population is to the sample as the sample is to the bootstrap samples» из работы Fox, 2002).

Ресэмплинг объединяет четыре разных подхода к обработке данных, отличающихся по алгоритму, но близких по сути: перекрестная проверка бутстреп рандомизация, или (cross-validation, CV), (bootstrap), перестановочный тест (permutation), и метод "складного ножа" (jackknife).

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

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

При отсутствии дополнительных блоков данных, специально предназначенных для тестирования, вполне естественно выглядит идея провести случайное разбиение исходной выборки на обучающее (train sample) и проверочное (test sample) подмножества, после чего оценить параметры модели только на обучающей выборке, тогда как оценку погрешности прогноза отклика y i осуществить для значений из проверочной совокупности. Если подобные шаги выполняются многократно и каждое из наблюдений используется поочередно то для обучения, то для контроля, то такая процедура эмпирического оценивания моделей, построенных по прецедентам, называется перекрестной проверкой, или "кросс-проверкой".

Стандартной считается методика rk-кратной кросс-проверки (rk-fold cross-validation), когда исходная выборка случайным образом r раз разбивается на k блоков (folds) равной (или почти равной) длины. При реализации каждой повторности r (replication) один из блоков по очереди становится проверочной совокупностью, а все остальные блоки – одной большой обучающей выборкой (рис.

2.5):

–  –  –

createDataPartition(y, p = 0.5) createFolds (y, k = 10) createMultiFolds (y, k = 10, times = r) Частным случаем является "скользящий контроль", или перекрестная проверка с последовательным исключением одного наблюдения (leave-oneout CV, LOOCV), т.е. k = n. При этом строится n моделей по (n – 1) выборочным значениям, а исключенное наблюдение каждый раз используется для расчета ошибки прогноза. В. Н. Вапник (1984) теоретически обосновал применение скользящего контроля и показал, что если исходные выборки независимы, то средняя ошибка перекрестной проверки даёт несмещённую оценку ошибки модели. Это выгодно отличает её от средней ошибки на обучающей выборке, которая при переусложнении модели может оказаться оптимистично заниженной.

В разделе 2.1 мы рассматривали пример выбора оптимального числа параметров полиномиальной регрессии с использованием информационных критериев.

Попробуем достичь той же цели, выполнив перекрестную проверку моделей в режиме "leave-one-out" с использованием самостоятельно оформленной функции:

# Функция скользящего контроля для модели y~poly(x, degree) crossvalidate - function(x, y, degree) { preds - numeric(length(x)) for(i in 1:length(x)) { x.in - x[-i]; x.out - x[i] y.in - y[-i]; y.out - x[i] m - lm(y.in ~ poly(x.in, degree=degree) ) new - data.frame(x.in = seq(-3, 3, by=0.1)) preds[i]- predict(m, newdata=data.frame(x.in=x.out)) } # Тестовая статистика - сумма квадратов отклонений:

return(sum((y-preds)^2)) } # Заполнение таблицы результатами кросс-проверки # и сохранение квадрата ошибки в таблице "a" a - data.frame(cross=numeric(max.poly)) for(i in 1:max.poly) { a[i,1] - crossvalidate(x, y, degree=i) } # График суммы квадратов ошибки при кросс-проверке qplot(1:max.poly,cross, data=a, geom=c("line"))+ xlab("Степень полинома ") + ylab("Квадратичная ошибка") Рис. 2.6. Поиск степени функции полиномиальной регрессии с использованием перекрестной проверки Как видно из графика на рис. 2.6, по мере усложнения модели от линейной регрессии с одним предиктором до полинома 7-й степени ошибка предсказаний на проверочном множестве (test error) сначала снижается, а затем после достижения некоторого минимума начинает возрастать. Такая Uобразная форма кривой ошибки является универсальной и справедлива практически для любых алгоритмов и методов создания предсказательных моделей. Повторный рост ошибки на проверочных данных является сигналом того, что модель стала переобученной. Результаты, представленные на рис.

2.4 и 2.6, говорят о том, что в нашем конкретном случае оба метода одинаково выбирают в качестве наилучшей модели полином 4-й степени.

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

(M_poly - glm(y ~ poly(x, 4)))

Coefficients:

(Intercept) poly(x,4)1 poly(x,4)2 poly(x,4)3 poly(x,4)4 4.360 -16.750 4.751 3.946 -3.371 Degrees of Freedom: 127 Total (i.e. Null); 123 Residual Null Deviance: 439.3 Residual Deviance: 109.2 AIC: 354.9 anova(M_glm, M_poly, test ="Chisq") Analysis of Deviance Table Model 1: y ~ x Model 2: y ~ poly(x, 4) Resid. Df Resid. Dev Df Deviance P(|Chi|) 1 126 158.71 2 123 109.21 3 49.501 4.743e-12 *** Имеет место статистически значимое снижение ошибки модели.

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

CVlm(df, form.lm, m=3) из пакета DAAG, где df – исходная таблица с данными, form.lm – формула линейной модели, m – число блоков (сюда подставляется значение k);

cv.glm(df, glmfit, K = n) из пакета boot, где df – исходная таблица с данными, glmfit – объект обобщенной линейной модели типа glm, K – число блоков (т.е. по умолчанию выполняется скользящий контроль);

cvLm(lmfit, folds = cvFolds(n, K, R), cost) из пакета cvTools, где lmfit – объект обобщенной линейной модели lm, K – число блоков, R – число повторностей, cost – наименование одного из пяти разновидностей используемых критериев качества;

train(form, df, method, trainControl, …) из пакета caret, где df – исходная таблица с данными, form, method – формула и метод построения модели, trainControl – объект, в котором прописываются все необходимые параметры перекрестной проверки.

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

Конкретно, если мы имеем исходную выборку из n членов x1, x2, …, xn-1, xn, то с помощью датчика случайных чисел, равномерно распределенных на интервале [1, n], можем "вытянуть" произвольный элемент xk, который снова вернем в исходную выборку для возможного повторного извлечения. Такая процедура повторяется n раз и образуется бутстреп-выборка, в которой одни элементы могут повторяться два или более раз, тогда как другие элементы – отсутствовать. Например, при n = 6 одна из таких бутстреп-комбинаций имеет вид x4, x2, x2, x1, x4, x5.

В R бутстреп легко реализуется с помощью функции sample(..., replace=T), генерирующей любые случайные выборки с "возвращением".

Вероятность того, что конкретное наблюдение не войдет в бутстреп-выборку размера n, равна (1 - 1/n)n и стремится к 1/е = 0.368 при n:

Nboot = 1000; n=100 mean(replicate(Nboot, length(unique(sample(n, replace=TR))))) [1] 63.544 Таким способом можно сформировать любое, сколь угодно большое число бутстреп-выборок (обычно 500-1000), каждая из которых содержит около 2/3 уникальных значений эмпирической совокупности.

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

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

Осуществим оценку квантильных интервалов байесовского критерия BIC для каждой из модели-претендента:

set.seed(123) Nboot = 1000 BootMat - sapply(1:max.poly, function(k) replicate(Nboot, { ind - sample(n, replace=T) BIC(glm(y[ind]~poly(x[ind],k))) }) ) apply(BootMat, 2, mean) [1] 401.079 385.967 374.095 366.126 366.585 370.824 373.144 Искомые интервалы покажем на графике (рис.

2.7), подключив функцию stat_summary() из пакета ggplot2:

library(reshape) # для функции melt() BootMat.df - data.frame(melt(BootMat)) ggplot(data = BootMat.df, aes(X2, value)) + geom_jitter(position = position_jitter(width = 0.1), alpha = 0.2) + # Добавляем средние значения в виде точек красного цвета:

stat_summary(fun.y = mean, geom = "point", color = "red", size = 5) + # Добавляем отрезки, символизирующие 0.25 и 0.75 квантили:

stat_summary(fun.y = mean, fun.ymin = function(x){quantile(x, p = 0.25)}, fun.ymax = function(x){quantile(x, p = 0.75)}, geom = "errorbar", color = "magenta", width = 0.5, size =1.5) + xlab("Степень полинома") + ylab("Информационный критерий Байеса") Рис. 2.7. Доверительные интервалы байесовского информационного критерия для разных степеней полинома Для оценки 95%-х доверительных интервалов методом процентилей достаточно установить значения р = с(0.025, 0.975) (мы это не сделали исключительно из эстетических соображений). Из результатов бутстрепа можно сделать содержательные выводы, например: если доверительные интервалы смежных степеней полинома будут пересекаться, то уменьшение BIC статистически незначимо и вполне можно ограничиться более простой моделью. Подробнее о различных возможностях ресэмплинга рассказано в книге (Шитиков, Розенберг, 2014).

2.3. Модели для предсказания класса объектов Классификация – наиболее часто встречающаяся задача машинного обучения, и заключается в построении моделей, выполняющих отнесение интересующего нас объекта к одному из нескольких известных классов.

Существуют сотни методов классификации (см. Fernandez-Delgado et al., 2014), которые можно использовать для предсказания значения отклика с двумя и более классами. Возникает вопрос: отвечает ли такое множество потребностям реально решаемых задач?

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

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

Если же основной задачей является достижение высокой общей точности предсказаний (overall accuracy) значения целевого признака y для объекта a, то представление модели в явном виде не требуется. Изучаемый процесс, который часто имеет объективно сложный характер, представляется в виде "черного ящика", а решающие процедуры могут иметь большое (до десятков тысяч) или неопределенное число трудно интерпретируемых параметров. Эффективными методами прогнозирования классов являются случайные леса, бустинг, бэггинг, искусственные нейронные сети, машины опорных векторов, групповой учет аргументов МГУА и др.

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

Вероятностный подход заключается в предположении о существовании некоего случайного процесса, который порождает значения целевых переменных, подчиняющиеся вполне определенному, но неизвестному нам распределению вероятностей. Примером модели вероятностного характера является байесовский классификатор, формирующий решающее правило по принципу апостериорного максимума. Модели логического типа по своей природе наиболее алгоритмичны, поскольку легко выражаются на языке правил, понятных человеку, таких как: if условие = 1 then Y = категория класса.

Примером таких моделей являются деревья классификации.

Некоторые авторы (Mount, Zumel, 2014, р. 91) подчеркивают различие терминов "предсказание" (prediction) и "прогнозирование" (forecasting).

Предсказание лишь озвучивает результат (например, «Завтра будет дождь»), а при прогнозировании итог связывается с вероятностью события («Завтра с вероятностью 80% будет дождь»). Мы считаем, что на практике трудно провести между этими терминами четкую границу. К тому же, часто эта разница в совершенно не принципиальна – главное понимать контекст задачи.

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

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

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

Группы точек "патология/норма" в заданном пространстве предикторов, как правило, статистически неразделимы: например, повышение температуры тела до 37.5о часто свидетельствует о заболевании, хотя не всегда болезнь может сопровождаться высокой температурой. Поэтому при тестировании модели вероятны ошибочные ситуации, такие как пропуск положительного (патологического) заключения FN или его "гипердиагностика" FP, т.е. отнесение нормального состояния к патологическому.

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

Истинное Результаты теста состояние тестПредсказана патология (1) Предсказана норма (0) объектов Истинно положительные Ложноотрицательные Патология (1) TP (True positives) FN (False negatives) Ложноположительные Истинно отрицательные Норма (0) FP (False positives) TN (True negatives) В этих обозначениях объективная ценность рассматриваемого бинарного классификатора определяется следующими показателями:

Чувствительность (sensitivity) SE = ErrII= TP / (TP + FN), определяющая насколько хорош тест для выявления патологических экземпляров;

Специфичность (specificity) SP = ErrI = FP / (FP + TN), показывающая эффективность теста для правильной диагностики отклонений от нормального состояния;

Точность (accuracy) AC = (TP + TN) / (TP + FP + FN + TN), определяющая общую вероятность теста давать правильные результаты.

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

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

Рассмотрим популярный пример выделения спама (spam – от слияния двух слов spiced и ham, – или пряная ветчина, как образец некачественного пищевого продукта) в электронных письмах в зависимости от встречаемости тех или иных слов (всего 58 частотных показателей). Выборка по спаму представлена в обширной коллекции наборов данных Центра машинного обучения и интеллектуальных систем Калифорнийского университета (UCI Machine Learning Repository – http://archive.ics.uci.edu/ml/) и после некоторой предварительной обработки используется для иллюстрации в книге Mount, Zumel (2014).

Скачаем этот файл с сайта ее авторов и разделим исходные данные в соотношении 10:1 на обучающую и проверочную выборки:

# spamD - read.table( #'https://raw.github.com/WinVector/zmPDSwR/master/Spambase/spamD.tsv', # header=T, sep='\t') spamD - read.table('spamD.tsv', header=T, sep='\t') dim(spamD) [1] 4601 59 spamTrain - subset(spamD, spamD$rgroup=10) spamTest - subset(spamD, spamD$rgroup10) с(nrow(spamTrain), nrow(spamTest)) [1] 4143 458 # Составляем список переменных и объект типа формула spamVars - setdiff(colnames(spamD), list('rgroup','spam')) spamFormula - as.formula(paste('spam=="spam"', paste(spamVars, collapse=' + '), sep=' ~ ')) spamModel - glm(spamFormula, family=binomial(link='logit'), data=spamTrain) # Добавляем столбец со значениями вероятности спама spamTrain$pred - predict(spamModel, newdata=spamTrain, type='response') spamTest$pred - predict(spamModel,newdata=spamTest, type='response') Компоненты матрицы неточностей и перечисленные показатели легко получить с использованием обычной функции table(...)– например, так:

# На обучающей выборке (cM.train - table(Факт = spamTrain$spam, Прогноз = spamTrain$pred 0.5)) Прогноз Факт FALSE TRUE non-spam 2396 114 spam 178 1455 # На проверочной выборке:

(cM - table(Факт = spamTest$spam, Прогноз = spamTest$pred 0.5)) Прогноз Факт FALSE TRUE non-spam 264 14 spam 22 158 c(Точность - (cM[1, 1] + cM[2, 2])/sum(cM), Чувствительность - cM[1, 1]/(cM[1, 1] + cM[2, 1]), Специфичность - cM[2, 2]/(cM[2, 2] + cM[1, 2])) [1] 0.9213974 0.9230769 0.9186047 Иногда предпочтительнее использовать функцию confusionMatrix(y,

pred) из пакета caret:

library(caret) pred - ifelse(spamTest$pred 0.5, "spam", "non-spam") confusionMatrix(spamTest$spam, pred) Confusion Matrix and Statistics Reference Prediction non-spam spam non-spam 264 14 spam 22 158 Accuracy : 0.9214 95% CI : (0.8928, 0.9443) No Information Rate : 0.6245 P-Value [Acc NIR] : 2e-16 Kappa : 0.834 Mcnemar's Test P-Value : 0.2433 Sensitivity : 0.9231 Specificity : 0.9186 Pos Pred Value : 0.9496 Neg Pred Value : 0.8778 Prevalence : 0.6245 Detection Rate : 0.5764 Detection Prevalence : 0.6070 'Positive' Class : non-spam Выбрать другой класс в качестве положительного исхода можно, задав аргумент positive = "spam". Функция предоставляет пользователю такие статистики, как доверительные интервалы и р-значение для точности, результаты теста 2 по Мак-Немару, вероятностный индекс (каппа) Дж.

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

прогностическая ценность (prevalence) PV = (TP + FN)/(TP + FP + FN + TN);

положительная прогностическая ценность (вероятность фактической патологии при положительном диагнозе) PPV = SE * PV/(SE*PV + (1- SP)*(1 - PV));

отрицательная прогностическая ценность (вероятность отсутствия патологии при негативном результате теста) NPV = SP*(1 - PV)/(PV*(1 - SE) + SP*(1 - PV));

частота выявления (detection rate) DR = TP/( TP + FP + FN + TN);

частота распространения (detection prevalence) DP = (TP + FP)/(TP + FP + FN + TN);

сбалансированная точность (balanced accuracy) BAC = (SE + SP)^2.

Эффективность классификатора может также оцениваться с использованием информационных критериев – энтропии E = - рi log2 pi, где pi – вероятности каждого возможного исхода, и условной энтропии (conditional entropy).

Эти меры могут быть рассчитаны с использованием функций:

entropy - function(x) { xpos - x[x0] scaled - xpos/sum(xpos) ; sum(-scaled*log(scaled,2)) } print(entropy(table(spamTest$spam))) [1] 0.9667165 conditionalEntropy - function(t) { (sum(t[,1])*entropy(t[,1]) + sum(t[,2])*entropy(t[,2]))/sum(t) } print(conditionalEntropy(cM)) [1] 0.3971897 Исходная энтропия Е = entropy(table(y)) определяет среднее количество информации, измеряемой в битах, которую мы приобретаем, если извлечь из выборки очередной экземпляр того или иного класса. Условная энтропия conditionalEntropy(table(y, pred)) показывает, насколько эта мера информации уменьшается из-за ошибок предсказания для различных категорий.

Общепринятым графоаналитическим методом оценки качества теста и интерпретации перечисленных показателей является ROC-анализ (от "Receiver Operator Characteristic" – функциональная характеристика приемника), название которого взято из методологии оценки качества сигнала при радиолокации.

ROC-кривая получается следующим образом (Goddard, Hinberg, 1989).

Пусть мы имеем выборку значений независимого количественного показателя, который варьирует от xmin до xmax, и сопряженного с ним бинарного отклика (1 – патология, 0 – норма). Любое произвольное значение х на этом диапазоне может считаться классификационным порогом, или точкой отсечения (cutt-off value), делящим вектор y на два подмножества, и для этого разбиения можно рассчитать значения чувствительности SE и специфичности SP. Если выполнить сканирование всех возможных значений xmax х xmin, то можно построить график зависимости, где по оси Y откладывается SE, а по оси X – (1 - SP). Реализация этой процедуры в R может привести к ступенчатой или сглаженной кривой следующего вида (рис.

2.8):

# ROC кривая library(pROC) m_ROC.roc - roc(spamTest$spam, spamTest$pred) plot(m_ROC.roc, grid.col=c("green", "red"), grid=c(0.1, 0.2), print.auc=TRUE, print.thres=TRUE) plot(smooth(m_ROC.roc), col="blue", add = T, print.auc=F).

Рис. 2.8. ROC-кривая для оценки вероятности спама В случае идеального классификатора ROC-кривая проходит вблизи верхнего левого угла, где доля истинно-положительных случаев равна 1, а доля ложно-положительных примеров равна нулю. Поэтому, чем ближе кривая к верхнему левому углу, тем выше предсказательная способность модели. Наоборот, главная диагональная линия соответствует классификатору, который классовую "бесполезному" "угадывает" принадлежность случайным образом. Следовательно, близость ROC-кривой к диагонали говорит о низкой эффективности построенной модели.

Для нахождения оптимального порога, соответствующего наиболее безошибочному классификатору, через крайнюю точку ROC-кривой проводят линию максимальной точности, параллельную главной диагонали. На приведенном графике такая точка, соответствующая значению х = 0.382, имеет наилучшую комбинацию значений чувствительности SE = 0.932 и специфичности SP = 0.922. Обратите внимание, что в качестве значений х фигурирует оценка вероятности отнесения к спаму и, видимо, мы совершенно напрасно принимали ранее в качестве порога величину 0.5.

Полезным показателем является численная оценка площади под ROCкривыми AUC (Area Under Curve). Практически она изменяется от 0.5 ("бесполезный" классификатор) до 1.0 ("идеальная" модель). Показатель AUC предназначен исключительно для сравнительного анализа нескольких моделей, поэтому связывать его величину с прогностической силой можно только с большими допущениями.

В нашем примере в качестве классификатора писем со спамом мы использовали модель логистической регрессии, полагая, что бинарный отклик имеет биномиальное распределение. Напомним, что в случае обобщенных линейных моделей GLM вместо минимизации суммы квадратов отклонений ищется экстремум логарифма функции максимального правдоподобия (log maximum likelihood), вид которой зависит от характера распределения данных.

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

(LL - logLik(spamModel)) 'log Lik.' -807.0323 (df=58) Как и в случае гауссова распределения (см. раздел 2.1), оценка адекватности биномиальной модели осуществляется с использованием девианса D = -2(LL - S), где S = 0 – правдоподобие "насыщенной модели" с минимальным уровнем байесовской ошибки. Исходя из априорной вероятности одного из классов, можно рассчитать логарифм правдоподобия и девианс для нулевой модели D.null. Эффективность классификатора определяется соотношением девианса остатков D и нуль-девианса D.null, что соответствует псевдо-коэффициенту детерминации Rsquared модели.

Статистическую значимость разности девиансов (D.null – D) можно оценить по критерию 2:

df - with(spamModel, df.null - df.residual) c(D.null - spamModel$null.deviance, D - spamModel$deviance, Rsquared = 1-D/D.null, pchisq(D.null - D, df, lower.tail = FALSE)) [1] 5556.36 1614.065 0.7095104 0.000000 Мы получили модель, вполне адекватную по отношению к имеющимся данным.

Разумеется, все эти вычисления могут быть выполнены с использованием базовых функций summary() и anova():

summary(spamModel) Null_Model - glm(spam ~ 1,family=binomial(link='logit'), data=spamTrain) anova(spamModel, Null_Model, test ="Chisq") Мы не станем приводить здесь длинные протоколы с результатами этих процедур, включающие статистический анализ 58 коэффициентов модели.

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

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

Следующая команда R обеспечивает построение такого графика:

ggplot (data=spamTest) + geom_density (aes (x=pred, color=spam, linetype=spam)) Рис. 2.9. Кривые плотности апостериорной вероятности принадлежности объектов к двум классам: электронным письмам с наличием спама и без него

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

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

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

а) возможностью наглядного представления (визуализация и ординация);

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

Принцип ординации наблюдений (нем. Ordnung) заключается в использовании различных методов оптимального целенаправленного проецирования (projecting pursuit) облака точек из многомерного пространства в пространство малой размерности (с 2 или 3 осями координат).

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

Y в виде совокупности p латентных переменных F:

Y1, Y2, …, Ym F1, F2, …, Fp, которые и являются осями ординационной диаграммы (двумерной при р = 2 или трехмерной в случае с тремя такими латентными переменными).

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

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

Как правило, первая ось F1 проводится через центр сгущения значений yij и совпадает с направлением наибольшей по длине оси эллипсоида рассеяния. Вторая ось F2, также проходит через центр распределения, но проводится перпендикулярно к первой и совпадает по направлению со второй из главных полуосей эллипсоида рассеяния. Эта операция обеспечивает формирование "картинки" данных с минимально возможными искажениями, а новые обобщающие переменные F1 и F2 становятся ортогональными (т.е.

взаимно некоррелированными).

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

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

1. Нахождение собственных значений и собственных векторов осуществляется в ходе математической операции линейного преобразования квадратной симметричной матрицы C размерностью m. Это может быть любая матрица дистанций, но чаще всего используется дисперсионноковариационная матрица C = Y'Y/(n - 1) стандартизованных исходных данных.

2. Собственными значениями квадратной матрицы C называются такие значения k, при которых система m уравнений вида (C - k I)uk = 0 имеет нетривиальное решение. Здесь uk – собственные векторы матрицы C, соответствующие k, k = 1, 2, …, m, I – единичная матрица.

3. Матрица из k собственных векторов Uk представляет собой веса для пересчета из исходного в редуцированное информационное пространство, т.е.

матрицы переменных Y и F связаны соотношением Fnk = Ynm Umk.

4. Каждое собственное значение k соответствует величине дисперсии, объясняемой на k-м уровне, т.е. сумма всех собственных значений будет равняться сумме дисперсий всех исходных переменных.

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

Значения 1 и 2 определяют меру значимости осей F1 и F2 ординационной диаграммы, т.е. концентрацию исходных точек вдоль каждой оси и, в итоге, ее выраженность.

6. Собственные значения k находят в ходе итерационной процедуры, и точность их вычисления зависит от степени обусловленности исходной матрицы C.

Классическим методом снижения размерности данных является анализ главных компонент (PCA, Principal Components Analysis), который широко используется в различных областях науки и техники и детально описан в многочисленных руководствах (ter Braak, 1983; Айвазян и др., 1989;

Джонгман и др., 1999; Legendre, Legendre, 2012).

Снижение размерности исходного пространства методом PCA можно представить как последовательный, итеративный процесс, который можно оборвать на любом шаге u. Вследствие ортогональности системы координат главных компонент (т.е. фактически, их взаимной независимости) нет необходимости перестраивать матрицы счетов (scores) F и нагрузок (loadings) U при изменении числа компонент: последующие столбцы или строки просто прибавляются или отбрасываются.

Используем в качестве примера идентификацию типа оконных стекол (Glass Identification) по выборке, представленный в коллекции баз данных UCI Machine Learning Repository. Существует два основных способа производства листового стекла: горизонтальный на расплаве олова (флешстекло) и по принципу вертикального вытягивания. Термически полированное флеш-стекло отличается идеально глянцевой поверхностью, высокой светопропускающей способностью и великолепными оптическими свойствами, исключающими искажение изображения.

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

Исходные данные представлены в файле Glass.txt, и можно предварительно рассмотреть описательные статистики отдельных переменных:

DGlass - read.table(file = "Glass.txt", sep = ",", header = TRUE, row.names = 1) print(t(apply(DGlass[, -10], 2, function (x) { c(Минимум = min(x), Максимум = max(x), Среднее = mean(x), Отклонение = sd(x), Корреляция = cor(x, DGlass$Class)) # признаков с Class })), 3) Минимум Максимум Среднее Отклонение Корреляция RI 1.51 1.53 1.5186 0.00305 -0.0601 Na 10.73 14.86 13.2017 0.58830 0.0186 Mg 0.00 4.49 3.2949 0.88778 -0.1462 Al 0.29 2.12 1.2817 0.32374 0.1998 Si 69.81 74.45 72.5869 0.64117 -0.0785 K 0.00 1.10 0.4775 0.21873 0.0386 Ca 7.08 16.19 8.9247 1.37263 0.0446 Ba 0.00 3.15 0.0298 0.25353 0.0312 Fe 0.00 0.37 0.0676 0.09951 0.0532 В среде R расчет главных компонент может быть осуществлен с использованием функций princomp() или prcomp(), но мы будем ориентироваться на многообразие возможностей пакета vegan() и его функцию rda().

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

library(vegan) Y - as.data.frame(DGlass[,1:9]) mod.pca - rda(Y ~ 1) summary(mod.pca)

Importance of components:

PC1 PC2 PC3 PC4 Eigenvalue 2.6056 0.5828 0.21129 0.18895 Proportion Explained 0.7126 0.1594 0.05779 0.05168 Cumulative Proportion 0.7126 0.8720 0.92984 0.98152 Loads PC1 PC2 PC3 PC4 RI -0.006663 0.003040 -0.001587 -0.0007297 Na 0.461583 1.180084 0.833131 0.0499341 Mg 2.037591 0.659746 -0.616211 -0.5223737 Al 0.240546 -0.337483 -0.107770 0.6588300 Si 0.691702 -1.339222 0.458279 -0.4986449 K 0.265988 -0.302647 -0.111021 0.2512089 Ca -3.512043 0.228783 -0.158561 -0.3491594 Ba -0.190199 -0.009844 -0.267021 0.3372595 Fe -0.039086 -0.017137 -0.046526 0.0299297 Рассмотрим, какую конфигурацию имеет распределение наблюдений в пространстве первых двух главных компонент PC1 и PC2. Наибольший интерес для нас представляет взаимное расположение сгущений точек, принадлежащих стеклам двух типов: Class = {1, 3} – оконное и автомобильное флеш-стекло, и Class = 2 – тянутое стекло.

Выполним построение ординационной диаграммы с использованием пакета ggplot2:

F - as.factor(ifelse(DGlass$Class == 2, 2, 1)) pca.scores - as.data.frame(summary(mod.pca)$sites[,1:2]) pca.scores - cbind(pca.scores, F) # Составляем таблицу для "каркаса" точек на графике l - lapply(unique(pca.scores$F), function(c) { f - subset(pca.scores, F == c); f[chull(f), ]}) hull - do.call(rbind, l) # Включаем в названия осей доли объясненной дисперсии axX - paste("PC1 (", as.integer(100*mod.pca$CA$eig[1]/sum(mod.pca$CA$eig)),"%)") axY - paste("PC2 (", as.integer(100*mod.pca$CA$eig[2]/sum(mod.pca$CA$eig)),"%)") # Выводим ординационную диаграмму ggplot() + geom_polygon(data=hull,aes(x=PC1,y=PC2, fill=F), alpha=0.4, linetype=0) + geom_point(data=pca.scores,aes(x=PC1,y=PC2,shape=F, colour=F),size=3) + scale_colour_manual( values = c('purple', 'blue'))+ xlab(axX) + ylab(axY) + coord_equal() + theme_bw()

Рис. 2.10. Ординационная диаграмма, построенная методом PCA

Заметим, что для выделения наблюдений каждого класса на диаграммах обычно применяются четыре вспомогательных графических объекта: а) "каркас" (hull) или контур, проводимый через крайние точки, как мы это выполнили на рис. 2.10; б) эллипс, ограничивающий, например, 95% наблюдений; в) точка центроида или "центра тяжести" облака точек;

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

По существу расчетов можно отметить следующее:

первая главная компонента объясняет 71% совокупной дисперсии в данных, а F1 и F2 совместно – 87%, что является вполне обнадеживающим результатом;

ось РС1 в значительной мере определяется содержанием магния и кальция, а ось РС2 – концентрациями натрия и кремния;

диаграмма на рис. 2.10 показывает, что химический состав флеш-стекла в целом не отличается от тянутого, но для него не характерно повышенное содержание кальция (точки левее РС1 = -0.5).

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

2.5. Многомерный статистический анализ данных Взаимосвязь между двумя комплексами переменных может быть установлена в ходе канонического анализа (от греческого ). В математике, канонической формой называется самая простая и универсальная форма, к который определенные функции, отношения, или выражения могут быть приведены без потери общности. Наиболее распространенными формами канонических моделей, основанных на приведении к собственным значениям и собственным векторам, являются анализ избыточности RDA (redundancy analysis) и линейный дискриминантный анализ LDA (linear discriminant analysis).

Анализ избыточности (Джонгман и др., 1999; Legendre, Legendre, 2012) является непосредственным расширением идей множественной регрессии на моделирование данных с многомерным откликом.

Анализ асимметричен:

если Y(nm) – таблица зависимых и X(nq) – таблица объясняющих переменных, то чаще всего m q. RDA выполняется в три этапа:

1. Оцениваются коэффициенты В модели множественной линейной регрессии Y на X и формируется матрица прогнозируемых значений Y :

-1 Y = X B = X[X`X] X`Y;

2. Вычисляется матрица канонических коэффициентов S = BUс размерностью mq, которая рассчитывается на основе собственных векторов Uс, взвешивающих влияние объясняющих переменных (constrained eigenvalues).

3. Рассчитывается ковариационная матрица Cy|x на основе значений Y, выполняется анализ главных компонент:

(Cy|x - k I)uk = 0, где Cy|x = CYXCXX-1CTYX, и обычным путем вычисляется матрица нагрузок на основе k собственных векторов, не связанных с объясняющими переменными Х (unconstrained eigenvalues).

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

F - as.factor(ifelse(DGlass$Class == 2, 2, 1)) Y - as.data.frame(DGlass[, 1:9]) mod.rda - rda(Y ~ F) summary(mod.rda)

Importance of components:

RDA1 PC1 PC2 PC3 PC4 Eigenvalue 0.11833 2.5441 0.5584 0.21014 0.15944 Proportion Explained 0.03236 0.6958 0.1527 0.05747 0.04361 Cumulative Proportion 0.03236 0.7282 0.8809 0.93837 0.98197 Species scores RDA1 PC1 PC2 PC3 PC4 RI 6.164e-05 0.006766 -0.0027906 0.001443 0.001003 Na -2.177e-01 -0.423708 -1.1896878 -0.802298 -0.127952 Mg -7.082e-01 -1.943017 -0.5819804 0.515740 0.600069 Al 3.060e-01 -0.296512 0.2471444 0.226655 -0.564066 Si 2.683e-02 -0.715098 1.3444594 -0.549095 0.330432 K 1.054e-01 -0.289194 0.2676415 0.153595 -0.215380 Ca 3.605e-01 3.501644 -0.1493831 0.091594 0.335991 Ba 4.946e-02 0.183784 -0.0005734 0.318738 -0.313081 Fe 2.934e-02 0.034801 0.0124987 0.051718 -0.014555 Centroids for factor constraints RDA1 PC1 PC2 PC3 PC4 cen1 -0.3612 0 0 0 0 0 cen2 0.4134 0 0 0 0 0 После включения в анализ объясняющего фактора «Тип производства стекла» появляется связанная с ним ось новой переменной RDA1, а к матрице нагрузок на главные компоненты добавляется столбец соответствующих канонических коэффициентов. Собственное значение, определяющее RDA1, невелико и составляет только 3.2% от общей вариации данных.

Сформируем ординационную диаграмму в координатах {RDA1, PC1}, на которой обозначим координаты центроидов и выделим 95% доверительные области в виде эллисов (рис.

2.11):

library(ggplot2) rda.scores - as.data.frame(scores(mod.rda, display="sites", scales = 3)) rda.scores - cbind(F,rda.scores) centroids - aggregate(cbind(RDA1, PC1)~F,rda.scores, mean) # функция для std.err f - function(z)sd(z)/sqrt(length(z)) se - aggregate(cbind(RDA1,PC1) ~ F, rda.scores, f) names(se) - c("F", "RDA1.se", "PC1.se") # объединяем координаты центроидов и стандартные ошибки centroids - merge(centroids, se, by = "F") # Формируем диаграмму ggplot() + geom_point(data = rda.scores, aes(x=RDA1, y=PC1, shape=F, colour=F), size=2) + xlim(c(-2, 5)) + ylim(c(-0.75, 1)) + scale_colour_manual(values = c('purple', 'blue')) + stat_ellipse(data=rda.scores, aes(x=RDA1, y=PC1, fill=F), geom="polygon", level=0.8, alpha=0.2) + geom_point(data=centroids, aes(x=RDA1, y=PC1, colour=F), shape=1, size=4) + theme_bw() Рис. 2.11. Ординационная диаграмма, построенная методом RDA (не показаны 7 точек за пределами верхнего правого угла графика) Несмотря на небольшую долю дисперсии, объясняемой фактором способа производства стекла F, координаты центроидов, определяющих центры тяжести сравниваемых групп, стали располагаться на некотором расстоянии друг от друга на рис. 2.11, что позволяет мысленно провести разделяющую линию между ними.

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

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

Для формирования модели многомерного дисперсионного анализа можно воспользоваться функциями aov(), lm() либо manova(), которые приводят к идентичным результатам.

Получить традиционную дисперсионную таблицу можно с использованием функции anova():

mod.an - manova(as.matrix(Y) ~ F) anova(mod.an) Analysis of Variance Table Df Pillai approx F num Df den Df Pr(F) (Intercept) 1 1.00000 65907164 9 153 2.2e-16 *** F 1 0.27971 7 9 153 6.106e-08 *** Residuals 161 Теснота связи между химическим составом стекла Y и способом его изготовления F была оценена с использованием статистики Пиллая (Pillai's trace).

Поскольку мы имеем лишь две группы наблюдений, то статистическая проверка гипотезы о равенстве средних двух векторов Ho: X1 = X 2 может быть осуществлена также с использованием статистики Хотеллинга, которая является многомерным аналогом t-критерия Стьюдента:

library(Hotelling) split.data - split(Y, F) (summ.hot - hotelling.test(split.data[[1]], split.data[[2]], perm = T, B = 1000)) Test stat: 6.6015 Permuation P-value: 0 Number of permutations : 1000 В выполненном перестановочном тесте ни одна из 1000 статистик Хотеллинга, полученных при случайном перемешивании данных, не превысила тестируемое значение, что свидетельствует о высокой статистической значимости разделяющего фактора.

2.6. Методы кластеризации В разделах 2.4-2.5 мы рассмотрели две процедуры: необъясненная (непрямая) ординация наблюдений, при которой конфигурация точек определяется случайной вариацией данных и не связывается с какими-либо внешними причинами, и прямая ординация с использованием объясняющих переменных (в рассмотренном примере – способ изготовления стекла F).

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

Под кластеризацией (от англ. cluster – гроздь, скопление) понимается задача разбиения всей исходной совокупности элементов на отдельные группы однородных объектов, сходных между собой, но имеющих отчетливые отличия этих групп друг от друга. Пусть d(хi, хj) – некоторая мера близости между каждой парой классифицируемых объектов i и j. В качестве таковой может использоваться любая полезная функция: евклидово или манхэттенское расстояние, коэффициент корреляции Пирсона, расстояние 2, коэффициенты сходства Жаккара, Съеренсена, Ренконена и многие другие.

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

Рассмотрим в качестве примера набор данных USArrests о криминогенной обстановке по штатам США (из пакета cluster). С помощью тепловой карты, представленной на рис. 2.12, можно: а) разбить 50 штатов США на группы по криминальной напряженности, б) установить характер взаимосвязи между числом арестованных за убийство (Murder), изнасилование (Rape), разбой (Assault), а также долей городских жителей (UrbanPop) и, наконец, в) убедиться в том, что число изнасилований в Аляске много меньше, чем в Род-Айленде.

library(cluster) data("USArrests") x = as.matrix(USArrests) rc - rainbow(nrow(x), start = 0, end =.3) cc - rainbow(ncol(x), start = 0, end =.3) hv - heatmap(x, scale = "col", RowSideColors = rc, ColSideColors = cc, margins = c(10,10), cexCol = 1.5, cexRow = 1) Рис. 2.12. Пример тепловой карты с дендрограммами Однако существенным недостатком иерархических методов является их детерминированность: объединение классов на более высоком уровне полностью определяется результатами агломерации на нижних уровнях. Для оценки того, какой из вариантов кластеризации в наибольшей степени отражает близость объектов в исходном признаковом пространстве, могут быть использованы кофенетическая корреляция (cophenetic correlation) или различные ранговые индексы.

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

Выполним разбиение 50 штатов на 5 классов по степени их криминализации:

df.stand - as.data.frame(scale(USArrests)) (clus - kmeans(df.stand, centers = 5)) K-means clustering with 5 clusters of sizes 13, 11, 12, 7, 7

Cluster means:

Murder Assault UrbanPop Rape 1 -0.2162425 -0.2611064 -0.04793489 -0.06172647 2 -1.1034717 -1.1654231 -0.99194587 -1.04874074 3 0.7298036 1.1188219 0.75717991 1.32135653 4 1.5803956 0.9662584 -0.77751086 0.04844071 5 -0.6958674 -0.5679476 1.12728218 -0.55096728

Within cluster sum of squares by cluster:

[1] 10.860162 8.499862 18.257332 6.128432 5.244931 (between_SS / total_SS = 75.0 %) Мы видим, что к каждому из пяти кластеров относится от 7 до 13 штатов, причем наибольший криминальный риск имеет место в 3 и 4 группах.

Любой алгоритм кластеризации может считаться результативным, если выполняется гипотеза компактности (Загоруйко, 1999), т.е. можно найти такое разбиение объектов на группы, что расстояния между объектами из одной группы (intra-cluster distances) будут меньше некоторого значения 0, а между объектами из разных групп (cross-cluster distance) – больше.

Можно сформировать таблицу всех этих расстояний как показано ниже:

–  –  –

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

c.pca - prcomp(USArrests, center = TRUE, scale = TRUE) d - data.frame(x=c.pca$x[, 1], y=c.pca$x[, 2]) d$cluster - clus$cluster library('ggplot2'); library('grDevices') h - do.call(rbind, lapply(unique(clus$cluster), function(c) { f - subset(d,cluster==c); f[chull(f),]})) ggplot() + geom_text(data=d, aes(label=cluster, x=x, y=y, color=cluster), size=3) + geom_polygon(data=h, aes(x=x, y=y, group=cluster, fill=as.factor(cluster)), alpha=0.4, linetype=0) + theme(legend.position = "none")

Рис. 2.13. Пример кластеризации методом k средних

Другим современным подходом к кластеризации объектов являются алгоритмы типа нечетких C-средних и Гюстафсона–Кесселя (Барсегян и др., 2004), которые ищут кластеры в пространстве нечетких множеств в форме эллипсоидов, что делает эти методы более гибкими при решении различных практических задач. В литературе также описывается множество других методов кластеризации, не использующих матрицы сходств и основанных на оценивании функций плотности статистического распределения, эвристических алгоритмах перебора, идеях математического программирования. Каждый из них имеет специфическую идеологическую основу и подходы к построению критериев качества моделей, которые мы рассмотрим далее в главе 10.

3. ПАКЕТ СARET - ИНСТРУМЕНТ ПОСТРОЕНИЯ СТАТИСТИЧЕСКИХ

МОДЕЛЕЙ В R

3.1. Универсальный интерфейс доступа к функциям машинного обучения в пакете caret Использование сложных алгоритмов построения статистических моделей становится все более распространенной тенденцией в самых разных областях от академических исследований до всевозможных бизнесприложений. Среда статистических вычислений R отличается особенно богатым набором пакетов для создания различных моделей классификации и регрессии. Однако такое высокое разнообразие реализованных алгоритмов создает и ряд некоторых проблем.

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

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

Разработчики поставили перед собой следующие задачи:

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

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

создание пакета с постоянно расширяющимся набором методов;

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

Результатом этой работы стал пакет caret (сокращение от Classication and Regression Training), который сегодня стал одним из наиболее популярных инструментов среди пользователей R, занимающихся разработкой предсказательных моделей. Основные возможности пакета caret достаточно полно представлены в публикациях разработчиков (Kuhn, 2008, 2012, 2013; Kuhn, Johnson, 2013).

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

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

1. Определение наборов данных и их предобработка (при необходимости)

2. Определение спецификации параметров модели

3. Цикл для каждого параметра модели:

4. | Цикл для каждой итерации ресэмплинга:

5. | | Выделение тестовой выборки 6. | | Подгонка модели по объектам обучающей выборки 7. | | Прогнозирование отклика для тестовых объектов 8. | end 9. | Вычисление показателей средней эффективности прогноза

10. end

11. Установление оптимальных параметров модели

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

1. Разбиение исходных данных на обучающую и контрольную выборки.

Формирование индексов любой такой последовательности для вектора x случайным образом в заданном соотношении p осуществляется функцией FALSE), что, в целом, createDataPartition(x, p=3/4, list = эквивалентно использованию функции sample(length(x), size = 0.75*length(x)). C помощью аналогичных функций createResample(), createFolds() и createMultiFolds() можно создать бутстреп-выборки или произвольные разбиения на части в любом количестве повторностей.

2. Разведочный анализ исходных переменных. Графический анализ характера распределения данных и их взаимной зависимости, наличия выбросов и других аномальных ситуаций может быть выполнен с помощью функции featurePlot(). Функция nearZeroVar() позволяет исключить предикторы, дисперсия значений которых близка к нулю, а findCorrelation() позволяет выявить переменные, которые в значительной мере коррелируют с другими признаками.

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

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

5. Оценка "важности" предикторов (variable importance) и селекция оптимального их набора. Функция varImp() рассчитывает для каждой переменной количественные показатели, отражающих их вклад в получение точных предсказаний в рамках построенной модели. Отбор информативного комплекса переменных может быть осуществлен также с помощью функции rfe(), выполняющей рекурсивное исключение переменных, или функции gafs (), в которой реализован генетический алгоритм.

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

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

В качестве примера используем набор данных GermanCredit, входящий в состав пакета caret. Он содержит информацию по 1000 клиентам одного из немецких банков, (подробное описание см. на сайте UCI Machine Learning Repository).

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

отклика Class, принимающего два значения Good (хороший) и Bad (плохой).

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

Подобные предикторы с околонулевой дисперсией (near-zero variance) рекомендуется удалять из дальнейшего анализа (Kuhn, Johnson, 2013).

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

library(caret) data(GermanCredit) (u -unique(GermanCredit$ResidenceDuration)) # Доля уникальных значений:

length(u)/nrow(GermanCredit) [1] 4 2 3 1 [1] 0.004 Как видим, 1000 клиентов представлены только четырьмя уникальными значениями, доля которых составляет лишь 0.4% от их общего числа.

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

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

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

(t- sort(table(GermanCredit$ResidenceDuration), decreasing = TRUE)) t[1]/t[2] [1] 1.340909 На практике предлагается придерживаться следующих эмпирических правил для заключения о том, что некоторый предиктор обладает околонулевой дисперсией:

доля его уникальных значений от общего числа наблюдений составляет не более 10%;

отношение частот первых двух наиболее обычных его значений превышает 20.

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

# Создадим копию данных без столбца с откликом Class:

gcred = GermanCredit[, -10] # Функция nearZeroVar() возращает вектор номеров переменных, # обладающих околонулевой дисперсией:

(nz = nearZeroVar(gcred)) print("Имена этих переменных:"); names(gcred)[nz] # Удаляем предикторы с околонулевой дисперсией:

gcred.clean = gcred[, -nz] [1] 9 14 15 23 24 26 27 29 33 44 46 53 58 [1] "Имена этих переменных:" [1] "ForeignWorker" "CreditHistory.NoCredit.AllPaid" [3] "CreditHistory.ThisBank.AllPaid" "Purpose.DomesticAppliance" [5] "Purpose.Repairs" "Purpose.Vacation" [7] "Purpose.Retraining" "Purpose.Other" [9] "SavingsAccountBonds.gt.1000" "Personal.Female.Single" [11] "OtherDebtorsGuarantors.CoApplicant" "OtherInstallmentPlans.Stores" [13] "Job.UnemployedUnskilled" Описанная процедура вызывает, вероятно, некоторые сомнения у читателя, знакомого с основами метрологии. Во-первых, сама по себе близость дисперсии к нулевому значению ни о чем не говорит, поскольку все зависит от шкалы измерений (для субмолекулярных конструкций изменения размера в миллимикроны могут быть существенными, в то время, как для астрономических наблюдений ошибка в несколько километров является ничтожной). Поэтому дисперсия измерений считается недопустимо малой, если она не превосходит оценки погрешности измерений (или ошибки воспроизведения опыта). Во-вторых, следует уточнить, что эта процедура является эвристикой, полезной лишь для наблюдений, измеренных в порядковых или счетных шкалах с небольшим интервалом значений.

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

Выберем наиболее коррелирующие пары переменных из таблицы

GermanCredit. Составим для этого специальную функцию (А. Шипунов):

# Наибольшие значения треугольной матрицы top.mat - function(X, level=0.45, N=12, values=TRUE) { X.nam - row.names(X) X.tri - as.vector(lower.tri(X)) X.rep.g - rep(X.nam, length(X.nam))[X.tri] X.rep.e - rep(X.nam, each=length(X.nam))[X.tri] X.vec - as.vector(X)[X.tri] X.df - data.frame(Var1=X.rep.g, Var2=X.rep.e, Value=X.vec) {if (values) {X.df - X.df[abs(X.df$Value) = level, ] X.df - X.df[order(-abs(X.df$Value)), ]} else {X.df - X.df[order(-abs(X.df$Value)), ] X.df - X.df[1:N, ]}} row.names(X.df) - seq(1, along=X.df$Value) return(X.df) } top.mat(cor(gcred.clean)) Var1 Var2 Value 1 OtherInstallmentPlans.None OtherInstallmentPlans.Bank -0.8405461 2 Housing.ForFree Property.Unknown 0.7798526 3 Personal.Male.Single Personal.Female.NotSingle -0.7380357 4 Housing.Own Housing.Rent -0.7359677 5 OtherDebtorsGuarantors.Guarantor OtherDebtorsGuarantors.None -0.7314079 6 CreditHistory.Critical CreditHistory.PaidDuly -0.6836174

–  –  –

В состав пакета caret входит функция findCorrelation(), которая, как следует из ее названия, находит предикторы, чей уровень корреляции с другими предикторами в среднем превышает некоторый заданный пользователем порог (аргумент cutoff):

# Функция findCorrelation() возвращает вектор # номеров переменных с высокой корреляцией:

(highCor = findCorrelation(cor(gcred.clean), cutoff = 0.75)) print("Имена этих переменных:") names( gcred.clean)[highCor] # Удаляем эти переменные:

gcred.clean = gcred.clean[, -highCor] [1] 40 42 [1] Имена этих переменных:

[1] "Property.Unknown" [2] "OtherInstallmentPlans.None" Наконец, специальная функция findLinearCombos() предназначена для нахождения и исключения переменных, связанных линейными зависимостями. Если в выборке, например, есть переменная, которая может быть выражена через сумму нескольких других переменных, то такой предиктор можно рассматривать как малоинформативный (несущий дублирующую информацию) и его можно исключить из построения модели.

Результатом функции findLinearCombos() является список из двух элементов: $linearCombos – список найденных линейных комбинаций и $remove – вектор индексов переменных, которые можно выразить через линейную комбинацию остальных переменных.

(linCombo - findLinearCombos(gcred.clean)) # Удаляем эти переменные:

gcred.clean = gcred.clean[, -linCombo$remove] dim(gcred.clean) $linearCombos[[1]] [1] 30 9 10 11 12 26 27 28 29 $linearCombos[[2]] [1] 34 9 10 11 12 31 32 33 $linearCombos[[3]] [1] 43 9 10 11 12 41 42 $remove [1] 30 34 43 В результате всех этих операций число предикторов сократилось с 61 до 43.

Отметим, что очистка исходных данных от "ненужных" предикторов с использованием представленных функций может показаться излишне радикальной. Во-первых, например, непонятно, почему рекомендуется удалить признак OtherInstallmentPlans.None, а не составляющий с ним корреляционную пару OtherInstallmentPlans.Bank.

Во-вторых, разработаны более продвинутые методы минимизации "корреляционных плеяд":

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

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

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

caret. Ее основными аргументами являются следующие:

preProcess(x, method = c("center", "scale"), na.remove = TRUE, verbose = FALSE), где x – таблица или матрица с исходными данными (все переменные должны быть количественными);

method – список с названиями методов трансформации;

verbose – флаг для указания необходимости выводить протокол обработки.

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

преобразование отдельных предикторов: "center", "scale", "range","expoTrans", "BoxCox", "YeoJohnson";

групповая трансформация подмножества переменных: "spatialSign", "pca" и "ica";

заполнение пропущенных значений: "knnImpute", "bagImpute", "medianImpute".

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

Она заключается в вычитании из исходного значения xi некоторой переменной X соответствующего среднего значения x ("центрирование", или "center") и последующего деления полученной разницы на стандартное отклонение этой переменной x ("нормализация", или "scale"):

xi= (xi x )/x.

В результате стандартизации все количественные переменные будут иметь среднее значение, равное 0, и стандартное отклонение, равное 1.

Например, в таблице GermanCredit содержатся переменные, измеренные в разных шкалах: размер кредита Amount измеряется в немецких марках, возраст клиента Age – в годах, и т.д. Как следствие, размах значений переменных также существенно разнится: индикаторные переменные по определению варьируют от 0 до 1, тогда как размер кредита изменяется от 250 до 18424 марок.

Выполним стандартизацию трех метрических переменных:

data(GermanCredit) TransPred - c("Duration","Amount","Age") preVar - preProcess(GermanCredit[, TransPred]) TransVar = predict(preVar, GermanCredit[, TransPred]) print("До преобразования:") summary(GermanCredit[,TransPred]) print("После преобразования:") summary(TransVar) # при необходимости, обновление переменных GermanCredit[,TransPred] - TransVar [1] "До преобразования:" Duration Amount Age Min. : 4.0 Min. : 250 Min. :19.00 1st Qu.:12.0 1st Qu.: 1366 1st Qu.:27.00 Median :18.0 Median : 2320 Median :33.00 Mean :20.9 Mean : 3271 Mean :35.55 3rd Qu.:24.0 3rd Qu.: 3972 3rd Qu.:42.00 Max. :72.0 Max. :18424 Max. :75.00 [1] "После преобразования:" Duration Amount Age Min. :-1.402e+00 Min. :-1.070e+00 Min. :-1.455e+00 1st Qu.:-7.383e-01 1st Qu.:-6.751e-01 1st Qu.:-7.513e-01 Median :-2.407e-01 Median :-3.372e-01 Median :-2.238e-01 Mean : 1.260e-16 Mean : 6.628e-17 Mean : 5.330e-17 3rd Qu.: 2.568e-01 3rd Qu.: 2.483e-01 3rd Qu.: 5.674e-01 Max. : 4.237e+00 Max. : 5.368e+00 Max. : 3.468e+00 Опция method функции preProcess() может принимать и другие значения. В частности, значение "range" приводит значения переменных к диапазону [0, 1], а "expoTrans" выполняет вычисление экспоненциальной функции.

Некоторые из переменных имеют также явно выраженные асимметричные распределения (например, Amount и Age), что может представлять проблему для классических статистических методов. Часто решить эту проблему позволяют такие простые преобразования исходных значений, как извлечение квадратного корня, логарифмирование или расчет обратных значений. Если истинное нормализующее преобразование неизвестно, лучшим считается преобразование Бокса-Кокса (Box-Cox transformation – Box, Cox, 1964), которое позволяет найти оптимальное решение, в первую очередь, для нормализации дисперсии.

Универсальное семейство преобразований Бокса-Кокса (БК) случайной величины x зависит от значения параметра :

при 0 мы имеем дело со степенным преобразованием вида x 1 где показатель степени может принимать любые y ( ) =, положительные или отрицательные значения;

в частных случаях после подстановки параметра в основную формулу будем иметь: y = 1/x при = -1; y = 1 / x при = -0.5; y = x при = 0.5;

y = x2 при = 2 и т.д.;

при = 0 используется логарифмическое преобразование y() = ln(x) (поскольку деление на нуль приводит к неопределенности).

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

Применим к метрическим показателям таблицы БК-преобразование preBox - preProcess(GermanCredit[, TransPred], method="BoxCox") Pre-processing: Box-Cox transformation

Lambda estimates for Box-Cox transformation:

0.1, -0.1, -0.7 Для предиктора Amount мы получили расчетное значение = -0.1. Результат преобразования с этим параметром представим графически (рис. 3.1):

BoxVar - predict(preBox,GermanCredit[,TransPred]) y - factor(GermanCredit$Class) trellis.par.set(theme = col.whitebg(), warn = FALSE) featurePlot(GermanCredit$Amount, y, "density", labels = c("Amount","")) featurePlot(BoxVar[,2], y, "density", labels = c("Box-Amount","")) Рис. 3.1. Результат преобразования Бокса-Кокса для размера кредита в немецком банке При значении параметра method="YeoJohnson" выполняется трансформация Йео-Джонсона (Yeo-Johnson), которая весьма похожа на БКпреобразование, но учитывает в расчетах нулевые и отрицательные значения обрабатываемых переменных. Другой функцией, также позволяющей оценить оптимальное значение для набора исходных данных, является функция BoxCoxTrans() из пакета caret.

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



Pages:   || 2 | 3 | 4 | 5 |

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

«Спіс крыніц і літаратуры асноўная 1. Хан-Пира Э. И. Архивоведческое терминоведение: пособие по спецкурсу.– М.: МГИАИ, 1990.– 136 с.2. Делопроизводство и документоведение: Краткий словарь современной терминологии / авт.-...»

«64 3. ТЕХНОЛОГИЧЕСКАЯ ЧАСТЬ.3.1. Расчет производственной программы 3.1.1. Планирование раскроя и баланс перерабатываемого сырья Рациональное использование сырья и получение пилопродукции определенных размеров и...»

«ОКРУЖНАЯ ДЬОКУУС КАП КУОРАТ АДМИНИСТРАЦИЯ У 0 К У Р У Г У I ! ДЬАЬЛДТАТА ГОРОДА ЯКУТСКА П РО !О КО I заседания комиссии по рассмотрению предложений о присвоении имен государственных, муниципальных, общественных деятелей му и и ци и а 1ы 1ы м и ред п р и я ти я м и у ч режде н ия м городского округа "город...»

«Оценка грамматической стороны речи и связной речи Оценка грамматической стороны речи и связной речи Оценка связной речи. Постепенно у ребенка развивается способность описывать, объяснять, пересказывать все то, что он видит кругом и понимает. Таким образом, формируется связная речь, кот...»

«ПРИМЕНЕНИЕ ПРОИЗВОДНОЙ ФУНКЦИИ Уравнение касательной Рассмотрим следующую задачу: требуется составить уравнение касательной l, проведенной к графику функции в точке. Согласно геометрическому смыслу производн...»

«Pilipenko N. "A humble man": the poetic world-view of F.Tyutchev. The principles of artistic consciousness by F. Tyutchev are analyzed from positions of Christian anthropology through its main concepts – Love, Humility, Soul; the comp...»

«Особенности подготовки и реализации проектов мобильного маркетинга в сфере В2В Когда возникает необходимость в реализации проектов мобильного маркетинга на B2B рынке? Необходимость задействовать в кампании широкую аудиторию Адресное взаимодействие с потребителем Контроль за акцией Повышение лояльности к бренду: актуальный ти...»

«АРБИТРАЖНЫЙ СУД КЕМЕРОВСКОЙ ОБЛАСТИ Красная ул., д.8, Кемерово, 650000 http://www. kemerovo. arbitr. ru/ тел. (384-2) 58-17-59, факс (384-2) 58-37-05 e-mail: info@kemerovo. arbitr. ru. ИМЕНЕМ РОССИЙСКОЙ ФЕДЕРА...»

«1938 УСПЕХИ ФИЗИЧЕСКИХ НАУК Т. XIX, вып. 3 УСТАНАВЛИВАЮЩИЕСЯ ПРОЦЕССЫ В АКУСТИКЕ 1) Г. Бакгауз СОДЕРЖА И Б IV. Н е к о т о р ы е в о п р о с ы э л е к т р о а к у с т и к и.9. Методы и;слеювания нестационарных акустических процессов.10. Устанавливаю диеся процессы в электроакустических аппаратах. V. У с т а н...»

«Little Greene Характеристика: Intelligent Exterior Eggshell (Водная полуматовая краска для внешних работ) ОПИСАНИЕ Качественная водная полуматовая краска, образует прочное покрытие, стойкое к воздействию UV-лучей; подходит для работ по...»

«Муниципальное автономное образовательное учреждение для обучающихся с ограниченными возможностями здоровья средняя общеобразовательная школа-интернат № 6 города Ялуторовска "Согласовано" "Согласовано" "Утверждено" Протокол...»

«Проект ПРАВИТЕЛЬСТВО РОССИЙСКОЙ ФЕДЕРАЦИИ РАСПОРЯЖЕНИЕ от "" 2011 г. № МОСКВА 1. Утвердить прилагаемую Стратегию развития страховой деятельности в Российской Федерации до 2020 года (далее Стратегия).2. Федеральным органам исполнительной власти руководствоваться положения...»

«ПРАВИТЕЛЬСТВО НОВГОРОДСКОЙ ОБЛАСТИ ПОСТАНОВЛЕНИЕ 27.05.2016 № 194 Великий Новгород Об утверждении Положения о порядке предоставления в 2016-2017 годах субсидий социально ориентированным некоммерческим организациям Новгородской области В целях реализации мероприятий государственной пр...»

«Монтаж лепнины из полиуретана, установка лепнины Для монтажа рекомендуется использовать специальные клеи от компании Европласт: Клей Европласт Экстра Клей Европласт Стандарт Специальный клей для Специальный клей для соединения карнизов и монтажа карнизов, молдингов меж...»

«Утверждено “ ” г. ЗАО "ФБ ММВБ" (подпись уполномоченного лица) (печать) ИЗМЕНЕНИЯ В РЕШЕНИЕ О ВЫПУСКЕ ЦЕННЫХ БУМАГ Акционерный коммерческий банк АК БАРС (открытое акционерное общество) облигации документарные на предъявителя с обязательным централизованным хранением серии БО-05 биржевые неконвертируемые процентные с возможностью доср...»

«DOI: 10.15393/j9.art.2008.272 Н. А. Криничная Петрозаводск "Я БЫЛ В ДУХЕ В ДЕНЬ ВОСКРЕСНЫЙ." (Стихотворение Н. А. Клюева в контексте христианских и народно-христианских традиций) Эпиграфом к названному с...»

«Инструкция по эксплуатации Стоматологический интернет-магазин www.stomamart.ru toptopkr@gmail.com STA TM Компьютеризированная Система локальной анестезии с наконечником STA Wand STA-5220 220 Вольт Caution: Federal law restricts thi...»

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

«Б. В. А Н А II Ь И Ч УЧЕТНО-ССУДНЫЙ БА Н К ПЕРСИИ В 1894-1907 гг. Русское самодержавие втянулось в империалистическое сопер­ ничество за сферы влияния на огромных территориях полуколо­ ниальных и зависимых стран Дальнего и Сред...»

«Навчальна програма з дисципліни “ Управління в системах зв’язку” 1. Введение 1.1. Объект изучения. Объект изучения – автоматические системы и устройства в технике связи, методы анализа и синтеза разомкнутых и замкнутых автоматических устройс...»

«Судовые энергетические установки 2015 – № 35 142 УДК 629.129.359.1 Репетей В.Д., Костенко П.А. ОНМА ОСОБЕННОСТИ БУКСИРОВКИ СУДОВ В ПОРТОВЫХ ВОДАХ Буксирная система "буксировщик – буксирная связь – буксируемый объект" представляет весьма сложную в управлении структуру. Не все этапы буксирной операции достаточно изучены, а потому требуют определённ...»

«Теоретична і дидактична філологія. Випуск 20, 2015.16. Niedzielski Nancy A. Folk Linguistics / Nancy A. Niedzielski, Dennis Preston // Trends in Linguistics (Studies and Monofraphs 122). – Berlin : Walter de Gruyter, 2000. – 375 p.17. Radden Gnter. The folk model of language / Gnter Radden. – [Electronic resource]. – Режим дос...»

«http://vmireskazki.ru vmireskazki.ru › Европейские сказки › Польские сказки Свинья не коза, а дуб не берёза Польские сказки Вернулся домой Простачок, отдохнул от долгого пути. А тут и праздник подошел. Запряг он свою лошадёнку, взвалил на тел...»

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

«Письма протоіерея Арсенія ЛеОединцева, б. благочиннаго церквей берега Крыма, к ъ преосвященному Иннокенті", архіепископу Херсон Таврическому, съ донесеніями о ю д военныіъ дйствШ и состояніи и духовенства во время 11-тиісячной осады Севастополях). П и с ь м о 47-е. Мы...»

«Русский язык. 8 класс. Тематическая работа "Обособленные члены предложения. Обращения, вводные и вставные конструкции, междометия, оформление чужой речи". 1 Демонстрационный вариант Инструкция по выполнению работы На выполнение тематической работы по русскому языку даётся 45 минут....»

«Инструкция о бюджетному учету 148н 2-04-2016 1 Непоправимо проломленная это обожествление. Самовосстанавливающееся обыгрывание силою не заебывает пляжный карбюратор докрасна клявшим транслятором. Сохраняют ли правда не высунувшие подкомиссии? Отборные коллективизации кощунствуют...»

«Российская Академия Наук Институт философии ЧЕЛОВЕЧЕСКИЙ ПОТЕНЦИАЛ КАК КРИТИЧЕСКИЙ РЕСУРС РОССИИ Москва УДК 308+300-31 ББК 60.59(2)+15.56 Ч-39 Ответственный редактор доктор филос. наук Б.Г. Юдин Рецензен...»

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

«№3(2012) (2012) Учредитель, редакция и издатель журнала : ООО "Витпостер" Главный редактор: АВРУТИН Анатолий Юрьевич Редакционная коллегия: Анатолий АНДРЕЕВ Глеб АРТХАНОВ протоиерей Павел БОЯНКОВ Алексей ВАРАКСИН Иван ГОЛУБНИЧИЙ (Москва) Светлана ЕВСЕЕВА Николай КОЗЛОВ Николай КОНЯЕВ (Санкт-Петербург) Владим...»










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

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