Regresión con datos de NFL

Dataset

Utilizaremos los datos de #Tidytuesday, la cual es una iniciativa de la comunidad de R, en la cual se publica semanalmente un set de datos con el objetivo de practicar las habilidades de procesamiento, visualización y modelado de datos. El dataset elegido contiene datos de la NFL, liga de fútbol americano de USA.

Lectura de datos

attendance %>% glimpse()

## Rows: 10,846
## Columns: 8
## $ team              <chr> "Arizona", "Arizona", "Arizona", "Arizona", "Ariz...
## $ team_name         <chr> "Cardinals", "Cardinals", "Cardinals", "Cardinals...
## $ year              <dbl> 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2...
## $ total             <dbl> 893926, 893926, 893926, 893926, 893926, 893926, 8...
## $ home              <dbl> 387475, 387475, 387475, 387475, 387475, 387475, 3...
## $ away              <dbl> 506451, 506451, 506451, 506451, 506451, 506451, 5...
## $ week              <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15...
## $ weekly_attendance <dbl> 77434, 66009, NA, 71801, 66985, 44296, 38293, 629...

standings %>%  glimpse()

## Rows: 638
## Columns: 15
## $ team                 <chr> "Miami", "Indianapolis", "New York", "Buffalo"...
## $ team_name            <chr> "Dolphins", "Colts", "Jets", "Bills", "Patriot...
## $ year                 <dbl> 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000...
## $ wins                 <dbl> 11, 10, 9, 8, 5, 13, 12, 9, 7, 4, 3, 12, 11, 7...
## $ loss                 <dbl> 5, 6, 7, 8, 11, 3, 4, 7, 9, 12, 13, 4, 5, 9, 1...
## $ points_for           <dbl> 323, 429, 321, 315, 276, 346, 333, 321, 367, 1...
## $ points_against       <dbl> 226, 326, 321, 350, 338, 191, 165, 255, 327, 3...
## $ points_differential  <dbl> 97, 103, 0, -35, -62, 155, 168, 66, 40, -174, ...
## $ margin_of_victory    <dbl> 6.1, 6.4, 0.0, -2.2, -3.9, 9.7, 10.5, 4.1, 2.5...
## $ strength_of_schedule <dbl> 1.0, 1.5, 3.5, 2.2, 1.4, -1.3, -2.5, -0.2, -1....
## $ simple_rating        <dbl> 7.1, 7.9, 3.5, 0.0, -2.5, 8.3, 8.0, 3.9, 1.1, ...
## $ offensive_ranking    <dbl> 0.0, 7.1, 1.4, 0.5, -2.7, 1.5, 0.0, 0.6, 3.2, ...
## $ defensive_ranking    <dbl> 7.1, 0.8, 2.2, -0.5, 0.2, 6.8, 8.0, 3.3, -2.1,...
## $ playoffs             <chr> "Playoffs", "Playoffs", "No Playoffs", "No Pla...
## $ sb_winner            <chr> "No Superbowl", "No Superbowl", "No Superbowl"...

Unión de los datasets

Ya leímos el dataset attendance y standings:

attendance: datos de la asistencia semanal de los equipos.
standings: datos de la clasificación de los equipos con sus respectivas marcas.

A continuación vamos a unir ambos datasets utilizando la sentencia left_join y mediante las siguientes variables: c(“year”, “team_name”, “team”))

Objetivo

Nuestro objetivo es obtener un modelo capaz predecir la asistencia semanal o ‘weekly_attendance’ a la NFL a partir del conjunto de datos #TidyTuesday. Hasta ahora, sabemos que es necesario utilizar modelos de regresión.

Análisis exploratorio inicial

Tengamos en cuenta que para los 32 equipos de la NFL, existen años en que lo dieron todo y aún así no llegaron a las eliminatorias o playoffs, lo que será bueno para modelar.

¿Cuánto influye el ´margin_of_victory´, una medida de los puntos anotados en relación con los puntos pérdidos, para llegar a los playoffs o eliminatorias?

attendance_joined %>%
  distinct(team_name, year, margin_of_victory, playoffs) %>%
  ggplot(aes(margin_of_victory, fill = playoffs)) +
  geom_histogram(position = "identity", alpha = 0.7) +
  labs(
    x = "Margen de victoria",
    y = "Númmero de equipos",
    fill = NULL
  )

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

