---
title: 'How to enter the "Madness"? Analysis on the picking and seeding of the tournament teams.'
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')
```

```{r, message = FALSE}
library(tidyverse)
library(ggplot2)
library(cowplot)
library(gridExtra)
library(knitr)
library(magrittr)
library(GGally)
library(reshape2)
library(reshape)
library(MASS)
library(rattle)
library(heplots)
library(InformationValue)
library(klaR)
```

# The "Big" Picture

The NCAA March Madness tournament structure comes from a selection and seeding process, which is conducted by a selection committee. This process is composed of three major phases, (1) selecting of the 36 "at-large" teams, (2) seeding the 68 teams and (3) generating a tournament bracket. The committee adheres to the principles of selecting the *36 best teams* (other than the 32 conference champions) and achieving *reasonable competitive balance* in each region of the bracket. Even though the whole process has clear guidelines, the decision-making process and the criteria on which each committee member relies on is unknown and varies from year to year.

I seek to infer the committee's criteria for picking the at-large teams entering in the tournament by analyzing the regular-season data from 2003 and later. Along the way, I will also explore the relations between the seedings with other iconic ranking methods and explore their ability on forcasting the match results.

(More detailed and "official" guideline about the selection process: https://www.ncaa.com/news/basketball-men/article/2018-10-19/how-field-68-teams-picked-march-madness)

# Seeding and Its Effect on Quantifying the Teams

## Introduction

Before investigating the committee's hidden formula for selecting and seeding the teams, it is interesting to see whether they are making sensible picks. Specifically, are these top teams truly strong?

There are two intuitive perspective evaluating these seedings: *hindsight* and *foresight*. From the perspective of hindsight, we want to explore how well does the seeds reflect these teams' regular-season performance. However, the perspective of foresight seeks to examine evaluation of the seeds by examining their predictability for the tournament.

I will use both box scores and sabermetrics to quantify each team's performance. Kaggle provided detailed per-game box score data since 2002-03 season. For the sabermetrics, I used the definations and calculation from Jason Zivkovic's 2019 Kaggle notebook [1] and Kubatko's paper. [2]


## Hindsight Analysis


```{r}
teams = read.csv("../input/march-madness-analytics-2020/2020DataFiles/2020DataFiles/2020-Mens-Data/MDataFiles_Stage1/MTeams.csv") %>% filter(FirstD1Season<=2019 & LastD1Season>=2003)
#reg_season_compact = read.csv("../input/march-madness-analytics-2020/2020DataFiles/2020DataFiles/2020-Mens-Data/MDataFiles_Stage1/MRegularSeasonCompactResults.csv") %>% filter(Season >= 2003)
reg_season_detail = reg_season_compact = read.csv("../input/march-madness-analytics-2020/2020DataFiles/2020DataFiles/2020-Mens-Data/MDataFiles_Stage1/MRegularSeasonDetailedResults.csv") %>% filter(Season >= 2003)

#tourney_compact = read.csv("../input/march-madness-analytics-2020/2020DataFiles/2020DataFiles/2020-Mens-Data/MDataFiles_Stage1/MNCAATourneyCompactResults.csv") %>% filter(Season >= 2003)

seeds = read.csv("../input/march-madness-analytics-2020/2020DataFiles/2020DataFiles/2020-Mens-Data/MDataFiles_Stage1/MNCAATourneySeeds.csv", stringsAsFactors = F)  %>% filter(Season >= 2003)
seeds$Seed = seeds$Seed %>% str_replace_all("[:alpha:]","")
seeds$Seed = as.integer(seeds$Seed)

tourney_team = seeds %>% dplyr::select(Season, TeamID, Seed)
```

### Seeds vs. Box Scores

Each year, teams in the March Madness tournament come from different conferences. Therefore, it is not garaunteed that these teams had played against each other in the regular season. Also, to consider as many match histories as possible, it is desirable to look, it would be a good idea to look at each teams' average performance during the season.

The questions that I intend to explore are whether higher seeds perform better, using box scores as a reference and determining which game statistics have more effect on the final seedings.

The plots below all have seeding positions on the horizontal axis and the vertical axis shows the teams' seasonal means in each feature. The red line is a fitted linear regression between the independent and dependent variable.

Most of the statistics have a weak linear correlation with the seedings. Since lower number in seeding represents a team being strong, most of the correlation coefficients are negative. 

Among all the box scores, the mean score, field goals made, assists and blocks have higher coorelation with the seedings and their slopes for the linear regression line are also steeper. Therefore, the committee may regard teams achieve higher scores and cooperative moves such as, assist and blocks to be stronger teams. 


```{r}
teamID = teams$TeamID
win_feature = grep("W", colnames(reg_season_detail), value=T)
lose_feature = grep("L", colnames(reg_season_detail), value=T)
feature = str_replace_all(win_feature, "W", "")
team_wins = reg_season_detail %>%  dplyr::select(Season, all_of(win_feature))
team_loses = reg_season_detail %>% dplyr::select(Season, all_of(lose_feature))
colnames(team_wins) = c("Season",feature)
colnames(team_loses) = c("Season",feature)
all_team_reg_season_avg_stats = dplyr::bind_rows(team_wins,team_loses) %>% dplyr::select(-Loc) %>% group_by(Season, TeamID) %>% summarise_all(mean)
tourney_team_reg_season_avg_performance = all_team_reg_season_avg_stats%>% right_join(tourney_team,by=c('TeamID','Season'))
```

```{r, fig.width = 10, fig.height = 10,  message=FALSE}
average_per_seed_per_year = tourney_team_reg_season_avg_performance %>% group_by(Seed, Season) %>% summarise_all(mean)
average_per_seed = tourney_team_reg_season_avg_performance %>% group_by(Seed) %>% summarise_all(mean)

