tiers.Rmd

---
title: "tiers"
author: "Huey Kwik"
date: "May 30, 2016"
output: html_document
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)

library(readr)
library(dplyr)
library(ggplot2)
library(mclust)
library(gsubfn)
library(XML)
library(tictoc)
library(ggrepel)

source("util.R")
source("tier-functions.R")

# Initialize constants
pos.list = c('qb','rb','wr','te','k','dst')

kYear = 2015

# Plotting
font = 3.5
barsize = 1.5  
dotsize = 2  
```

# Introduction

This builds upon Boris Chen's work in applying Gaussian Mixture Models to fantasy
football data. Instead of applying GMM to rankings, we apply them to projected 
scores. 

## Grabbing Data

I needed the following data for the 2015 season:

* Projected and actual scores from [Fantasy Football Analytics](apps.fantasyfootballanalytics.net/projections)
* Expert Consensus Rankings from FantasyPros.

The projected scores for each week are the average of the scores from CBS, ESPN,
FantasyFootball Nerd, FantasySharks, FFToday, NFL NumberFire, and YahooSports. 
I used Yahoo! standard scoring.

I wrote a modified version of [Chen's code](https://github.com/borisachen/fftiers) 
to download the ECR data:

```{r download ECR}
download.data = function(week, pos.list) {
  for (mp in pos.list) {
 	 	# Remove old data files.
    rmold1 = paste('rm data/2015/week-', week, '-',mp,'-raw.xls', sep='')
  	system(rmold1)
  	
  	urlStr = fn$identity("http://www.fantasypros.com/nfl/rankings/`mp`.php?week=`week`&export=xls")
  	dlPath = fn$identity("data/2015/week-`week`-`mp`-raw.xls")
  	
	  rmold2 = paste('rm ~ data/2015/week_', week, '_', mp, '.tsv', sep='')
  	system(rmold2)
    download.file(urlStr, destfile=dlPath)  
    sedstr = paste("sed '1,4d' data/2015/week-", week, '-',mp,'-raw.xls', 
  			  ' > data/2015/week_', week, '_', mp, '.tsv',sep="")
    system(sedstr)
  }	  
}

#  Error when trying to download using knitr, but works otherwise.
# sapply(1:17, function(week) download.data(week))
```

## Rank Tier example

Since Chen’s blog doesn’t have published tiers for Weeks 1–16, I used Chen’s code to generate Rank Tiers myself.
I stuck to the default number of tiers used in his code, but it’s possible that the tiers I generated differ slightly from the ones that were published on his blog during the season.

```{r computeTiers util function}
kCutoffs = c(QB = 24, RB = 40, WR = 60, TE = 24, K = 24, DST = 24)
kTiers = c(QB = 8, RB = 9, WR = 12, TE = 8, K = 5, DST = 6)

computeTiers = function(scores, k, reverse = FALSE) {
  # reverse: Useful for clustering on projected points. We want highest projected points to be rank 1, so 
  # we reverse the levels.
  clusters = NULL
  while (is.null(clusters)) {
    tryCatch({ 
      clusters = Mclust(scores, G = k)
      n_clusters = length(unique(clusters$classification))
    }, 
    warning = function(w) { warning(w); return(NULL) },
    error = function(e) { warning(e); return(NULL) })
    
    if (!is.null(clusters)) break 
    k = k - 1
  }
  
  n_clusters = length(unique(clusters$classification))
  tiers = factor(clusters$classification)
  if (reverse) {
    levels(tiers) = rev(levels(tiers))
    levels(tiers) = n_clusters:1
  } else {
    levels(tiers) = 1:n_clusters
  }
    
  return(tiers)
}
```

The tiers we generate should look pretty similar to what we see on Chen's site.
```{r rankTier example, warning = FALSE}
week = 17
pos = "QB"

ecr_df = read_ecr_data(fn$identity("data/2015/week_`week`_`tolower(pos)`.tsv"))
ecr_df = ecr_df[1:kCutoffs[pos], ]  # The data is ordered from best rank to worst rank  
ecr_df$TierRank = computeTiers(ecr_df$Avg.Rank, kTiers[pos])

ecr_df$nchar = nchar(as.character(ecr_df$Player.Name))  # For formatting later

# Calculate position rank, negative so lowest rank will be on top in the plot
# below
ecr_df$position.rank = -seq(nrow(ecr_df))

# We put Avg.Rank as y because geom_errorbar requires ymin/ymax. We then flip the 
# coordinates.
p = ggplot(ecr_df, aes(x = position.rank, y = Avg.Rank))
p = p + geom_errorbar(aes(ymin = Avg.Rank - Std.Dev/2, ymax = Avg.Rank + Std.Dev/2, width=0.2, colour=TierRank), size=barsize*0.8, alpha=0.4)
p = p + coord_flip()
p = p + geom_text(aes(label=Player.Name, colour=TierRank, y = Avg.Rank - nchar/6 - Std.Dev/1.4), size=font)
p = p + scale_x_continuous("Expert Consensus Rank")
p = p + ylab("Average Expert Rank")
p
```

## Projection Tier Example

Similarly, here is an example of computing tiers based on projections.
```{r, warning = FALSE}