¿Hay cambios en cada semana de la temporada?

attendance_joined %>%
  mutate(week = factor(week)) %>%
  ggplot(aes(week, weekly_attendance, fill = week)) +
  geom_boxplot(show.legend = FALSE, outlier.alpha = 0.5) +
  labs(
    x = "Semana de la temporada de la NFL",
    y = "Asistencia semanal de juegos de la NFL"
  )

Hemos realizado un breve análisis exploratorio de datos, el cual siempre es una parte importante de la tarea de modelado. El siguiente paso es generar un conjunto de datos para modelar.

Eliminemos las semanas que cada equipo no jugó (es decir, donde la asistencia semanal es NA). Conservemos únicamente las columnas que queremos usar para el modelado. Por ejemplo, mantendremos margin_of_victory y Strength_of_schedule, pero no simple_rating, que es la suma de esas dos primeras cantidades.

attendance_df <- attendance_joined %>%
  filter(!is.na(weekly_attendance)) %>%
  select(
    weekly_attendance, team_name, year, week,
    margin_of_victory, strength_of_schedule, playoffs
  )
attendance_df

## # A tibble: 10,208 x 7
##    weekly_attendan~ team_name  year  week margin_of_victo~ strength_of_sch~
##               <dbl> <chr>     <dbl> <dbl>            <dbl>            <dbl>
##  1            77434 Cardinals  2000     1            -14.6             -0.7
##  2            66009 Cardinals  2000     2            -14.6             -0.7
##  3            71801 Cardinals  2000     4            -14.6             -0.7
##  4            66985 Cardinals  2000     5            -14.6             -0.7
##  5            44296 Cardinals  2000     6            -14.6             -0.7
##  6            38293 Cardinals  2000     7            -14.6             -0.7
##  7            62981 Cardinals  2000     8            -14.6             -0.7
##  8            35286 Cardinals  2000     9            -14.6             -0.7
##  9            52244 Cardinals  2000    10            -14.6             -0.7
## 10            64223 Cardinals  2000    11            -14.6             -0.7
## # ... with 10,198 more rows, and 1 more variable: playoffs <chr>

attendance_df_new <- attendance_df %>% drop_na()
attendance_df_new

## # A tibble: 10,208 x 7
##    weekly_attendan~ team_name  year  week margin_of_victo~ strength_of_sch~
##               <dbl> <chr>     <dbl> <dbl>            <dbl>            <dbl>
##  1            77434 Cardinals  2000     1            -14.6             -0.7
##  2            66009 Cardinals  2000     2            -14.6             -0.7
##  3            71801 Cardinals  2000     4            -14.6             -0.7
##  4            66985 Cardinals  2000     5            -14.6             -0.7
##  5            44296 Cardinals  2000     6            -14.6             -0.7
##  6            38293 Cardinals  2000     7            -14.6             -0.7
##  7            62981 Cardinals  2000     8            -14.6             -0.7
##  8            35286 Cardinals  2000     9            -14.6             -0.7
##  9            52244 Cardinals  2000    10            -14.6             -0.7
## 10            64223 Cardinals  2000    11            -14.6             -0.7
## # ... with 10,198 more rows, and 1 more variable: playoffs <chr>

Modelos simples con tidymodels

¡Ahora es el momento de cargar el metapaquete tidymodels! 💪 El primer paso aquí es dividir nuestros datos en pruebas de entrenamiento y pruebas. Podemos usar initial_split () para crear estos conjuntos de datos, divididos para que cada uno tenga aproximadamente la misma cantidad de ejemplos de equipos que pasaron a los playoffs.

library(tidymodels)  #rsample

set.seed(1234)
attendance_split <- attendance_df_new %>%
  initial_split(strata = playoffs)

nfl_train <- training(attendance_split)
nfl_test <- testing(attendance_split)

nfl_test

## # A tibble: 2,552 x 7
##    weekly_attendan~ team_name  year  week margin_of_victo~ strength_of_sch~
##               <dbl> <chr>     <dbl> <dbl>            <dbl>            <dbl>
##  1            65356 Cardinals  2000    12            -14.6             -0.7
##  2            50289 Cardinals  2000    14            -14.6             -0.7
##  3            37452 Cardinals  2000    16            -14.6             -0.7
##  4            65711 Cardinals  2000    17            -14.6             -0.7
##  5            73025 Falcons    2000     3            -10.1              1.5
##  6            66019 Falcons    2000     7            -10.1              1.5
##  7            44680 Falcons    2000    14            -10.1              1.5
##  8            41017 Falcons    2000    17            -10.1              1.5
##  9            68481 Ravens     2000     4             10.5             -2.5
## 10            83252 Ravens     2000     7             10.5             -2.5
## # ... with 2,542 more rows, and 1 more variable: playoffs <chr>