box_score_plot <- function(seed,variable, plot_name){

        tourney_team_reg_season_avg_performance %>% 
        ggplot(aes_string(seed, y=variable, color=variable))+
        geom_point()+
        geom_smooth(method = "lm", colour="red",se=F)+
        ggtitle(plot_name, subtitle=paste0("Correlation", round(cor(tourney_team_reg_season_avg_performance[,seed],tourney_team_reg_season_avg_performance[,variable]),2)))+
        theme_minimal()+
        scale_color_gradient(low="Cyan1", high="darkblue")     
    
}

b0 <- box_score_plot('Seed','Score',"average score")

b1 <- box_score_plot('Seed','FGM',"field goals made")

b2 <- box_score_plot('Seed',"FGA","field goals attempted")

b3 <- box_score_plot('Seed',"FGM3","three pointers made")

b4 <- box_score_plot('Seed',"FGA3","three pointers attempted")

b5 <- box_score_plot('Seed',"FTM","free throws made")

b6 <- box_score_plot('Seed',"FTA","free throws attempted")

b7 <- box_score_plot('Seed',"OR","offensive rebounds")

b8 <- box_score_plot('Seed',"DR","defensive rebounds")

b9 <- box_score_plot('Seed',"Ast","assists")

b10 <- box_score_plot('Seed',"TO","turnovers committed")

b11 <- box_score_plot('Seed',"Stl","steals")

b12 <- box_score_plot('Seed',"Blk","blocks")

b13 <- box_score_plot('Seed',"PF","personal fouls")



grid.arrange(b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,ncol=4)

```

Each year, each seed position is shared by 4–6 teams. Therefore, to better examine the linear relationship between seedings and box scores, I calculated the average box scores per seed position. The heat map below shows the correlation between these averaged box scores.

The heat map shows that three-point goals and attempts, free throws, and steals have low correlation with the seedings. Team's seedings, average scorings and successful field goals have high correlation with the two offensive and defensive rebounds (OR, DR), assist (Ast) and blocks (Blk). Based on these high correlations, we interpret that aggressive teams are more like to attain higher score and have higher seed position.

```{r, fig.width = 8, fig.height = 8,  message=FALSE}
average_per_seed_per_year = tourney_team_reg_season_avg_performance %>% group_by(Seed, Season) %>% summarise_all(mean)
average_per_seed = tourney_team_reg_season_avg_performance %>% group_by(Seed) %>% summarise_all(mean)
ggcorr(average_per_seed[,-c(2,3)], palette = "RdBu", label = TRUE)
```

### Seeds vs. Sabermatrics

Possession is calculated by the weighted sum of field goal attempts, turnovers, and free throw attempts minus the offensive rebounds. The offensive and defensive ratings and strength of schedule is calculated based on the team's and opponent score divided by their possession. The graph shows that possession has almost no impact on the seeding positions. 

However, the linear correlations between seedings and the two ratings and the strength of schedule are quite strong. Notice, a team's average seasonal offensive rating equals the average score divided by average possession, which has a -0.48 correlation with the seeding comparing to the -0.32 correlation between average score alone with the seedings (shown in the graph above). Therefore, by dividing the possession, the average score of a team has more significant linear relationship with its seeding. After subtracting the offensive rating by the defensive rating (which equals the opponent's offensive rating), the strength of schedule has shown even stronger linear relationship with the seeding. With this observation, it is a reasonable guess that, the committee will look more into the relative performance of the teams compare to their component during the seeding process. 

The winning ratio of the teams shows a strong linear correlation with the seeding. To other to be selected as the top seed, the teams has to at least maintain an 80% winning ratio in the regular season. It is interesting to observe a surge of winning ration at seed 12. This sudden increase does not appear significantly on the other sabermetrics.

```{r}
reg_season_stats <- read.csv("../input/march-madness-analytics-2020/2020DataFiles/2020DataFiles/2020-Mens-Data/MDataFiles_Stage1/MRegularSeasonDetailedResults.csv", stringsAsFactors = FALSE)
# regular season
reg_season_stats <- reg_season_stats %>%
  mutate(WPoss = WFGA + (WFTA * 0.475) + WTO - WOR,
         LPoss = LFGA + (LFTA * 0.475) + LTO - LOR) %>% 
  filter(Season>=2003)

