Используйте исторические данные о бейсболе и логистические регрессии с R, чтобы предсказать количество ранов, забитых в играх Мировой серии.

Введение

Вы когда-нибудь задумывались, что общего между американским бейсболом, машинным обучением и статистикой? Что ж, тебе сегодня повезло!

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

Есть несколько предположений относительно этого проекта. Во-первых, у вас есть настройка среды R. Я собираюсь использовать файл R Markdown (RMD) в RStudio на Mac для этого проекта, но код будет отлично работать в других операционных системах и средах. Во-вторых, у вас есть хотя бы элементарное представление о бейсболе. Я ни в коем случае не специалист по бейсболу, но особых знаний не потребуется. Пока вы знаете, что такое хоумран, у нас все будет хорошо!

Для вашего удобства файл RMD доступен на моем GitHub здесь.

Итак, давайте перейдем к целям и приступим к написанию кода!

Цели

  1. Узнайте об импорте данных без использования локальных файлов или загрузки данных в среду вручную.
  2. Разделите данные на наборы для обучения и тестирования для целей машинного обучения.
  3. Используйте ненормальные распределения данных для машинного обучения
  4. Создание моделей логистической регрессии с использованием метода пошаговой регрессии
  5. Используйте модель, чтобы делать прогнозы и проверять точность нашей модели.

Данные

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

Вот сайт Шона: http://www.seanlahman.com/baseball-archive/statistics/

Теперь пришло время кодировать!

Код

Загрузить библиотеки

В этом проекте мы будем использовать пакеты tidyverse и Lahman.

Вот код:

# Uncomment the commands below if you have not installed either the tidyverse or Lahman packages
# install.packages("tidyverse")
# install.packages("Lahman")
# Load Libraries
require(tidyverse)
require(Lahman)

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

Вот код:

# Imports the BattingPost data set from the Lahman database. Notice how there is not a file to read in. The data is part of the package!
data(BattingPost)
# Check the data import
# Notice again that we did not have create and object and assign it data like we normally do with R. It just knows to create a data frame called "BattingPost"
head(BattingPost)

Вот результат:

Фильтровать данные

Теперь у нас есть данные в R. Нам нужно немного почистить его. Во-первых, поскольку это проект о Мировой серии, нам нужно убедиться, что в нашем столбце round есть только значения WS. Во-вторых, нам нужно убедиться, что столбец yearID больше или равен 1920 году. Я не историк бейсбола, но именно тогда началась современная эра живых мячей [2]. Если вам интересна такая история, взгляните на статью Брэдфорда Дулиттла на ESPN об этом здесь [2].

Вот код:

# Filter down the data a bit. We went the year to be 1920 and later because that's when the "Live Ball Era" starts. We also want only WS in the round column because that's for the World Series.
SlimBattingPost <- BattingPost %>%
  filter(yearID >= 1920) %>%
  filter(round == "WS")
# Check the data
head(SlimBattingPost)

Вот результат:

Давайте немного проясним это для тех, кто следит за нами. Пара %›%, которую вы видите пару раз, — это то, что известно как символ трубы [3]. Когда я использую tidyverse, в частности его часть dplyr, для работы с данными, мне нравится думать об этом как о том, что каждая функция или условие, следующие за этим символом вертикальной черты, делают вашу гипотетическую воронку данных немного меньше.

Нам редко нужны все данные в наборе данных, поэтому «заливать» их через «воронку» все меньшего и меньшего диаметра — хороший способ визуализировать то, что здесь происходит. Просто применив эти два фильтра, мы увеличили количество строк с 14 750 в исходных данных до 4 280 строк в наших новых. Вы можете получить эти числа, просмотрев свои объекты в RStudio или используя функцию nrow().

Вот код:

# Print number of rows in each data set
nrow(BattingPost)
nrow(SlimBattingPost)

Введение в машинное обучение

Давайте демистифицируем, что означает машинное обучение на простом английском языке, используя наш пример. То, что мы собираемся сделать, известно как «контролируемое машинное обучение», потому что мы, люди, «контролируем», какие данные передаются в машину [3].

