Создание предсказательных моделей: основные шаги

Как отмечалось в одном из предыдущих сообщений, статистические модели создаются либо с целью получения предсказаний, либо для понимания взаимоотношений между переменной-откликом и предикторами. Создание предсказательных моделей (англ. "predictive models", в некоторых областях также "forecasting models") включает ряд стандартных шагов, обзор которых приведен ниже. Этим сообщением я начинаю новую серию публикаций по методам машинного обучения ("machine learning"), или "статистического обучения" ("statistical learning"; James et al. 2013), и построению предсказательных моделей с их помощью.

1. Разведочный анализ данных

Для простоты изложения здесь мы предположим, что данные уже получены из соответствующих ситуации источников (хранилища данных, рынки данных, базы неструктурированных данных, протоколы лабораторных или полевых экспериментов и т.п.), "очищены" от ошибочных наблюдений и оформлены как "опрятные данные" ("tidy data", Wickham 2014; см. также R-пакеты tidyr и dplyr). Разведочный анализ данных (РАД) является первым и очень важным шагом в процессе создания предсказательных моделей. Главная его цель - понимание свойств имеющихся в наличии переменных (распределение, наличие выбросов, необходимость трансформации и др.) и выявлении характера взаимоотношений между откликом и предикторами. Конкретные задачи, решаемые в ходе РАД, будут несколько варьировать от случая к случаю. Некоторые часто распространенные этапы РАД при построении статистических моделей были описаны в предыдущих сообщениях (см. ссылки в конце этой статьи).

В качестве простого примера используем данные по электрическому сопротивлению (Ом) мякоти фруктов киви в зависимости от процентного содержания в ней сока (Harker & Maindonald 1994). Таблица с этими данными (fruitohms) входит в состав пакета DAAG, который является приложением к книге Maindonald (2010) Data Analysis and Graphics Using R. Чаще всего наиболее эффективным способом выявления характера связи между подлежащей предсказанию переменной и предикторами является графический. В рассматриваемом примере у нас есть лишь один предиктор - содержание сока в мякоти фруктов, и поэтому достаточно изобразить данные в виде диаграммы рассеяния:

install.packages("DAAG")
library(DAAG)
data("fruitohms")
 
plot(ohms ~ juice, xlab = "Содержание сока, %", 
     ylab = "Сопротивление, Ом", data = fruitohms)

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

При наличии нескольких предикторов необходимо было бы исследовать их связь как с переменной-откликом, так и между собой. Полученная информация могла бы оказаться важной с точки зрения того, какие предикторы стоит включать в модель, необходимо ли некоторые из них преобразовать, как именно сделать преобразования, и т.д. Как правило, чем больше предикторов мы включаем в модель, тем более точными становятся предсказания. Однако при избыточном усложнении модели качество ее предсказаний начнет снижаться в связи с явлением, известным как переобучение модели ("model overfitting"). Для борьбы с переобучением используют разнообразные методы отбора информативных предикторов ("variable selection", или "feature selection") и перекрестную проверку ("cross-validation").

2. Разбиение исходных данных на обучающую и проверочную выборки

Одним из наиболее распространенных подходов при построении предсказательных моделей является выполненное случайным образом разбиение исходных данных на обучающую ("training set") и проверочную ("test set", или "validation set") выборки. Как следует из названия, обучающая выборка используется для "обучения" той или иной модели, т.е. для построения математических отношений между переменной-откликом и предикторами, тогда как проверочная выборка применяется для оценки предсказательных свойств модели на данных, которые она (модель) раньше "не видела". Как правило, обучающая выборка составляет 75-80% от объема исходных данных, хотя каких-то строгих правил в этом отношении не существует. Более того, при наличии небольшого числа наблюдений разбиение исходных данных на обучающую и проверочную выборки не всегда является возможным. В таких случаях используется перекрестная проверка с последовательным исключением одного наблюдения ("leave one out cross-validation", LOOC), когда строится n моделей по (n - 1) выборочным значениям, а исключенное значение каждый раз используется для оценки качества прогноза.

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

В нашем пример с фруктами киви мы не собираемся прогнозировать будущее - мы просто хотим иметь возможность предсказывать значения электрического сопротивления для той же генеральной совокупности, из которой были взяты исходные данные. Поэтому достаточно выполнить простое случайно разбиение данных на обучающую (80%) и проверочную (20%) выборки. Ниже приведен один из нескольких возможных способов того, как это можно сделать в R:

set.seed(101) # для воспроизводимости результата
 
# случайный выбор индексных номеров наблюдений, которые
# войдут в обучающую выборку:
tr.index = sample(1:nrow(fruitohms), nrow(fruitohms)*0.8) 
 
# обучающие данные:
trSet = fruitohms[tr.index, ]
 
# проверочные данные:
testSet = fruitohms[-tr.index, ]

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

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

Даже для нашего простого примера с фруктами киви можно подогнать несколько моделей, различающихся математически. Простейшей будет модель линейной регрессии (см. рис. ниже слева): \(y_i = \beta_0 + \beta_{1} x_i + \epsilon_{i}\). Однако, как мы выяснили в ходе РАД, зависимость между содержанием сока в мякоти фруктов и электрическим сопротивлением не совсем линейна, что можно выразить, например, при помощи полинома второй степени (см. рис. ниже справа): \(y_i = \beta_0 + \beta_{1} x_i + \beta_{2} x^2_{i} + \epsilon_{i}\).

Для выбора оптимальной модели аналитику необходимо определиться с критерием, который позволит выполнить объективное сравнение качества альтернативных моделей. При решении регрессионных задач, когда зависимой является количественная переменная, такого рода критерии, как правило, основаны на остатках ("residuals"), т.е. разницах между наблюдаемыми и предсказанными значениями отклика: \(y_i - \hat{y_i}\). В частности, распространенным критерием является квадратный корень из среднеквадратичной ошибки ("root mean squared error", RMSE), который рассчитывается как

\[RMSE=\sqrt{\frac{1}{n}\sum_{i=1}^{n}(y_i - \hat{y_i})^2}.\]

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

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

M1 = lm(ohms ~ juice, trSet) # простая линейная модель
M2 = lm(ohms ~ juice + I(juice^2), trSet) # полином 2-й степени
 
RMSE1 = sqrt(sum((trSet$ohms - predict(M1))^2)/nrow(trSet))
[1] 1123.313
 
RMSE2 = sqrt(sum((trSet$ohms - predict(M2))^2)/nrow(trSet))
[1] 1038.543

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

Для решения указанной проблемы (т.е. для получения несмещенных оценок предсказательных свойств моделей) применяются методы формирования повторных выборок ("resampling"), при помощи которых обучающие данные случайным образом разбиваются на несколько частей. Эти части затем попеременно используются для подгонки разных версий одной и той же модели. Каждый раз данные, не участвующие в подгонке конкретной версии модели, используются для расчета соответствующего критерия качества (например, RMSE). Полученные таким образом значения этого критерия далее усредняются, что позволяет получить несмещенную оценку предсказательных свойств модели. Наиболее распространенным способом получения подобных несмещенных оценок критериев качества является упомянутая выше перекрестная проверка. Обычно применяют 10-кратную перекрестную проверку ("10-fold cross-validation"), при которой обучающие данные случайным образом разбиваются на 10 блоков ("folds"). Одну и ту же модель затем подгоняют 10 раз по данным из 9 разных блоков, используя десятый блок для расчета критерия качества.

Выполним 10-кратную перекрестную проверку для нашего примера. Для этого воспользуемся функцией cvLm() из пакета cvTools, которая по умолчанию использует RMSE как критерий качества модели:

install.packages("cvTools")
library(cvTools)
 
cvLm(M1, K = 10)
10-fold CV results:
  CV 
1141 
 
cvLm(M2, K = 10)
10-fold CV results:
  CV 
1072

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

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


4. Выбор оптимальной модели на основе проверочных данных

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

sqrt(sum((testSet$ohms - predict(M2, testSet))^2)/nrow(testSet))
[1] 1022

что очень неплохо, особенно учитывая, что RMSE, полученный для этой модели в ходе 10-кратной перекрестной проверки был равен 1072. Для сравнения,  RMSE по проверочным данным для модели M1 составил

sqrt(sum((testSet$ohms - predict(M1, testSet))^2)/nrow(testSet))
[1] 1090

Таким образом, окончательный выбор следует сделать в пользу M2, т.к. эта модель с большой вероятностью будет давать более точные предсказания, чем M1.


Последующие сообщения будут посвящены разнообразным аспектам описанных выше шагов процедуры создания предсказательных моделей. Упор будет сделан на возможности, реализованные в пакете caret (Kuhn 2008; Kuhn & Johnson 2013).