advanced_stats_per_game <- rbind(
  reg_season_stats %>%
    dplyr::select(Season, DayNum, TeamID=WTeamID, Score=WScore, OScore=LScore, Poss=WPoss, OPoss=LPoss, WLoc, NumOT, FGM=WFGM, FGA=WFGA, FGM3=WFGM3,  FGA3=WFGA3, FTM=WFTM, FTA=WFTA, OR=WOR, DR=WDR, Ast=WAst, TO=WTO, Stl=WStl, Blk=WBlk, PF=WPF, OFGM=LFGM, OFGA=LFGA, OFGM3=LFGM3, OFGA3=LFGA3, OFTM=LFTM, OFTA=LFTA, O_OR=LOR, ODR=LDR, OAst=LAst, OTO=LTO, OStl=LStl, OBlk=LBlk, OPF=LPF) %>%
    mutate(Winner=1),

  reg_season_stats %>%
    dplyr::select(Season, DayNum, TeamID=LTeamID, Score=LScore, OScore=WScore, Poss=LPoss, OPoss=WPoss, WLoc, NumOT, FGM=LFGM, FGA=LFGA, FGM3=LFGM3,  FGA3=LFGA3, FTM=LFTM, FTA=LFTA, OR=LOR, DR=LDR, Ast=LAst, TO=LTO, Stl=LStl, Blk=LBlk, PF=LPF, OFGM=WFGM, OFGA=WFGA, OFGM3=WFGM3, OFGA3=WFGA3, OFTM=WFTM, OFTA=WFTA, O_OR=WOR, ODR=WDR, OAst=WAst, OTO=WTO, OStl=WStl, OBlk=WBlk, OPF=WPF) %>%
    mutate(Winner=0)) %>%
  mutate(PtsDiff = Score - OScore)


advanced_stats_per_game <- advanced_stats_per_game %>%
  mutate(OffRtg = 100 * (Score / Poss),
         DefRtg = 100 * (OScore / OPoss),
         SoS = OffRtg - DefRtg,
         Pie = Score + FGM + FTM - FGA - FTA + DR + (0.5 * OR) + Ast + Stl + (0.5 * Blk) - PF - TO,
         OPie = OScore + OFGM + OFTM - OFGA - OFTA + ODR + (0.5 * O_OR) + OAst + OStl + (0.5 * OBlk) - OPF - OTO,
         Tie = Pie / (Pie + OPie) * 100,
         AstRatio = 100 * Ast / (FGA + (0.475 * FTA) + Ast + TO),
         TORatio = 100 * TO / (FGA + (0.475 * FTA) + Ast + TO),
         TSPerc = 100 * Score / (2 * (FGA + (0.475 * FTA))),
         EFGPerc = 100* (FGM + 0.5 * FGM3) / FGA,
         FTR = FTA / FGA,
         ORbP = OR / (OR + ODR),
         DRbP = DR / (DR + O_OR),
         TRbP = (DR + OR) / (DR + OR + ODR + O_OR),
         WinRate = Winner) %>%
  dplyr::select(Season, TeamID, WinRate, PtsDiff, Poss, OffRtg, DefRtg, SoS, Pie, Tie, AstRatio, TORatio, TSPerc, EFGPerc, FTR, ORbP, DRbP, TRbP)

all_team_reg_season_adv_stats = advanced_stats_per_game %>% group_by(Season,TeamID) %>% summarise_all(mean)
#all_team_reg_season_adv_stats_rank = all_team_reg_season_adv_stats
team_count = nrow(all_team_reg_season_adv_stats)
#for (i in seq(2, ncol(all_team_reg_season_adv_stats))){
#    temp = cbind(all_team_reg_season_adv_stats[order(data.frame(all_team_reg_season_adv_stats[,i]),decreasing = T),], Rank=seq(team_count))
#    all_team_reg_season_adv_stats_rank[,i]=temp[order(temp$TeamID),ncol(temp)] 
#}
tourney_team_reg_season_adv = all_team_reg_season_adv_stats %>% right_join(tourney_team,by=c('Season','TeamID'))
#ggcorr(tourney_team_reg_season_adv[,-c(1,2)], palette = "RdBu", label = TRUE)
```

```{r, fig.width = 10, fig.height = 10,  message=FALSE}
adv_score_plot <- function(seed,variable, plot_name){
    tourney_team_reg_season_adv %>% 
        ggplot(aes_string(seed, y=variable, color=variable))+
        geom_point()+
        geom_smooth(method = "lm", colour="red",se=F)+
        ggtitle(plot_name,subtitle=paste0("Correlation", round(cor(tourney_team_reg_season_adv[,seed],tourney_team_reg_season_adv[,variable]),2)))+
        theme_minimal()+
        scale_color_gradient(low="cyan1", high="darkblue")     
}

b0 <- adv_score_plot('Seed','WinRate',"Winning Rate")

b1 <- adv_score_plot('Seed','PtsDiff',"Margin against Opponents")

b2 <- adv_score_plot('Seed',"Poss","Possession")

b3 <- adv_score_plot('Seed',"OffRtg","Offensive Rating")

b4 <- adv_score_plot('Seed',"DefRtg","Deffensive Rating")

b5 <- adv_score_plot('Seed',"SoS","Strength of Schedule")

b6 <- adv_score_plot('Seed',"Pie","Performance Impact Estimator")