nfl_train

## # A tibble: 7,656 x 7
##    weekly_attendan~ team_name  year  week margin_of_victo~ strength_of_sch~
##               <dbl> <chr>     <dbl> <dbl>            <dbl>            <dbl>
##  1            77434 Cardinals  2000     1            -14.6             -0.7
##  2            66009 Cardinals  2000     2            -14.6             -0.7
##  3            71801 Cardinals  2000     4            -14.6             -0.7
##  4            66985 Cardinals  2000     5            -14.6             -0.7
##  5            44296 Cardinals  2000     6            -14.6             -0.7
##  6            38293 Cardinals  2000     7            -14.6             -0.7
##  7            62981 Cardinals  2000     8            -14.6             -0.7
##  8            35286 Cardinals  2000     9            -14.6             -0.7
##  9            52244 Cardinals  2000    10            -14.6             -0.7
## 10            64223 Cardinals  2000    11            -14.6             -0.7
## # ... with 7,646 more rows, and 1 more variable: playoffs <chr>

Especificación del modelo lm

lm_spec <- linear_reg() %>%   #parsnip
  set_engine(engine = "lm")

Ajuste del modelo utilizando lm

lm_fit <- lm_spec %>%
  fit(weekly_attendance ~ .,
    data = nfl_train)

lm_fit

## parsnip model object
## 
## Fit time:  51ms 
## 
## Call:
## stats::lm(formula = weekly_attendance ~ ., data = data)
## 
## Coefficients:
##          (Intercept)        team_nameBears      team_nameBengals  
##            -81107.86              -2879.80              -4875.47  
##       team_nameBills      team_nameBroncos       team_nameBrowns  
##              -361.08               2805.24               -324.11  
##  team_nameBuccaneers    team_nameCardinals     team_nameChargers  
##             -3063.65              -6139.80              -6489.31  
##      team_nameChiefs        team_nameColts      team_nameCowboys  
##              1974.33              -3392.79               6068.70  
##    team_nameDolphins       team_nameEagles      team_nameFalcons  
##               139.68               1259.16               -204.17  
##      team_nameGiants      team_nameJaguars         team_nameJets  
##              5447.10              -3095.46               4044.23  
##       team_nameLions      team_namePackers     team_namePanthers  
##             -3480.69               1114.11               1227.32  
##    team_namePatriots      team_nameRaiders         team_nameRams  
##              -214.20              -6324.74              -2884.85  
##      team_nameRavens     team_nameRedskins       team_nameSaints  
##              -398.90               6447.07                380.98  
##    team_nameSeahawks     team_nameSteelers       team_nameTexans  
##             -1405.89              -3567.81                264.07  
##      team_nameTitans      team_nameVikings                  year  
##             -1118.23              -3183.08                 74.73  
##                 week     margin_of_victory  strength_of_schedule  
##               -72.83                137.58                230.74  
##     playoffsPlayoffs  
##              -427.94

Ya fiteamos el primer modelo, ahora vamos por el segundo.

rf_spec <- rand_forest(mode = "regression") %>%
  set_engine("ranger")

rf_spec

## Random Forest Model Specification (regression)
## 
## Computational engine: ranger

Ajuste del modelo utilizando random forest

library('ranger')
rf_fit <- rf_spec %>%
  fit(weekly_attendance ~ ., data = nfl_train)

rf_fit

## parsnip model object
## 
## Fit time:  4.2s 
## Ranger result
## 
## Call:
##  ranger::ranger(formula = weekly_attendance ~ ., data = data,      num.threads = 1, verbose = FALSE, seed = sample.int(10^5,          1)) 
## 
## Type:                             Regression 
## Number of trees:                  500 
## Sample size:                      7656 
## Number of independent variables:  6 
## Mtry:                             2 
## Target node size:                 5 
## Variable importance mode:         none 
## Splitrule:                        variance 
## OOB prediction error (MSE):       74734819 
## R squared (OOB):                  0.08221351