sdf = read_score_data(fn$identity("data/2015/FFA-CustomRankings-Week-`week`.csv"), week)
pos_df = sdf %>% filter(position == pos, positionRank <= kCutoffs[pos]) %>% arrange(desc(points))

pos_df = sdf %>% filter(position == pos, positionRank <= kCutoffs[pos]) %>% arrange(desc(points))
pos_df$TierProj = computeTiers(pos_df$points, kTiers[pos], reverse = TRUE)

players_to_label = c("Cam Newton", "Ben Roethlisberger", "Drew Brees", "Ryan Fitzpatrick", "Derek Carr", "Brock Osweiler", "Ryan Tannehill", "Zach Mettenberger")

names(pos_df)[25] = "Tier"
pos_df$Tier = factor(pos_df$Tier, levels = rev(levels(pos_df$Tier)))
p = ggplot(pos_df, aes(x = positionRank, y = points))
p = p + geom_point(aes(colour=Tier)) + geom_text_repel(aes(label=playername, colour=Tier), data = subset(pos_df, playername %in% players_to_label))
p = p + scale_x_continuous("Position Rank")
p = p + ylab("Projected Points")
p 
```

## Evaluating Accuracy

[Tier accuracy](http://www.borischen.co/2013/11/week-10-retrospective.html) 
is defined as the percentage of time a higher tier results in a higher median score 
than a lower tier.

```{r computeAccuracy}
computeAccuracy = function(scores, tiers) {
  n_tiers = length(unique(tiers))
  
  # It's possible the clustering algorith does not find tiers, especially
  # on bootstrapped data.
  if (n_tiers == 1) {
    return(0)
  }
  n_correct = 0
  n_comparisons = 0
  
  medians = c()
  for (i in seq(n_tiers)) {
    median = median(scores[tiers == i])
    medians = c(medians, median)
  }
  
  for (i in seq(n_tiers-1)) {
    curr_median = medians[i]
    below_medians = medians[(i+1):n_tiers]
  
    comparisons = sapply(below_medians, function(x) curr_median > x)
  
    n_correct = n_correct + sum(comparisons)
    n_comparisons = n_comparisons + length(comparisons)
  }
  
  return(n_correct / n_comparisons)
}
```

# Computing accuracy for Rank Tiers and Projection Tiers
Here we compute the accuracy for all weeks of the 2015 season. In the code below,
we also try many random partitions.

```{r 2015, warning = FALSE}

joinECRScores = function(ecr_df, all_pos_df, pos) {
  if (pos == "DST") {
    ecr_scores = inner_join(ecr_df, all_pos_df, by = c("Team" = "team"))
  } else {
    ecr_scores = inner_join(ecr_df, all_pos_df, by = c("Player.Name" = "playername"))
  }
  return(ecr_scores)
}

tier_ecr_accuracy = data.frame()
tier_proj_accuracy = data.frame()
tier_proj_random_accuracy = data.frame()

tic(quiet = FALSE)
set.seed(1) 
kTries = 1 #200000  #200000
for (week in seq(17)) {
  # Parser warnings are due to trailing comma at the end of each row
  sdf = read_score_data(fn$identity("data/2015/FFA-CustomRankings-Week-`week`.csv"), week)
  
  week_tier_ecr_accuracy = c()
  week_tier_proj_accuracy = c()
  week_tier_proj_random_accuracy = c()
  
  for (pos in toupper(pos.list)) {
    ## Calculate TierProj
    pos_df = sdf %>% filter(position == pos, positionRank <= kCutoffs[pos]) %>% arrange(desc(points))

    pos_df$TierProj = computeTiers(pos_df$points, kTiers[pos], reverse = TRUE)

    ## Calculate TierProjRandomPartition
    best_indexes = bestPartition(kTiers[pos], pos_df$points, kTries)
    pos_df$TierProjRandom = expand(1:kTiers[pos], best_indexes, nrow(pos_df))

    ## Calculate TierRank
    ecr_df = read_ecr_data(fn$identity("data/2015/week_`week`_`tolower(pos)`.tsv"))
    ecr_df = ecr_df[1:kCutoffs[pos], ]  # The data is ordered from best rank to worst rank  
    ecr_df$TierRank = computeTiers(ecr_df$Avg.Rank, kTiers[pos])
    
    all_pos_df = sdf %>% filter(position == pos)
    ecr_scores = joinECRScores(ecr_df, all_pos_df, pos)
    
    acc_proj = computeAccuracy(pos_df$actualPoints, pos_df$TierProj)
    acc_proj_random = computeAccuracy(pos_df$actualPoints, pos_df$TierProjRandom)
    acc_ecr = computeAccuracy(ecr_scores$actualPoints, ecr_scores$TierRank)
    
    week_tier_proj_accuracy  = c(week_tier_proj_accuracy, acc_proj)
    week_tier_ecr_accuracy = c(week_tier_ecr_accuracy, acc_ecr)
    week_tier_proj_random_accuracy = c(week_tier_proj_random_accuracy, acc_proj_random)
  }

  tier_proj_accuracy = rbind(tier_proj_accuracy, week_tier_proj_accuracy)
  tier_ecr_accuracy = rbind(tier_ecr_accuracy, week_tier_ecr_accuracy)
  tier_proj_random_accuracy = rbind(tier_proj_random_accuracy, week_tier_proj_random_accuracy)

  names(tier_proj_accuracy) = toupper(pos.list)
  names(tier_proj_random_accuracy) = toupper(pos.list)
  names(tier_ecr_accuracy) = toupper(pos.list)
}
toc()
```

# Results

```{r}
  round(tier_proj_accuracy * 100, 1)
  round(tier_ecr_accuracy * 100, 1)
  round(tier_proj_random_accuracy * 100, 1)

  mean_proj = sapply(tier_proj_accuracy, mean)
  mean_tier = sapply(tier_ecr_accuracy, mean)
  mean_proj_random = sapply(tier_proj_random_accuracy, mean)
  
  round(mean_proj * 100, 1)
  round(mean_tier * 100, 1)
  round(mean_proj_random * 100, 1)