b7 <- adv_score_plot('Seed',"Tie","Team Impact Estimator")

b8 <- adv_score_plot('Seed',"AstRatio","Assist Ratio")

b9 <- adv_score_plot('Seed',"TORatio","Turnover Ratio")

b10 <- adv_score_plot('Seed',"TSPerc","True Shooting %")

b11 <- adv_score_plot('Seed',"EFGPerc","Effective Field Goal %")

b12 <- adv_score_plot('Seed',"FTR","Free Throw Rate")

b13 <- adv_score_plot('Seed',"ORbP","Offensive Rebound %")

b14 <- adv_score_plot('Seed',"DRbP","Defensive Rebound %")

b15 <- adv_score_plot('Seed',"TRbP","Total Rebound %")



grid.arrange(b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14,b15,ncol=4)

```

The heat map shows the correlation between the seasonal average of the sabermetrics and the seeding positions.

In comparison to the heat map for basic box scores, the sabermetrics show much stronger correlation between each features. 

Among all of these sabermetrics, possession, free-throw ration and defensive rebound percentage seem to be largely independent of the other quantities and the seedings.

```{r, fig.width = 10, fig.height = 10,  message=FALSE}
adv_average_per_seed_per_year = tourney_team_reg_season_adv %>% group_by(Seed, Season) %>% summarise_all(mean)
adv_average_per_seed = tourney_team_reg_season_adv %>% group_by(Seed) %>% summarise_all(mean)

ggcorr(adv_average_per_seed[,-c(2,3)], low = "steelblue", mid = "white", high = "red", label = TRUE)
```

## Foresight Analysis

Using the seeding to evaluate the team performances in the tournament offers a hindsight perspective for evaluating the committee's selections. I will first explore the relation between the difference between box scores and sabermetrics and the difference between seedings for each tournament game. I will then build a naive linear model to predict the margin of victory of the tournament games using the difference of seed between each pair of teams.

```{r}
tourney_compact = read.csv("../input/march-madness-analytics-2020/2020DataFiles/2020DataFiles/2020-Mens-Data/MDataFiles_Stage1/MNCAATourneyCompactResults.csv") %>% filter(Season >= 2003)
tourney_detailed = read.csv("../input/march-madness-analytics-2020/2020DataFiles/2020DataFiles/2020-Mens-Data/MDataFiles_Stage1/MNCAATourneyDetailedResults.csv") %>% filter(Season >= 2003)
seeds = read.csv("../input/march-madness-analytics-2020/2020DataFiles/2020DataFiles/2020-Mens-Data/MDataFiles_Stage1/MNCAATourneySeeds.csv", stringsAsFactors = F)  %>% filter(Season >= 2003)
seeds$Seed = seeds$Seed %>% str_replace_all("[:alpha:]","")
seeds$Seed = as.integer(seeds$Seed)

tourney_team = seeds %>% dplyr::select(Season, TeamID, Seed)
tourney_detailed = tourney_detailed  %>% mutate(TeamID=WTeamID) %>% left_join(tourney_team,by=c("TeamID", "Season")) %>% 
    mutate(WSeed=Seed, TeamID=LTeamID) %>% dplyr::select(-c(Seed))  %>% left_join(tourney_team,by=c("TeamID", "Season")) %>% 
    mutate(LSeed=Seed) %>% dplyr::select(-c(Seed))  %>% filter(DayNum >= 136)

touney_diff_extend <- rbind(
    tourney_detailed  %>% dplyr::select(-c(WLoc, NumOT)) %>% mutate(DScore = WScore-LScore,
                                                                 DSeed = WSeed - LSeed,
                                                                 DFGM = WFGM - LFGM,
                                                                 DFGA = WFGA - LFGA,
                                                                 DFGM3 = WFGM3 - LFGM3,
                                                                 DFGA3 = WFGA3 - LFGA3,
                                                                 DFTM = WFTM - LFTM,
                                                                 DFTA = WFTA - LFTA,
                                                                 DOR = WOR - LOR,
                                                                 DDR = WDR - LDR,
                                                                 DAst = WAst - LAst,
                                                                 DTO = WTO - LTO,
                                                                 DStl = WStl - LStl,
                                                                 DBlk = WBlk - LBlk,
                                                                 DPF = WPF - LPF),
    tourney_detailed  %>% dplyr::select(-c(WLoc, NumOT)) %>% mutate(DScore = -WScore+LScore,
                                                                 DSeed = -WSeed + LSeed,
                                                                 DFGM = -WFGM + LFGM,
                                                                 DFGA = -WFGA + LFGA,
                                                                 DFGM3 = -WFGM3 + LFGM3,
                                                                 DFGA3 = -WFGA3 + LFGA3,
                                                                 DFTM = -WFTM + LFTM,
                                                                 DFTA = -WFTA + LFTA,
                                                                 DOR = -WOR + LOR,
                                                                 DDR = -WDR + LDR,
                                                                 DAst = -WAst + LAst,
                                                                 DTO = -WTO + LTO,
                                                                 DStl = -WStl + LStl,
                                                                 DBlk = -WBlk + LBlk,   
                                                                 DPF = WPF - LPF)
)
tourney_diff_only <- touney_diff_extend  %>% dplyr::select(Season, WTeamID, LTeamID,DSeed,DScore, DFGM, DFGA, DFGM3, DFGA3, DFTM, DFTA, DOR,DDR, DAst, DTO, DStl, DBlk, DPF)
```
### Tournament Box Score and Seed Difference

Let the feature "DSeed" denote the seed difference between the winner and loser of a game. Therefore, a negative DSeed means that the stronger team (team with higher seed position) won.

From my plot DSeed against DFGA3, observe that stronger teams are attempting both fewer two-point shots and fewer three-point shots than their weaker opponents. However, the stronger teams succeed more often in their goal attempts.

Higher-seeded teams are also have more defensive rebounds, assists and blocks, observed from their steeply-sloped negative correlations.


```{r,fig.width = 10, fig.height = 10,  message=FALSE}
box_score_plot <- function(seed,variable, plot_name){
    tourney_diff_only %>% 
        ggplot(aes_string(seed, y=variable, color=variable))+
        geom_point()+
        geom_smooth(method = "lm", colour="red", se=F)+
        ggtitle(plot_name,subtitle=paste0("Correlation", round(cor(tourney_diff_only[,seed],tourney_diff_only[,variable]),2)))+
        theme_minimal()+
        scale_color_gradient(low="cyan1", high="darkblue")     
}
b0 <- box_score_plot('DSeed','DScore',"Average Score")

