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.

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 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 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}} } ]\]

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)

Install package from CRAN or development version from github.

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

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

`gpheats`

- results SGP heats. Column`rank`

is a numeric version of column`position`

- rider position in race.`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
```

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
```

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
#> [1] "dbl"
#>
#> $`Overall Accuracy`
#> [1] 0.6167273
#>
#> $`Number of pairs`
#> [1] 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)
#> [1] "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.

`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
```

`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
```

`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)
#> [1] FALSE
```

`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)
```

`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
```