# sport an R package for online update algorithms

#### 2019-01-07

Name sport is an abbreviation for Sequential Pairwise Online Rating Techniques. Package contains functions calculating ratings for two-player or multi-player matchups. Methods included in package are able to estimate ratings (players strengths) and their evolution in time, also able to predict output of challenge. Algorithms are based on Bayesian Approximation Method, and they don’t involve any matrix inversions nor likelihood estimation. sport incorporates methods such glicko, glicko2, bayesian Bradley-Terry, dynamic logistic regression. Parameters are updated sequentially, and computation doesn’t require any additional RAM to make estimation feasible. Additionally, base of the package is written in C++ what makes sport computation even faster.

# Theory

Problem of sport matchups falls into subject of paired comparison modeling and choice modeling. Estimating player skills is equivalent to estimating preference of choice between two alternatives. Just as one product is more preferred over another to buy, similarly better player is more preferred to win over worst. As player/event and alternative/experiment can be used interchangeably, for ease of use sport nomenclature is adapted (player/event).

Algorithms implemented in a sport package works similarly, as all using Bayesian Approximation Method. Algorithms works as follows: At the moment player i competes with player j while both have initial $$R_i$$ and $$R_j$$ ratings. Prior to event, probability that player i win over player j is $$\hat{Y_i}$$. After event is finished when true result $$Y_{ij}$$ is observed, initial believe about rating is changed $$R_i^{'} \leftarrow R_i$$ according to the prediction error $$( Y_{ij} - \hat{Y_{ij}} )$$ and some constant $$K$$. Updates are summed as player can compete with more than one player in particular event.

$\large R_i^{'} \leftarrow R_i + \sum_{j \neq i}{ K * ( Y_{ij} - \hat{Y_{ij}}} )$ Where: $\large \hat{Y} = P(X_i > X_j)$ $K - learning rate$

Outcome probability function is based on Bradley-Terry model designed to predict outcome of pairwise comparison. For multi-player matchups where output is a ranking, sport package uses the same data transformation as in exploded logit - ranking is then presented as combination of all possible pairs competing within same event.

## Glicko rating system

Glicko is the first bayesian online update algorithm incorporating rating volatility to rating and outcome computation. Glicko system is not balanced, and sum of rating rewards of all players are not zero. In one 2-players event, reward of player i differs from reward of player q as it depends on their individual ratings deviation. Rating values oscillates around r=1500 with max deviation rd<=350.

For deeper knowledge read Mark E. Glickman (1999).

Update Rules:

$\hat{Y_{ij}} = P(X_i>X_j) = \frac{1}{ 1 + 10^{-g(RD_{ij}) * (R_i-R_j)/400}}$

${R'}_i = R_i + \frac{1}{\frac{1}{{RD}^2_{i}} + \frac{1}{d^2_i}} * \sum_j{g(RD_j) * (Y_{ij} - \hat{Y_{ij}}) }$

${RD'}_i = \sqrt{(\frac{1}{{RD}^2_{i}} + \frac{1}{d^2_i}})^{-1}$

## Glicko2 rating system

Glicko2 improved predecessor by adding volatile parameter $$\sigma_i$$ which increase/decrease rating deviation in periods when player performance differs from expected. Sigma is estimated iteratively using Illinois algorithm, which converges quickly not affecting computation time. Rating values oscillates around r=1500 with max deviation rd<=350.

For further knowledge read Mark E. Glickman (2013)

$\hat{Y_{ij}} = \frac{1}{1 + e^{-g(\phi_{ij})*(\mu_i - \mu_j)} }$

${\phi'}_i = \frac{1}{\sqrt{ \frac{1}{ { {\phi_i}^2 + {\sigma'_i}^2}} + \frac{1}{v} }}$

${\mu'_i} = \mu_i + {\phi'}_i * \sum_j{g(\phi_j)*(Y_{ij} - \hat{Y_{ij}})}$

The fastest algorithm with simple formula. Original BT formula lacks variance parameter, and this method incorporates rating deviation into model. BBT also prevents against fast rd decline to zero using gamma and kappa.

For further knowledge read Ruby C. Weng and Chih-Jen Lin (2011)

$\hat{Y_{ij}} = P(X_i>X_j) = \frac{e^{R_i/c_{i_j}}}{e^{R_i/c_{ij}} + e^{R_j/c_{ij}}}$