b1 <- box_score_plot('DSeed','DFGM',"field goals made")

b2 <- box_score_plot('DSeed',"DFGA","field goals attempted")

b3 <- box_score_plot('DSeed',"DFGM3","three pointers made")

b4 <- box_score_plot('DSeed',"DFGA3","three pointers attempted")

b5 <- box_score_plot('DSeed',"DFTM","free throws made")

b6 <- box_score_plot('DSeed',"DFTA","free throws attempted")

b7 <- box_score_plot('DSeed',"DOR","offensive rebounds")

b8 <- box_score_plot('DSeed',"DDR","defensive rebounds")

b9 <- box_score_plot('DSeed',"DAst","assists")

b10 <- box_score_plot('DSeed',"DTO","turnovers committed")

b11 <- box_score_plot('DSeed',"DStl","steals")

b12 <- box_score_plot('DSeed',"DBlk","blocks")

b13 <- box_score_plot('DSeed',"DPF","personal fouls")

grid.arrange(b0,b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,ncol=4)
```

### Sabermetric and Seed Difference



```{r}

tourney_stats <- tourney_detailed %>%
  mutate(WPoss = WFGA + (WFTA * 0.475) + WTO - WOR,
         LPoss = LFGA + (LFTA * 0.475) + LTO - LOR) %>% 
  filter(Season>=2003)


tourney_advanced_stats_per_game <- rbind(
  tourney_stats %>%
    dplyr::select(Season, DayNum,Seed=WSeed,OSeed=LSeed,TeamID=WTeamID, Score=WScore, OScore=LScore, Poss=WPoss, OPoss=LPoss, WLoc, NumOT, FGM=WFGM, FGA=WFGA, FGM3=WFGM3,  FGA3=WFGA3, FTM=WFTM, FTA=WFTA, OR=WOR, DR=WDR, Ast=WAst, TO=WTO, Stl=WStl, Blk=WBlk, PF=WPF, OFGM=LFGM, OFGA=LFGA, OFGM3=LFGM3, OFGA3=LFGA3, OFTM=LFTM, OFTA=LFTA, O_OR=LOR, ODR=LDR, OAst=LAst, OTO=LTO, OStl=LStl, OBlk=LBlk, OPF=LPF) %>%
    mutate(PtsDiff = Score - OScore) %>% 
    mutate(SeedDiff = Seed-OSeed),

  tourney_stats %>%
    dplyr::select(Season, DayNum,Seed=LSeed,OSeed=WSeed,TeamID=LTeamID, Score=LScore, OScore=WScore, Poss=LPoss, OPoss=WPoss, WLoc, NumOT, FGM=LFGM, FGA=LFGA, FGM3=LFGM3,  FGA3=LFGA3, FTM=LFTM, FTA=LFTA, OR=LOR, DR=LDR, Ast=LAst, TO=LTO, Stl=LStl, Blk=LBlk, PF=LPF, OFGM=WFGM, OFGA=WFGA, OFGM3=WFGM3, OFGA3=WFGA3, OFTM=WFTM, OFTA=WFTA, O_OR=WOR, ODR=WDR, OAst=WAst, OTO=WTO, OStl=WStl, OBlk=WBlk, OPF=WPF) %>%
    mutate(PtsDiff = Score - OScore) %>% 
    mutate(SeedDiff = Seed-OSeed))


