Как бы я рассказал про линейную регрессию (если б меня кто-то спросил)

Про мотивацию ^[1]. Если любопытно почитать мысли автора, раскрывайте этот абзац

Пиша Расписывая первые строчки этой заметки на 100+ страниц я не могу не ответить на важный вопрос, во первых, самому себе, во вторых, читателю – “Зачем опять писать про линейную регрессию”. Во первых, потому что так хочется. Во вторых, потому что я считаю что могу по новому и непохоже на других авторов рассказать про машинное обучение ^[2] на примере линейной регрессии. Чем же хочется выделится на этой сцене:

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

• Использовать визуальный ряд как основную повествовательную компоненту. Такая статья в виде комикса – когда текст читать полезно, но в целом не обязательно, потому что даже пробежавшись по картинками и анимациям можно будет получить полное представление о том как оно работает. Такой концепт меня также привлекает своей полезностью когда времени мало, а напомнить себе как оно там делается, хочется (например когда через 40 минут завтра собеседование на ML инженера);

• Я в пункте выше писал про “комикс” – забудьте. Компьютерная наука ^[3] – это среда, раскрывающаяся в динамике, поэтому я считаю, что важно показать ключевые элементы в виде анимаций;

• (Почти) каждый график и алгоритм должен быть описан в коде, чтобы можно было запустить самому из любопытства. Или использовать в своих докладах;

• Фокус на практике по принципу “от боли ^[4] к решению” – при продвижении по вниз по тексту, каждый последующий шаг совершается для решения возникшей проблемы на текущем шаге. Я надеюсь что такая последовательность создает необходимую связность между блоками.

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

Для кого эта статья

Да не читайте вы этот абзац, – лучше проскролльте статью в течение двух минут до конца. Акцентируя внимание только на визуализациях вы сможете понять стоит ли вчитываться (поскольку все ключевые темы обозначены не только в тексте, но и в графиках и анимациях). От себя же добавлю, что данный пост может быть полезен как и начинающим дата саентистам и всем специалистам кто так или иначе работает с данными, так и коллегам с опытом ^[5], которым хочется освежить в памяти ^[6] некоторые аспекты.

О чем речь

В данном посте три повествовательных акта:

• Что такое линейная регрессия, зачем нужна и как её построить

• Как измерить качество модели

• Как улучшить качество модели, если мы недовольны 

Глобально, мы поговорим о:

1. моделировании на основе данных;

2. аналитическом решении линейной регрессии и почему иногда им воспользоваться не получится; 

3. способах измерения качества моделирования – визуально и при помощи метрик;

4. многомерной (многофакторной) линейной регрессии когда для предсказания используется много данных;

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

6. методах улучшения качества моделирования: от усложнения модели до её упрощения (регуляризации).

В частности, мы рассмотрим: 

1. метод наименьших квадратов для парной регрессии;

2. регрессионные метрики R², RMSE, MAE, MAPE, SMAPE, коэффициент корреляции и коэффициент детерминации а также методы визуальной оценки, например графики остатков;

3. метод максимального правдоподобия и предсказательный интервал;

4. разбиение на тренировочную и тестовую выборки: зачем нужно и как делать;

5. методы фильтрации выбросов: Консенсус случайной выборки, Махаланобисово расстояние, Локальный фактор выброса, Расстояние Кука;

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

7. линейную алгебру стоящую за методом наименьших квадратов и то как она работает для многомерной регрессии;

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

9. L1 и L2 регуляризацию линейных моделей;

10. кросс валидацию и оптимизацию гиперпараметров. 

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

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

Все анимации и картинки в статье оригинальные и выполнены автором.

Маленький литературный обзор

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

• Линейная регрессия на Python: объясняем на пальцах (2019) ^[7] – материал с наглядным примером и кодом на Python для инициализации модели линейной регрессии;

• Как я покупал квартиру ^[8], классика из хабра 2013 – статья представляет собой больше мини-историю по сбору и анализу данных с построением линейной модели в конце, чем пост полностью посвященный регрессии. Зато для анализа используется R – добавляет разнообразия;

• Линейные модели. Учебник по машинному обучению от преподавателей и выпускников Школы анализа данных ^[9] – раздел хендбука удачно объясняющий основную и смежные (например, численные методы) темы, хотя и сжато;

• Линейная регрессия. Основная идея, модификации и реализация с нуля на Python ^[10], хабр 2024 года – хорошо продемонстрированы разные вариации линейной регрессии, в том числе, и код как их обучать и визуализировать;

• Основы линейной регрессии ^[11], статья на хабре от 2020 года – на мой взгляд хороший заход с объяснения интерполяции и аппроксимации, а также формальными выкладками. Особенно полезным, считаю, может быть залипание на гифки;

• Открытый курс машинного обучения. Тема 4. Линейные модели классификации и регрессии ^[12], Open Data Science на хабре от 2017 – статья может являться неплохой входной точкой в курс, плюс хорошо расставлены акценты на регуляризации

Отдельно хочу выделить “Курс Линейные модели, дисперсионный и регрессионный анализ с использованием R” (см. материалы курса в гитхаб репозитории ^[13]): курс охватывает большое количество тем и при этом соблюдается (на мой взгляд) баланс между простотой и формальностью изложения. 

И конечно, не обходите стороной классические работы. В данном разделе они не приведены в виде списка литературы, однако ниже на них будут ссылки. Эти статьи приведены сразу после фрагмента который на них ссылается в квадратных скобках в формате [Авторы(ы). Название. Год публикации. Ссылка на первоисточник]

Хорошая модель начинается с данных

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

• Количество комнат в квартире

• Стоимость апартаментов (для удобства, стоимость будет указана в долларах $, иначе бы много нулей мелькало на экране)

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

Зачем нужна модель 

Как говорил британский статистик Джордж Бокс, “Все модели неточны, но некоторые из них полезны”. И полезны они для выявления закономерностей в данных. Как только закономерности описаны в виде математического выражения (модели), мы можем использовать это для, например, генерации предсказаний (Рисунок 2). 

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

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

