---
title: 'Jump Shot to Conclusions - March Madness EDA'
date: '`r Sys.Date()`'
output:
  html_document:
    number_sections: true
    fig_caption: true
    toc: true
    fig_width: 7
    fig_height: 4.5
    theme: cosmo
    highlight: tango
    code_folding: hide
---

```{r setup, include=FALSE, echo=FALSE}
knitr::opts_chunk$set(echo=TRUE, error=FALSE)
knitr::opts_chunk$set(out.width="100%", fig.height = 4.5, split=FALSE, fig.align = 'default')
```

<div 
    style="background-image: url('https://images.unsplash.com/photo-1538405505827-e519f0efb68e?ixlib=rb-1.2.1&ixid=eyJhcHBfaWQiOjEyMDd9&auto=format&fit=crop&w=1350&q=80'); 
    width:100%; 
    height:600px; 
    background-position:center;">&nbsp;
</div>

Photo by DiAnte Squire on Unsplash.


# Introduction

Welcome to an initial Exploratory Data Analysis for the [March Madness 2020](https://www.kaggle.com/c/google-cloud-ncaa-march-madness-2020-division-1-mens-tournament) competition with [tidyverse](http://tidyverse.org/) tools. This is an exploration of the Mens Tournament data, but I'm sure you can apply similar ideas to the [Womens Tournament](https://www.kaggle.com/c/google-cloud-ncaa-march-madness-2020-division-1-womens-tournament/) challenge.

**The goal** is to predict the detailed outcome of this year's final bracket in the US [National Collegiate Athletic Association (NCAA)](https://en.wikipedia.org/wiki/National_Collegiate_Athletic_Association) basketball tournament. We are tasked to predict a probability for which team will win in every possible matchup between two teams in the tournament. More about [March Madness on Wikipedia](https://en.wikipedia.org/wiki/NCAA_Division_I_Men%27s_Basketball_Tournament). This competition starts with a month of evaluating models on historical data. Thereafter, a second stage will run in which the model is scored against the live results of the bracket.

**The metric** is [log loss](https://www.kaggle.com/c/google-cloud-ncaa-march-madness-2020-division-1-mens-tournament/overview/evaluation), which means that the more certain your algorithm is (i.e. winning probability close to 0 or 1) the more severe is the penalty if it's wrong.

**The data:** This competition comes with tons of [data](https://www.kaggle.com/c/google-cloud-ncaa-march-madness-2020-division-1-mens-tournament/data), so make sure to check out the data overview below.

**Related Competitions:** This is far from the first time that Kaggle is hosting a March Madness challenge. You can find lots inspiration in competitions from previous years: [2019](https://www.kaggle.com/c/mens-machine-learning-competition-2019), [2018](https://www.kaggle.com/c/mens-machine-learning-competition-2018).

*Structure:* This notebook is intended to be read from top to bottom, but the first part will contain a comprehensive exploratory analysis and in the second part I focus on a specific aspect of the data. The notebook will also evolve non-linearly: I will analyse the Basic Data first, and look at simple correlations and insights. Then I will gradually explore the other data sets and add new aspect to the multi-feature analysis. (This notebook also uses tabs to keep its subsections tidy and organised).

Alright, enough warm up. Here's our jump ball and let's go!



# Preparations {.tabset .tabset-fade .tabset-pills}

## Load libraries

We load a range of libraries for general data wrangling and general visualisation together with more specialised tools.

```{r, message = FALSE}
# general visualisation
library('ggplot2') # visualisation
library('scales') # visualisation
library('patchwork') # visualisation
library('RColorBrewer') # visualisation
library('corrplot') # visualisation
library('ggthemes') # visualisation
library('ggrepel') # visualisation

# general data manipulation
library('dplyr') # data manipulation
library('readr') # input/output
library('vroom') # input/output
library('skimr') # overview
library('tibble') # data wrangling
library('tidyr') # data wrangling
library('stringr') # string manipulation
library('forcats') # factor manipulation

# specific visualisation
library('alluvial') # visualisation
library('ggrepel') # visualisation
library('ggforce') # visualisation
library('ggridges') # visualisation
library('gganimate') # animations
library('GGally') # visualisation
library('ggExtra') # visualisation
library('viridis') # visualisation
library('usmap') # geo

# specific data manipulation
library('lazyeval') # data wrangling
library('broom') # data wrangling
library('purrr') # data wrangling
library('reshape2') # data wrangling
library('rlang') # encoding
library('kableExtra') # display

# modelling
library(xgboost) # model
library(parsnip) # model
library(yardstick) # metrics
```


## Helper functions

We make use of a brief helper function to compute binomial confidence intervals.

```{r}
# function to extract binomial confidence levels
get_binCI <- function(x,n) as.list(setNames(binom.test(x,n)$conf.int, c("lwr", "upr")))
```


## Load data

We use the [vroom](https://cran.r-project.org/web/packages/vroom/vignettes/vroom.html) package for reading data:

```{r echo=FALSE}
if (dir.exists("/kaggle")){
  path <- "/kaggle/input/google-cloud-ncaa-march-madness-2020-division-1-mens-tournament/"
  cpath <- "/kaggle/input/world-cities-database/"
} else {
  path <- ""
  cpath <- ""
}

subpath <- "MDataFiles_Stage1/"
```

```{r warning=FALSE, results=FALSE}
# basic data
teams <- vroom(str_c(path, subpath, "MTeams.csv"), col_types = cols())
seasons <- vroom(str_c(path, subpath, "MSeasons.csv"), col_types = cols())
seeds <- vroom(str_c(path, subpath, "MNCAATourneySeeds.csv"), col_types = cols())
regular_res <- vroom(str_c(path, subpath, "MRegularSeasonCompactResults.csv"), col_types = cols())
bracket_res <- vroom(str_c(path, subpath, "MNCAATourneyCompactResults.csv"), col_types = cols())
sample_submit <- vroom(str_c(path, "MSampleSubmissionStage1_2020.csv"), col_types = cols())

# beyond section 1
regular_detail <- vroom(str_c(path, subpath, "MRegularSeasonDetailedResults.csv"), col_types = cols())
bracket_detail <- vroom(str_c(path, subpath, "MNCAATourneyDetailedResults.csv"), col_types = cols())
cities <- vroom(str_c(path, subpath, "Cities.csv"), col_types = cols())
cities_games <- vroom(str_c(path, subpath, "MGameCities.csv"), col_types = cols())
ranks <- vroom(str_c(path, subpath, "MMasseyOrdinals.csv"), col_types = cols())


# events15 <- vroom(str_c(path, "MEvents2015.csv"), col_types = cols())
# events16 <- vroom(str_c(path, "MEvents2016.csv"), col_types = cols())
# events17 <- vroom(str_c(path, "MEvents2017.csv"), col_types = cols())
# events18 <- vroom(str_c(path, "MEvents2018.csv"), col_types = cols())
#events19 <- vroom(str_c(path, "MEvents2019.csv"), col_types = cols())

#players <- vroom(str_c(path, "MPlayers.csv"), col_types = cols())


# auxil data:
world_cities <- vroom(str_c(cpath, "worldcitiespop.csv"), col_types = cols())
```


# The 2019 bracket - A visual aid

Before we look into the data, let's get an idea of what "the bracket" means.

Here is the kind tournament structure which the teams will have to navigate. This is the **2019 bracket** which can be found on the [official NCAA site](https://www.ncaa.com/news/basketball-men/ncaa-bracket-march-madness):

<br>
<center><img src="https://www.ncaa.com/sites/default/files/public/styles/original/public-s3/images/2019/04/09/ncaa-tournament-bracket-2019-scores-games-virginia-texas-tech.png?itok=0E3VNWmI"></center>
<br>

Last year, Virginia made its way from the South region all the way to beat Texas Tech and become national champion. The little numbers on the outside of the team names at the start of the bracket are the **seeds**; those determine the path of the teams. These seeds correspond to a ranking based on the team's performance in the regular season: the number 1 performing team will get seeded 1st.

Typically, strong teams have higher seeds (i.e. lower numbers) and will meet weaker teams early on, which gives them a good chance to make it far before they meet similarly high-seeded teams from the other corners of the bracket. But of course sometimes there are suprises and everybody loves an underdog (or a "cinderella story"). In 2019 the biggest upset was probably Oregon (seed 12) beating Wisconsin (seed 5) and advancing further before being knocked out by the eventual champions Virginia. Other than that, you will notice that all the teams that made it to the last 8 had pretty high seeds. There's your baseline prediction :-)


# Overview: Data file structure and content {.tabset .tabset-fade .tabset-pills}

To start our exploration, we'll get a quick overview of the data sets using the *summary* and *glimpse* tools. We will follow the characterisation of data in the [data description](https://www.kaggle.com/c/google-cloud-ncaa-march-madness-2020-division-1-mens-tournament/data). Some initial insights can be derived from the data. All background info was condensed from the excellent data description.

**Note, that seasons always last from November to March and the years in the data refer to the year the season ended.** E.g. the 2018/19 season would be identified as `2019`.


## Basic data - Data Section 1 {.tabset .tabset-fade .tabset-pills}

### Teams

Glimpse:

```{r}
glimpse(teams)
```

- These are the `TeamID`s and corresponding compact `TeamName`s, together with the first and last season that this team was playing in Divison 1. From season to season, certain teams can drop from or ascend to Division 1 and will thus appear and disappear in our data.

- We can already see that most, but not all, teams have been in Division 1 since at least 1985 and are still present there today.


### Seasons

Summary:

```{r}
summary(seasons %>% 
          mutate_if(is.character, as.factor))
```


Last 5 rows:

```{r}
seasons %>% 
  tail(5) %>% 
  kable() %>% 
  kable_styling()
```


- Here we have the historical seasons since 1985. **The `Season` variable indicates the year the season ended.** The `DayZero` feature allows us to relate the date to the relative day numbering used in the game data files. **The final game of the season ("national championship") is always played on day 154 (a Monday).**

- The purpose of the `Region` features is to encode the bracket setup for the 4 regions "East", "West", "Midwest", and "Southeast". `RegionW` will always play `RegionX` in the semi-finals (i.e. they are on the same side of the bracket). Same with `RegionY` vs `RegionZ` on the other side of the bracket. For 2020, the bracket structure will only be known on Sunday March 15 ("Selection Sunday") since it depends on the rankings of the 4 top seeds from each region.


### Seeds

Glimpse:

```{r}
glimpse(seeds)
```

- The seeds data contains the historical starting `Seed`s for each `TeamID` and `Season`. This tells us how the bracket is populated and which teams are facing off against each other.

- The `Seed` feature starts with a letter that encodes the region (W/X/Y/Z) and then includes the seed number within the region. Some seeds contain characters at the end to indicate a team that needed to go through a pre-bracket play-in to qualify for participation in the bracket. Higher seeds (i.e. smaller numbers) are stronger teams that are more likely to make it far in the tournament.


### Regular Season Results

First 5 rows:

```{r}
head(regular_res, 5) %>% 
  kable() %>% 
  kable_styling()
```


Summary:

```{r}
summary(regular_res %>% 
          mutate(WLoc = as.factor(WLoc)))
```


- The compact results contain the results of all games in the regular part of each `Season`. This regular part are all games leading up to but not including the bracket. They always occur on `DayNum <= 132` (which is Selection Sunday).

- In `WTeamID` we find the ID of the winning team and `WScore` is their score for this game. The `L` prefix columns contain the corresponding information for the losing team. The `WLoc` feature indicates whether the winning team played at home ('H'), away ('A'), or on a neutral court ('N').

- `NumOT` is the number of overtimes in the game. In basketball, 5-minute overtimes occur when the score is tied at the end of the final quarter. There can be multiple overtimes to determine a winner (here the maximum is 6; which I imagine to be rather exhausting).


### Bracket Results


First 5 rows:

```{r}
head(bracket_res, 5) %>% 
  kable() %>% 
  kable_styling()
```

- The bracket results data have the same structure as the regular season results above. Location is always neutral: `WLoc = N`.

- The days of the game always start on `DayNum = 136` or `DayNum = 134` if there were play-in games (= "first four"; played on days 134 and 135).

- Through the bracket structure we can link the stages of the tournament to specific days: the "first round" (64 teams) is played on days 136 and 137, the "second round" (32 teams) on days 138/139, the "regional semifinals" (= "Sweet Sixteen", last 16 teams) happen on days 143/144, "regional finals" on days 145/146 (last 8; i.e national quarterfinals aka "Elite Eight"), the "national semifinals" on day 152 (those are also "Final Four" teams), and the final ("National Championship") always happens on day 154. (See the bracket in Sect. 3 for the 2019 setup.)


###  Missing data

Good news: There are no missing data in any of the Section 1 data files:

```{r}
missing <- tibble(
  table_name = c("Teams", "Seasons", "Seeds", "Regular Season Results", "Bracket Results"),
  NA_count = c(sum(is.na(teams)), sum(is.na(seasons)), sum(is.na(seeds)), sum(is.na(regular_res)), sum(is.na(bracket_res)))
)

missing %>% 
  kable() %>% 
  kable_styling()
```


### Reformating features

This is for visualisation purposes. We change some column types and add a few convenience features.

```{r}
#seasons <- seasons %>% 
#  mutate_if(is.character, as.factor)

regular_res <- regular_res %>% 
  mutate(WLoc = as.factor(WLoc))

bracket_res <- bracket_res %>% 
  select(-WLoc)

```


## Team Box Scores - Data Section 2 {.tabset .tabset-fade .tabset-pills}

### Regular Season - Detailed Results

Here we use the *skim* tool from the [skimr](https://github.com/ropensci/skimr) package for a quick overview:

```{r}
regular_detail %>% 
  skim()
```

We find:

- The regular season detailed results are essentially an extension of the compact regular season results (here: `regular_res`). They contain the same information (i.e. Season, TeamIDs, Scores) but add a whole bunch of additional features such as field goals made, 3-pointers attempted, rebounds, fouls, steals, blocks, and more. Just like in the regular season data, all features are available for winning ("W") and losing ("L") teams.

- There are no missing values. The data description notes that "field goals" refers to both 2-pointers and 3-pointers.

- However, there is a big caveat: all this rich data is only available for Season 2003 and after. From 2003 onward, though, it can essentially replace the regular season data completely. Here is a comparison of the first 6 rows of regular (>= 2003) data compared to the first 6 rows of detailed data:


```{r}
head(regular_res %>% filter(Season >= 2003))
```

```{r}
head(regular_detail %>% 
       select(1:8))
```



### Bracket - Detailed Results

Here is the same `skim` view for the detailed bracket data:

```{r}
bracket_detail %>% 
  skim()
```

We find:

- The statements from the previous subsection also apply here. The detailed bracket data starts in 2003, too. Also here there are no missing values.


## Geography - Data Section 3 {.tabset .tabset-fade .tabset-pills}

Matches are being played in certain cities. Almost all cities are in the US.

### Master list of cities

```{r}
glimpse(cities)
```

This list provides the names and states of the cities, together with a unique ID. Those are all the cities that ever have hosted an NCAA game relevant to this competition.


### Game Cities

This file links all the regular season games, the bracket matches, and other post-season tournaments to the data of the games themselves. It contains only matches during or after the 2010 season. Those are the first rows:

```{r}
cities_games %>% 
  head(5) %>% 
  kable() %>% 
  kable_styling()
```

We see that the Season, (DayNum) and Team IDs allow us to uniquely join this table to the `Compact` or `Detailed` regular/bracket results. The `CityID` joins to the "Cities" table. The `CRType` describes which part of the season the game was played in:

```{r}
cities_games %>% 
  distinct(CRType) %>% 
  kable() %>% 
  kable_styling()
```

Here, "Regular" stands of course for the regular season, "NCAA" is the bracket, and "Secondary" are post-season matches.


## Public Rankings - Data Section 4 {.tabset .tabset-fade .tabset-pills}

Going back to the 2003 Season, this dataset provides team rankings from a large number of different methodologies and sources. The data description notes that all systems are expressed here as ordinal ranks, rather than numerical ratings which would include information about the specific gap between two teams.

Here are 6 random rows:

```{r}
set.seed(4321)
ranks %>% 
  sample_n(6) %>% 
  kable() %>% 
  kable_styling()
```


We find:

- `Season` and `TeamID` link to our game stats. `OrdinalRank` is obviously the ranking itself.

- Ranking is done at multiple times per season, usually on a weekly basis, and `RankingDayNum` gives the day of the ranking in the same relative system as our match data. The [data description](https://www.kaggle.com/c/google-cloud-ncaa-march-madness-2020-division-1-mens-tournament/data) tells us that the final pre-tournament rankings have `RankingDayNum = 133`.

- There are a total of 174(!) ranking systems (`SystemName`) in this dataset, and those are the 5 with the most entries:

```{r}
ranks %>% 
  count(SystemName) %>% 
  arrange(desc(n)) %>% 
  head(5) %>% 
  kable() %>% 
  kable_styling()
```





# Basic Data: Visual EDA + Baseline Model

After familiarising ourselves with the basic structure of the data, we will now use visualisation to explore the various facets of our data sets. Here we will start with Data Section 1 and build a baseline model. The other data files will be the subject of the next section.

Data section 1 covers teams, seasons, and match overview statistics for the regular seasons and the brackets. All data starts in 1985.


## Teams

We first want to get an idea from our `teams` data frame about how common it is that teams drop out or appear during our coverage time frame. We plot histograms of all first seasons next to the corresponding view of all last seasons. Then we complement this view with a count of all instances where a team had its first 1st Division season in 1985 or a team is currently in the 1st division

```{r fig.cap ="Fig. 1", fig.height = 6}
p1 <- teams %>% 
  select(-TeamName) %>% 
  pivot_longer(ends_with("Season"), names_to = "first_last", values_to = "year") %>% 
  ggplot(aes(year, fill = first_last)) +
  geom_histogram(bins = 50) +
  scale_fill_manual(values = c("orange", "darkgreen")) +
  theme_fivethirtyeight() +
  theme(legend.position = "none") +
  facet_wrap(~ first_last) +
  labs(x = "Year", y = "", title = "Count of teams with first/last season in year")

p2 <- teams %>% 
  select(-TeamName) %>% 
  pivot_longer(ends_with("Season"), names_to = "first_last", values_to = "year") %>% 
  mutate(check = (first_last == "FirstD1Season" & year == 1985) | (first_last == "LastD1Season" & year == 2020)) %>% 
  group_by(first_last, check) %>% 
  count() %>% 
  ungroup() %>% 
  mutate(first_last = if_else(first_last == "FirstD1Season", "FirstD1Season = 1985", "LastD1Season = 2020")) %>% 
  ggplot(aes(check, n, fill = check)) +
  geom_col() +
  geom_label(aes(label = n)) +
  theme_fivethirtyeight() +
  theme(legend.position = "none") +
  facet_wrap(~ first_last) +
  labs(x = "Year", y = "", title = "Number of Teams with 1st/last season 1985/2020")

p1 / p2
```

We find:

- It's uncommon but not rare that 1st division teams appear in our data after 1985. The last season, however, looks much more focused on 2020.

- The aggregate view confirms this impression: 85 teams had their initial 1st division season after 1985, but only 14 are not present in the current season. As the data description states, there are 353 teams in the current season 2019/20.



## Seasons and regions

Now we take the `seasons` data to determine whether there are any typical starting positions on the bracket for the various regions. Remember that this encodes which regions were facing off against each other from the start and which could only meet in the final.

We use a facetted barplot to count the occurences of each region within our 4 region categories. The bars are roughly ordered by size and the regions have the same colours independent of in which facet they occur:

```{r fig.cap ="Fig. 2"}
seasons %>% 
  select(Season, starts_with("Region")) %>% 
  pivot_longer(starts_with("Region"), names_to = "reg_code", values_to = "regs") %>% 
  filter(!str_detect(regs, "TBD")) %>% 
  count(reg_code, regs) %>% 
  ggplot(aes(reorder(regs, n, FUN = min), n, fill = regs)) +
  geom_col() +
  geom_label(aes(label = n)) +
  coord_flip() +
  theme_fivethirtyeight() +
  theme(legend.position = "none") +
  labs(x = "", y = "", title = "Region codes over all seasons",
       subtitle = "Regions W & X and regions Y & Z are always on the same side of the bracket") +
  facet_wrap(~ reg_code, scales = "free_y")
```

We find:

- Starting out with "East" as "RegionW" is the alphabetic convention. The matching "RegionX" is pretty equally distributed between "Midwest", "West", and "South". Conversely, regions "Y" and "Z" are most often "Midwest" and "West".

Note, that **the region codes "W, X, Y, Z" are a convention of this Kaggle dataset**, with the aim to homogenise the changing region names over time. There are always 4 regions, but their names changed over time (e.g. "South" was "Southeast" before 2011).


So what are the most common pairings in our data? Let's combine the entries in regions "W vs X" and "Y vs Z" to figure out who's most likely to meet from the outset:

```{r fig.cap ="Fig. 3"}
seasons %>% 
  unite("Region W vs X", RegionW, RegionX, sep = " vs ") %>% 
  unite("Region Y vs Z", RegionY, RegionZ, sep = " vs ") %>% 
  pivot_longer(starts_with("Region"), names_to = "bracket", values_to = "pairing") %>% 
  count(bracket, pairing) %>% 
  filter(!str_detect(pairing, "TBD")) %>% 
  ggplot(aes(reorder(pairing, n, FUN = min), n, fill = n)) +
  geom_col() +
  coord_flip() +
  scale_y_continuous(breaks = pretty_breaks()) +
  scale_fill_viridis_c(option = "viridis") +
  theme_fivethirtyeight() +
  theme(legend.position = "none") +
  #facet_wrap(~bracket) +
  labs(x = "", y = "", title = "Most frequent regional pairings")
  
```

We find:

- "Midwest vs West" is historically most common, whereas "East vs Southeast" was rare.

- The various small regions illustrate the changing names of the (always) 4 regions over the years. According to [Wikipedia](https://en.wikipedia.org/wiki/NCAA_Division_I_Men%27s_Basketball_Tournament#Seeding_history_and_statistics), while East and West had always the same names, "Midwest" was "Southwest" in 2011 and "South" was alternating with "Southeast" (used in 1985–1999 and 2011).



## Seeds

The tournament seeding is based on the team's performance in the regular season just preceding the bracket. The teams are ranked from 1 to (now) 68 and seeds are assigned accordingly. Therefore, the seed is probably the strongest predictor for a good tournament performance. Here we take a closer look at the previous seeds.

First off: The most common teams in each seed over the complete time frame. Keep in mind, that there are 4 regions per bracket and therefore 4 top seeds (and 4 of every of the 16 seeds, of course). Here the facet plot shows the count for each seed, with seeds coded by facet & colour simultaneously:

```{r fig.cap ="Fig. 4", fig.height = 6}
foo <- seeds %>% 
  mutate(seed = str_sub(Seed, 2,3)) %>% 
  group_by(TeamID, seed) %>% 
  count() %>% 
  ungroup() %>% 
  arrange(seed, desc(n)) %>% 
  inner_join(teams %>% select(TeamID, TeamName), by = "TeamID") %>% 
  group_by(seed) %>% 
  dplyr::slice(c(1:4)) %>% 
  ungroup()

foo %>% 
  ggplot(aes(reorder(TeamName, n, FUN = min), n, fill = seed)) +
  geom_col() +
  coord_flip() +
  facet_wrap(~seed, scales = "free") +
  theme_fivethirtyeight() +
  theme(legend.position = "none") +
  labs(title = "Top teams for each seed (all time)")
```

We find:

- Duke has show a strong all-time performance: it is tied with "North Carolina" and "Kansas" for 1st place in top seed, has the strongest standing in 2nd seed, and is even present in the 3rd seed top list. Kansas and Kentucky are almost as strong. (Note, that because we restrict the list for each seed at 4 teams there might sometimes be other teams that share the same count as the lowest team count, but have been cut off. This is due to the overall small numbers.)

- Last year's winner "Virginia" turns up in the top teams for seed 5, but not higher. The lower seeds are pretty diverse with almost no repetitions.

Note, that the `TeamName` here is a standardised version that is consistent throughout the data files provided by Kaggle. In any external files those names will likely be different. A helpful look-up table with alternative or historical names is included in the `MTeamSpellings` file.


So this is the aggregated picture, but how did the seedings evolve over the years? There are too many teams in our data to plot all of them, but we can have a closer look at the history of some of the top teams:


```{r fig.cap ="Fig. 5", fig.height = 5}
foo <- seeds %>% 
  mutate(seed = str_sub(Seed, 2,3)) %>% 
  group_by(TeamID, seed) %>% 
  count() %>% 
  ungroup() %>% 
  arrange(seed, desc(n))

top_seeds <- foo %>% 
  head(12) %>% 
  distinct(TeamID)

seeds %>% 
  mutate(seed = as.integer(str_sub(Seed, 2,3))) %>% 
  inner_join(teams %>% select(TeamID, TeamName), by = "TeamID") %>% 
  inner_join(top_seeds, by = "TeamID")  %>% 
  #filter(TeamName == "Duke") %>% 
  ggplot(aes(Season, seed, col = TeamName, group = TeamName)) +
  geom_line() +
  #geom_point() +
  scale_y_reverse(breaks = seq(1,15,2)) +
  theme_fivethirtyeight() +
  theme(legend.position = "none") +
  facet_wrap(~TeamName) +
  labs(x = "Season", y = "Seed", title = "Variations in seeding of top teams",
       subtitle = "Even strong teams show occasional dips in regular season performance")

```

We find:

- While there is a certain consistency for teams like "Duke" or "Kansas", all of the strong teams show quite a bit of variation in their seeds over the last 35 years. This is perhaps not surprising, since college athletes often look for a placement in the NBA, especially the strongest players.

- For some teams like "Oklahoma", "Connecticut", and "Ohio State" there is a downturn in performance over the last year, the fates smile upon Villanove, (recently) Arizona, and of course Virginia.


Let's linger a bit on last years finalists: Virginia and Texas Tech. Note, that in the time series below both have notable gaps in years where they didn't reach the final bracket at all:

```{r fig.cap ="Fig. 6"}
seeds %>% 
  mutate(seed = as.integer(str_sub(Seed, 2,3))) %>% 
  inner_join(teams %>% select(TeamID, TeamName), by = "TeamID") %>% 
  filter(TeamName %in% c("Virginia", "Texas Tech")) %>% 
  mutate(label = if_else(Season == 2007, TeamName, NA_character_)) %>% 
  ggplot(aes(Season, seed, col = TeamName, group = TeamName)) +
  geom_line() +
  geom_point() +
  geom_label_repel(aes(label = label), nudge_x = 0, na.rm = TRUE, alpha = 0.8) +
  scale_y_reverse(breaks = seq(1,15,2)) +
  theme_fivethirtyeight() +
  theme(legend.position = "none") +
  labs(x = "Season", y = "Seed", title = "2019 finalists seeding history",
       subtitle = "Recent upwards trends for both Virginia and Texas Tech")
```

We find:

- While both colleges have been doing well in the last few year, they have previously shown average seeding performances (and in the case of Texas Tech rather long gaps).



## Regular Season Results

There is quite a bit of useful information in the history of the regular season games, because it tracks the teams performances up to the March Madness bracket. We start with a couple of overview visuals:

```{r fig.cap ="Fig. 7", fig.height = 6}
p1 <- regular_res %>% 
  count(WLoc) %>% 
  add_tally(n, name = "total") %>% 
  mutate(percent = n/total) %>% 
  mutate(foo = if_else(WLoc == "H", "red", "grey70")) %>% 
  ggplot(aes(reorder(WLoc, -percent, FUN = min), percent, fill = foo)) +
  geom_col() +
  scale_fill_identity() +
  scale_y_continuous(labels = scales::percent) +
  theme_hc() +
  labs(x = "", y = "", title = "Winning advantage: Home", 
       subtitle = "VS % of wins Away / Neutral Court")

p2 <- regular_res %>% 
  select(ends_with("Score")) %>% 
  rownames_to_column() %>% 
  pivot_longer(ends_with("Score"), names_to = "pos", values_to = "score") %>%
  mutate(pos = if_else(pos == "WScore", "Winner", "Loser")) %>% 
  ggplot(aes(score, fill = pos)) +
  geom_density(alpha = 0.5) +
  coord_cartesian(xlim = c(25, 130)) +
  theme_hc() +
  theme(legend.position = "top") +
  labs(x = "Score", y = "", fill = "", title = "Global Scores of\nwinning & losing teams")

p3 <- regular_res %>% 
  select(Season, ends_with("Score")) %>% 
  rownames_to_column() %>% 
  pivot_longer(ends_with("Score"), names_to = "pos", values_to = "score") %>% 
  mutate(pos = if_else(pos == "WScore", "Winner", "Loser")) %>% 
  ggplot(aes(x = score, y = as.factor(Season), fill = pos)) +
  geom_density_ridges(alpha = 0.5, bandwidth = 2) +
  coord_cartesian(xlim = c(25, 130)) +
  theme_hc() +
  theme(legend.position = "none") +
  labs(x = "Scores", y = "", title = "Scores per Season")

((p1 / p2) | p3)

```

We find:

- There is a definite home advantage: Almost 60% of a time a team is winning, they do so at home. This assumes that the home vs away games are somewhat equally distributed, but over all teams they should be; unless already stronger teams would somehow be allowed to play more frequently at home. **Of course, all the bracket games happen on neutral ground, so the winning location is not a direct predictor here but rather an adjustment factor for encoding past performance.**

- Looking at the overall distributions of winning and losing scores we see a lot of overlap. Scoring more than 100 points virtually guarantees a win, but already the average winning score of about 75 points has a signficant overlap with the tail of the losing distribution.

- We use a ridgeline plot to compare those winning and losing score distributions for each season over time. There is a trend towards somewhat higher winning scores over time, although this fluctuates a bit throughout the years. Overall, there's no big difference between average winning and losing scores for either Season


But what's more important than the overall scores is of course the difference between the winning and the losing score. A win is a win, sure, but a difference of 20 points should inspire more confidence as a prediction input than a difference of 2 points. We will look at the winning margins, i.e. the difference between winning and losing score, aggregated by decade and by time in the season which we split into 3 roughly equal parts: (note the x-axis log scale)

```{r fig.cap ="Fig. 8"}
regular_res %>% 
  select(DayNum, Season, ends_with("Score")) %>% 
  mutate(decade = str_c(as.character(Season %/% 10 * 10), "s")) %>% 
  mutate(day_bin = case_when(
    DayNum < 45 ~ "early season",
    between(DayNum, 45, 90) ~ "mid season",
    DayNum > 90 ~ "late season"
  )) %>% 
  mutate(day_bin = fct_relevel(as.factor(day_bin), "late season", after = Inf)) %>% 
  mutate(win = WScore - LScore) %>% 
  ggplot(aes(win, fill = decade)) +
  geom_density(bw = 0.08) +
  facet_grid(day_bin ~ decade) +
  scale_x_log10(breaks = c(1, 2, 5, 10, 20, 50)) +
  theme_fivethirtyeight() +
  theme(legend.position = "none") +
  labs(title = "Consistent winning point margins in Regular Season",
       subtitle = "Early Season < day 45; Late Season > day 90. Log Scale.")
```

We find that there is really no big difference between the winning point margins in the early season (before day 45), mid season (days 45-90), and late season (after day 90). This picture holds true for the four decades that our data covers.


Now let's move on to the teams themselves. The performance data for the regular season also tells us which teams were consistently winning and which ones were losing. Since the a stronger performance in the regular season results in a better seed, this will be not only a predictor not only for how well a team might be doing in the bracket but also which path they might take.

Here we look at the six top performing teams per decade based on the number of games they won. For our visualisation we choose ordered barplots with decades as facets, plus a Tufte-style theme:

```{r fig.cap ="Fig. 9"}
regular_res %>% 
  mutate(decade = str_c(as.character(Season %/% 10 * 10), "s")) %>% 
  group_by(decade, WTeamID) %>%
  summarise(n = n()) %>% 
  top_n(6, n) %>% 
  arrange(decade, desc(n)) %>% 
  mutate(ord = seq(6,1)) %>% 
  left_join(teams %>% select(TeamName, TeamID), by = c("WTeamID" = "TeamID")) %>% 
  ggplot(aes(ord, n, col = decade, fill = decade)) +
  #ggplot(aes(fct_reorder(TeamName, ord, .desc = TRUE), n, col = decade, fill = decade)) +
  geom_col() +
  geom_text(aes(label = TeamName), nudge_y = 10, hjust = "left") +
  #geom_label_repel(aes(label = TeamName), nudge_y = 20, na.rm = TRUE, alpha = 0.8) +
  coord_flip(ylim = c(0, 350)) +
  facet_wrap(~ decade, scales = "free") +
  theme_tufte() +
  theme(legend.position = "none", axis.text.y = element_blank(), axis.ticks.y = element_blank()) +
  labs(x = "", y = "", title = "Top 6 teams per decade (regular season)",
       subtitle = "Based on the total number of games won.")
```

We find:

- Back in the 80s, "UNLV" was in a very famous position indeed. "North Carolina" didn't really manage to translate its performances into the new millenium.

- The presence of "Duke" in each of the facets, and of "Kansas" in the 3 recent ones, underlines the historical consistency of these teams. "Gozanga" has been doing pretty well after 2000, too.


We also have information about overtimes in our data. Let's see which teams played the most (and won the most) overtime games since 1985. For the winning percentage we only consider teams that played at least 10 games with overtime. We're not counting more than 1 overtime per game; those are all counted as 1:

```{r fig.cap ="Fig. 10", fig.height=2.5}
prep <- regular_res %>% 
  filter(NumOT > 0) %>% 
  select(WTeamID, LTeamID) %>% 
  pivot_longer(ends_with("ID"), names_to = "cat", values_to = "TeamID") %>% 
  left_join(teams %>% select(starts_with("Team")), by = "TeamID") %>% 
  mutate(cat = if_else(cat == "WTeamID", "win", "lose")) %>% 
  group_by(TeamName) %>% 
  summarise(total = n(),
            win_perc = mean(cat == "win")*100) %>% 
  filter(total > 10) %>% 
  pivot_longer(cols = c(total, win_perc), names_to = "type", values_to = "val")

prep %>% 
  group_by(type) %>%
  arrange(desc(val)) %>% 
  dplyr::slice(c(1:5)) %>% 
  ungroup() %>% 
  mutate(type = if_else(type == "total", "Total Wins", "Winning Percentage")) %>% 
  ggplot(aes(reorder(TeamName, val, FUN = min), val, fill = type)) +
  geom_col() +
  coord_flip() +
  theme_fivethirtyeight() +
  theme(legend.position = "none") +
  facet_wrap(~ type, scales = "free") +
  labs(x = "", y = "", title = "Highest Overtime stats (all time)")
```

We find:

- "West Virginia" and "Austin Peay" have both played more than 60 games of overtime since 1985. (Note, that the ranking does not take into account how long a team has been in the 1st division during that time.)

- And if you're going into overtime, you really don't want to go up against "Temple". "Baylor" and "Gonzaga" are also among the teams that are most likely to convert an overtime to their advantage. (Note, that this ranking does not correct for a possible Home court advantage.)


After getting a better feel for this data, we will check out some match results. Let's start with very one-sided matchups, where one team almost always triumphs over the other. From a visual point of view, we choose a facetted barplot with a special colour coding of orange for the winners, and with labels printed on the bars instead of at the axes. See how the `space == "free"` parameter in `facet_grid` allows us to squish the facet with the losses on one side (the bars & axes are still to scale horizontally):

```{r fig.cap ="Fig. 11", fig.height = 6}
foo <- regular_res %>% 
  count(WTeamID, LTeamID) %>% 
  arrange(desc(n)) %>% 
  left_join(teams %>% select(TeamID, TeamName), by = c("WTeamID" = "TeamID")) %>% 
  rename(WTeamName = TeamName) %>% 
  left_join(teams %>% select(TeamID, TeamName), by = c("LTeamID" = "TeamID")) %>% 
  rename(LTeamName = TeamName)

bar <- foo %>% 
  select(-ends_with("ID")) %>% 
  rename(wins = n) %>% 
  left_join(foo %>% select(-ends_with("ID")), by = c("WTeamName" = "LTeamName", "LTeamName" = "WTeamName")) %>% 
  rename(losses = n) %>% 
  mutate(win_perc = wins/(wins + losses) * 100)

p1 <- bar %>% 
  group_by(WTeamName) %>% 
  arrange(desc(wins)) %>% 
  dplyr::slice(1) %>% 
  ungroup() %>% 
  arrange(desc(wins)) %>% 
  head(5) %>%
  mutate(rank = seq(5,1)) %>%  
  pivot_longer(cols = c(wins, losses), names_to = "cat", values_to = "n") %>% 
  mutate(cat = fct_relevel(as.factor(cat), "losses", after = Inf)) %>% 
  # mutate(WTeamName = if_else(cat == "wins", str_c(WTeamName, " vs"), LTeamName)) %>% 
  mutate(TeamName = if_else(cat == "wins", str_c(WTeamName, " vs ", LTeamName), "")) %>% 
  # rename(TeamName = WTeamName) %>% 
  ggplot(aes(rank, n, fill = cat)) +
  geom_col() +
  geom_text(aes(x = rank, y = 0, label = TeamName), nudge_y = 1, hjust = "left", col = "black") +
  coord_flip() +
  scale_fill_manual(values = c("orange", "grey30")) +
  facet_grid(~cat, space = "free", scales = "free") +
  theme_hc() +
  theme(legend.position = "none", axis.text.y = element_blank(), axis.ticks.y = element_blank()) +
  labs(x = "", y = "", title = "Top One-sided matchups (of all time)",
       subtitle = "By number of wins. One pairing per winning team only.")

p2 <- bar %>% 
  filter(wins + losses > 50) %>% 
  group_by(WTeamName) %>% 
  arrange(desc(win_perc)) %>% 
  dplyr::slice(1) %>% 
  ungroup() %>% 
  arrange(desc(win_perc)) %>%
  head(5) %>%
  mutate(rank = seq(5,1)) %>%  
  pivot_longer(cols = c(wins, losses), names_to = "cat", values_to = "n") %>% 
  mutate(cat = fct_relevel(as.factor(cat), "losses", after = Inf)) %>% 
  mutate(TeamName = if_else(cat == "wins", str_c(WTeamName, " vs ", LTeamName, sprintf("  (%.1f %%)", win_perc)), "")) %>% 
  ggplot(aes(rank, n, fill = cat)) +
  geom_col() +
  geom_text(aes(x = rank, y = 0, label = TeamName), nudge_y = 1, hjust = "left", col = "black") +
  coord_flip() +
  scale_fill_manual(values = c("orange", "grey30")) +
  facet_grid(~cat, space = "free", scales = "free") +
  theme_hc() +
  theme(legend.position = "none", axis.text.y = element_blank(), axis.ticks.y = element_blank()) +
  labs(x = "", y = "", title = "", subtitle = "By winning percentage (> 50 games played)")

p1 / p2
```

We find:

- There's a rather one-sided local rivalry between "Kansas" and "Kansas St". The "Kansas" team also enjoys travelling to Colorado.

- On the flip side, "Gonzaga" is probably not too welcome in "San Diego". The neighbours matchup of "Boston Univ" vs "New Hampshire" makes it onto both lists.


Of course, a lot can change in 35 years, as we have seen in the time series of seeds above. Let's have a look at the same stats over the last 10 years only. Since we've spent quite a bit of thought on designing this layout let's get some more mileage out of it. This time we'll be looking at the winning percentage only, and will also be changing the colours to avoid confusion:


```{r fig.cap ="Fig. 12", fig.height = 7}
foo <- regular_res %>% 
  filter(Season >= 2009) %>% 
  count(WTeamID, LTeamID) %>% 
  arrange(desc(n)) %>% 
  left_join(teams %>% select(TeamID, TeamName), by = c("WTeamID" = "TeamID")) %>% 
  rename(WTeamName = TeamName) %>% 
  left_join(teams %>% select(TeamID, TeamName), by = c("LTeamID" = "TeamID")) %>% 
  rename(LTeamName = TeamName)

bar <- foo %>% 
  select(-ends_with("ID")) %>% 
  rename(wins = n) %>% 
  left_join(foo %>% select(-ends_with("ID")), by = c("WTeamName" = "LTeamName", "LTeamName" = "WTeamName")) %>% 
  rename(losses = n) %>% 
  replace_na(list(wins = 0, losses = 0)) %>% 
  mutate(win_perc = wins/(wins + losses) * 100)

bar %>% 
  filter(wins + losses > 20) %>% 
  group_by(WTeamName) %>% 
  arrange(desc(win_perc)) %>% 
  dplyr::slice(1) %>% 
  ungroup() %>% 
  arrange(desc(win_perc)) %>%
  head(20) %>%
  mutate(rank = seq(20,1)) %>%  
  pivot_longer(cols = c(wins, losses), names_to = "cat", values_to = "n") %>% 
  mutate(cat = fct_relevel(as.factor(cat), "losses", after = Inf)) %>% 
  mutate(TeamName = if_else(cat == "wins", str_c(WTeamName, " vs ", LTeamName, sprintf("  (%.1f %%)", win_perc)), "")) %>% 
  ggplot(aes(rank, n, fill = cat)) +
  geom_col() +
  geom_text(aes(x = rank, y = 0, label = TeamName), nudge_y = 1, hjust = "left", col = "black") +
  coord_flip() +
  scale_y_continuous(breaks = pretty_breaks(4)) +
  scale_fill_manual(values = c("lightblue", "grey30")) +
  facet_grid(~cat, space = "free", scales = "free") +
  theme_hc() +
  theme(legend.position = "none", axis.text.y = element_blank(), axis.ticks.y = element_blank()) +
  labs(x = "", y = "", title = "Top 20 One-sided matchups (last 10 years)",
       subtitle = "By winning percentage (> 20 games played)")
```

We find:

- "Gonzaga" is still in a prominent position when visiting "Pepperdine", and "Kansas" vs "Texas" is further down the list at 86%.

- We also see interesting local rivalries such as "North Carolina" vs "NC State" or "Princeton" vs "Columbia".

- With odds so heavily stacked, there's a significant psychological pressure on the "losing" side to avoid another defeat.


After all this intricate ggplot2 wrangling, let's have a simple plot to close this sub-section. Here is the number of unique teams per `Season` over the last 35 years:

```{r fig.cap ="Fig. 13", fig.height = 2}
regular_res %>% 
  select(Season, ends_with("ID")) %>% 
  pivot_longer(ends_with("ID"), names_to = "foo", values_to = "TeamID") %>% 
  distinct(Season, TeamID) %>% 
  count(Season) %>% 
  ggplot(aes(Season, n)) +
  geom_line(col = "blue") +
  theme_fivethirtyeight() +
  labs(title = "Number of Teams per Regular Season")
```

We find:

- the number of teams has increased quite a lot over time. From 282 in 1985 to 353 in 2019 is a 25% increase. Consequently, more games will be played during the regular season.

- **This tells us that we have to correct for this change when looking at raw numbers from `Seasons` that are further in the past.** Or use relative rates only.



## Historical Bracket results

After looking at the recorded results during the regular season, let's now check out the historical bracket results. This table has the same structure as the previous one; only I've already removed the winning location feature since it is always `WLoc = N` (neutral court).

Now the `DayNum` feature reflects how far each team managed to progress into the tournament. We're going to look at this aspect first. Here is a lookup table based on the tournament schedule described in the 2019 bracket visual above and in the data description. I'm adding a `rank` column to simplify plotting:

```{r}
stage_fct <- c("First Four", "First Round", "Second Round",
               "Regional Semifinals", "Regional Finals",
               "National Semifinals", "National Championship")

stage_fct_plus <- c("First Four", "First Round", "Second Round",
               "Regional Semifinals", "Regional Finals",
               "National Semifinals", "National Championship", "Champion")

bracket_stage <- tibble(
  DayNum = c(seq(134,139), seq(143, 146), 152, 154),
  stage = c(rep("First Four", 2), rep("First Round", 2), rep("Second Round", 2),
            rep("Regional Semifinals", 2), rep("Regional Finals", 2), "National Semifinals",
            "National Championship"),
  rank = c(7, 7, 6, 6, 5, 5, 4, 4, 3, 3, 2, 1)
) %>% 
  mutate(stage = fct_relevel(as.factor(stage), stage_fct))

bracket_stage %>% 
  head(3) %>% 
  kable() %>% 
  kable_styling()
```


Let's visualise the teams with the most appearances in the National Championship (i.e. the final game of the tournament):



```{r fig.cap ="Fig. 14", fig.height = 7}
foo <- bracket_res %>% 
  select(Season, DayNum, WTeamID, LTeamID) %>% 
  pivot_longer(ends_with("ID"), names_to = "win_lose", values_to = "TeamID") %>% 
  group_by(Season, TeamID) %>% 
  summarise(DayNum = max(DayNum)) %>% 
  ungroup() %>% 
  left_join(bracket_stage, by = "DayNum") %>% 
  left_join(teams %>% select(starts_with("Team")), by = "TeamID") 

top <- foo %>% 
  filter(stage == "National Championship") %>% 
  count(TeamName) %>% 
  arrange(desc(n)) %>% 
  filter(n > 2) %>% 
  mutate(top = n > 6)
  
p1 <- top %>% 
  ggplot(aes(reorder(TeamName, -n, FUN = min), n, fill = top)) +
  geom_col() +
  scale_x_discrete(labels = function(x) lapply(str_wrap(x, width = 10), paste, collapse="\n")) +
  scale_y_continuous(breaks = pretty_breaks()) +
  scale_fill_manual(values = c("orange", "red")) +
  theme_fivethirtyeight() +
  theme(legend.position = "none", title = element_text(size = 8)) +
  labs(x = "", y = "", title = "Teams with > 6 Championship appearances (all time)")
  
p2 <- foo %>% 
  filter(TeamName %in% c(top$TeamName, "Virginia", "Texas Tech") & Season >= 2010) %>% 
  filter(!(TeamName %in% c("Connecticut", "Syracuse"))) %>% 
  ggplot(aes(Season, stage, col = TeamName, group = TeamName)) +
  geom_line() +
  geom_point() +
  facet_wrap(~ TeamName) +
  scale_x_continuous(breaks = pretty_breaks()) +
  theme_fivethirtyeight() +
  theme(legend.position = "none", title = element_text(size = 8)) +
  labs(x = "", y = "", title = "Best bracket stage for top teams + 2019 finalists (last 10 years)")

p1 / p2 + plot_layout(heights = c(2, 5))
```

We find:

- "Duke" is holds the all time record of 9 Championship appearances. No other team has more than 5, with "Kansas", "Kentucky", and "Michigan" tied for 2nd place.

- Over the last 10 year, "Duke" has been present in every bracket. There is a lot of variance in their performance, with 2 first round exits and 2 championship appearances.

- The 2019 finalists did step up their game significantly compared to past years. This is especially true for "Texas Tech", who in the 2010s only made it past the 1st round in 2018 (regional finals exit).

*Iterating the tournament structure and nomenclature:* the "first round" (64 teams) is played on days 136 and 137, the "second round" (32 teams) on days 138/139, the "regional semifinals" (= "Sweet Sixteen", last 16 teams) happen on days 143/144, "regional finals" on days 145/146 (last 8; i.e national quarterfinals aka "Elite Eight"), the "national semifinals" on day 152 (those are also "Final Four" teams), and the final ("National Championship") always happens on day 154.


Could there be any quasi-cyclical patterns over a longer time range? Hear me out: There is a realistic scenario in sports where a coach builds a great team, they win a lot, then the best players leave, and success only comes back after the team has been re-built. Especially in college basketball, where the NBA is extremly tempting to any young prodigy, such a scenario is not a-priori ridiculous.

Let's look at the bracket performance of some of the strongest teams over a larger time range:

```{r fig.cap ="Fig. 15", fig.height = 5.5}
foo %>% 
  filter(TeamName %in% c("Duke", "Kansas", "Florida", "Kentucky")) %>% 
  ggplot(aes(Season, stage, col = TeamName, group = TeamName)) +
  geom_smooth(method = "loess", span = 1/3, formula = "y~x") +
  geom_line() +
  geom_point() +
  facet_wrap(~ TeamName, nrow = 2) +
  scale_x_continuous(breaks = pretty_breaks()) +
  theme_tufte() +
  theme(legend.position = "none") +
  labs(x = "", y = "", title = "Long-term trends for top teams (w/ loess smoother)")
```

We find:

- There are certainly some year-to-year trends in the data: "Florida" and "Kentucky" seem to be on a downward path. "Kansas" is pretty stable in reaching the last 16 or 8. So is "Duke", and both are always able to go all the way to the final.

- Could these patterns be the primary source of prediction? No, because there's lots of variability in the data. Could they be used as one of many input features? I don't see why not, as long as there's not too much overlap with other predictors.

- Are those patterns cyclical? There's not enough data to tell, unfortunately. "Kentucky" or "Florida" look tempting, but we have at best 2 cycles and that's way to little to detect periodicity. And for "Kansas": those high-frequency oscillations after 2000 might look tempting at first, but there's no telling whether they're artefacts.

- One could try to link the variability to other features, or maybe apply domain knowledge. But other than that I'd be very careful in extrapolating from them.


Similar to the regular season overtime stats above, let's check out how teams perform in overtime during the bracket:

```{r fig.cap ="Fig. 16", fig.height=2.5}
prep <- bracket_res %>% 
  filter(NumOT > 0) %>% 
  select(WTeamID, LTeamID) %>% 
  pivot_longer(ends_with("ID"), names_to = "cat", values_to = "TeamID") %>% 
  left_join(teams %>% select(starts_with("Team")), by = "TeamID") %>% 
  mutate(cat = if_else(cat == "WTeamID", "win", "lose")) %>% 
  group_by(TeamName) %>% 
  summarise(total = n(),
            win_perc = mean(cat == "win")*100) %>% 
  filter(total > 5) %>% 
  pivot_longer(cols = c(total, win_perc), names_to = "type", values_to = "val")

prep %>% 
  group_by(type) %>%
  arrange(desc(val)) %>% 
  dplyr::slice(c(1:5)) %>% 
  ungroup() %>% 
  mutate(type = if_else(type == "total", "Total Wins", "Winning Percentage")) %>% 
  ggplot(aes(reorder(TeamName, val, FUN = min), val, fill = type)) +
  geom_col() +
  coord_flip() +
  theme_fivethirtyeight() +
  theme(legend.position = "none") +
  facet_wrap(~ type, scales = "free") +
  labs(x = "", y = "", title = "Highest Overtime stats during bracket (all time)")
```

We find:

- "Kentucky" and "Florida" have both needed the stamina to play 8 overtimes during their bracket appearances, but of the two only "Florida" managed to win more than 50% of those.

- For the winning percentage stats we only considered teams with more than 5 overtime appearances. "Michigan" had 6 and managed to win them all - an impressive performance under the added pressure of the elimination rounds.


## Combining Basic Data features

Ok, we got ourselves a good overview of the datasets in our Basic Data group, and explored their various features. Now we will put this information together and study some multi-feature correlations.

First off, let's link the performance during the regular season to the initial bracket seeds. The complete [seeding process](https://en.wikipedia.org/wiki/NCAA_basketball_tournament_selection_process) is a rather complex mix of performance indicators and scheduling constraints; for instance the selection comittee will try to avoid repeating match-ups from the regular season or last year's bracket in the early rounds ([source](https://en.wikipedia.org/wiki/NCAA_Division_I_Men%27s_Basketball_Tournament#Seeding_and_bracket)). Let's see whether we can boil this down to a few, robust main criteria.

Expand this code cell to see how I prepare the regular season data to extract win/loss score margins, win margins, and join the different data frames together.

```{r}
# preparation df with score margin per game + win/loss indicator
prep <- regular_res %>% 
  mutate(score_margin = WScore - LScore) %>% 
  select(Season, ends_with("TeamID"), score_margin) %>% 
  pivot_longer(ends_with("ID"), names_to = "cat", values_to = "TeamID") 

# win/loss margin per Season and TeamID
foo <- prep %>% 
  group_by(Season, TeamID, cat) %>% 
  summarise(mean_margin = mean(score_margin)) %>% 
  ungroup() %>% 
  mutate(wl = if_else(cat == "WTeamID", "wins", "losses")) %>% 
  mutate(mean_margin = if_else(wl == "losses", -mean_margin, mean_margin)) %>% 
  select(-cat)

# extract win percentage, join with prep + foo + seeds df
wins_seeds <- prep %>%
  group_by(Season, TeamID) %>% 
  summarise(win_perc = mean(cat == "WTeamID")) %>% 
  ungroup() %>% 
  left_join(teams %>% select(starts_with("Team")), by = "TeamID") %>% 
  left_join(seeds, by = c("TeamID", "Season")) %>% 
  filter(!is.na(Seed)) %>% 
  mutate(Seed = as.integer(str_sub(Seed, 2, 3))) %>%
  mutate(decade = str_c(as.character(Season %/% 10 * 10), "s")) %>% 
  left_join(foo, by = c("Season", "TeamID"))
```


Here we plot the winning percentage during the regular Season on the x-axis vs the initial bracket seeds on the y-axis. We break up the plot by decade (facet + colour) to avoid overcrowding:

```{r fig.cap ="Fig. 17", fig.height = 5.5}
 wins_seeds %>% 
  ggplot(aes(win_perc, Seed, col = decade)) +
  geom_point() +
  scale_x_continuous(labels = scales::percent) +
  scale_y_reverse(breaks = seq(1, 15, 2)) +
  facet_wrap(~ decade) +
  theme_fivethirtyeight() +
  theme(legend.position = "none") +
  labs(x = "", y = "", col = "Decade", title = "Winning percentage vs Seeds",
       subtitle = "Bracket seeding broadly correlated with % regular season wins.")
```

We find:

- There's a broad correlation of seeding with winning percentage: generally you need more than 80% wins to become a 1st seed team, and no seed above 6 or 7 had less than 60% wins. This picture is pretty consistent throughout the decades. (Remember that there are 4 teams for each seed; corresponding to the 4 regions).

- However, there is also a lot of scatter. Some teams with > 80% wins where placed as low as seed 15. The reasons for this might have something to do with scheduling constraints and winning margins (we get to the latter in a moment), but also with the fact that teams performances will be ranked relative to their competitors.

This means that for instance if everyone else had a mediocre year with 70% wins, then 80% can get you far. On the other hand, if the season was more polarised, with a lot of teams approaching 85-90% winning stats, then your 80% will be ranked lower.

Let's investigate this thought further and look at performance within a season. In the following plot, the left-side panels are related to each other, and the same holds for the right-side panels. The visual shows the results of applying a simple linear model to predict the seed from the winning percentage (left) and winning % together with winning margin (right) for each individual season. The top panels always show some examples for the 3 more recent seasons, and the bottom panels show the distribution of R-square values for the linear models.

(From a coding perspective, we use a combination of `group_by`, `nest`, and `purrr`'s `map_dbl` to extract model R-squares per Season in a vectorised way.)

```{r fig.cap ="Fig. 18", fig.height = 5.5}
p1 <- wins_seeds %>% 
  filter(Season >= 2017) %>% 
  ggplot(aes(win_perc, Seed), col = "blue") +
  geom_smooth(method = "lm", se = FALSE, linetype = 2, col = "black", formula = 'y ~ x') +
  geom_point(col = "blue") +
  scale_x_continuous(labels = scales::percent, breaks = c(0.5, 0.8)) +
  scale_y_reverse(breaks = seq(1, 16, 3)) +
  facet_wrap(~ Season, nrow = 1) +
  theme_hc() +
  theme(legend.position = "none", title = element_text(size = 10)) +
  labs(x = "", y = "", title = "Win % vs Seeds: marginal correlation",
       subtitle = str_c("Lower panel: Linear model ", as.character(expression(R^2)), " all seasons."))
  
p2 <- wins_seeds %>% 
  group_by(Season) %>% 
  nest() %>% 
  mutate(r_squared = map_dbl(data, ~summary(lm(.x$Seed ~ .x$win_perc))$r.squared)) %>% 
  ggplot(aes(r_squared)) +
  geom_density(fill = "blue", bw = 0.01) +
  theme_tufte() +
  labs(x = expression(R^2), y = "")

binsize <- 10

p3 <- wins_seeds %>% 
  filter(Season >= 2017 & wl == "wins") %>% 
  arrange(mean_margin) %>%
  mutate(margin_group = case_when(
    mean_margin < 12 ~ "< 12",
    between(mean_margin, 12, 14) ~ "< 12 - 14",
    between(mean_margin, 14, 18) ~ "< 14 - 18",
    mean_margin > 18 ~ "> 18"
  )) %>% 
  # mutate(margin_bin = mean_margin %/% binsize,
  #        margin_group = str_c(sprintf("%02i", binsize * margin_bin + 1), ' to ', binsize * (margin_bin + 1))) %>% 
  ggplot(aes(win_perc, Seed, col = margin_group)) +
  geom_point(size = 2.5) +
  scale_x_continuous(labels = scales::percent, breaks = c(0.5, 0.8)) +
  scale_y_reverse(breaks = seq(1, 16, 3)) +
  facet_wrap(~ Season, nrow = 1) +
  theme_hc() +
  theme(legend.position = "top", title = element_text(size = 10),
        legend.text = element_text(size = 7)) +
  labs(x = "", y = "", col = "", title = "Winning margin doesn't help",
       subtitle = "Margin binned into 3 classes: ")
  
p4 <- wins_seeds %>% 
  group_by(Season) %>% 
  nest() %>% 
  mutate(r_squared = map_dbl(data, ~summary(lm(.x$Seed ~ .x$win_perc + .$mean_margin))$r.squared)) %>% 
  ggplot(aes(r_squared)) +
  geom_density(fill = "orange", bw = 0.01) +
  theme_tufte() +
  labs(x = expression(R^2), y = "")

((p1 / p2)) | ((p3 / p4))
```

We find:

- There is a marginal relation between regular season winning percentage and seeding rank. The black dashed line show a naive linear model, but even those are primarily affected by the most extreme percentages for seeds 1 and 16. The R-squared values are between 0.1 and 0.4, which is not much to write home about.

- Adding the winning margin as an additional parameter doesn't really help. If you can't see a difference between the R-squared values then that's because there isn't really any. In the scatter plot, I'm binning the mean winning margin into 4 hand-designed classes to avoid having too many colours on the plot. (The linear model uses the un-binned mean margin). You see that between ranks 6 and 13 there is a lot of scatter of the different bins.


Our prediction task, however, is more related to the teams performances during the bracket than to the seeds themselves. Still, it's good to have a feeling for how the seeds come to be. Let's look at the bracket performance by using the ranking proxies we built above. Here we plot the highest bracket stage that the teams reached for each Season since 2000:


```{r}
prep <- bracket_res %>% 
  select(Season, DayNum, WTeamID, LTeamID) %>% 
  pivot_longer(ends_with("ID"), names_to = "win_lose", values_to = "TeamID") %>% 
  left_join(bracket_stage, by = "DayNum") %>% 
  mutate(stage = as.character(stage)) %>% 
  mutate(rank = if_else( (stage == "National Championship" & win_lose == "WTeamID"), 0, rank ),
         stage = if_else(rank == 0, "Champion", stage)) %>% 
  mutate(stage = fct_relevel(as.factor(stage), stage_fct_plus)) %>% 
  group_by(Season, TeamID) %>% 
  summarise(rank = min(rank)) %>% 
  ungroup() %>% 
  left_join(bracket_stage, by = "rank") %>% 
  mutate(stage = as.character(stage)) %>%
  replace_na(list(stage = "Champion")) %>% 
  mutate(stage = fct_relevel(as.factor(stage), stage_fct_plus))
```


```{r fig.cap ="Fig. 19", fig.height = 6.5}
wins_seeds %>% 
  filter(wl == "wins" & Season >= 2000) %>% 
  left_join(prep, by = c("TeamID", "Season")) %>% 
  filter(!is.na(rank)) %>% 
  ggplot(aes(win_perc, stage, col = as.factor(Season))) +
  geom_point() +
  scale_x_continuous(labels = scales::percent, breaks = c(0.5, 0.7, 0.9)) +
  facet_wrap(~ Season) +
  theme_hc() +
  theme(legend.position = "none") +
  labs(x = "", y = "", title = "Best Bracket performance correlates weakly with win %")
```

We find:

- There is not much correlation here either. Finalists and Winners don't even always have a strong winning percentage (see e.g. 2011 or 2014).

- This will always be the problem with elimination style tournaments, where a single bad performance can break your run.


## Head to Head predictions - Case Study

So let's get closer to our prediction goal by looking at individual head-to-head comparisons and how the features we have found so far reflect the outcome. This is a simple case study that is intended to illustrate what we are aiming at; not an exhaustive analysis.

We're randomly picking 2 teams with a long and somewhat imbalanced game history:

```{r}
sel_teams <- regular_res %>% 
  group_by(WTeamID, LTeamID) %>% 
  count() %>% 
  ungroup() %>% 
  arrange(desc(n)) %>% 
  dplyr::slice(15) %>% 
  select(-n) %>% 
  pivot_longer(everything(), names_to = "foo", values_to = "teams") %>% 
  pull(teams)

foo <- regular_res %>% 
  filter( (WTeamID == sel_teams[1] & LTeamID == sel_teams[2]) | (WTeamID == sel_teams[2] & LTeamID == sel_teams[1]) ) %>% 
  count(WTeamID) %>% 
  rename(wins = n, TeamID = WTeamID)

h2h_teams <- teams %>% 
  filter(TeamID %in% sel_teams) %>% 
  left_join(foo, by = "TeamID") 

h2h_teams %>% 
  kable() %>% 
  kable_styling()
```

Here is an overview of their encounters: We plot the Season against the score differences. We choose a simple scatter plot and augment it with a facetted and flipped histogram plot to emphasise the spread of scores. This is similar to a marginal histogram, but the facets here add some clarity (instead of having overlapping distributions) and the facet heads provide the colour legend:

```{r fig.cap ="Fig. 20", fig.height = 3.5}
foo <- regular_res %>% 
  filter( (WTeamID == sel_teams[1] & LTeamID == sel_teams[2]) | (WTeamID == sel_teams[2] & LTeamID == sel_teams[1]) ) %>% 
  mutate(h2h_diff = WScore - LScore) %>% 
  select(-DayNum, -NumOT, -WScore, -LScore) %>%
  pivot_longer(ends_with("ID"), names_to = "win_lose", values_to = "TeamID") %>% 
  filter(win_lose == "WTeamID") %>% 
  left_join(h2h_teams %>% select(TeamID, TeamName), by = "TeamID") 

p1 <- foo %>% 
  ggplot(aes(Season, h2h_diff, col = TeamName)) +
  #geom_line(col = "grey50") +
  geom_point(size = 3) +
  theme_hc() +
  theme(legend.position = "none", title = element_text(size = 10)) +
  labs(x = "", y = "", title = "Case study: Gonzaga beats Santa Clara since 2000",
       subtitle = "Both y-axes show winner score - loser score.")

p2 <- foo %>% 
  ggplot(aes(h2h_diff, fill = TeamName)) +
  geom_histogram(bins = 30, alpha = 1) +
  coord_flip() +
  facet_grid(~ TeamName) +
  theme_hc() +
  scale_y_continuous(breaks = seq(0,15,3)) +
  theme(legend.position = "none") +
  labs(x = "", y = "")
  

p1 + p2 + plot_layout(widths = c(7, 3))
```

We find:

- While in the 80s and 90s "Santa Clara" was winning most of their games, and the score differences were relative similar, this all changed around the year 2000. "Gonzaga" started dominating and their winning leads grew. Since 2010 a "Santa Clara" win has become rare, and I would advice betting against "Gonzaga" in their current match ups.


And now it comes down to a classification problem. We label the games that "Gonzaga" won with `1` and the ones they lost with `0`. Then we compare the distributions of some of the features we found to be interesting in our global analysis of the dataset. So we're collecting the average statistics for score, winning difference, and losing difference, together with the winning percentage for both of our teams.

If you unfold the code block you'll find a little trick on running `pivot_longer` for more than 2 pairs of variables using the the `names_pattern` parameter.

From a plotting perspective, here we reach deeper into the rich feature set of the [patchwork](https://github.com/thomasp85/patchwork) package, and give the shared legend its own quadrant in the layout. We also add a global title and subtitle:

```{r fig.cap ="Fig. 21", fig.height = 5.5}
gu_sc <- regular_res %>% 
  filter( (WTeamID == sel_teams[1] & LTeamID == sel_teams[2]) | (WTeamID == sel_teams[2] & LTeamID == sel_teams[1]) ) %>% 
  mutate(target = if_else(WTeamID == 1211, 1, 0)) %>% 
  select(Season, target)

gu <- regular_res %>% 
  filter(WTeamID == sel_teams[1] | LTeamID == sel_teams[1]) %>% 
  select(Season, ends_with("ID"), ends_with("Score")) %>% 
  mutate(diff_score = WScore - LScore) %>% 
  rename(W_TeamID = WTeamID, L_TeamID = LTeamID, W_Score = WScore, L_Score = LScore) %>% 
  pivot_longer(matches("^[WL]."), names_to = c("win_lose", ".value"), names_pattern = "(.)_(.+)") %>% 
  filter(TeamID == sel_teams[1]) %>% 
  mutate(diff_score = if_else(win_lose == "L", -diff_score, diff_score)) %>% 
  group_by(Season) %>% 
  summarise(diff_score = mean(diff_score),
            win_lose = mean(win_lose == "W"),
            mean_score = mean(Score))

p1 <- gu_sc %>% 
  left_join(gu, by = "Season") %>% 
  ggplot(aes(diff_score, fill = as.factor(target))) +
  geom_density(alpha = 0.5) +
  labs(x = "", y = "", title = "Score Difference", fill = "Gonzaga Won") +
  theme_hc()

p2 <- gu_sc %>% 
  left_join(gu, by = "Season") %>% 
  ggplot(aes(win_lose, fill = as.factor(target))) +
  geom_density(alpha = 0.5) +
  scale_x_continuous(labels = scales::percent) +
  labs(x = "", y = "", title = "Winning Percentage", fill = "Gonzaga Won") +
  theme_hc()

p3 <- gu_sc %>% 
  left_join(gu, by = "Season") %>% 
  ggplot(aes(mean_score, fill = as.factor(target))) +
  geom_density(alpha = 0.5) +
  labs(x = "", y = "", title = "Mean Score", fill = "Gonzaga Won") +
  theme_hc()

p3 + p1 + p2 + guide_area() + 
  plot_layout(guides = 'collect') +
  plot_annotation(title = 'Gonzaga vs Santa Clara results reflect global Gonzaga stats',
                  subtitle = 'Showing distributions of 3 average stats per season for Gonzaga win/loss vs Santa Clara')
```

We find:

- Now the distributions make more sense: in seasons where "Gonzaga" won against "Santa Clara" they had significantly higher average scores, winning percentages, and score differences. Thus the season's performance as measured by those 3 parameters allows for a reasonably good estimate on whether "Gonzaga" would win. (Those statistics also include the games played against "Santa Clara", but that's just 2 data points among many.)

- Although the two distributions (win vs loss) are different, they also have a notable overlap. Part of this is the remaining uncertainty that we would need to reduce using other features; e.g. ones that take into account the specific characteristics of the "Gozanga" and "Santa Clara" teams. We'll get more into this in the Additional Data sections below.

- Note, that another part of the overlap is due to there being more than 1 game played per season between two teams. As we saw in the scatter plot prior, the results of these games may well be different (e.g. 1 win & 1 loss) but the season's statistics are the same. This also can be addressed using additional features that change from game to game within a season.


## Understanding the Submission Format

Before moving on to a baseline model, let's take a closer look at the submission format and make sure that we know exactly what to predict in the 1st and 2nd stage of this competition. **Only the stage 2 results matter for the final leaderboard.**

It is important to note that the stage 1 predictions use a test set from the 2015-19 bracket stages. **These results are part of the training data** in the `MNCAATourneyCompactResults.csv` file (called "Bracket Results" in this notebook). It is therefore trivial to get a "perfect score" on the stage 1 leaderboard by simply submitting the correct results; but it is also completely pointless. The function of the test data set in the 1st stage is to help evaluate our model performances. The models will ultimately have to predict the results of the unseen 2020 bracket (stage 2) and should be designed with this goal in mind.

Apart from the obvious, willful misuse of 2015-19 bracket data to predict the 1st stage test set, there is also the more subtle danger of stage 1 groundtruth leaking into our models. This would lead to overfitting and bad performance on the important stage 2 data. As has been pointed out in the [discussions](https://www.kaggle.com/c/google-cloud-ncaa-march-madness-2020-division-1-mens-tournament/discussion/131028), the winning scores from previous years were in the range of just below 0.5 (same metric). Any score much lower than that is likely due to overfitting.

The format of the stage 1 and stage 2 submission files will be the same, so let's see what's expected from us. In this competition, there is no separate test data set, but the sample submission file contains all the info we need. Here are the first 5 lines of that sample submission file (`MSampleSubmissionStage1_2020.csv`):


```{r}
sample_submit %>% 
  head(5) %>% 
  kable() %>% 
  kable_styling()
```

There are only 2 columns:

- The prediction column `Pred` contains only dummy entries of 0.5 (i.e. 50% winning chance). These are the probabilities we will need to predict.

- The `ID` column tells us what these probabilities correspond to. It contains 3 values: Season, TeamID1, TeamID2. Here, in stage 1, Season is in 2015-19. The entries are sorted by those 3 values in that order; i.e. **TeamID1 < TeamID2 and there are no duplicates** (i.e. only `1107_1112` but not `1112_1107`). It's not necessary to have both combinations in the data, because in a 2 team match the winning probability for one is always 1 - winning probability for the other.

- The combination of TeamIDs cover every possible game between teams that made the bracket for that year. For each season in 2015-19, and also for 2020, this means 68 teams in total.

For more details and some examples see the [data description](https://www.kaggle.com/c/google-cloud-ncaa-march-madness-2020-division-1-mens-tournament/data). 

So, is the stage 1 leaderboard useful at all? Not really for model validation because the true results are part of the training data and you can just check them locally. But it is useful for trying out the submission format and making sure that your predictions can be properly read by the Kaggle scoring system. You don't want to carefully build a model but then get the submission format wrong and error out of the competition.

The final stage 2 data, including an updated sample submission file, will be available after Selection Sunday (March 15).


PS: I can now happily show you the effect of overfitting on the leaderboard, because I accidentally submitted a toy prediction that was wildly overfit; instead of the non overfit toy prediction I wanted to submit. Not even gonna pretend this is some weird strategic move ... . This is what happens if you give cryptic names to your submission files. Still, there might be another teaching moment here (not just on proper file names for me): this 0.3 score used all the regular season + bracket data but only very few features. It's trivial to overfit if you're training on the literal ground truth. And all the results below about 0.49 are similarly overfit. Just in case you didn't believe me before.

Anyway, moving on ...


## Baseline Model

The modelling approach we're taking will be very similar to the Python [starter Kernel](https://www.kaggle.com/hiromoon166/2020-basic-starter-kernel) contributed this year by [hiromu](https://www.kaggle.com/hiromoon166). A Kernel which, by the way, [traces](https://www.kaggle.com/addisonhoward/basic-starter-kernel-ncaa-men-s-dataset-2019) its [lineage](https://www.kaggle.com/kplauritzen/notebookde27b18258) all the way back to a [2016 R script](https://www.kaggle.com/jaredcross/getting-started) by [Jared Cross](https://www.kaggle.com/jaredcross).

The basic design has stayed the same: we compute our features for each team, then join them to the historical bracket data on team ID and season. To express the win & loss conditions, we design 2 data frames: one from the perspective of the winning team and the other from the losing team. Then we stack them together to get our training data.

The features we will use are:

- winning percentages from the regular season data, and the differences with respect to the same value for the other team (i.e. win percentage team 1 - win percentage team 2)

- avarage scores from the regular season data, and difference wrt other team

- bracket seeds and seed difference

Since we want even this baseline model to be reasonably realistic, we will not include bracket results or regular season data after 2014 (the former would be the ground truth). (Regular season results in 2015 would not have the ground truth, but there's no bracket season to join them to, so we couldn't use them in this setting. More on this below.)

```{r}
base_regular <- regular_res %>% 
  filter(Season <= 2014) %>% 
  select(Season, ends_with("ID"), ends_with("Score"))

base_bracket <- bracket_res %>% 
  filter(Season <= 2014) %>% 
  select(Season, WTeamID, LTeamID) %>% 
  left_join(seeds, by = c("Season", "WTeamID" = "TeamID")) %>% 
  rename(Seed1 = Seed) %>% 
  left_join(seeds, by = c("Season", "LTeamID" = "TeamID")) %>% 
  rename(Seed2 = Seed) %>% 
  mutate(Seed1 = as.numeric(str_sub(Seed1, 2, 3)),
         Seed2 = as.numeric(str_sub(Seed2, 2, 3)))
```


We extract the regular season's mean score and winning percentage for each team in one fell swoop, then join them onto the bracket data:

```{r}
foo <- base_regular %>% 
  rename(W_TeamID = WTeamID, L_TeamID = LTeamID, W_Score = WScore, L_Score = LScore) %>% 
  pivot_longer(matches("^[WL]."), names_to = c("win_lose", ".value"), names_pattern = "(.)_(.+)") %>% 
  group_by(Season, TeamID) %>% 
  summarise(win_percentage = mean(win_lose == "W"),
            mean_score = mean(Score)) %>% 
  ungroup()

train_prep <- base_bracket %>% 
  left_join(foo, by = c("Season", "WTeamID" = "TeamID")) %>% 
  rename(win_percent_1 = win_percentage, mean_score_1 = mean_score) %>% 
  left_join(foo, by = c("Season", "LTeamID" = "TeamID")) %>% 
  rename(win_percent_2 = win_percentage, mean_score_2 = mean_score) %>% 
  mutate(win_diff = win_percent_1 - win_percent_2,
         score_diff = mean_score_1 - mean_score_2,
         seed_diff = Seed1 - Seed2) %>% 
  select(ends_with("1"), ends_with("2"), ends_with("diff"))

train_prep %>% 
  head(5) %>% 
  kable() %>% 
  kable_styling()
```


This leaves us with a data frame from the point of view from the winning team; i.e. all the `target` results are `1`. We then flip those numbers to include the perspective of the losing team (where `target = 0`). This creates the following 2 data frames (first 5 rows shown):

```{r}
win <- train_prep %>% 
  select(-ends_with("2")) %>% 
  # mutate(target = "win")
  mutate(target = 1)

lose <- train_prep %>% 
      select(-ends_with("1")) %>% 
      rename_all( ~ str_replace(., "2", "1")) %>% 
      mutate_at(vars(ends_with("diff")), function(x) {-x}) %>% 
      # mutate(target = "lose")
      mutate(target = 0)

win %>% head(5) %>% 
  kable() %>% 
  kable_styling()
```

```{r}
lose %>% head(5) %>% 
  kable() %>% 
  kable_styling()
```

As you can see, they are mirror images of each other, in that the 1st has all the wins and the second has all the losses. The `diff` columns are negative, and the `_1` columns are relative to the winning vs losing team. Now we bind those two tables together to create our training data. Here is a summary:

```{r}
train <- win %>% 
  bind_rows(lose)
#  mutate(target = as.factor(target))

train %>% 
  summary()
```


For modelling we will be using the trusty [XGBoost](https://xgboost.ai/), which this week announced its [version 1.0](https://twitter.com/TerryTangYuan/status/1230501041365561344)! 

The XGBoost parameters are defined in this code chunk here. We're keeping it very basic, given the limited number of features so far. This is just a starting point for you to add features and tinker with hyper-parameters. Importantly, it only includes data from before 2015 and is therefore free of leakage.

Note, that because we're using trees we're not scaling or normalising the input features. If you want to use a different model then that's something you need to keep in mind.


```{r}
dtrain <- xgb.DMatrix(as.matrix(train %>% select(-target)),label = train$target)

xgb_params <- list(colsample_bytree = 1., #variables per tree 
                   subsample = 1., #data subset per tree 
                   booster = "gbtree",
                   max_depth = 4, #tree levels
                   eta = 0.01, #learning rate
                   eval_metric = "logloss", 
                   objective = "binary:logistic",
                   early_stopping_rounds = 10,
                   seed = 4321
                   )

watchlist <- list(train=dtrain)
```


We're training for 500 rounds. I'm setting a seed to keep the results stable over future Kernel runs.

```{r}
set.seed(4321)
xgb_fit <- xgb.train(params = xgb_params,
                     data = dtrain,
                     print_every_n = 100,
                     watchlist = watchlist,
                     nrounds = 500)
```

Training log loss doesn't look super overfit, which is reassuring.

Those are the resulting feature importances:

```{r fig.cap ="Fig. 22", fig.height = 3.5}
imp_matrix <- as_tibble(xgb.importance(feature_names = colnames(train %>% select(-target)), model = xgb_fit))

imp_matrix %>%
  ggplot(aes(reorder(Feature, Gain, FUN = max), Gain, fill = Feature)) +
  geom_col() +
  coord_flip() +
  theme_fivethirtyeight() +
  theme(legend.position = "none") +
  labs(x = "", y = "Importance", title = "Feature Importance for Baseline Model")
```

We find:

- The seed difference is by far the most important feature. This reflects the general trend that the top seeds are more likely to go far.

- Of the other features, the difference in mean score for regular season games has the most impact. The absolute features (e.g. mean score or winning percentage) are always less important that the corresponding relative differences; but they appear to have an additional small impact nonetheless.


In order to make predictions on the test file, we need to prepare the same features for the test period. This is not entirely trivial, and it is worth spending a couple of thoughts on how it can be scaled to a larger feature set. Here we first extract the features for the regular season & bracket stage separately (in long format) and then join them together in a shape that has unique entries for each Season and TeamID:

```{r}
prep_test_regular <- regular_res %>% 
  filter(Season > 2014) %>% 
  select(Season, ends_with("ID"), ends_with("Score")) %>% 
  rename(W_TeamID = WTeamID, L_TeamID = LTeamID, W_Score = WScore, L_Score = LScore) %>% 
  pivot_longer(matches("^[WL]."), names_to = c("win_lose", ".value"), names_pattern = "(.)_(.+)") %>% 
  group_by(Season, TeamID) %>% 
  summarise(win_percentage = mean(win_lose == "W"),
            mean_score = mean(Score)) %>% 
  ungroup()

prep_test_bracket <- bracket_res %>% 
  filter(Season > 2014) %>% 
  select(Season, WTeamID, LTeamID) %>% 
  left_join(seeds, by = c("Season", "WTeamID" = "TeamID")) %>% 
  rename(W_Seed = Seed) %>% 
  left_join(seeds, by = c("Season", "LTeamID" = "TeamID")) %>% 
  rename(L_Seed = Seed, W_TeamID = WTeamID, L_TeamID = LTeamID) %>% 
  pivot_longer(matches("^[WL]."), names_to = c("win_lose", ".value"), names_pattern = "(.)_(.+)") %>% 
  mutate(Seed = as.numeric(str_sub(Seed, 2, 3))) %>% 
  distinct(Season, TeamID, Seed) %>% 
  arrange(Season, TeamID)


prep_test <- prep_test_regular %>% 
  left_join(prep_test_bracket, by = c("Season", "TeamID")) %>% 
  filter(!is.na(Seed))

prep_test %>% 
  head(5) %>% 
  kable() %>% 
  kable_styling()
```


Then we extract the corresponding `Season` and `TeamID1`, `TeamID2` from the `ID` column of the sample submission table, and join our features. The names need to be the same as for the training data.

```{r}
test_df <- sample_submit %>% 
  separate(ID, into = c("Season", "TeamID1", "TeamID2")) %>% 
  select(-Pred) %>% 
  mutate_if(is.character, as.numeric) %>% 
  left_join(prep_test, by = c("Season", "TeamID1" = "TeamID")) %>% 
  rename(win_percent_1 = win_percentage, mean_score_1 = mean_score, Seed1 = Seed) %>% 
  left_join(prep_test, by = c("Season", "TeamID2" = "TeamID")) %>% 
  rename(win_percent_2 = win_percentage, mean_score_2 = mean_score, Seed2 = Seed) %>% 
  mutate(win_diff = win_percent_1 - win_percent_2,
         score_diff = mean_score_1 - mean_score_2,
         seed_diff = Seed1 - Seed2) %>% 
  select(-Season, -starts_with("Team"), -ends_with("2"))

summary(test_df)
```


And then we predict on the test set and put the results into our submission file format. Remember that XGBoost expects matrices as input for train and predict, and that those matrices should have the same column names in the same order. Here are the first 5 rows of the submission file:

```{r}
dtest <- xgb.DMatrix(as.matrix(test_df %>% select(colnames(dtrain))))

submit_df <- sample_submit %>% 
  mutate(Pred = predict(xgb_fit, dtest))
  
submit_df %>% 
  head(5) %>% 
  kable() %>% 
  kable_styling()
```

We then write this table to an output file that we can directly submit from the Kernel:

```{r}
submit_df %>% 
  write_csv("submission.csv")
```



Given all that's been said about the Stage 1 ground truth being contained in our training data, it won't be a surprise that we can build our own Stage 1 testing setup, without needing the leaderboard. See also the discussion [here](https://www.kaggle.com/c/google-cloud-ncaa-march-madness-2020-division-1-mens-tournament/discussion/131539). For this we take our submitted probabilities and compare them with the match results from the bracket results data.

To compute the log loss we use a function from the [tidymodels](https://github.com/tidymodels) package [yardstick](https://tidymodels.github.io/yardstick/):


```{r}
foo <- bracket_res %>% 
  filter(Season >= 2015) %>% 
  select(Season, WTeamID, LTeamID, DayNum)

bar <- sample_submit %>% 
  separate(ID, into = c("Season", "TeamID1", "TeamID2")) %>% 
  select(-Pred) %>% 
  mutate_if(is.character, as.numeric) %>%
  left_join(foo, by = c("Season", "TeamID1" = "WTeamID", "TeamID2" = "LTeamID")) %>% 
  rename(win = DayNum) %>% 
  left_join(foo, by = c("Season", "TeamID1" = "LTeamID", "TeamID2" = "WTeamID")) %>% 
  rename(lose = DayNum) %>% 
  mutate(truth = as.factor(case_when(
    !is.na(win) ~ 1,
    !is.na(lose) ~ 0,
    TRUE ~ NA_real_
  )))

submit_df %>% 
  mutate(truth = bar %>% pull(truth)) %>% 
  yardstick::mn_log_loss(truth = fct_rev(truth), Pred) %>% 
  kable() %>% 
  kable_styling()
```

Here, the column `.estimate` holds our result. There's certainly room for improvement, but that's not a big surprise given that we only use a few features and a single model with very basic hyperparameters.

If you want to compute the loss yourself, then that's pretty straight-forward too: `loss = truth * log(Pred) + (1 - truth) * log(1 - Pred)`. Here is a code chunk that does just that and of course leads to the same result:


```{r}
foobar <- submit_df %>% 
  mutate(truth = as.integer(bar$truth)-1) %>% 
  filter(!is.na(truth)) 

foobar %>% 
  mutate(loss = truth * log(Pred) + (1 - truth) * log(1 - Pred)) %>% 
  summarise(loss = -sum(loss, na.rm = TRUE)/nrow(foobar)) %>% 
  kable() %>% 
  kable_styling()
```

So this is a baseline, without leakage, that you can use to measure more sophisticated models against.

Before we move on with further exploration, here's a couple of thoughts on the modelling:

- Could we have used seeds and seed differences from the 2015-2019 bracket data? Yes and no. We will have the 2020 seeds for our stage 2 predictions; so using the 2015 seeds would have been fine. However, including the 2016-19 seeds as well would give us some information about the future performance of teams. For instance 2017 seeds could have been used to not only predict the 2017 results but also the 2015 and 2016 results, thus leading to leakage. It's generally bad practice to include information about the future when predicting the future.

- With that in mind, it's worth thinking about the eventual goal of the competition and whether we're currently making the best use of the stage 1 predictions.

- I'm not using a validation set / validation loss when fitting. I'm doing that on purpose to encourage you to think about validation strategies, especially in the current setup where we're essentially counting each match twice. It's worth thinking about whether there's a danger of leakage between random sets or folds.


# Additional Data Visuals

Here we explore all the remaining data files that come bundled with this competition, and think about how to incorporate them in our analysis. In the [data description](https://www.kaggle.com/c/google-cloud-ncaa-march-madness-2020-division-1-mens-tournament/data) those files are categorised into Data Sections. We will follow these categories.

*Note, that I will simultaneously add material to the Overview section above.*


## Team Box Scores - Detailed Results

From 2003 onwards we have much more detailed information on game-level stats for both regular season and bracket. Let's start with visualising some regular season features.

We are given the number of shots made ("WFG"), as well as attempted ("WFGA"), for all field goals & for 3-pointers only (suffix "3"). Note, that the "WFGM" and "WFGA" features include both 2-pointer & 3-pointers. To measure the 2-pointers only, we subtract "WFGM3" from "WFGM" (same for attempts). Then we divide shots made by shots attempted to arrive at a shot percentage. Here we look at the distributions over all seasons using boxplots:

```{r fig.cap ="Fig. 23", fig.height = 4.5}
foo <- regular_detail %>% 
  select(Season, starts_with("WFG")) %>% 
  mutate(WFGM = WFGM - WFGM3,
         WFGA = WFGA - WFGA3) %>% 
  mutate(WFGR = WFGM/WFGA,
         WFGR3 = WFGM3/WFGA3) %>% 
  select(-matches("GM")) %>% 
  pivot_longer(starts_with("WFG"), names_to = "type", values_to = "points") %>% 
  mutate(shot = if_else(str_detect(type, "3"), "3-pointers", "2-pointers"),
         rate = if_else(str_detect(type, "R"), "success rate", "shot attempts")) 

bar <- regular_detail %>% 
  select(Season, starts_with("LFG")) %>% 
  mutate(LFGM = LFGM - LFGM3,
         LFGA = LFGA - LFGA3) %>% 
  mutate(LFGR = LFGM/LFGA,
         LFGR3 = LFGM3/LFGA3) %>% 
  select(-matches("GM")) %>% 
  pivot_longer(starts_with("LFG"), names_to = "type", values_to = "points") %>% 
  mutate(shot = if_else(str_detect(type, "3"), "3-pointers", "2-pointers"),
         rate = if_else(str_detect(type, "R"), "success rate", "shot attempts"))

foobar <- foo %>% 
  bind_rows(bar) %>% 
  mutate(team = fct_rev(as.factor(if_else(str_detect(type, "W"), "winner", "loser"))))

p1 <- foobar %>% 
  filter(rate == "success rate") %>% 
  ggplot(aes(shot, points, fill = team)) +
  geom_boxplot() +
  coord_flip() +
  scale_y_continuous(labels = scales::percent) +
  facet_wrap(~ shot, scales = "free_y", nrow = 2) +
  labs(x = "", y = "", fill = "", title = "Shot percentage") +
  theme_fivethirtyeight() +
  theme(legend.position = "none") +
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank())

p2 <- foobar %>% 
  filter(rate == "shot attempts") %>% 
  ggplot(aes(shot, points, fill = team)) +
  geom_boxplot() +
  coord_flip() +
  facet_wrap(~ shot, scales = "free", nrow = 2) +
  labs(x = "", y = "", fill = "", title = "Number of shot attempts") +
  theme_fivethirtyeight() +
  theme(legend.position = "none") +
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank())

p1 + p2 +
  plot_annotation(title = 'Regular season shot statistics: winning team (red) vs losing team (blue)')
```


We find:

- Winning teams (red) usually have 2-point shot percentages that are on average just above 50%; the median for the losing teams (blue) is a bit below that threshold. A similar observation can be made for the 3-pointers, where the winners score more reliably than the losers. And 2-point attempts overall are easier to score than 3-pointers. So far, not very suprising.

- However, if we look at the number of shot attempts we see that the winners (red) are, if anything, slightly behind the losing teams (blue). So it seems that what gives the winning teams the edge is the higher conversion rate of their shots, rather than creating a higher number of shot opportunities. As somebody who knows not a whole lot about basketball, this was a surprising find for me.


Let's look at the shot percentages since 2003:

```{r fig.cap ="Fig. 24", fig.height = 5.5}
foobar %>% 
  filter(rate == "success rate") %>% 
  group_by(Season, shot, team) %>% 
  summarise(mean_rate = mean(points),
            sd_rate = sd(points)) %>% 
  ungroup() %>% 
  mutate(Season = if_else(team == "winner", Season + 0.2, Season)) %>% 
  ggplot(aes(Season, mean_rate, col = team)) +
  geom_point() +
  scale_y_continuous(labels = scales::percent) +
  geom_errorbar(aes(ymin = mean_rate - sd_rate, ymax = mean_rate + sd_rate)) +
  facet_wrap(~ shot, nrow = 2, scales = "free") +
  theme_fivethirtyeight() +
  labs(x = "", y = "", title = "Shot percentages over time",
       subtitle = "Mean percentages with standard deviation error bars")
```

We find:

- The 2-pointers might have become slightly more efficient in the last years for the winning teams. The 3-point percentages show no movement at all.

- This stability over time might be an advantage when using shot percentage features in our prediction.


We will collect the remaining stats in 1 comprehensive overview layout. This plot will get a little bit busy, so take your time digesting it. We're also using a short helper function to build overlapping density plots more efficiently:

```{r fig.cap ="Fig. 25", fig.height = 5}
p1 <- regular_detail %>% 
  select(Season, matches("FT")) %>% 
  mutate(WFTR = WFTM/WFTA,
         LFTR = LFTM/LFTA) %>% 
  select(Season, ends_with("TR")) %>% 
  pivot_longer(ends_with("TR"), names_to = "team", values_to = "ft") %>% 
  mutate(team = fct_rev(as.factor(if_else(str_detect(team, "W"), "winner", "loser")))) %>% 
  filter(ft > 0) %>% 
  ggplot(aes(ft, fill = team)) +
  geom_density(bw = .05, alpha = 0.5) +
  scale_x_continuous(labels = scales::percent) +
  theme_tufte() +
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank()) +
  labs(x = "", y = "", title = "Free Throw percentage")

p2 <- regular_detail %>% 
  select(Season, matches("OR"), matches("DR"), -matches("Score")) %>% 
  rename(W_OR = WOR, L_OR = LOR, W_DR = WDR, L_DR = LDR) %>% 
  pivot_longer(matches("^[WL]."), names_to = c("win_lose", ".value"), names_pattern = "(.)_(.+)") %>% 
  mutate(team = fct_rev(as.factor(if_else(str_detect(win_lose, "W"), "winner", "loser")))) %>% 
  pivot_longer(cols = c(OR, DR), names_to = "type", values_to = "rebounds") %>% 
  mutate(type = if_else(type == "OR", "Offensive Rebounds", "Defensive Rebounds")) %>% 
  ggplot(aes(rebounds, fill = team)) +
  geom_density(bw = 1, alpha = 0.5) +
  facet_wrap(~type) +
  theme_tufte() +
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank()) +
  labs(x = "", y = "", title = "# Rebounds")

win_vs_lose <- function(df, feat, title_name, bw){
  df %>% 
    select(Season, contains(feat)) %>% 
    pivot_longer(contains(feat), names_to = "team", values_to = "feat") %>% 
    mutate(team = fct_rev(as.factor(if_else(str_detect(team, "W"), "winner", "loser")))) %>% 
    ggplot(aes(feat, fill = team)) +
    geom_density(bw = bw, alpha = 0.5) +
    theme_tufte() +
    theme(axis.text.y = element_blank(), axis.ticks.y = element_blank()) +
    labs(x = "", y = "", title = title_name)
}

p3 <- win_vs_lose(regular_detail, "Ast", "# Assists", 1)

p4 <- win_vs_lose(regular_detail, "TO", "# Turnovers", 1)

p5 <- win_vs_lose(regular_detail, "Stl", "# Steals", 1)

p6 <- win_vs_lose(regular_detail, "Blk", "# Blocks", 1)

p7 <- win_vs_lose(regular_detail, "PF", "# Fouls", 1)

layout <- "
ABB
CDE
FGH
"

p1 + p2 + p3 + p4 + p5 + p6 + p7 + guide_area() +
  plot_layout(design = layout, guides = 'collect') +
  plot_annotation(title = 'Regular Season stats - collected overview')
  
```

We find:

- Winning teams have visibly more Defensive Rebounds (but not Offensive ones!) and Assists; as well as slightly more Steals, Blocks, and a higher Free Throw percentage. They also commit fewer fouls and sligthly fewer Turnovers (aka instances where the winning team lost the ball to the other team; i.e. having fewer of those is better). 

- This difference between Defensive and Offensive Rebounds is really interesting! The offensive ones don't seem to matter at all, whereas a good defence after a missed shot has a large impact.


Those rebounds and turnovers are interesting; let's look at them in more detail. Here we'll be using another set of ridgeline plots from the [ggridges](https://cran.r-project.org/web/packages/ggridges/index.html) package to visualise their distributions from season to season:

```{r fig.cap ="Fig. 26", fig.height = 5}
p1 <- regular_detail %>% 
  select(Season, matches("OR"), matches("DR"), -matches("Score")) %>% 
  rename(W_OR = WOR, L_OR = LOR, W_DR = WDR, L_DR = LDR) %>% 
  pivot_longer(matches("^[WL]."), names_to = c("win_lose", ".value"), names_pattern = "(.)_(.+)") %>% 
  mutate(team = fct_rev(as.factor(if_else(str_detect(win_lose, "W"), "winner", "loser")))) %>% 
  pivot_longer(cols = c(OR, DR), names_to = "type", values_to = "rebounds") %>% 
  mutate(type = if_else(type == "OR", "Offensive Rebounds", "Defensive Rebounds")) %>% 
  filter(type == "Defensive Rebounds") %>% 
  ggplot(aes(x = rebounds, y = as.factor(Season), fill = team)) +
  geom_density_ridges(bandwidth = 1, alpha = 0.5) +
  geom_vline(xintercept = 25, linetype = 2) +
  theme_hc() +
  theme(legend.position = "bottom") +
  labs(x = "", y = "", title = "# Defensive Rebounds")

p2 <- regular_detail %>% 
  select(Season, matches("TO")) %>% 
  pivot_longer(matches("TO"), names_to = "team", values_to = "feat") %>%
  mutate(team = fct_rev(as.factor(if_else(str_detect(team, "W"), "winner", "loser")))) %>% 
  ggplot(aes(x = feat, y = as.factor(Season), fill = team)) +
  geom_density_ridges(bandwidth = 1, alpha = 0.5) +
  geom_vline(xintercept = 13, linetype = 2) +
  coord_cartesian(xlim = c(0, 35)) +
  theme_hc() +
  theme(legend.position = "bottom") +
  labs(x = "", y = "", title = "# Turnovers")

p1 + p2 +
  plot_annotation(title = 'Regular Seasons 2003 vs 2019: More Defensive Rebounds, fewer turnovers')
```

We find:

- Over the last few years, the average number of defensive rebounds has increased slightly. This appears to be true for the winning and losing teams, meaning that the gap between the two teams appears to remain approximately constant.

- Similarly, the number of Turnovers has somewhat decreased over the years, so that the 2003 average for the winning team is practically the same as the 2019 average for the losing team. The winning team's very slight advantage persists over the years.


Those are the detailed regular season stats, and we have exactly the same features for the brackets. Instead of repeating the above analysis step by step, I will try to fit everything into a single comprehensive overview plot; expanding the idea of Fig. 25. Feel free to repeat some of my other analysis steps in a forked version of this notebook. This code block will reuse the plotting function for overlapping density plots from above.


```{r fig.cap ="Fig. 27", fig.height = 6}
foo <- bracket_detail %>% 
  select(Season, starts_with("WFG")) %>% 
  mutate(WFGM = WFGM - WFGM3,
         WFGA = WFGA - WFGA3) %>% 
  mutate(WFGR = WFGM/WFGA,
         WFGR3 = WFGM3/WFGA3) %>% 
  select(-matches("GM")) %>% 
  pivot_longer(starts_with("WFG"), names_to = "type", values_to = "points") %>% 
  mutate(shot = if_else(str_detect(type, "3"), "3-pointers", "2-pointers"),
         rate = if_else(str_detect(type, "R"), "success rate", "shot attempts")) 

bar <- bracket_detail %>% 
  select(Season, starts_with("LFG")) %>% 
  mutate(LFGM = LFGM - LFGM3,
         LFGA = LFGA - LFGA3) %>% 
  mutate(LFGR = LFGM/LFGA,
         LFGR3 = LFGM3/LFGA3) %>% 
  select(-matches("GM")) %>% 
  pivot_longer(starts_with("LFG"), names_to = "type", values_to = "points") %>% 
  mutate(shot = if_else(str_detect(type, "3"), "3-pointers", "2-pointers"),
         rate = if_else(str_detect(type, "R"), "success rate", "shot attempts"))

foobar <- foo %>% 
  bind_rows(bar) %>% 
  mutate(team = fct_rev(as.factor(if_else(str_detect(type, "W"), "winner", "loser"))))

p1 <- foobar %>% 
  filter(rate == "success rate") %>% 
  ggplot(aes(points, fill = team)) +
  geom_density(bw = 0.05, alpha = 0.5) +
  scale_x_continuous(labels = scales::percent) +
  facet_wrap(~ shot, scales = "free", nrow = 1) +
  labs(x = "", y = "", fill = "", title = "Shot percentage") +
  theme_tufte() +
  theme(legend.position = "none") +
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank())

p2 <- foobar %>% 
  filter(rate == "shot attempts") %>% 
  ggplot(aes(points, fill = team)) +
  geom_density(bw = 2, alpha = 0.5) +
  facet_wrap(~ shot, scales = "free", nrow = 1) +
  labs(x = "", y = "", fill = "", title = "Shots attempted") +
  theme_tufte() +
  theme(legend.position = "none") +
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank())


p3 <- bracket_detail %>% 
  select(Season, matches("FT")) %>% 
  mutate(WFTR = WFTM/WFTA,
         LFTR = LFTM/LFTA) %>% 
  select(Season, ends_with("TR")) %>% 
  pivot_longer(ends_with("TR"), names_to = "team", values_to = "ft") %>% 
  mutate(team = fct_rev(as.factor(if_else(str_detect(team, "W"), "winner", "loser")))) %>% 
  filter(ft > 0) %>% 
  ggplot(aes(ft, fill = team)) +
  geom_density(bw = .05, alpha = 0.5) +
  scale_x_continuous(labels = scales::percent) +
  theme_tufte() +
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank()) +
  labs(x = "", y = "", title = "Free Throw percentage")

p4 <- bracket_detail %>% 
  select(Season, matches("OR"), matches("DR"), -matches("Score")) %>% 
  rename(W_OR = WOR, L_OR = LOR, W_DR = WDR, L_DR = LDR) %>% 
  pivot_longer(matches("^[WL]."), names_to = c("win_lose", ".value"), names_pattern = "(.)_(.+)") %>% 
  mutate(team = fct_rev(as.factor(if_else(str_detect(win_lose, "W"), "winner", "loser")))) %>% 
  pivot_longer(cols = c(OR, DR), names_to = "type", values_to = "rebounds") %>% 
  mutate(type = if_else(type == "OR", "Offensive Rebounds", "Defensive Rebounds")) %>% 
  ggplot(aes(rebounds, fill = team)) +
  geom_density(bw = 2, alpha = 0.5) +
  facet_wrap(~type) +
  theme_tufte() +
  theme(axis.text.y = element_blank(), axis.ticks.y = element_blank()) +
  labs(x = "", y = "", title = "# Rebounds")

p5 <- win_vs_lose(bracket_detail, "Ast", "# Assists", 2)

p6 <- win_vs_lose(bracket_detail, "TO", "# Turnovers", 2)

p7 <- win_vs_lose(bracket_detail, "Stl", "# Steals", 2)

p8 <- win_vs_lose(bracket_detail, "Blk", "# Blocks", 2)

p9 <- win_vs_lose(bracket_detail, "PF", "# Fouls", 2)

layout <- "
AAABBB
CCDDDD
EEFFGG
HHIIJJ
"

p1 + p2 + p3 + p4 + p5 + p6 + p7 + p8 + p9 + guide_area() +
  plot_layout(design = layout, guides = 'collect') +
  plot_annotation(title = 'Bracket stats - all detailed features')
  
```

We find:

- Many aspects look similar to the regular season stats: winners have higher shot percentages, more defensive rebounds and assists, slightly more blocks, and fewer fouls.

- The differences wrt the regular season are slight: the distributions for steals appears overlap a bit more, the 2-pointers for win vs lose are virtually indistinguishable as are the turnovers.

- Most of the interesting features in the regular season stats should be just as usable for the bracket stats.


To close this subsection, here's a brief look at some multi-feature connections. We take the difference in defensive rebounds between winning team - losing team, as well as the difference in personal fouls. Then we group both features, independently, into bins of 5 (e.g. 5 to 10, 10 to 15, etc). And the plot shows the aggregate counts for each combination of two feature levels. Here, the colours and the labels both indicate counts per bin:

```{r fig.cap ="Fig. 28", fig.height = 4}
binsize <- 5

foo <- bracket_detail %>% 
  select(Season, ends_with("TeamID"), contains("DR"), contains("PF")) %>% 
  mutate(dr_diff = WDR - LDR,
         pf_diff = WPF - LPF) %>% 
  select(Season, ends_with("TeamID"), ends_with("diff")) %>% 
  mutate(dr_diff_bin = dr_diff %/% binsize,
         dr_diff_group = as.factor(str_c(sprintf("%i", binsize * dr_diff_bin), ' to ', binsize * (dr_diff_bin + 1)))) %>% 
  mutate(pf_diff_bin = pf_diff %/% binsize,
         pf_diff_group = as.factor(str_c(sprintf("%i", binsize * pf_diff_bin), ' to ', binsize * (pf_diff_bin + 1)))) %>% 
  mutate(dr_diff_group = fct_relevel(dr_diff_group,
                                     c("-15 to -10", "-10 to -5", "-5 to 0", "0 to 5", "5 to 10",
                                       "10 to 15", "15 to 20", "20 to 25", "25 to 30"))) %>% 
  mutate(pf_diff_group = fct_relevel(pf_diff_group,
                                     c("-20 to -15", "-15 to -10", "-10 to -5", "-5 to 0", "0 to 5", "5 to 10")))

foo %>% 
  count(dr_diff_group, pf_diff_group) %>% 
  ggplot(aes(dr_diff_group, pf_diff_group, fill = n)) +
  geom_tile() +
  scale_fill_gradient(low = "#ffdf3f", high = "#5c46ff", trans = "log10") +
  geom_text(aes(label = n), color = "black", size = 3) +
  theme_fivethirtyeight() +
  theme(axis.title = element_text(), axis.title.y = element_text(angle = 0)) +
  theme(legend.position = "none") +
  labs(x = "Defensive Rebounds", y = "Fouls", title = "Bracket: Win - Lose for Rebounds and Fouls",
       subtitle = "Differences (win - lose) in DR and Fouls are grouped each in bins of 5")
```

We find:

- There is not that much of a correlation between the two features. The binned counts decay pretty smoothly from the center around -5 to 0 fouls and 0 - 10 defensive rebounds.

- Among the edge cases, the upper left corner is probably most noteworthy when compared to the upper right and lower left corners.


## Geography

We have geography information on all games from 2010 onward: namely US city and US state. I'm not going to spend a great deal of time on the geography features, but I want at least to provide some overview maps.

The `cities` data frame contains the US `State` names for some aggregate counts. We are also joining a Kaggle Dataset of [world cities](https://www.kaggle.com/max-mind/world-cities-database) to quickly get coordinates for 2/3 of our cities. To get all of them you would have to build a more specialised table that includes all the missing city names and their coordinates.

```{r}
city_loc <- world_cities %>% 
  filter(Country == "us") %>% 
  rename(State = Region) %>% 
  select(City, State, Latitude, Longitude) %>% 
  mutate(City = str_to_sentence(City)) %>% 
  inner_join(cities, by = c("City", "State")) %>% 
  left_join(cities_games, by = "CityID") %>% 
  count(City, Latitude, Longitude) %>% 
  select(longitude = Longitude, latitude = Latitude, City, n) %>% 
  usmap_transform() %>% 
  filter(longitude < -68)
  
state_loc <- cities %>% 
  left_join(cities_games, by = "CityID") %>% 
  count(State) %>% 
  rename(state = State)
```


This is a map with the aggregate number of games per US state visualised as the colour of the state. Those are all games for which we have cities (i.e. since 2010). The lighter the colour, the more games there were. On top of that, I'm adding total counts of games per city as blue circles, the size of which corresponds to the number of games. Those are only 2/3 of the cities in our dataset. The size and colour scales are logarithmic:

```{r fig.cap ="Fig. 29", fig.height = 6, message=FALSE, warning=FALSE}
plot_usmap(data = state_loc, values = "n",  color = "grey20",
           labels=FALSE, label_color = "grey20", ) +
  geom_point(data = city_loc, aes(x = longitude.1, y = latitude.1, size = n), color = "darkblue", alpha = 0.6) +
  scale_fill_viridis(option = "viridis", begin = 0.4, end = 1, trans = "log10") + #, breaks = c(100, 1000, 3000)) +
  #scale_size(trans = "log10") +
  scale_size(breaks = c(10, 100, 200, 500)) +
  theme(legend.position = "bottom") + #, title = element_text(size = 20),legend.text = element_text(size = 11)) 
  labs(size = "# Games: per city", fill = "per state", title = "Number of Games per State and for Selected Cities")
```

We find:

- Lots of games in California, Texas, and New York state. The northern states have the fewest games.

- The city map is essentially a population map, even without the missing 1/3 of cities. As with a population map, it becomes obvious that the distances between the mid-western cities are much larger than in the east.

One possible feature here, which has been suggested before, is to take into account how far (i.e. long) teams would travel to their next bracket game. This would require getting all the coordinates, but otherwise is pretty straightforward. It's worth testing, I guess, but I have my doubts whether travel times really matter anymore on this quasi-professional level.


## Team Rankings

Derived from Kenneth Massey's [College Basketball Ranking Composite page](https://www.masseyratings.com/cb/compare.htm), this dataset features team rankings, starting in 2003, from a large number of independent sources. There might be quite a lot of potential for additional predictive power here.

Let's get an overview first. There are in total 174 different ranking systems in this dataset. We would want to find out which ones are most useful for our analysis. In order to be useful, a rating system should:

1. Cover many Seasons to provide sufficient training data.

1. Provide pre-bracket ratings (day 133) that we can expect to have available for 2020 before the competition deadline.

1. Include the relevant teams (i.e. the ones that made it to the bracket), to provide relevant training data. We don't care so much about rating inaccuricies for teams that never made it to the bracket in the first place.

1. Be different from other ratings. As it is generally the case for any predictor features, if two features are very closely related then one of them can generally be discarded because it doesn't provide new information.


With that in mind, here are the 11 ranking systems that go back to 2003 or 2004 (BIH, CNG, WOB). This is a visualisation of a correlation matrix for their pre-bracket ranking (day 133) in the 2019 Season. The other seasons look similar. We're using the fantastic `ggcorr` tool from the [GGally package](https://cran.r-project.org/web/packages/GGally/index.html):

```{r fig.cap ="Fig. 30", fig.height = 4, message=FALSE}
foo <- bracket_res %>% 
  select(Season, ends_with("TeamID")) %>% 
  pivot_longer(ends_with("TeamID"), names_to = "foo", values_to = "TeamID")

bar <- ranks %>% 
  semi_join(foo, by = c("Season", "TeamID")) %>% 
  filter(RankingDayNum == 133) %>% 
  left_join(teams %>% select(TeamID, TeamName), by = "TeamID")
  
top_sys <- bar %>% 
  count(SystemName) %>% 
  arrange(desc(n)) %>% 
  filter(n >= 1067) %>% 
#  filter(n == max(n)) %>% 
  pull(SystemName)

bar %>% 
  filter(SystemName %in% top_sys) %>% 
  filter(Season == 2019) %>% 
  select(-Season, -RankingDayNum) %>% 
  pivot_wider(names_from = SystemName, values_from = OrdinalRank) %>% 
  #filter(POM <= 68) %>% 
  select(-TeamID, -TeamName) %>% 
  ggcorr(method = c("pairwise","spearman"), label = TRUE, angle = -0, hjust = 0.2) +
  coord_flip() +
  theme_fivethirtyeight() +
  theme(legend.position = "right") +
  ggtitle("Strong correlations between historical Ranking Systems",
          subtitle = "BIH, CNG, WOB ranking since 2004; rest since 2003. Here: 2019 Season.")
```

We find:

- There are very strong similarities between the systems. Some are virtually identical from the point of view of a correlation matrix. From this big picture perspective, using more than one system would provide only marginal improvements.

- Note, that this plot tells us something about the strength of a (linear) relationship. It doesn't really touch much on outliers and the reasons that could be causing them.


To have a more detailed view on how individual rating systems compare, here we plot the POM ratings vs COL for the 2019 season bracket stage. All these teams made it to the list of 68 on Selection Sunday. We use colour coding to show the difference in ranking between the two systems:


```{r fig.cap ="Fig. 31", fig.height = 5}
bar %>% 
  filter(Season == 2019 & SystemName %in% c("POM", "COL")) %>% 
  pivot_wider(names_from = SystemName, values_from = OrdinalRank) %>%
  mutate(rank_diff = POM - COL) %>% 
  mutate(lab = if_else(abs(rank_diff) > 20, TeamName, "")) %>% 
  ggplot(aes(POM, COL, col = rank_diff)) +
  geom_abline(slope = 1, linetype = 2, col = "grey30") +
  geom_point() +
  geom_label_repel(aes(label = lab), alpha = 0.8) +
  scale_color_viridis(option = "viridis", begin = 0, end = 1) +
  theme_fivethirtyeight() +
  labs(col = "Rank Difference (POM - COL)", title = "Ratings comparison: POM vs COL for 2019 Season")
```

We find:

- As expected from the previous plot, the correlation is strong. However, there are also notable outliers and more cases where the POM ranks were larger numbers than the COL ones (i.e. the yellow and light green data points vs the dark blue ones).

- We've marked the largest outliers with 20 ranking points difference in either direction. There is only 1 case (St Mary's CA) where POM ranked a team higher than COL by more than 20 ranks (i.e. the numerical value of the POM rank is lower than the COL rank).

- This gives us another way to distinguish two rating systems: the direction and magnitude of their relative ranking differences.


---

I stopped working on this Notebook after the 2020 bracket was cancelled due to Covid-19. However, I believe that there is quite a bit of interesting visuals and methodology here, also for the 2021 competition and beyond.

Best of success in future March Madness challenges!