tourney_advanced_stats_per_game <- tourney_advanced_stats_per_game %>%
  mutate(OffRtg = 100 * (Score / Poss),
         DefRtg = 100 * (OScore / OPoss),
         SoS = OffRtg - DefRtg,
         Pie = Score + FGM + FTM - FGA - FTA + DR + (0.5 * OR) + Ast + Stl + (0.5 * Blk) - PF - TO,
         OPie = OScore + OFGM + OFTM - OFGA - OFTA + ODR + (0.5 * O_OR) + OAst + OStl + (0.5 * OBlk) - OPF - OTO,
         Tie = Pie / (Pie + OPie) * 100,
         AstRatio = 100 * Ast / (FGA + (0.475 * FTA) + Ast + TO),
         TORatio = 100 * TO / (FGA + (0.475 * FTA) + Ast + TO),
         TSPerc = 100 * Score / (2 * (FGA + (0.475 * FTA))),
         EFGPerc = 100* (FGM + 0.5 * FGM3) / FGA,
         FTR = FTA / FGA,
         ORbP = OR / (OR + ODR),
         DRbP = DR / (DR + O_OR),
         TRbP = (DR + OR) / (DR + OR + ODR + O_OR),) %>%
  dplyr::select(Season, TeamID, SeedDiff, PtsDiff, Poss, OffRtg, DefRtg, SoS, Pie, Tie, AstRatio, TORatio, TSPerc, EFGPerc, FTR, ORbP, DRbP, TRbP)

```
There is no significant linear correlation between the possession, turnover ratio, and free-throw rate against the seeding difference in each game. 

Among all of the sabermetrics, strength of schedule and the Team Impact Estimator have the strongest correlation against the seeding differences. These two quantities are all calculated by combinations of multiple box scores; the strong linear relation between the seed difference and these two metrics represent the seed indeed reflect a team's strength.

```{r, fig.width = 10, fig.height = 10,  message=FALSE}
adv_score_plot <- function(seed,variable, plot_name){
    tourney_advanced_stats_per_game %>% 
        ggplot(aes_string(seed, y=variable, color=variable))+
        geom_point()+
        geom_smooth(method = "lm", colour="red",se=F)+
        ggtitle(plot_name,subtitle=paste0("Correlation", round(cor(tourney_advanced_stats_per_game[,seed],tourney_advanced_stats_per_game[,variable]),2)))+
        theme_minimal()+
        scale_color_gradient(low="cyan1", high="darkblue")     
}



b2 <- adv_score_plot('SeedDiff',"Poss","Possession")

b3 <- adv_score_plot('SeedDiff',"OffRtg","Offensive Rating")

b4 <- adv_score_plot('SeedDiff',"DefRtg","Deffensive Rating")

b5 <- adv_score_plot('SeedDiff',"SoS","Strength of Schedule")

b6 <- adv_score_plot('SeedDiff',"Pie","Performance Impact Estimator")

b7 <- adv_score_plot('SeedDiff',"Tie","Team Impact Estimator")

b8 <- adv_score_plot('SeedDiff',"AstRatio","Assist Ratio")

b9 <- adv_score_plot('SeedDiff',"TORatio","Turnover Ratio")

b10 <- adv_score_plot('SeedDiff',"TSPerc","True Shooting %")

b11 <- adv_score_plot('SeedDiff',"EFGPerc","Effective Field Goal %")

b12 <- adv_score_plot('SeedDiff',"FTR","Free Throw Rate")

b13 <- adv_score_plot('SeedDiff',"ORbP","Offensive Rebound %")

b14 <- adv_score_plot('SeedDiff',"DRbP","Defensive Rebound %")

b15 <- adv_score_plot('SeedDiff',"TRbP","Total Rebound %")



grid.arrange(b1,b2,b3,b4,b5,b6,b7,b8,b9,b10,b11,b12,b13,b14,b15,ncol=4)
```

### Forcasting Margin of Victory Using Seeds

Because of the significant linear correlation between the seed and the margin of victory [3], it is possible to generate an accurate forecast of the margin of victory of the games solely using the seeds.

The horizontal axis is the seed difference using the higher seed minus the lower seed and the margin of victory is calculated using lower-seeded team's score to minus the score of the higher-seeded team. Therefore, points with a negative value represet the upsets (meaning weaker team beating the stronger one, in terms of the seeding)

This simple linear model is able to achieve an a value of the sum of residuals that is close to 0 and a value of the sum of squared residuals of about 130.


```{r, message=F}

tourney_compact = read.csv("../input/march-madness-analytics-2020/2020DataFiles/2020DataFiles/2020-Mens-Data/MDataFiles_Stage1/MNCAATourneyCompactResults.csv") %>% filter(Season >= 2003)