${R'}_i = R_i + \sum_j{\frac{RD_i^2}{c_{ij}}*(Y_{ij} - \hat{Y_{ij}})}$

${RD'}_i = RD_i * [ 1 - \frac{RD_{ij}^2}{RD_i^2}\sum_j{ \gamma_j * (\frac{RD_i}{c_{ij}})^2* \hat{Y_{ij}}\hat{Y_{ji}} } ]$

## Dynamic Bayesian Logit

Following algorithm gives some advantages over mentioned rating systems, adding other important factors to estimation process making final ratings unbiased. Algorithm perform better in disciples where other variables can make a difference in result eg. home field advantage. DBL implements Extended Kalman Filter learning rule, and allows to estimate multiple parameters in addition to player ratings. DBL is a Dynamic Logit extended to usage in pairwise comparisons by modeling differences in players characteristics. Classic Bradley-Terry model is enriched by moderation element $$K(s_t)$$ which adds prior uncertainty to output prediction.

$\hat{Y_{ij}} = \frac{ e^{-K(s_t)w _t^T(x_{it}-x_{jt})} }{1+e^{-K(s_t)w _t^T(x_{it}-x_{jt})}}$ Parameters for player i competing with player j are estimated using EKF update rule. $\hat{\omega}_{it} = \hat{\omega}_{i(t-1)} + \frac{RD^2_{i(t-1)}}{1+\hat{Y_{ij}} (1-\hat{Y_{ij}})} x_t (Y_{ij} - \hat{Y_{ij}})$ $RD^2_{i t} = RD^2_{i(t-1)} - \frac{\hat{Y_{ij}}(1-\hat{Y_{ij}})}{1+\hat{Y_{ij}} (1-\hat{Y_{ij}})s_t^2}(RD^2_{i(t-1)}x_i)(RD^2_{i(t-1)}x_i)^T$

For further knowledge read Stephen J. Roberts, William Penny (2011)

# Package Usage

## Installation

Install package from CRAN or development version from github.

# devtools::install_github("gogonzo/sport")
# install.packages("sport")
library(sport)

## Available Data

Package contains actual data from Speedway Grand-Prix. There are two data.frames:

1. gpheats - results SGP heats. Column rank is a numeric version of column position - rider position in race.
2. gpsquads - summarized results of the events, with sum of point and final position.
str(gpheats)
#> 'data.frame':    20649 obs. of  11 variables:
#>  $id : num 1 1 1 1 2 2 2 2 3 3 ... #>$ season  : int  1995 1995 1995 1995 1995 1995 1995 1995 1995 1995 ...
#>  $date : POSIXct, format: "1995-05-20 18:00:00" "1995-05-20 18:00:00" ... #>$ round   : int  1 1 1 1 1 1 1 1 1 1 ...
#>  $name : chr "Speedway Grand Prix of Poland" "Speedway Grand Prix of Poland" "Speedway Grand Prix of Poland" "Speedway Grand Prix of Poland" ... #>$ heat    : int  1 1 1 1 2 2 2 2 3 3 ...
#>  $field : int 1 2 3 4 1 2 3 4 1 2 ... #>$ rider   : chr  "Tomasz Gollob" "Gary Havelock" "Chris Louis" "Tony Rickardsson" ...
#>  $points : int 2 0 3 1 3 0 1 2 0 2 ... #>$ position: chr  "2" "4" "1" "3" ...
#>  $rank : num 2 4 1 3 1 4 3 2 4 2 ... Data used in sport package must be in so called long format. Typically data.frame contains at least id, name and rank, with one row for one player within specific match. Package allows for any number of players within event and allows ties also. For all games, output needs to be a rank/position in event. Don’t mix up rank output with typical 1-win, 0-lost. In sport package output for two player game is 1-winner 2-looser. Below example of two matches with 4 players each. #> id rider rank #> 1 1 Tomasz Gollob 2 #> 2 1 Gary Havelock 4 #> 3 1 Chris Louis 1 #> 4 1 Tony Rickardsson 3 #> 5 2 Sam Ermolenko 1 #> 6 2 Jan Staechmann 4 #> 7 2 Tommy Knudsen 3 #> 8 2 Henrik Gustafsson 2 ## Estimate dynamic ratings To compute ratings using each algorithms one has to specify formula. Form rank | id ~ name is required, which estimates name - rating of a player, by observing outputs - rank, nested within particular event - id. Variable names in formula are unrestricted, but model structure remains the same. All methods are named method_run. formula = rank|id ~ name glicko <- glicko_run( formula = rank|id ~ rider, data = gpheats ) glicko2 <- glicko2_run( formula = rank|id ~ rider, data = gpheats ) bbt <- bbt_run( formula = rank|id ~ rider, data = gpheats ) dbl <- dbl_run( formula = rank|id ~ rider, data = gpheats ) print(dbl) #> #> Call: rank | id ~ rider #> #> Number of unique pairs: 31081 #> #> Accuracy of the model: 0.62 #> #> True probabilities and Accuracy in predicted intervals: #> Interval Model probability True probability Accuracy n #> 1: [0,0.1] 0.076 0.164 0.835 436 #> 2: (0.1,0.2] 0.158 0.252 0.747 2217 #> 3: (0.2,0.3] 0.257 0.300 0.699 5447 #> 4: (0.3,0.4] 0.353 0.372 0.626 10113 #> 5: (0.4,0.5] 0.450 0.454 0.545 12911 #> 6: (0.5,0.6] 0.550 0.546 0.545 12825 #> 7: (0.6,0.7] 0.647 0.628 0.626 10113 #> 8: (0.7,0.8] 0.743 0.700 0.699 5447 #> 9: (0.8,0.9] 0.842 0.748 0.747 2217 #> 10: (0.9,1] 0.924 0.836 0.835 436 ## Output Objects returned by method_run are of class rating and have their own print summary which provides most important informations. -print.sport shows condensed informations about model performance like accuracy and consistency of model predictions with observed probabilities. More profound summarization are given by summary by showing ratings, ratings deviations and comparing model win probabilities with observed. summary(dbl) #>$formula
#> rank | id ~ rider
#>
#> $method #>  "dbl" #> #>$Overall Accuracy
#>  0.6167273
#>
#> $Number of pairs #>  62162 #> #>$r
#>                  name      r    rd Model probability True probability
#>   1:    Tomasz Gollob  1.329 0.005             0.572            0.587
#>   2:    Gary Havelock  0.696 0.111             0.410            0.476
#>   3:      Chris Louis  0.205 0.023             0.523            0.505
#>   4: Tony Rickardsson  1.793 0.012             0.684            0.703
#>   5:    Sam Ermolenko  0.552 0.064             0.534            0.546
#>  ---
#> 209:   Justin Sedgmen -0.508 0.511             0.347            0.267
#> 210:    Rohan Tungate  2.880 0.712             0.537            1.000
#> 211:    Maksym Drabik  0.693 0.731             0.297            0.429
#> 212:       Dan Bewley -0.271 0.837             0.307            0.250
#> 213:       Joel Kling  0.644 0.831             0.305            0.500
#>      Accuracy pairings
#>   1:    0.591     2718
#>   2:    0.659      126
#>   3:    0.580      605
#>   4:    0.699     1313
#>   5:    0.523      216
#>  ---
#> 209:    0.800       15
#> 210:    0.625        8
#> 211:    0.714        7
#> 212:    0.750        4
#> 213:    0.500        4

To visualize top n ratings with their 95% confidence interval one can use dedicated plot.rating function. For “bdl” method top coefficients are presented not necessarily ratings. It’s also possible to examine ratings evolution in time, by specifying players argument.

plot(glicko, n=15) plot(glicko, players = c("Greg Hancock","Nicki Pedersen","Jason Crump")) Except dedicated print,summary and plot there is possibility to extract more detailed information to be analyzed. rating object contains following elements:

names(glicko)
#>  "final_r"  "final_rd" "r"        "pairs"
• rating$final_r and rating$final_rd contains ratings and ratings deviations estimations.
• r contains data.frame with sequential ratings estimations from first event to the last. Number of rows in r equals number of rows in input data.
• pairs pairwise combinations of players in analyzed events with prior probability and result of a challenge.
tail(glicko$r) #> id name r rd p_win #> 1: 5154 Martin Vaculik 1556.186 15.598332 0.1134510 #> 2: 5154 Patryk Dudek 1668.852 19.472240 0.4256283 #> 3: 5155 Matej Žagar 1587.954 9.007093 0.2132507 #> 4: 5155 Fredrik Lindgren 1588.903 8.318461 0.2167572 #> 5: 5155 Martin Vaculik 1554.414 15.552804 0.1448214 #> 6: 5155 Nicki Pedersen 1652.299 6.677019 0.4251707 tail(glicko$pairs)
#>      id           name         opponent         P Y
#> 1: 5155 Martin Vaculik      Matej Žagar 0.4548475 0
#> 2: 5155 Martin Vaculik Fredrik Lindgren 0.4528676 0
#> 3: 5155 Martin Vaculik   Nicki Pedersen 0.3656971 0
#> 4: 5155 Nicki Pedersen      Matej Žagar 0.5914313 1
#> 5: 5155 Nicki Pedersen Fredrik Lindgren 0.5895059 1
#> 6: 5155 Nicki Pedersen   Martin Vaculik 0.6343029 1

Examples presented in package overview might be sufficient in most cases, but sometimes it is necessary to adjust algorithms to fit data better. One characteristic of the online update algorithms is that variance other parameters quickly drops to zero. Especially, when the number of events for the player is big ($n_i>100$), after hundreds iterations rating parameters are very difficult to change, and output probabilities are so extreme. To avoid these mistakes some additional controls should be applied, which is explained in this section with easy to learn examples.

## Prior beliefs about r and rd

Main functionality which is common between all algorithms is to specify prior r and rd. Both parameters can be set by creating named vectors. Let’s suppose we have 4 players c("A","B","C","D") competing in an event, and we have players prior r and rd estimates. It’s important to have r and rd names corresponding with levels of name variable. One can run algorithm, to obtain new estimates.

library(dplyr); library(magrittr) # for examples purpose
data <- data.frame( id = 1, name = c( "A", "B", "C", "D" ),  rank  = c( 3, 4, 1, 2 ))
r     <- setNames( c(1500, 1400, 1550, 1700), c("A","B","C","D") )
rd    <- setNames( c(200, 30, 100, 300), c("A","B","C","D") )

model <- glicko_run(rank|id ~ name, data=data, r=r, rd=rd)
print(model$final_r) #> A B C D #> 1464.297 1396.039 1606.521 1674.836 We can also run models using previously estimated parameters from model$final_r and model$final_rd in next event. data2 <- data.frame( id = 2, name = c( "A", "B", "C", "D" ), rank = 1:4 ) r <- model$final_r
rd    <- model$final_rd glicko_run(rank|id ~ name, data, r=r, rd=rd)$final_r
#>        A        B        C        D
#> 1451.417 1392.132 1650.762 1661.553

## Tuning update size by weight

All algorithms have a weight argument which increases or decreases update change. Higher weight increasing impact of corresponding event. Effect of the weight on update size can be expressed directly by following formula - $$\small R_i^{'} \leftarrow R_i + weight \sum_{j \neq i}{ K ( Y_{ij} - \hat{Y_{ij}}} )$$. To specify weight one needs to add create additional column in input data, and name a column in weight argument. For example weight could depend on importance of competition. In speedway Grand-Prix last three heats determine event winner and thus have higher weight.

gpheats %<>% mutate(weight = ifelse(heat >= (max(heat)-3),2,1) )
glicko  <- glicko_run(rank|id ~ rider, gpheats, weight="weight")
#> r is missing and will set to default=1500
#> rd is missing and will set to default=350

## Avoiding excessive RD shringe with kappa

In situation when games are playing very frequently by a player, rd can quickly decrease to zero, making further changes limited. Setting kappa disallows deviation decrease to be lower than specified fraction of rd. In other words final rd can’t be lower than initial RD times kappa $$\small RD' \geq RD * kappa$$.

bbt1 <- bbt_run(rank|id~rider, gpheats,kappa=0.9)
#> r is missing and will set to default=25
#> rd is missing and will set to default=8.33333333333333
bbt2 <- bbt_run(rank|id~rider, gpheats,kappa=0.5)
#> r is missing and will set to default=25
#> rd is missing and will set to default=8.33333333333333

all(bbt1$final_rd > bbt2$final_rd)
#>  FALSE

## Control output uncertainty by beta

To increase/decrease uncertainty of challenge output one can define beta. Increasing beta flatten probabilities of winning challenge equalizing competitors chances. This element differs from sigma by not affecting rd directly but to adjust prediction. This argument can be used in bbt and dbl by identifying column in input data where beta is stored. It is worth to specify beta in events when level of competition is higher than usual e.g. derby, knock-out phase, decisive matches, exhibitions and friendlies. Let’s assume that beta follows the same rule as previously set weight.

gpheats %<>% mutate(beta = ifelse(heat >= (max(heat)-3),2,1)  )
dbl <- dbl_run(rank|id~rider, beta="beta", data=gpheats) 

## Influence ratings volatility by sigma

In some periods player ratings tends to fluctuate more, when players form is more unstable. Sigma is na additional parameter controlling volatility of the ratings, making update bigger for higher values of sigma. Strictly, sigma can be expressed by $$\small rd_i' = \sqrt{rd_i^2+\sigma_i^2}$$. Sigma can be specified in two ways: - In glicko2 sigma is another parameter per player estimated simultaneously with r and rd. One can specify prior sigma in the same manner as r and rd by creating named vector of positive values. Sigma values depends on a particular application, and priors should be chosen by user.

sigma   <- unique(gpheats$rider) %>% setNames( runif(0.1,0.5, n = length(.)) , . ) glicko2 <- glicko2_run(rank|id~rider,gpheats, sigma=sigma) #> r is missing and will set to default=1500 #> rd is missing and will set to default=350 In glicko and bbt sigma needs to be fixed before iteration by adding new column to input data. For example sigma can be column of time difference between events in which player participated - The longer break in competing, the higher sigma is. It is recommended to keep values close to data$sigma.

# bbt example
gpheats %<>%
group_by(rider) %>%
mutate( days_scaled = as.integer(date - lag(date))/90,
days_scaled = if_else(days_scaled>1,1.0, days_scaled)) %>%
filter( !is.na(days_scaled) )

bbt <- bbt_run( rank|id ~ rider, data=gpheats, sigma="days_scaled")
#> r is missing and will set to default=25
#> rd is missing and will set to default=8.33333333333333