Aclaración: hemos ajustado ambos modelos utilizando ´nfl_train´, es decir, los datos de entrenamiento. No hemos tocado los datos de las pruebas durante el entrenamiento.

Evaluación de los modelos

Cuando sea el momento de evaluar nuestros modelos (para estimar qué tan bien funcionarán nuestros modelos con datos nuevos), analizaremos nfl_test. Podemos predecir cuál será la asistencia semanal a los partidos de la NFL tanto para los datos de entrenamiento como para los datos de prueba utilizando los modelos lm y de random forest. Uno de los objetivos de tidymodels es poder utilizar código como el siguiente de formas predecibles y coherentes para muchos tipos de modelos, y utilizar las herramientas de tidyverse adecuadas para este tipo de tareas.

results_train <- lm_fit %>%
  predict(new_data = nfl_train) %>%
  mutate(
    truth = nfl_train$weekly_attendance,
    model = "lm"
  ) %>%
  bind_rows(rf_fit %>%
    predict(new_data = nfl_train) %>%
    mutate(
      truth = nfl_train$weekly_attendance,
      model = "rf"
    ))

results_test <- lm_fit %>%
  predict(new_data = nfl_test) %>%
  mutate(
    truth = nfl_test$weekly_attendance,
    model = "lm"
  ) %>%
  bind_rows(rf_fit %>%
    predict(new_data = nfl_test) %>%
    mutate(
      truth = nfl_test$weekly_attendance,
      model = "rf"
    ))
results_test

## # A tibble: 5,104 x 3
##     .pred truth model
##     <dbl> <dbl> <chr>
##  1 59162. 65356 lm   
##  2 59017. 50289 lm   
##  3 58871. 37452 lm   
##  4 58798. 65711 lm   
##  5 66880. 73025 lm   
##  6 66589. 66019 lm   
##  7 66079. 44680 lm   
##  8 65860. 41017 lm   
##  9 68096. 68481 lm   
## 10 67877. 83252 lm   
## # ... with 5,094 more rows

Para este modelo de regresión, la métrica que tendremos en cuenta para evaluar el modelo es el rmse o la raíz del error cuadrático medio de lo que hemos hecho hasta ahora.

results_train %>%
  group_by(model) %>%
  rmse(truth = truth, estimate = .pred) %>% view

results_test %>%
  group_by(model) %>%
  rmse(truth = truth, estimate = .pred) #%>% view

## # A tibble: 2 x 4
##   model .metric .estimator .estimate
##   <chr> <chr>   <chr>          <dbl>
## 1 lm    rmse    standard       8351.
## 2 rf    rmse    standard       8573.

results_test

## # A tibble: 5,104 x 3
##     .pred truth model
##     <dbl> <dbl> <chr>
##  1 59162. 65356 lm   
##  2 59017. 50289 lm   
##  3 58871. 37452 lm   
##  4 58798. 65711 lm   
##  5 66880. 73025 lm   
##  6 66589. 66019 lm   
##  7 66079. 44680 lm   
##  8 65860. 41017 lm   
##  9 68096. 68481 lm   
## 10 67877. 83252 lm   
## # ... with 5,094 more rows

¿Malas decisiones?

Si miramos los datos de entrenamiento, el modelo de bosque aleatorio funcionó mucho mejor que el modelo lineal; el rmse es mucho menor. Sin embargo, no se puede decir lo mismo de los datos de prueba. 😭 La métrica para el entrenamiento y las pruebas para el modelo lineal es aproximadamente la misma, lo que significa que no hemos sobreajustado. Para el modelo de bosque aleatorio, la rmse es mucho más alta para los datos de prueba que para los datos de entrenamiento. Nuestros datos de entrenamiento no nos dan una buena idea de cómo funcionará nuestro modelo, y este poderoso algoritmo ML se ha sobreajustado a este conjunto de datos.

results_test %>%
  mutate(train = "testing") %>%
  bind_rows(results_train %>%
    mutate(train = "training")) %>%
  ggplot(aes(truth, .pred, color = model)) +
  geom_abline(lty = 2, color = "gray80", size = 1.5) +
  geom_point(alpha = 0.5) +
  facet_wrap(~train) +
  labs(
    x = "Valor real",
    y = "Asitencia predicha",
    color = "Tipo de modelo"
  )

Regresión con datos de NFL_attendance

Basado en el blog de Julia Silge

29/9/2020