seeds_performance <- tourney_compact %>% dplyr::select(Season, DayNum, WTeamID, WScore, LTeamID, LScore) %>% 
  filter(Season >= 2003, DayNum >=136) %>%
  left_join(tourney_team, by = c("Season", "WTeamID" = "TeamID")) %>%
  left_join(tourney_team, by = c("Season", "LTeamID" = "TeamID")) %>%
  mutate(WinnerSeed = Seed.x, LoserSeed = Seed.y) %>%
  mutate(winner_higher_seed = ifelse(WinnerSeed < LoserSeed, "Higher Seed Wins", ifelse(LoserSeed < WinnerSeed, "Lower Seed Wins", "Same Seed"))) %>%
  mutate(winner_higher_seed = factor(winner_higher_seed, levels = c("Lower Seed Wins", "Same Seed", "Higher Seed Wins"))) %>% 
  mutate(higher_seed = ifelse(WinnerSeed < LoserSeed, LoserSeed, ifelse(LoserSeed < WinnerSeed, WinnerSeed, WinnerSeed))) %>%
  mutate(lower_seed = ifelse(WinnerSeed < LoserSeed, WinnerSeed, ifelse(LoserSeed < WinnerSeed, LoserSeed, LoserSeed))) %>%
  mutate(higher_seed_score = ifelse(WinnerSeed < LoserSeed, LScore, ifelse(LoserSeed < WinnerSeed, WScore, WScore))) %>%
  mutate(lower_seed_score = ifelse(WinnerSeed < LoserSeed, WScore, ifelse(LoserSeed < WinnerSeed, LScore, LScore))) %>%
  mutate(seed_diff = higher_seed-lower_seed)  %>% 
  mutate(margin_of_vic = lower_seed_score - higher_seed_score)


seeds_performance  %>% ggplot(aes(x= seed_diff, y= margin_of_vic, colour = margin_of_vic)) + 
                        geom_point() +
                        geom_smooth(method = "lm", se=F, colour= "red")+
                        geom_hline(yintercept = 0)
first_round_fit <- lm(seeds_performance$margin_of_vic ~ seeds_performance$seed_diff, data=seeds_performance)
summary(first_round_fit)

```

# Clustering on Teams Entering The Tournament

In this section, I am using linear dicriminant analysis and logistic regression to build a classification model to generate my own expert pick for teams entering into the March Madness. In this process, I will try to identify the statistical features that are most likely considered by the selection committee while choosing teams entering in the tournament.

## Linear Discriminant Analysis on Box Scores And Sabermetrics

### Selection by Box Scores

The naivest model is to use solely the average box score per season for each team to determine whether it can enter the tournament. 

The model shows that the average number of offensive rebounds, and steals perform best determining the tournament team, with a coefficient larger than 0.3 in magnitude.

```{r, message=FALSE}
avg_box_scores = all_team_reg_season_avg_stats
success_team = data.frame(Season = tourney_team$Season, TeamID = tourney_team$TeamID, Success = rep(1,nrow(tourney_team)))
avg_box_scores = avg_box_scores %>% left_join(success_team, by = c('Season','TeamID'))
avg_box_scores$Success=avg_box_scores$Success %>% replace_na(0) %>% factor()

#avg_box_scores.lda = lda(Success~Score+FGM+FGA+FGM3+FGA3+FTM+FTA+OR+DR+Ast+TO+Stl+Blk+PF, data = avg_box_scores)
avg_box_scores.lda = lda(Success~FGM+FGA+FGM3+FGA3+FTM+FTA+OR+DR+Ast+TO+Stl+Blk+PF, data = avg_box_scores)
#avg_box_scores.logistic = glm(formula = Success~Score+FGM+FGA+FGM3+FGA3+FTM+FTA+OR+DR+Ast+TO+Stl+Blk+PF, family = binomial, data = avg_box_scores)
avg_box_scores.lda
```
### Selection by Sabermetrics

For the per-season average of sabermetrics, the winning ratio and offensive rebound percentage perform significantly well in discriminate the tournament teams (having coefficient with high magnitude). The turnover ratio, true shoot percentage, free throw rate are also important in the discrimination.

```{r, message=FALSE}
avg_adv_stat = all_team_reg_season_adv_stats
avg_adv_stat = avg_box_scores %>% dplyr::select(Season,TeamID, Success) %>% right_join(avg_adv_stat, by=c('Season','TeamID'))