```

Projection Tiers outperform Rank Tiers for all the positions listed above for
the 2015 season.

# Bootstrapping

To get a better sense of how these two algorithms compete on more data, there 
are two approaches we can take. 

The most obvious approach is to look at historical data. Unfortunately, FantasyPros
does not keep a record of ECR from past seasons. I tried using the 
[Wayback Machine](https://archive.org/web/), but there simply were not enough 
snapshots to make this worthwhile. 

Chen has published tiers for the 2013-2015 seasons on his blog, but the data isn't
always available. For instance, some of the links to the tiers in the 2014 season
are dead. However, the 2013 season is available so I might try that next.

Instead, we create bootstrap samples for the 2015 season and then see how well
each method performs. 

```{r bootstrapping, warning = FALSE}

bootstrapCompare = function(n) {
  results_df = data.frame()
  
  for (week in seq(17)) {
    sdf = read_score_data(fn$identity("data/2015/FFA-CustomRankings-Week-`week`.csv"), week)
    row = c()
    for (pos in toupper(pos.list)) {
      n_proj_wins = 0
      pos_df = sdf %>% filter(position == pos)
      ecr_df = read_ecr_data(fn$identity("data/2015/week_`week`_`tolower(pos)`.tsv"))
      ecr_scores = joinECRScores(ecr_df, pos_df, pos) 
      for (i in seq(n)) {
        sample = ecr_scores[sample(nrow(ecr_scores), nrow(ecr_scores), replace = TRUE), ]
    
        #pos_df = sample %>% filter(positionRank <= kCutoffs[pos])
        pos_df = sample %>% arrange(desc(points))
        pos_df = pos_df[1:kCutoffs[pos], ]
        # select(pos_df, Player.Name, points, actualPoints)  # For debugging
        pos_df$TierProj = computeTiers(pos_df$points, kTiers[pos], reverse = TRUE)
        select(pos_df, Player.Name, points, TierProj)
        
        ## CalculateTierRank
        ecr_df = sample %>% arrange(Rank)
        ecr_df = ecr_df[1:kCutoffs[pos], ]
        #select(ecr_df, Player.Name, points, actualPoints, Avg.Rank)  # For debugging
        ecr_df$TierRank = computeTiers(ecr_df$Avg.Rank, kTiers[pos])
        select(ecr_df, Player.Name, Avg.Rank, TierRank)
        
        acc = computeAccuracy(pos_df$actualPoints, pos_df$TierProj)
        
        acc_ecr = computeAccuracy(ecr_df$actualPoints, ecr_df$TierRank)
        
        n_proj_wins = n_proj_wins + (acc > acc_ecr) 
      }
      row = c(row, n_proj_wins / n)
    }
    results_df = rbind(results_df, row)
  }
  names(results_df) = toupper(pos.list)
  return(results_df)
}
```

```{r bootstrap results, warning = FALSE}
set.seed(1)
sapply(bootstrapCompare(100), mean)
```

The numbers above represent when Tier Projections was more accurate than Tier Ranks.
Interestingly enough, Tier Projections seems to do better for QBs, RBs, and WRs, but worse for the other positions. 
In particular, it performs significantly worse for Defense, but that seems to be corroborated by the fact that [projections perform worse than rankings for Defense (and Kickers)](http://fantasyfootballanalytics.net/2016/04/accuracy-of-rankings-vs-projections.html).