При использовании подобных реальных данных мы хотим разделить некоторую часть всего набора данных на «обучающие» данные, а остаток — на «тестовые» данные [3]. Это сделано по вполне рациональной причине. Если вы хотите, чтобы машина могла делать прогнозы, вам нужно проверить ее на практике, чтобы убедиться, что она действительно работает правильно [3].

Традиционно и по уважительной причине я могу добавить, что разделение составляет 80% данных обучения и 20% данных тестирования [3]. Я видел различные другие комбинации, такие как 70/30 и 75/25, но мне нравится разделение 80/20, потому что я обычно получаю наилучшие результаты, и это звучит так же, как принцип Парето [4]. Давай, я не могу с собой поделать. Моя степень бакалавра в области экономики.

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

Установите семя

Прежде чем мы разделим данные, мы действительно должны установить начальное число [3]. Мы собираемся использовать случайное (ну, если быть точным, псевдослучайное) разбиение. Это означает, что код будет разбивать данные на разные наборы каждый раз, если мы не используем функцию set.seed() с числом, чтобы гарантировать воспроизводимость [3]. Для тех, кто следит за нами, это означает, что ваш код будет соответствовать моему, если мы используем тот же номер.

Вот код:

# Set Seed
set.seed(1337)

Установка начального значения не имеет никакого вывода. Это просто то, что обеспечивает правильную установку скрытых параметров [3].

Данные обучения и тестирования

Есть масса способов сделать это, но это самый простой способ, который я когда-либо придумывал или видел на практике. Вот секрет. Мы собираемся создать столбец с произвольным идентификационным номером и добавить его к имеющимся у нас данным. Это позволяет нам легко использовать одну из самых хитрых функций в R — функцию anti_join() из пакета dplyr внутри tidyverse. Что это делает, так это получает все строки из одной таблицы, которые не имеют соответствия в другой таблице. Это подло. Мне это нравится, и вы можете прочитать все об этом в документации здесь.

Вот код:

# Creates an ID column so we can more easily sort train from test
SlimBattingPost$id <- 1:nrow(SlimBattingPost)
# Creates set of data randomly sampling 80% of total data
train <- SlimBattingPost %>% dplyr::sample_frac(.80)
# Creates set of data with other 20%
test <- dplyr::anti_join(SlimBattingPost, train, by = 'id')
# Check the data
head(test)
paste("The test data has this many rows:", nrow(test))
head(train)
paste("The train data has this many rows:",nrow(train))

Вот результат:

Определить распределение

Большую часть времени машинное обучение имеет дело с нормально распределенными данными [3]. Когда вы думаете о нормальной кривой нормального распределения, которую мы все видели много раз, она называется «гауссовым» распределением, потому что так сказали математики [3]. Это выглядит так:

В нормальном распределении есть среднее (мю, синяя линия) со стандартными отклонениями (сигма, зеленые линии) от среднего. Не очень важно вдаваться в подробности этого, но это означает, что ~ 95% данных находятся в пределах 2 стандартных отклонений (2 сигма) от среднего значения в любом случае [3]. Есть множество причин, по которым машинному обучению нужны нормальные данные, но иногда мы получаем данные, которые не являются нормальными.

Другое распределение известно как распределение Пуассона [3]. Конечно, оно названо в честь другого математика, но дело в том, что это стандартный тип распределения, в котором большая часть данных сильно искажена [3]. Чаще всего это выглядит примерно так:

Для машинного обучения нас действительно заботит распределение переменной, которую мы пытаемся предсказать [3]. В нашем случае запуски (R) будут зависимой переменной. Давайте сделаем очень быструю гистограмму, чтобы проверить, имеем ли мы нормальное распределение или распределение Пуассона.

Вот код:

# Visually determine type of distribution of Runs (R) in dataset
hist(SlimBattingPost$R)
# Classic Poisson distribution. Almost all data is 0 or 1, heavily skewed

Вот результат:

Переменная run (R) выглядит как довольно классическое распределение Пуассона!

Построить модель

А вот и самое интересное! Теперь мы можем построить модель логистической регрессии! Давайте посмотрим на код, а затем поговорим о нем.

Вот код:

# Start with everything else as independent variables with runs (R) as the dependent variable
# For a data dictionary about what each abbreviation means, go to: http://www.seanlahman.com/files/database/readme2017.txt
# Search for "BattingPost" in that document and the third match should be the table where the abbreviations are defined. G = Games, AB = At Bats, etc.
# Create the logistic regression with everything in it, make sure to include the Poisson distribution as the family. If we had normal data, we would use Gaussian in its place.
fitAll <- glm(R ~ G + AB + H + X2B + X3B + HR + RBI + SB + CS + BB + SO + IBB + HBP + SH + SF + GIDP , data = train, family = "poisson")
summary(fitAll)

Вот результат:

Здесь есть что распаковать.

Во-первых, мы создали объект модели под названием «fitAll», который использует функцию glm() с R в качестве зависимой переменной, отделенной символом ~ от всех независимых переменных, разделенных знаком +, а также используя наше обучение данных (помните 80% из предыдущих?) и использование распределения Пуассона для получения более точных результатов [3].

Как только мы создадим модель, функция summary() выдаст нам указанные выше выходные данные в простой форме.

Первая часть сводки — «Вызов:», в которой просто указано, какая модель использовалась на самом деле. Здесь не о чем говорить.

«Остатки отклонения» пытаются рассказать нам о разнице между реальностью и ожиданиями при работе с ошибками [3]. На самом деле, это не все, что имеет отношение к изучению логистической регрессии на данный момент, поэтому мы продолжим двигаться вперед.

Раздел «Коэффициенты:» — это то, на чем мы должны сосредоточиться. Первый столбец — это все переменные плюс точка пересечения. Помните формулу по алгебре I класса «y = mx + b»? Это своего рода экстремальная версия значения b для перехвата [3].

Столбец «Оценка» сообщает нам, является ли влияние каждого значения положительным или отрицательным на количество засчитываемых запусков, а также насколько велик его коэффициент [3]. Это что-то вроде группы различных значений m, где значение x является числом в фактическом наборе данных. Уравнение может выглядеть примерно так: y = m₁x + m₂x + m₃x …. + b со всеми различными коэффициентами и значениями [3].

«Стд. Столбец «Ошибка» — это статистическая мера того, насколько мы ошибаемся [3]. Значение z используется в статистике для получения z-показателя на кривой распределения [3]. Все это не так актуально для нас сейчас.

Что наиболее важно для нас, так это столбец, элегантно помеченный как «Pr(›|z|)», который представляет собой действительно сложный способ сказать p-значение [3]. Для p-значений нас интересует определенный уровень уверенности в том, что то, что мы говорим, значимо [3]. Кроме того, p-значения находятся в диапазоне от 0 до 1, поэтому, если есть отрицательные числа или числа больше 1, математика где-то неверна [3].

Например, если мы должны быть уверены как минимум на 95% в том, что независимая переменная является статистически значимым предиктором зависимой переменной, значение p будет меньше 0,05 [3]. Если мы хотим быть на 90% уверены, что переменная значима, нам нужно увидеть значение p меньше 0,10 [3]. Чтобы быть уверенным на 99%, нам нужно увидеть значение p меньше 0,01 [3]. Видишь, к чему я клоню?

R выполняет для нас небольшую удобную функцию, помечая, на каком уровне значимости находится каждая переменная. Ставим один «.» рядом с переменной, как с переменными SB и SF, означает, что мы на 90–95 % уверены, что эти переменные статистически значимы [3]. «*» рядом с переменной означает, что мы на 95–99% уверены, что переменная является статистически значимой, а остальные сдвигаются вверх с точки зрения того, что они являются сильными предикторами, чем меньше получаются p-значения [3].

Цель состоит в том, чтобы в конечном итоге получить модель, которая имеет только статистически значимые независимые переменные. Итак, как мы это делаем? Мы используем методологию, называемую «ступенчатой ​​регрессией» [3].

Пошаговая регрессия

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

Вот код:

# Start step-wise regression.
# Looks like G is the least significant
# Make sure to create a new object name so we do not overwrite our models as we go along!
fit1 <- glm(R ~ AB + H + X2B + X3B + HR + RBI + SB + CS + BB + SO + IBB + HBP + SH + SF + GIDP , data = train, family = "poisson")
summary(fit1)

Вот результат:

Итак, что вы заметили? После того, как мы убрали переменную G как наименее значимую, многие p-значения изменились! Это нормально. Когда все числа обрабатываются, то, что считается статистически значимыми изменениями, поскольку переменные удаляются [3]. Это тоже имеет смысл. Когда вы думаете об этом, не имеет смысла то, что игра 1, 3 или 6 действительно влияет на количество забитых пробежек.

Примечание. Когда вы удаляете переменную, не забудьте также удалить «+» со стороны независимой переменной в формуле, чтобы избежать ненужных ошибок! Также не забудьте переименовать свои модели в другие имена объектов, чтобы избежать перезаписи прошлой работы во время обучения!

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

Вот код:

# Looks like SO is the least significant
# Make sure to create a new object name so we do not overwrite our models as we go along!
fit2 <- glm(R ~ AB + H + X2B + X3B + HR + RBI + SB + CS + BB + IBB + HBP + SH + SF + GIDP , data = train, family = "poisson")
summary(fit2)

Вот результат:

Надеюсь, на этом этапе вы уловили шаблон. Пока в модели не останутся только статистически значимые переменные до 95% доверительного порога (p-значение ‹ 0,05), мы будем наблюдать за этим процессом.

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

Вот код:

# Final Fit
fitFinal <- glm(R ~ AB + H + X2B + X3B + HR + CS + BB + IBB + HBP + GIDP , data = train, family = "poisson")
summary(fitFinal)

Вот результат:

Делать предсказания

Настало время, когда мы видим, насколько хорошо мы справились с нашей моделью. К счастью для нас, в R есть для этого простой в использовании инструмент — функция predict(). Теперь мы будем использовать нашу окончательную модель логистической регрессии, протестируем ее на наших тестовых данных (20%, которые мы зарезервировали ранее), а затем посмотрим, насколько мы точны!

Вот код:

# Create predictions
predictions <- predict(fitFinal, test, type = 'response')
# Check the output
head(predictions)

Вот результат:

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

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

Вот код:

# Add predictions to test data and create new data frame
predictDF <- data.frame(test, predictions)
# Create new data frame with less columns
SlimPredictDF <- select(predictDF, "yearID", "round", "playerID", "teamID", "R", "predictions")
# Add rounded predictions as a column
SlimPredictDF$roundedPredictions <- round(SlimPredictDF$predictions, 0)
# Create Boolean Column to see if real and predictions match
SlimPredictDF$TFmatch <- SlimPredictDF$R == SlimPredictDF$roundedPredictions
# Check data structure 
head(SlimPredictDF)

Вот результат:

Результаты

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

Вот код:

# Get the results!
results_table <- SlimPredictDF %>%
  group_by(TFmatch) %>%
  summarise(count = n())
results_table

Вот результат:

Что ж, ~66% — это не так уж и плохо, не так ли? Простое быстрое деление 564 на общее число 856 приводит нас к этому результату.

Результаты — бонусные баллы

Другой способ проверить наши результаты — создать быструю линейную модель с помощью функции lm() и посмотреть на наши значения R-квадрата. Мы получим два варианта R-квадрата, но они будут близкими. Это скажет нам, какая часть нашей зависимой переменной (R, run) объясняется нашей независимой переменной (roundedPredictions).

Вот код:

# Simple linear model to get p-vale for whether real Runs (R) are significantly prediction by predictions
fitLM <- lm(R ~ roundedPredictions, data = SlimPredictDF)
summary(fitLM)

Вот результат:

В итоге мы получили аналогичный ответ: при проверке таким образом мы были правы примерно на 63%.

Вывод

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

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

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

использованная литература

[1] С. Лахман, Загрузить базу данных бейсбола Лахмана (2020 г.), http://www.seanlahman.com/baseball-archive/statistics/

[2] Б. Дулиттл, Сто лет в эру живых мячей, лучше ли играть в бейсбол? (2019), https://www.espn.com/mlb/story/_/id /27546614/столетняя-живая-эра-бейсбола-лучшей-игры

[3] Р. Кабакофф, R в действии (2-е изд.) (2015 г.), Шелтер-Айленд, Нью-Йорк: Manning Publications Co.

[4]. К. Крузе, Правило 80/20 и как оно может изменить вашу жизнь (2016), https://www.forbes.com/sites/kevinkruse/2016/03/07/80- 20-правило/#230185973814