#avg_adv_stat.lda = lda(Success~WinRate+PtsDiff+Poss+OffRtg+DefRtg+SoS+Pie+Tie+AstRatio+TORatio+TSPerc+EFGPerc+FTR+ORbP+DRbP+TRbP, data = avg_adv_stat)
avg_adv_stat.lda = lda(Success~WinRate+Poss+OffRtg+DefRtg+Tie+AstRatio+TORatio+TSPerc+FTR+ORbP+DRbP, data = avg_adv_stat)
avg_adv_stat.logistic = glm(formula = Success~WinRate+PtsDiff+Poss+OffRtg+DefRtg+SoS+Pie+Tie+AstRatio+TORatio+TSPerc+EFGPerc+FTR+ORbP+DRbP+TRbP, family = binomial, data = avg_adv_stat)
avg_adv_stat.lda
```

```{r, message=FALSE}
box_adv_stat = avg_adv_stat %>% left_join(avg_box_scores, by = c('Season', 'TeamID', 'Success'))
box_adv_stat.lda = lda(Success~FGM+FGM3+FTM+FTA+OR+DR+Stl+Blk+PF+WinRate+Poss+OffRtg+DefRtg+Tie+AstRatio+TORatio+TSPerc+FTR+ORbP+DRbP, data = box_adv_stat)
#box_adv_stat.lda = lda(Success~FGM+FGA+FGM3+FGA3+FTM+FTA+OR+DR+Ast+TO+Stl+Blk+PF+WinRate+Poss+OffRtg+DefRtg+Tie+AstRatio+TORatio+TSPerc+FTR+ORbP+DRbP, data = box_adv_stat)
box_adv_stat.logistic = glm(formula = Success~FGM+FGA+FGM3+FGA3+FTM+FTA+OR+DR+Ast+TO+Stl+Blk+PF+WinRate+Poss+OffRtg+DefRtg+Tie+AstRatio+TORatio+TSPerc+FTR+ORbP+DRbP, family = binomial, data = box_adv_stat)
```
## Cross-validation Forcasting

To examine the performance of the linear discriminant model in picking the tournament teams, I use a cross-validation method to test the accuracy. I use the data from all seasons between the 2002–03 and 2018–19 seasons. I use one season as the validation and the rest as the training data each time.

formulate three different linear discriminant models: (1) using only box scores, (2) using only sabermetrics and two groups metrics together. Another logistic regression using both the box score and sabermetrics are also formulated. 

The final results shows that the third discriminant analysis models with both group of metrics and logistic regression model performed better, reaching a 90% accuracy picking the tournament team.


```{r, message=FALSE}
years = seq(2003,2019)
for (i in years){
    train = avg_box_scores %>% filter(Season != i)
    test = avg_box_scores %>% filter(Season == i)
    #avg_box_scores.lda = lda(Success~FGM+FGA+FGM3+FGA3+FTM+FTA+OR+DR+Ast+TO+Stl+Blk+PF, data = train)
    avg_box_scores.lda = lda(Success~FGM+FGM3+FTM+OR+DR+Ast+TO+Stl+Blk+PF, data = train)

    predicted.qda1 = predict(avg_box_scores.lda, newdata = test)
    table(test$Success, predicted.qda1$class, dnn = c('Actual Group','Predicted Group'))
    #cat("box_score Performance on year",i ,mean(predicted.qda$class == test$Success),"\n")
    
    train2 = avg_adv_stat %>% filter(Season != i)
    test2 = avg_adv_stat %>% filter(Season == i)
    #avg_adv_stat.lda = lda(Success~WinRate+Poss+OffRtg+DefRtg+Tie+AstRatio+TORatio+TSPerc+FTR+ORbP+DRbP, data = train2)
    avg_adv_stat.lda = lda(Success~WinRate+Poss+SoS+Tie+AstRatio+TORatio+TSPerc+FTR+TRbP, data = train2)
    predicted.qda2 = predict(avg_adv_stat.lda, newdata = test2)
    table(test$Success, predicted.qda2$class, dnn = c('Actual Group','Predicted Group'))
    #cat("adv_stats Performance on year",i ,mean(predicted.qda$class == test2$Success),"\n")
    
    train3 = box_adv_stat %>% filter(Season != i)
    test3 = box_adv_stat %>% filter(Season == i)
    #box_adv_stat.lda = lda(Success~FGM+FGM3+FTM+FTA+OR+DR+Stl+Blk+PF+WinRate+Poss+OffRtg+DefRtg+Tie+AstRatio+TORatio+TSPerc+FTR+ORbP+DRbP, data = train3)
    box_adv_stat.lda = lda(Success~FGM+FGM3+FTM+OR+DR+Ast+TO+Stl+Blk+PF+WinRate+Poss+SoS+Tie+AstRatio+TORatio+TSPerc+FTR+TRbP, data = train3)
    predicted.qda3 = predict(box_adv_stat.lda, newdata = test3)
    table(test$Success, predicted.qda3$class, dnn = c('Actual Group','Predicted Group'))
    #cat("box_adv_stat.lda Performance on year",i ,mean(predicted.qda$class == test3$Success),"\n")
    
    train4 = box_adv_stat %>% filter(Season != i)
    test4 = box_adv_stat %>% filter(Season == i)
    box_adv_stat.logistic = glm(formula = Success~FGM+FGM3+FTM+Stl+Blk+PF+WinRate+Poss+OffRtg+DefRtg+Tie+AstRatio+TORatio+TSPerc+FTR+ORbP+DRbP, family = binomial, data = train4)
    predicted.qda4 =  plogis(predict(box_adv_stat.logistic, test4))
    optCutOff <- optimalCutoff(test4$Success, predicted.qda4)[1] 
    predict_table = confusionMatrix(test4$Success, predicted.qda4, threshold = optCutOff)    
    cat("year", i, "box_score", mean(predicted.qda1$class == test$Success), "Adv_stats", mean(predicted.qda2$class == test2$Success), "mixed", mean(predicted.qda3$class == test3$Success), "logistic", (predict_table[1,1]+predict_table[2,2])/sum(predict_table), "\n")
}
```

# Reference
[1] Jason Zivkovic's notebook (https://www.kaggle.com/jaseziv83/a-recent-deep-look-at-the-men-s-ncaab#header)

[2] Kubatko, Justin, et al. "A starting point for analyzing basketball statistics." Journal of Quantitative Analysis in Sports 3.3 (2007).

[3] Smith, Tyler, and Neil C. Schwertman. "Can the NCAA basketball tournament seeding be used to predict margin of victory?." The American Statistician 53.2 (1999): 94-98.