где x – предиктор, y – целевая переменная [James, G., et al. Linear Regression. An Introduction to Statistical Learning, 2021. Открытая версия книги https://www.statlearning.com/ ^[16]]

То есть выражение это модель линейной регрессии. И – тоже. Разница только в коэффициентах. Поскольку коэффициенты в уравнении являются важнейшими параметрами, они имеют свои названия:

• – свободный член, или параметр (коэффициент) смещения

• – коэффициент наклона

Итак, когда мы строим модель линейной регрессии мы делаем следующее предположение: 

Предположение 1. Зависимость между признаками (независимыми переменными) и откликом (зависимой переменной) линейна [Kim, Hae-Young. Statistical notes for clinical researchers: simple linear regression 1 – basic concepts, 2018. https://www.rde.ac/upload/pdf/rde-43-e21.pdf ^[17]]

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

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

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

Немного про интерполяцию, экстраполяцию и аппроксимацию

Как построить простую модель

Итак, нам требуется подобрать коэффициенты и для уравнения ниже таким образом, чтобы прямая линия как можно точнее описывала эмпирические наблюдения (реальные данные): , где – количество комнат в квартире, – цена апартаментов, $

Почему именно такое уравнение и почему два коэффициента

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

Аналитическое решение

Для нахождения оптимальных значений коэффициентов воспользуемся аналитическим решением: подставим эмпирические данные из предыдущего раздела в заранее известную и выведенную умными людьми (Карл Гаусс, Адриен Лежандр) формулу. Аналитическое решение можно разбить на четыре шага (Рисунок 6). [Hastie, T., et al. Linear Methods for Regression (Chapter 3 in The Elements of Statistical Learning: Data Mining, Inference, and Prediction). 2009. https://hastie.su.domains/ElemStatLearn ^[20]]

Ошибка – тоже часть модели (когда двух коэффициентов недостаточно)

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

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

Как измерить качество модели

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

• Визуальная оценка

• Расчет метрик

Перед тем как предметно рассуждать о каждом пункте, предлагаю определить “качественность” как нашу удовлетворенность моделью. Удовлетворены же мы будем (в рамках данного поста) в том случае, если компонента ошибки будет минимально возможной.

Рассмотрим на примере исходных данных (см. Рисунок 1 ^[23]) и подставим разные значения коэффициентов в уравнение линейной регрессии. Затем сделаем предсказание для заранее известных примеров и сравним насколько сильно отличаются предсказанные и реальные значения (Таблица 1). Очевидно, что из множества комбинаций коэффициентов сдвига и наклона есть такая пара, которая приводит к минимальной средней ошибке.

Таблица 1. Сравнение ошибки на одном наблюдаемом значении (при количестве комнат равным двум) при различных значениях коэффициентов и

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

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

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

Визуальная оценка 

Воспользуемся формулой из раздела про Аналитическое решение (Рисунок 6 ^[25]). Подставим данные и получим следующие модели для каждого из датасетов: 

1. A:

2. B:

3. C:

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

Первый график, который есть смысл привести здесь, – это классическая точечная диаграмма (диаграмма рассеяния), где на ось x наносят значения признака, а по оси y приводятся показанные разным цветом предсказанные значения и реальные наблюдения. Анализировать такой рисунок довольно просто, – чем ближе модельная линия к реальным данным, тем качественнее модель. Ещё на графике видно зависимости между откликом и независимой переменной, поскольку перед глазами имеется и признак. [Piñeiro, G., et al. How to evaluate models: Observed vs. predicted or predicted vs. observed? 2008. https://doi.org/10.1016/j.ecolmodel.2008.05.006 ^[26]]

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

Решить проблему сложности визуализации, когда признаков становится много, можно с помощью диаграммы рассеяния предсказания vs наблюдения. График строится так: по оси x откладывают реальные значения, а по оси y – предсказанные (Рисунок 10) [Moriasi, D. N., et al. Hydrologic and Water Quality Models: Performance Measures and Evaluation Criteria. 2015. https://web.ics.purdue.edu/~mgitau/pdf/Moriasi%20et%20al%202015.pdf ^[28]]. Я встречал вариации, когда поступают наоборот: по оси x откладывают предсказанные значения. Тем не менее, свою функцию такие графики выполняют и так, и так – выбирайте на свой вкус ^[29].

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

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

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

Второй вариант, – классический Q-Q график двух выборок. В нём визуализируются не сырые упорядоченные значения а квантили, – данные разбиваются таким образом на конечное (я обычно использую около 100) число уровней. Такой график полезен, когда важно сравнить не отдельные пары “предсказание – отклик”, а общий рисунок данных (форму распределений): где находится медиана и насколько часто встречаются большие или очень маленькие значения.

Напоминание про то что такое квантиль и перцентиль

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

Временно вынося за скобки вероятность из определения выше, можно сказать, что квантиль – это такое значение для выборки, которое делит данные на части. Например, квантиль уровня 0.25 – это число, ниже которого (меньше) лежит 25% всех значений выборки. А квантиль уровня 0.9 – это значение, ниже которого 90% данных.

Для выборки [1, 3, 5, 7, 9] квантиль 0.5 (медиана) равен 5. Больше 5 есть только два значения: 7 и 9, и меньше тоже только два: 1 и 3.  

Квантиль уровня 0.25 ≈ 3. Квантиль уровня 0.75 ≈ 7. Смотрите пояснения на Рисунке ниже.

25й перцентиль (персентиль) также называется первым квартилем, 50й перцентиль – медиана или второй квартиль, 75й перцентиль – третий квартиль.

Во втором варианте, независимо от того насколько много объектов в датасете, количество отображаемых на графике точек всегда равно 99, а значит такой подход отлично масштабируется на сколь угодно большие выборки. На Рисунке 11 видно, что для датасета A квантили реальных и предсказанных значений лежат близко к диагонали – модель хорошая. В то время как  для датасета B правый хвост распределений (справа вверху) реальных и предсказанных значений расходится. Это значит, что модель хуже предсказывает цену дорогих квартир.

Для датасета C: 

• Ниже 25го перцентиля модельные данные лежат выше реальных

• В интерквартильном размахе (от 25го до 75го перцентиля) модельные значения лежат ниже 

• Выше 75го перцентиля опять хвост распределения предсказанных значений лежит выше реальных 

Если вам хочется углубиться в тему Q-Q графиков, то могу порекомендовать хабр статью Q-Q Plots. От чайника до профессионала за один гайд ^[32].

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

Именно по графику остатков удобнее всего анализировать следующие предположения для инициализации модели линейной регрессии (Предположение 1 ^[34] было выше по тексту – линейность): 

• Предположение 2. Нормальность остатков – ошибки модели (разности между фактическими и предсказанными значениями) должны быть нормально распределены (на вопрос “Что значит нормально распределены?” – ответ будет следующий. Большинство ошибок – маленькие, около нуля. А большие ошибки – редкость, и они встречаются примерно одинаково часто в обе стороны (в плюс и в минус)

• Предположение 3. Гомоскедастичность (постоянная дисперсия ошибок) – модель предсказывает с примерно одинаковой ошибкой и дешевые квартиры, и средний ценовой сегмент, и дорогие объекты

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

Из Рисунка 12 видно, что для датасета B не подтверждается второе предположение – с ростом количества комнат модель ошибается все больше и больше – дисперсия остатков увеличивается (разброс значений растет слева направо). Это значит, что ошибка для модели непостоянна и зависит от значения предиктора, – модель не учитывает какую-то скрытую (от нас и модели) закономерность и потому нестабильна.

Для датасета C не соблюдается условие нормальности остатков – модель иногда систематически завышает или систематически занижает значения уходя то выше, то ниже нуля в среднем. Кроме того на графике остатков вырисовываются “фигуры” (паттерны на графике) – это может быть признаком того что ошибки не независимы друг от друга (справедливости ради замечу, что может и не быть XD Но в любом случае, это сигнал, что что-то с моделью не ладно).

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

Я не знаю что такое частотное распределение

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

Частотное распределение – это упорядоченное представление, показывающее, сколько раз значения случайной величины попадают в определенные интервалы.

Для его построения:

1. Весь диапазон значений разбивают на k интервалов (интервальные частоты).

2. Подсчитывают, сколько наблюдений попало в каждый интервал – это называется абсолютной частотой.

3. Делят абсолютную частоту на общий объем выборки n – получают относительную частоту.

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

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

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

Существуют эмпирические формулы, которые позволяют подобрать адекватное значение количества интервалов в зависимости от объема выборки. Например формула Стерджеса или формула Райса (см. Рисунок ниже) [Sturges. The Choice of a Class Interval. 1926. DOI: 10.1080/01621459.1926.10502161], [Lane, David M., et. al. Histograms. https://onlinestatbook.com/2/graphing_distributions/histograms.html ^[38]].

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

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

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

Еще одним популярным способом представлять информацию о распределениях являются ящики с усами (Boxplot). Ящик с усами – это график, обобщённо описывающий распределение по квартилям. Он отображает:

• Медиану (второй квартиль)

• Первый Q1 и третий Q3 квартиль (25й и 75й перцентили соответственно) – границы «ящика»

• «Усы» – границы диапазона данных без выбросов

• Отдельные точки – выбросы

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

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

Частотное распределение является важным инструментом описания и изучения поведения ^[44] случайной величины по выборке. Разные визуальные методы, гистограмма, график ядерной оценки и ящик с усами, дополняют друг друга и позволяют исследователю получить полное представление о характере распределения, его симметричности, разбросе и выбросах. Далее по тексту мы будем пользоваться этими инструментами чтобы подсветить разные особенности в данных. А на данный момент раздел-напоминание закончен, – возвращаемся к главному нарративу про модели.

Продолжая тему работы с распределениями: распределение остатков, как мы заметили выше должно быть похоже на нормальное. И потому имеет смысл подготовить визуализации сравнения распределений – теоретического нормального и наблюдаемого у остатков модели. Удобно это делать через графики плотности распределения и квантильные диаграммы (квантили остатков Vs квантили нормального распределения). Оценка параметров нормального распределения производится по выборке остатков.  Поскольку такие графики работают лучше всего на больших объемах данных, то я предлагаю для наглядности искусственно увеличить размер выборок у остатков до 500 элементов в каждой, сохранив при этом основные черты поведения ^[45] остатков модели для каждого датасета (Рисунок 14). 

Как видно из рисунка 14, для датасета A* и B* распределение остатков хорошо аппроксимируется нормальным, хотя и для B* хвосты немного расходятся – модель допускает большие ошибки немного чаще, чем нам этого хотелось бы. В случае же с бимодальным C* различия особенно заметны – распределение остатков совсем не похоже на нормальное. 

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

Подытожим. Модель редко бывает идеальной – она будет ошибаться. А анализ ошибок при помощи графиков – удобный способ диагностики модели:

• Если предиктор всего один, то удобно построить график где показаны предсказанные значения вместе с реальными на оси y, а на оси x – признак. Сразу видно форму взаимосвязи между признаком и откликом;

• В дополнение строим график отклик (ось x) vs предсказания (ось y), на котором чем ближе данные к прямой линии, выходящей из левого нижнего угла, тем лучше. Еще эта визуализация удобна тем, что не зависит от количества признаков в модели;

• Хочется рассмотреть распределения предсказаний и откликов целиком, а не отдельными парами – выбираем Q-Q график;

• Если выборки очень большие, то можно снизить когнитивную нагрузку путем группировки  данных по уровням квантилей на Q-Q графике и тогда можно например, обойтись всего 100 точками;

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

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

Метрики

Дисклеймер про обозначения X и Y

В визуализациях этого раздела вы можете встретить не совсем привычные (в сравнении со смежной по теме литературой) обозначения. Так например, предсказанные значения будут обозначаться как X в то время как отклик – Y. Сделано это намеренно чтобы у читателя несмотря на привязку к конкретной теме анализа качества модели не сложилось впечатление ^[47] будто метрики можно рассчитывать только для пары “предсказания – наблюдения”. Вместо абстрактных X и Y могут быть абсолютно любые массивы – все зависит от решаемой вами задачи.

Другой причиной использовать именно эту пару является заметное отличие в написании между x и y, на графиках и анимациях легче отличить именно такие пары, чем например U и V, или классическую y и ŷ. Кроме того на анимациях ниже переменные X и Y размещены на соответствующих осях графиков.

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

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

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

Начнем с неисчерпывающего перечня: 

• Коэффициент корреляции (Correlation coefficient) [Pearson, Karl. VII. Note on regression and inheritance in the case of two parents. 1895. https://royalsocietypublishing.org/rspl/issue/58/347-352 ^[49]];

• Коэффициент детерминации R² – [Kvalseth, Tarald O. Cautionary Note about R². 1985. https://www.tandfonline.com/doi/abs/10.1080/00031305.1985.10479448 ^[50]];

• Смещение (Bias);

• Средняя абсолютная ошибка (Mean absolute error – MAE);

• Корень из средней квадратичной ошибки (Root mean square error – RMSE);

• Средняя абсолютная процентная ошибка (Mean absolute percentage error – MAPE);

• Симметричная средняя абсолютная процентная ошибка (Symmetric mean absolute percentage error – SMAPE);

• F тест (критерий Фишера, или F критерий) для проверки значима ли модель целиком;

• t критерий для проверки значимости предикторов и отклика;

• Критерий Дарбина – Уотсона для анализа остатков (Durbin-Watson test).

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

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

Вторая группа – зеленая. Эти метрики показывают ошибку в тех же величинах, что и отклик, то есть в данном случае в долларах $. Для всех трех метрик интерпретация звучит так – чем ближе значение метрики к нулю, тем лучше (Анимация 2).

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

Из Анимации 2 и Рисунка 16 можно заметить, что при увеличении расхождения между массивами X и Y, RMSE сильнее реагирует на большие ошибки, чем MAE. Это происходит потому что в RMSE ошибки возводятся в квадрат.

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

На Рисунке 17 видно, что разница, исчисляемая в единицах измерения величины между действительными и предсказанными значениями (абсолютное отклонение), не меняется в обоих случаях: она будет равна 0 в первой паре, 8 во второй и 47 для третьей. Для метрик, измеряемых в процентах, значения ошибок по понятным причинам уменьшатся поскольку абсолютные значения наблюдений растут. При этом для MAPE уменьшение при сдвиге более заметно потому что величины нормируются на значения выборки, а не на среднее суммы модулей как это реализовано в SMAPE. Эффект от различий в формулах в таком случае будет особенно заметен на малых значениях наблюдений и нивелироваться по мере удаления от нуля, что мы и видим (если хотите изучить этот момент подробнее, раскрывайте секцию “Особенности некоторых метрик”).

Особенности некоторых метрик

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

Начнем с небольшого эксперимента: используем данные из Рисунка 17, и исходный датасет наблюдения [1, 2, 3] – предсказания [1, 10, 50] будем двигать дальше и дальше от нуля путем прибавления 10 ко всем значениям на протяжении 10 итераций. И на каждом шаге будем рассчитывать три метрики: MAPE, SMAPE, MAE. Результаты можно увидеть на графике ниже: 

Как видно из Рисунка выше, чем большие значения включены в датасет, тем меньше различия между MAPE и SMAPE, и тем меньше ошибки, измеряемые в процентах. Выравнивание MAPE и SMAPЕ объясняется особенностями расчета, которые позволяют устранить эффект несимметричности MAPE особенно заметный именно на малых значениях наблюдений. MAE при этом неизменна, как и ожидалось.

Теперь разберем почему это называется “ассиметричностью”. Лучше всего показать на примере. Представим, что модель предсказала значение 110 когда наблюдение было равно 100. Тогда MAPE для данного наблюдения составит 10%. В случае же когда наблюдение равно 110, а предсказание 100, ошибка MAPE составит 9.1% несмотря на то что абсолютное отклонение в обоих случаях равно 10. Таким образом, метрика MAPE ассиметрична, потому что одинаковое абсолютное отклонение воспринимается по-разному, в зависимости от того, превышает ли прогноз истинное значение или наоборот.

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

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

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

Таблица 2. Некоторые метрики для задачи регрессии

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

Подведем промежуточный итог. Для оценки качества моделирования мы начинали с того, что составляли таблицу с предсказанными и реальными значениями (Таблица 1). Большие таблицы анализировать было бы тяжело, поэтому мы привели информацию в более удобоваримый вид при помощи графиков, то есть перешли к визуальной оценке (Рисунки 9-14). Затем упростили задачу еще немного – вместо проведения экспертной оценки по графикам мы рассчитали метрики (Рисунки 15, 16, 17 и Анимации 1-3). Но и тут есть подвох – пусть число(а) мы получили, все ещё наша ответственность – решить хорошая метрика или нет. Для Рисунка 15 использовался порог в 5% для MAPE. Однако эту эвристику никак нельзя обобщить на все линейные регрессии в мире: данные разные, (бизнес-)задачи разные и т.д. Для некоторых данных хорошей будет считаться такая модель, ошибка которой не превысила порог 7.5%, для других – 11.2%. 

F критерий

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

Если вы до этого никогда не встречались со статистическими тестами, имеет смысл начать с упрощенного определения. Статистический тест – это способ проверить, случайность ^[56] ли мы наблюдаем, или существует закономерность. Воспринимайте такой тест как черный ящик, который принимает на вход некоторые данные, и по определенным формулам рассчитывает ответ – несколько промежуточных значений (статистик и p-value) и вердикт (Рисунок 18) [Sureiman, Onchiri, et al. F-test of overall significance in regression analysis simplified. 2020. https://www.tandfonline.com/doi/full/10.1080/00031305.2016.1154108 ^[57]].

Как видно из Рисунка 18, важно перед тестированием определить порог (да, самое время поговорить про нюанс, нам и тут с порогом нужно повозиться, но тут значительно проще, – есть общепринятые наборы значений из которых можно выбрать), значение которого сильно влияет на заключительный ответ. Такой порог называется уровнем значимости. Значение 0.05 означает, что мы допускаем 5% шанс ошибочного отклонения нулевой гипотезы (гипотеза может звучать так: модель не лучше, чем наивное предсказание средним). Этот порог мы можем варьировать, например, в некоторых научных дисциплинах используют 0.01 или даже 0.001 (более строго), а в других 0.10 (менее строго). Если вам не до конца понятен “физический смысл” уровня значимости на этом этапе, – ничего страшного. В конце этого раздела есть подробное пояснение. На данный же момент будет достаточно зафиксировать следующее – статистические тесты которые мы будем обсуждать далее имеют параметр α, который мы, как исследователи-инженеры определяем под задачу. В нашем случае он равен 0.05.

Итак, статистический критерий позволяет взять данные и некоторые параметры, на основе чего рассчитывает статистики, которые используются для сравнения (больше ли статистика порогового значения или меньше). На основе этого сравнения делается вывод о значимости модели. Я рекомендую не переизобретать велосипед здесь и пользоваться статистическими пакетами для расчета критериев (отчасти поэтому я здесь не привожу формулы для расчета), – так надежнее. Что же касается вопроса какие значения сравнивать – F статистику с критическим значением F или p-value с уровнем значимости, лично я тяготею ко второму варианту в силу привычки. 

Мы можем использовать F тест для ответа на вопрос “Значима ли модель?”. Поскольку статистика – это математическая наука, предлагаю ниже формально описать две точки зрения ^[59] на получившуюся модель. Статистический тест поможет нам выбрать наиболее вероятную из этих гипотез. 

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

Альтернативная гипотеза () тогда будет звучать так: Хотя бы один коэффициент не равен нулю. Модель значима ведь она объясняет какую-то часть вариации целевой переменной.

Теперь проведем тесты на наших трех датасетах: A, B и C (Рисунок 19).

Как мы видим из Рисунка 19, во всех трех случаях p value меньше 0.05 (уровень значимости, мы его выбрали именно таким потому что это общепринятый порог по умолчанию и поскольку принятие неверной гипотезы в случае предсказания цены квартиры не влечет страшных последствий, как например в медицинских данных, предлагаю не делать его строже), значит мы отвергаем нулевую гипотезу H₀ для моделей A, B и C. После такой проверки мы можем говорить о статистической значимости моделей в целом – хотя бы один из признаков статистически связан с откликом.

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

К недостаткам такого метода оценки качества можно отнести его узкоспециализированность: F-критерий – это параметрический тест, разработанный специально для линейных моделей, то есть его не применить к, скажем, алгоритму случайного леса (это такой алгоритм в машинном обучении), в отличие от MAPE или MAE. Да и для линейных моделей этот статистический тест тоже требует выполнения стандартных предположений (см. Предположения 2-4 выше по тексту: независимость наблюдений, нормальность остатков и гомоскедастичность). Тем не менее, если вы заинтересовались темой, можете самостоятельно изучить материалы, например, по t-критерию для каждого предиктора, когда гипотезы проверяется для каждого коэффициента модели, а также критерий Дарбина – Уотсона. Или выберите любой другой критерий для изучения самостоятельно, мы здесь рассмотрели только лишь основной принцип. P.S. можете уделить особенное внимание формулам расчета статистик и той математике ^[61], которая за ними стоит.

Если не до конца понятно про уровень значимости a, пожалуйста посмотрите этот раздел

Каждый раз когда я пытался разобраться в том что значит уровень значимости, я сталкивался с невидимой стеной. Примеры посложнее приводили выкладки, которые я не понимал. Источники же “попроще” доносили – вот смотри пример, где все интуитивно понятно: 

• H₀ (нулевая гипотеза): У пациента нет рака

• Ошибка первого рода: Тест говорит “рак есть”, хотя его на самом деле нет

• Если уровень значимости α принимается равным 0.05 – В 5% случаев тест может ошибочно напугать здорового человека, сообщив, что у него рак.

• Поэтому в медицине часто выбирают низкий α (например, 0.01), чтобы свести к минимуму ложные тревоги

Но у нас то здесь данные и коэффициенты модели – всё фиксировано. Мы применяем F-критерий и получаем p-value < 0.05. Можно хоть 100 раз провести этот тест, – результат будет таким же, ведь модель одна и та же, и коэффициенты одни и те же. Вот пожалуйста – 100 раз получим подтверждение, что модель значима. Причем тут вообще порог в 5 процентов? Откуда берётся эта “вероятность”?

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

Подытожим, 0.05 p-value: даже если связи на самом деле нет, примерно в 5 случаях из 100 тест всё равно скажет “связь есть” из-за случайности в данных.

Чтобы разбавить текст, позвольте продемонстрировать это в виде анимации – мы сгенерируем 100 точек случайным образом и попытаемся из этой совокупности извлекать наборы данных по 30 объектов и строить по ним модель линейной регрессии. Повторим извлечение 20 раз. В таком случае при уровне значимости 5% мы допускаем, что 1 раз из 20 мы при помощи F критерия посчитаем что модель значима, хотя на самом деле взаимосвязи между величинами нет.

Правда, в случае когда никакой связи между x и y в действительности не было в 1 случае из 20 мы получили результат p-value при тестировании меньше 0.05. Выбери мы порог уровня значимости строже, например, 0.01 и мы бы избежали ошибки первого рода (случай когда мы отвергаем гипотезу H₀ (связи между x и y нет) и принимаем альтернативную (связь есть), хотя H₀ на самом деле верна).

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

Итогом данной обзорной главы про регрессионные метрики и F критерий обозначим:

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

• Метрики используются при оптимизации модели и потому важно понимать их свойства, например:

  • Метрики “зеленой группы” (RMSE, MAE, смещение) удобны тем что предоставляют число в исходных единицах измерения

  • Корень из средней квадратичной ошибки (RMSE) сильнее реагирует на большие ошибки и выбросы, чем средняя абсолютная ошибка (MAE)

  • “Синяя группа” (MAPE и SMAPE) – вычисляются в процентах и поэтому их часто удобно обсуждать в бизнес-контексте. При этом, когда значения отклика близки к нулю, эти метрики могут становиться нестабильными и давать вводящие в заблуждение оценки.

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

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

Неопределенность прогноза. Предсказательный интервал

На мой взгляд, интересной комбинацией визуальной оценки и формальных метрик является построение предсказательного интервала. Предсказательный интервал – это диапазон значений, внутри которого с заданной вероятностью окажется новое наблюдение. Он позволяет показать неопределённость прогноза, объединяя статистические метрики и наглядность визуализаций (Рисунок 20).

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

Вспомним. Для получения точечного прогноза в модель подставляют значения признаков (в случае парной регрессии – одного) и вычисляют предсказание. Однако прогноз крайне редко бывает абсолютно точным: как правило, присутствует случайная ошибка модели. 

При инициализации модели мы предполагаем, что делаем несерьезные (маленькие) ошибки чаще больших и при этом ошибки одинаково вероятны и в большую и в меньшую стороны. На основе этих двух предположений выстраивается статистико-вероятностный взгляд на тему инициализации модели линейной регрессии как неразрывной связки значений коэффициентов модели и распределения ошибок (Рисунок 21) [Fisher, R. A. On the Mathematical Foundations of Theoretical Statistics. 1922. https://doi.org/10.1098/rsta.1922.0009 ^[63]]. 

Как видно из Рисунка 21, разброс ошибок модели можно оценить путем расчета стандартного отклонения ошибок “сигма” (). Тут же можно говорить и о дисперсии ошибок – тоже подходящая метрика для измерения разброса, стандартное отклонение это квадратный корень из дисперсии.  Чем больше значение стандартного отклонения , тем выше неопределенность прогноза (см. секцию 2 на Рисунке 21).

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

• Шум в данных: чем больше шума, тем больше неопределенность

• Объем выборки: чем больше данных модель “видела” при обучении, тем надёжнее оценены коэффициенты и уже интервал

• Удаленность точки от центра данных: чем дальше новое значение признака от среднего, тем выше неопределённость.

Упрощенно, алгоритм построения предсказательного интервала выглядит так: 

1. Инициализируем модель (по формуле из раздела выше – Рисунок 6)

2. Считаем компоненту ошибки – остатки

3. Из остатков рассчитываем типичный размер ошибки

4. Получаем  точечный прогноз

5. Умножаем на поправки “сколько данных было в обучении”, “как далеко значения признака​ от центра данных”, “какой уровень доверия модели” (это вероятность того, что интервал накрывает значение параметра. Этот параметр мы выберем сами исходя из особенностей задачи, принцип здесь такой же как и уровня значимости при статистическом тестировании выше по тексту, классическое значение по умолчанию – 0.95)

Рассмотрим простейший пример, сгенерируем набор данных в 30 элементов с идеальной линейной зависимостью между признаком и откликом, построим модель и рассчитаем предсказательный интервал. А затем 1) добавим шум в данные, 2) увеличим объем выборки, 3) увеличим уровень доверия (от 90% до 99%, где при 99% ширина предсказательного интервала будет максимальной) – см. Анимация 4.

И отдельно рассмотрим как выглядит предсказательный интервал для датасетов A, B, C (Рисунок 22).

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

Разбиение на обучение и тест и оценка метрик 

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

Поэтому общепринятой практикой в машинном обучении является разбиение исходной выборки на части (подвыборки), на одной из которых модель “обучается” – обучающая, а на второй проверяется её обобщающая способность – тестовая (Рисунок 23).

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

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

• наконец то “ожила” метрика смещения – на тестовых данных она теперь не равна нулю, как для обучения, а изменяется и в большую (датасеты A и B) и в меньшую стороны (датасет C); 

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

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

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

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

• Сценарное тестирование. В ряде бизнес-задач важно заранее оценить, как модель поведёт себя в определённых критичных или редких ситуациях (например, при высоких или экстремальных значениях признаков). Такие случаи можно специально включить в тестовую выборку, даже если они не представлены или слабо представлены в тренировочной.

Представим что в мире всего 45 квартир…

Чтобы нам было проще идти дальше по тексту, примем важную условность для этой статьи – представим, что наш гипотетический мир в рамках которого мы строим модели очень маленький и в нем построили всего 45 квартир. Тогда все предыдущие попытки построить модель на выборках A, B и C были ничем иным как частными шагами к воспроизведению вот этой самой “первородной” зависимости из всех доступных данных. 

Тогда выборки A, B и C – это не совсем отдельные датасеты (представим, что это данные собранные из разных городов – A, B и C соответственно), а части генеральной совокупности D. Предположим, что мы можем объединить эти выборки и работать с ними как с одним целым (Рисунок 25).

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

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

В реальном мире  собрать данные по всем квартирам физически невозможно (слишком долго, дорого и трудоёмко), поэтому мы всегда берем только часть. Так и здесь мы собирали выборки и пытались оценить зависимость между переменными таким образом, чтобы быть максимально приближенными к зависимости из D. 

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

Улучшение качества модели

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

• Расширение выборки – увеличение количества объектов в выборке;

• Уменьшение выборки – исключение выбросов и других нежелательных строк из таблицы с данными;

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

• Упрощение модели – сокращение числа предикторов (иногда и это помогает улучшить метрики!);

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

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

Увеличение размеров выборки 

Помним, что значения из генеральной совокупности нам недоступны и мы можем извлекать их только частями. В рамках эксперимента мы будем формировать выборки по 10 и 20 квартир  случайным образом. Эксперименты повторим по 30 раз для каждого размера. Метрики будем замерять по 1) обучающей выборке 2) тестовой и 3) по всей совокупности значений – по 45 объектам. Попробуем оценить как увеличение объема выборок помогает построить более надежную модель для генеральной совокупности (Анимация 6).

Увеличение размеров выборки – занятие хорошее хотя бы потому, что математическая статистика лучше себя чувствует на больших числах. Потому и метрики будут получаться более устойчивыми, и статистические критерии будут более надежными (Рисунок 27).

Если вам более привычны ящики с усами (Boxplot) – раскрывайте раздел

Несмотря на то что мы здесь оперировали очень маленькими выборками (это отчасти сделано для наглядности), все же некоторые выводы, справедливые и для бОльших объемов, сделать мы можем на основе Анимации 6 и Рисунка 27. А именно:

• Среднее значение метрики ошибки RMSE на генеральной совокупности меньше в случае когда размер сэмплов составлял 20 элементов – 4088 против 4419 при 10 объектах. Это значит, что модель, построенная на большем объеме данных, имеет меньшую ошибку на генеральной совокупности.

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

Как видно, мы добились прироста в метрике на генеральной совокупности в случае когда использовали большие выборки (по 20 элементов). Так и Вы, уважаемый читатель, не забывайте измерять метрики после каких-либо изменений в данных или модели. Если изменение приводит к улучшению метрики – оставляем, наоборот – откатываем назад, – полагаемся на инженерный подход, а не на удачу! Хоть и в реальном имре замерять метрики по всем объектам у нас нет возможности, метрики на обучающих и тестовых выборках помогут выбрать правильный путь.

Уменьшение выборки (фильтрация выбросов)

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

Итак, мы попытались собрать побольше данных в выборку для обучения модели. Однако, представим, что нам не повезло: полученный семпл несмотря на свой размер в 20 элементов, не дает возможность построить модель которая была бы похожа на эталонную (Рисунок 28).

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

На этот раз попробуем решить проблему путем исключения некоторых объектов, наиболее непохожих на остальные (выбросы). Для этого существует большое количество методов, так или иначе оперирующих простым принципом отделения похожих объектов от непохожих по некоторому порогу (Рисунок 30) [Mandic-Rajcevic, et al. Methods for the Identification of Outliers and Their Influence on Exposure Assessment in Agricultural Pesticide Applicators: A Proposed Approach and Validation Using Biological Monitoring. 2019. https://doi.org/10.3390/toxics7030037 ^[74]]:

• Метод межквартильного размаха (IQR method) – непараметрический метод*

• Метод 3-сигм – параметрический* (предполагается распределение, чаще всего нормальное)

• Z-оценка (Z-score method) – параметрический*

• Модифицированная Z-оценка (на основе медианы) – параметрический*

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

В одномерных методах (Рисунок 30) не используется информация о предикторах – только об одной переменной – отклике. И потому, в частности, видно, что такие разбиения не учитывают тренд в данных. Другим ограничением является необходимость задавать порог, будь это 1.5 перед значениями разности квартилей, 3 в формуле трех сигм или порог для Z-оценки.

Следующая особенность: три из четырех приведенных методов фильтрации выбросов полагаются на предположение о форме распределения отклика. Если данные распределены нормально, или хотя бы с одной модой и без сильной асимметрии, то метод 3-сигм, Z оценка и модифицированная Z-оценка будут выдавать адекватный результат, однако если форма распределения в какой-то степени необычна, выбросы могут оказаться не выбросами вовсе. Поскольку на Рисунке 30 форма распределения не сильно отличается от нормального купола, в данном случае можно использовать именно эти стандартные методы.

Еще из интересностей здесь можно отметить, что метод трёх сигм – это частный случай Z-оценки при пороге 3.0, просто выраженный в исходных координатах y, а не в стандартизированном виде (в пространстве z-оценки). Можете обратить внимание на график: как там соотносятся линии 2σ для метода трех сигм и линии для z-оценки при пороге 2.0.

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

Рассматривая Рисунок 31 заметим, что “худшая” модель с точки зрения метрики RMSE на генеральной совокупности – именно та, которая была инициализирована на данных с выбросами. Лучшая по метрике получилась модель, которая обучалась на данных очищенных по Z оценке с порогом в 1.5.

По Рисунку 31 довольно легко делать выводы про эффективность того или иного метода исключения выбросов. Однако такое впечатление обманчиво по той причине, что мы сверяемся с метриками на генеральной совокупности D, что при реальной разработке недоступно. На закономерный вопрос “А что же тогда делать” – ответ один, экспериментировать. Иногда удобнее и быстрее всего “вычистить” тестовую выборку и замерять метрику по ней. В других случаях можно считать исключение выбросов успешным, если разница ошибок на обучении и тесте уменьшается. Серебрянной пули подходящей для всех случаев тут нет.

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

• Махаланобисово расстояние (Mahalanobis distance) [Etherington, Thomas R. Mahalanobis distances for ecological niche modelling and outlier detection: implications of sample size, error, and bias for selecting and parameterising a multivariate location and scatter method. 2021. https://doi.org/10.7717/peerj.11436 ^[77]]

• Локальный фактор выброса (LOF – Local Outlier Factor). Иногда называют – локальный уровень выброса [Breunig, Markus M., et al. LOF: Identifying Density-Based Local Outliers. 2000. https://webdocs.cs.ualberta.ca/~zaiane/pub/check/breunig.pdf ^[78]]

• Расстояние Кука (Cook’s Distance) [Cook, R. Dennis. Detection of Influential Observation in Linear Regression. 1977. https://experts.umn.edu/en/publications/detection-of-influential-observation-in-linear-regression-2 ^[79]]

• Консенсус случайной выборки (Random sample consensus – RANSAC) [Fischler, Martin A., et al. Random Sample Consensus: A Paradigm for Model Fitting with Applications to Image Analysis and Automated Cartography. 1981. – https://www.cs.ait.ac.th/~mdailey/cvreadings/Fischler-RANSAC.pdf ^[80]]

Каждый метод из Рисунка 32 можно обсуждать отдельно поскольку они уже значительно более продвинутые по сравнению с одномерными. Я же позволю здесь ограничиться визуализацией и не идти вглубь. Мы рассмотрим эти методы сугубо с практической точки зрения – посмотрим как их использование влияет на коэффициенты и метрики модели линейной регрессии (Рисунок 33).

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

Рассмотрим метод расчета рычага и Критерий Кука. Рычаг – это характеристика, показывающая “насколько объект странный по оси x”. Далеко ли расположен от центра. Если далеко, то рычаг у этого объекта будет большим. Подходящая метафора: качели – чем дальше сидишь от центра, тем сильнее влияешь на процесс качания. Критерий Кука определяет насколько выбранная точка способна изменить модель, если её убрать (зависит от величины рычага и остатка).

В примере который я привел выше для наглядности расчеты производятся итеративно. Однако в Python библиотеке sklearn например, расчет построен таким образом, что расстояние Кука определяется без реальной переоценки модели n раз.

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

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

• Выбор случайным образом подвыборки размера n объектов

• Обучение модели по выборке из выбранных n объектов

• Отсечение выбросов – исключение объектов, для которых ошибка модели превышает заданный порог

• Опциональный шаг: Обучение модели ещё раз на оставшихся не-выбросах и отсечение выбросов еще раз

• Подсчет количества не-выбросов (согласованных точек) – m

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

Результаты применения алгоритма RANSAC и метода расчета расстояния Кука показаны на Рисунке 34.

Исходя из результатов представленных на Рисунке 34, наиболее перспективной в дуэли оказалась модель построенная при помощи метода RANSAC.

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

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

Усложнение модели. Многомерная (многофакторная) линейная регрессия

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

Генерация новых признаков

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

Обратите внимание на то, что полученное уравнение теперь – модель полиномиальной регрессии, что позволяет моделировать нелинейные взаимосвязи в данных. Чем выше степень полинома, тем больше степеней свободы модель имеет (Рисунок 37).

Существует большое количество самых разнообразных преобразований, которые можно прикрутить к исходным данным. Однако в случае их применения, модель, получится не совсем линейная (что и можно заметить по форме получившихся модельных кривых на Рисунке 37). Поэтому я не буду останавливаться на них подробно в рамках этой статьи. Если стало любопытно, то можете прочитать подробнее про следующие преобразования которые можно делать над исходными признаками (и как пример литературы на которую можно опереться – Trevor Hastie, Robert Tibshirani, Jerome Friedman – The Elements of Statistical Learning ^[88]): 

• Функциональные преобразования

  • Логарифмы:

  • Обратные значения:

  • Корни:

  • Экспоненты:

  • Тригонометрические: особенно если признак имеет периодичность

  • Сигмоиды:

• Бинаризация и дискретизация

  • Биннинг: разбить признак X на интервалы, например,

  • Квантильный биннинг: разделение на равные по числу наблюдений группы

• Пороговые функции (привет, нейросети)

• Сплайны

• Вейвлет-преобразования

• Фурье-признаки

• и многие другие

Ставьте лайки – если наберется больше 500 плюсов, то обязательно выпущу отдельный пост посвященный генерации признаков, естественно, с большим количеством визуализаций XD

Сбор новых признаков

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

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

• Площадь квартиры, квадратные метры

• Удаленность от ближайшей станции метро, м

• Город

• Есть ли кондиционер

Обновленный датасет будет выглядеть следующим образом:

К вопросу визуализации 

Оглядываясь назад к Рисунку 1 (или к большинству рисунков встречавшихся в статье раньше) несложно заметить, что мы уже не охватываем все признаки двухмерной картинкой. Поэтому пора сделать новые визуализации, чтобы посмотреть на данные под новым углом (Рисунок 39 и Анимация 9). 

Рисунок лучше рассмотреть во всех деталях (Рисунок 40).

Из Анимации 9 видны две выразительных особенности набора данных:

• Чем ближе квартира к метро, тем выше стоимость; квартиры у метро обычно меньшей площади (наблюдение номер 2 из Рисунка 40);

• Наличие кондиционера – переменная, по которой заметно разделяется отклик (цена квартиры): когда в квартире установлен кондиционер, квартира обычно стоит дороже (наблюдение номер 6 из Рисунка 40).

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

Вспоминаем про Рисунок 5: почему цена все таки убывала

Вернемся к одному из первых рисунков в заметке (кликайте сюда чтобы перенестись -> figure_5 ^[91]) – там, где объяснялась концепция описания данных прямой линией, был пример из трех объектов, где цена убывала несмотря на увеличение количества комнат. Однако все становится очевидным если мы визуализируем данные с дополнительным признаком:

Можно хорошо увидеть причину снижения цены – несмотря на увеличение площади, квартиры значительно удалялась от станции метро. Пусть вас не смущает тривиальность примера, – за ним скрывается важная мысль, которая иногда пропадает из виду при анализе по настоящему больших и сложных данных: “мы не видим взаимосвязей между переменными за пределам данных, которые анализируем”. И делать выводы нужно осторожно, – всегда может открыться что-то новое когда в исходный набор данных добавится новая размерность.

Когда количество признаков увеличивается становится сложно строить такие попарные визуализации, которые были приведены на Рисунках 39 и 40. И если в вашем датасете много численных признаков, то довольно популярно отрисовывать матрицы корреляций (Рисунок 41) – их, я уверен, вы будете встречать очень часто, если продолжите интересоваться темой обработки данных.

Тут тот же принцип, что и с оценкой моделирования, – когнитивная нагрузка на инженера меньше при интерпретации чисел (по одному на пару) вместо разглядывания подграфиков. Итак, на Рисунке 41 видно, что цена положительно скоррелирована с признаками “Количество комнат” и “Площадь” и отрицательно с “Расстоянием до метро”, что логично ^[94], ведь обычно чем квартира ближе к метро и чем больше её площадь, тем она дороже. 

Отдельно хочется отметить почему часто визуализируется коэффициент корреляции, – всегда полезно смотреть есть ли в датасете сильно скоррелированные друг с другом предикторы (явление мультиколлинеарности). Именно это мы и видем для пары “Количество комнат” –  “Площадь”, где коэффициент пары равен единице. В таких случаях имеет смысл исключить один из признаков, поскольку он не несет полезной информации в модель в то время как ресурсы (например, вычислительные при подготовке данных и при оптимизации модели) забирает. Мультиколлинеарность приводит и к другим неприятным последствиям но поговорим об этом чуть позже.

Про важность предобработки (категориальных) предикторов

Как видно из Рисунка 39, теперь в таблице не только ровные числа количества комнат, но и не очень аккуратные дистанции до метро, так и совсем непонятные названия городов или ответ на вопрос “есть ли что-то в квартире” в виде текста.  

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

Начнем с признака попроще – “Есть ли кондиционер?”, поскольку он принимает только два значения – “да” или “нет”. Такие признаки принято кодировать (переводить из строкового представления) двумя числами, например (Рисунок 42): 

• да – 1

• нет – 0

Обратим внимание на то, что на Рисунке 42 не две разные модели, построенные каждая для своей подвыборки, а одна. Здесь коэффициент наклона b₁ фиксирован а сдвиг модельной линии будет отличаться для объектов когда бинарный признак равен 0 и когда равен 1. Так происходит потому что когда признак равен нулю, будет зануляться и произведение значения предиктора и коэффициента. Это работает хорошо, когда зависимость между признаками и откликом линейна и одинаково направлена для всех объектов выборки, но бинарный признак не поможет при моделировании более сложных разнонаправленных зависимостей (Рисунок 43).

Как видно из рисунка 43, в худшем случае модель с бинарным признаком сходится к модели с одним численным признаком. Чтобы эту “проблему” разрешить, можно вспомнить приемы из предыдущей главы, и сгенерировать новый признак – взаимодействие, или мы можем построить две разные модели для разных частей датасета (Рисунок 44).

Теперь когда с бинарным признаком разобрались, имеет смысл поговорить про более сложный случай, когда в столбце закодировано более двух уникальных значений. Существует большое количество методов кодировки категориальных значений, – некоторые из них показаны на Рисунке 45 [Unidata Blog, Encoding Categorical Variables: One-Hot vs. Label Encoding and Beyond. 2026. https://unidata.pro/blog/encoding-categorical-variables-one-hot-vs-label/ ^[96]], [Hancock, J. T., et al. Survey on categorical data for neural networks. 2020. https://doi.org/10.1186/s40537-020-00305-w ^[97]]. Однако я не буду приводить все из них здесь поскольку из своего личного опыта OneHotEncoding было достаточно для решения задач. Просто будьте в курсе что существуют разные механизмы кодировки.

Измерение важности признаков 

Теперь когда мы научились усложнять модель путем добавления новых признаков, имеет смысл поговорить о том, как комбинировать независимые переменные, оставлять только полезные. Безусловно, при раздувании признакового пространства путем генерации признаков или через сбор новых данных срабатывают такие ограничители как “здравый смысл” и “модель долго обучается”. Однако мы вполне можем положиться и на более эффективные эвристики, чтобы определить какие признаки имеет смысл в модели оставить. Начнем с самого простого, – рассмотрим поближе коэффициенты модели многомерной регрессии (Рисунок 46).

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

Нормализация и стандартизация предикторов

Нормализация, это трансформация данных, которая приводит значения массивов к одинаковому диапазону (Рисунок 47).

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

Какие именно формулы используются для нормализации и стандартизации, показаны на Рисунке 48.

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

Пространство коэффициентов модели когда признаки стандартизированы

После того как исходные признаки стандартизированы, а значит коэффициенты , , и т.д. “выровнены” (так удобнее варьировать коэффициенты), самое время более пристально посмотреть на то как значения этих коэффициентов влияют на ошибку модели. Для измерения ошибки воспользуемся метриками MAE и MAPE для парной регрессии и RMSE для многомерной.

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

Для наглядности, давайте посмотрим на то как по разному будут выглядеть ландшафты метрик, например для средней абсолютной процентной ошибки мы увидим следующее – поскольку MAPE чувствительна к ошибкам с небольшим значением отклика (“дешевые квартиры”), то минимум “растягивается” в вытянутую долину. И потому многие комбинации дают похожий MAPE, если линия хорошо попадает в область маленьких y, даже если на дорогих квартирах она заметно ошибается. (Анимация 11).

Далее увеличим количество признаков в модели, чтобы было необходимо найти оптимальное сочетание уже не двух, а трех коэффициентов (Анимации 12 и 13). 

На анимациях выше видно, что признаки имеют сильную линейную связь, например в Анимации 12, на проекции коэффициентов b₁ vs b₂ (плоскость слева на нижнем левом графике) виден линейный паттерн. Это говорит о том что 1) присутствует сильная обратная корреляция между признаками “количество комнат” и “расстояние до метро” и 2) несмотря на то что коэффициенты “плавают вдоль долины низких значений RMSE”, предсказания модели стабильны, – ошибка модели практически не изменяется. Значит, в том числе и о то что признаки дают похожую информацию. Такая же ситуация и с Анимацией 13, но здесь линейная связь между признаками еще сильнее, и уже прямая, а не обратная.

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

Обобщение аналитического решения на многомерный случай

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

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

Когда мы рассматривали формулу аналитического решения ранее, мы приводили расчеты в скалярном виде. Однако намного эффективнее перейти к векторной записи, для этого визуализируем исходные данные не в пространстве признаков, а в пространстве наблюдений (Рисунок 49).

Несмотря на то что поначалу такой взгляд на данные может казаться контринтуитивным, за ним не стоит никакой магии, данные все те же, изменилась лишь форма. А мы продолжаем, – из школьной программы нас учили (по крайней мере меня) называть вектором направленные отрезки. Эти “направленные отрезки” можно умножать на число и складывать друг с другом. Задача алгоритма линейной регрессии в векторном пространстве – найти такое преобразование вектора , чтобы получившийся вектор предсказаний (обозначется как с крышечкой наверху) был как можно больше похож на вектор откликов y. Давайте попробуем перебрать такие преобразования, начиная с умножения на число (Рисунок 50). 

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

Если из Рисунка 50 хорошо видно как производится умножение вектора на число, то вот на Рисунке 51 были замечены еще две операции над вектором, напомнить про которые имеет смысл отдельно (Рисунок 52).

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

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

Та область графика, в которую попадают возможные решения, можно очертить и получим плоскость – на анимации выше она для наглядности ограничена параллелограммом. Будем называть эту плоскость подпространством предсказаний и обозначать как . Как видно из Анимации 15, вектор откликов y не лежит в подпространстве решений: это значит что какое бы решение (вектор предсказаний) мы не нашли, он всегда будет немного отличаться от отклика. Наша задача – найти такой вектор предсказаний, чтобы он лежал как можно ближе к вектору и при этом принадлежал подпространству . 

Мы собрали это подпространство в визуализации выше путем сложения векторов и с разными коэффициентами перед ними. Однако существует и другая запись того же самого выражения, но более компактная – через матричное умножение. Для этого составим еще один вектор, но на этот раз, из коэффициентов и . Обозначим его как вектор B. Его можно преобразовать путем умножения на матрицу, чтобы повернуть вектор, растянуть или сжать его, а также отобразить в другое (под)пространство. И если взять матрицу составленную из векторов столбцов и и умножить на вектор , составленный из коэффициентов, то мы получим отображение y на подпространство предсказаний (Рисунок 53).

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

Присмотримся к картинке и попробуем найти отличия между векторами предсказаний разных моделей A, B и C: вектор модели C выглядит как тень вектора отклика на плоскости. Тогда процесс решения задачи линейной регрессии можно интерпретировать как проецирование вектора y на подпространство . Лучший прогноз из всех возможных – вектор заканчивающийся в ближайшей к отклику точка на плоскости. А ближайшая точка на плоскости, как можно вспомнить из базовых геометрических знаний, получается в основании перпендикуляра. Получившаяся вертикаль, это тоже вектор и его называет вектором остатков e, потому что он получается путем вычитания предсказаний из отклика (вспоминаем формулу остатков из главы про оценку качества модели с помощью визуализаций).

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

Название аналитического метода который мы только что рассмотрели – метод наименьших квадратов (МНК, или Ordinary Least Squares, OLS, по английски). Название такое, потому что мы подбирали такие коэффициенты, чтобы минимизировать сумму квадратов остатков модели (Рисунок 6, кликайте чтобы перенестись -> figure_6 ^[25]). В векторном пространстве размер остатков – это квадрат евклидова расстояния от точки отлика до подпространства (Рисунок 55). Наименьшие квадраты равно “наименьшее расстояние в квадрате”.

Вспомним цель текущего раздела – мы продирались через формулы и визуализации выше чтобы обобщить аналитическое решение на многомерный случай. И теперь самое время проверить как будет работать формула для случая когда признаков не один, а два! Рассмотрим на примере датасета с тремя объектами в который добавим один признак (Анимация 16).

На анимации 16 прокомментируем три важных наблюдения:

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

2. Смотрим направо – количество векторов не изменилось, поскольку в выборке как и прежде, три объекта

3. Теперь подпространство не плоскость на графике, а все пространство. Для наглядности значения были ограничены объемной фигурой, – параллелепипедом. И поскольку это подпространство полностью покрывает собой вектор отклик y, проекция отклика будет тривиальна. На Анимации это выглядит так: вектор отклика и вектор предсказаний совпадают. Остаток нулевой.

Когда аналитическое решение встречает трудности 

Представим, что нам не повезло, и новый признак новой информации не привносит. Пусть новый признак линейно выражается через два другие (сдвиг, и признак ) и тогда многоугольник опять сжимается до состояния плоскости – Анимация 17.

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

При таких исходных данных в расчетах проблема возникнет на этапе расчета обратной матрицы (Рисунок 56).

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

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

1. Используется не-квадратичная или негладкая функция потерь (L1, квантильная и т.д.) – тогда задача уже не сводится к методу наименьших квадратов;

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

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

Численные решения

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

Оптимизация полным перебором

Самый простой способ, – исследовать все комбинации коэффициентов с некоторым шагом. Метод полного перебора (exhaustive search) в данном случае – это перебор всех пар коэффициентов из заранее заданной дискретной сетки значений и выбор пары с минимальной ошибкой. Обычно выбирают метрику MSE (как RMSE но без квадратного корня) для измерения ошибки.

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

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

Оптимизация случайным поиском

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

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

Использование информации о направлении

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

Для наглядности упростим задачу – возьмем срез по одному параметру, одномерный массив, и рассмотрим на его примере. При движении по массиву воспользуемся тем, что мы на предыдущем шаге уже рассчитали значение ошибки (возьмем MSE для данного примера) и сравнивая текущее и предыдущее сделаем вывод куда имеет смысл сделать следующий шаг – Рисунок 57.

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

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

Явно отметитим, что в Анимации 20 шаг всегда равен одному интервалу (шагу сетки) и мы пока не используем никаких производных (предвкушая алгоритм градиентного спуска) – просто попарно сравниваем метрики.  

У описанного выше подхода есть один существенный недостаток, – он сильно зависит от размеров сетки. Так например, если сетка будет подробная, то придется сделать очень много шагов в направлении оптимума. И наоборот, если сетка слишком грубая – пропустим оптимум (Анимация 21).

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

Оптимизация методом градиентного спуска 

При уменьшении шага сетки при попарных сравнениях мы приходим к определению производной через предел (Рисунок 58).

Предлагаю ниже посерфить по ландшафту ошибки – рассмотрим анимацию где показаны величины градиента и антиградиента (Анимация 22). Как мы видим, шаг мы теперь можем выбирать произвольный – регулярной сеткой теперь не ограничены [Goh, Gabriel. Why Momentum Really Works. 2017. https://distill.pub/2017/momentum/ ^[114]].

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

Самое время посмотреть на градиентный спуск в действии (Анимация 23).

Посмотрите как сходится градиентный спуск при разных значениях темпа обучения (learning rate)

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

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

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

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

Регуляризация 

Вспомним эффект из Анимации 5 ^[121], где разные обучающие выборки приводили к разным оцененным коэффициентам несмотря на то что мы пытались восстановить одну и ту же исходную зависимость между признаком и откликом. Модель получалась нестабильна – чувствительна к разбиению. 

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

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

• Как регуляризация помогает при неоднородных разбиениях

• Как регуляризация помогает модели хорошо обобщать данные

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

Один из способов обеспечения численной стабилизации, это наложение ограничений на коэффициенты – регуляризация (Рисунок 63).

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

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

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

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

Переобучение

И раз мы заговорили о регуляризации, нельзя не упомянуть что степень регуляризации можно варьировать (Анимация 26).

Из Анимации 26 видно:

• Линия 1: Показаны значения коэффициентов перед признаками, метрики на обучении и тесте, а также график сравнения предсказаний и фактических значения для модели без регуляризации;

• Линия 2: Как себя ведет Lasso регрессия при разной силе регуляризации, ошибка сначала уменьшается на тестовых данных, но затем модель деградирует к предсказанию средним из-за чрезмерной регуляризации и значения коэффициентов перед признаками зануляются;

• Линия 3: Гребневая регрессия при увеличении силы регуляризации показывает всё меньшие и меньшие значения ошибки на тестовой выборке несмотря на то, что ошибка на обучении постепенно увеличивается.

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

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

Настройка гиперпараметров

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

Однако так мы рискуем настроить модель под конкретную тестовую выборку, поэтому в машинном обучении принято использовать кросс валидацию (Рисунок 68).

Линейная регрессия – это целый мир

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

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

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

Заключение

Комплимент дня Спасибо всем тем кто дочитал до этого момента.

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

Если вам понравились визуализации и примеры, и вы хотите использовать их в своих лекциях или выступлениях, пожалуйста делайте это. Все материалы и исходный код для их генерации есть в репозитории на GitHub – https://github.com/Dreamlone/linear-regression ^[127] 

Искренне Ваш, Михаил Сарафанов ^[128]