Learn R Programming

bpcs - A package for Bayesian Paired Comparison analysis with Stan

The bpcs package performs Bayesian estimation of Paired Comparison models utilizing Stan, such as variations of the Bradley-Terry (Bradley and Terry 1952) and the Davidson models (Davidson 1970).

Package documentation and vignette articles can be found at: https://davidissamattos.github.io/bpcs/

Installation

For the bpcs package to work, we rely upon the Stan software and the rstan package (Stan Development Team 2020).

To install the development version of the bpcs package, install directly from the Github repository.

remotes::install_github('davidissamattos/bpcs')

After installing, we load the package with:

library(bpcs)

Minimal example

The main function of the package is the bpc function. For the simple Bradley-Terry model, this function requires a specific type of data frame that contains:

  • Two columns containing the name of the contestants in the paired comparison
  • Two columns containing the score of each player OR one column containing the result of the match (0 if player0 won, 1 if player1 won, 2 if it was a tie)

We will utilize the tennis dataset available (Agresti 2003). The dataset can be seen below and is available as data(tennis_agresti):

dplyr::sample_n(tennis_agresti,10) %>% 
  knitr::kable()

player0

player1

y

id

Seles

Sanchez

0

14

Sabatini

Sanchez

1

42

Seles

Graf

0

1

Graf

Sabatini

1

22

Sabatini

Sanchez

1

41

Seles

Navratilova

1

12

Graf

Sanchez

0

32

Sabatini

Navratilova

0

35

Sabatini

Sanchez

0

39

Seles

Graf

1

4

Based on the scores of each contestant, the bpc function computes automatically who won the contest. Alternatively, you can provide a vector of who won if that is already available (for more information see ?bpc.

For the simple Bradley Terry Model we specify the model type as 'bt'. Here we hide the MCMC sampler chain messages for simplicity in the output.

m<-bpc(data = tennis_agresti, #datafrane
       player0 = 'player0', #name of the column for player 0
       player1 = 'player1', #name of the column for player 1
       result_column = 'y', #name of the column for the result of the match
       model_type = 'bt', #bt = Simple Bradley Terry model
       solve_ties = 'none' #there are no ties in the dataset so we can choose none here
       )
#> 
#> SAMPLING FOR MODEL 'bt' NOW (CHAIN 1).
#> Chain 1: 
#> Chain 1: Gradient evaluation took 8.7e-05 seconds
#> Chain 1: 1000 transitions using 10 leapfrog steps per transition would take 0.87 seconds.
#> Chain 1: Adjust your expectations accordingly!
#> Chain 1: 
#> Chain 1: 
#> Chain 1: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 1: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 1: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 1: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 1: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 1: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 1: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 1: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 1: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 1: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 1: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 1: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 1: 
#> Chain 1:  Elapsed Time: 0.299618 seconds (Warm-up)
#> Chain 1:                0.312838 seconds (Sampling)
#> Chain 1:                0.612456 seconds (Total)
#> Chain 1: 
#> 
#> SAMPLING FOR MODEL 'bt' NOW (CHAIN 2).
#> Chain 2: 
#> Chain 2: Gradient evaluation took 3.9e-05 seconds
#> Chain 2: 1000 transitions using 10 leapfrog steps per transition would take 0.39 seconds.
#> Chain 2: Adjust your expectations accordingly!
#> Chain 2: 
#> Chain 2: 
#> Chain 2: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 2: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 2: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 2: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 2: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 2: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 2: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 2: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 2: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 2: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 2: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 2: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 2: 
#> Chain 2:  Elapsed Time: 0.284584 seconds (Warm-up)
#> Chain 2:                0.344625 seconds (Sampling)
#> Chain 2:                0.629209 seconds (Total)
#> Chain 2: 
#> 
#> SAMPLING FOR MODEL 'bt' NOW (CHAIN 3).
#> Chain 3: 
#> Chain 3: Gradient evaluation took 4.1e-05 seconds
#> Chain 3: 1000 transitions using 10 leapfrog steps per transition would take 0.41 seconds.
#> Chain 3: Adjust your expectations accordingly!
#> Chain 3: 
#> Chain 3: 
#> Chain 3: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 3: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 3: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 3: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 3: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 3: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 3: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 3: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 3: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 3: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 3: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 3: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 3: 
#> Chain 3:  Elapsed Time: 0.305793 seconds (Warm-up)
#> Chain 3:                0.292243 seconds (Sampling)
#> Chain 3:                0.598036 seconds (Total)
#> Chain 3: 
#> 
#> SAMPLING FOR MODEL 'bt' NOW (CHAIN 4).
#> Chain 4: 
#> Chain 4: Gradient evaluation took 4.2e-05 seconds
#> Chain 4: 1000 transitions using 10 leapfrog steps per transition would take 0.42 seconds.
#> Chain 4: Adjust your expectations accordingly!
#> Chain 4: 
#> Chain 4: 
#> Chain 4: Iteration:    1 / 2000 [  0%]  (Warmup)
#> Chain 4: Iteration:  200 / 2000 [ 10%]  (Warmup)
#> Chain 4: Iteration:  400 / 2000 [ 20%]  (Warmup)
#> Chain 4: Iteration:  600 / 2000 [ 30%]  (Warmup)
#> Chain 4: Iteration:  800 / 2000 [ 40%]  (Warmup)
#> Chain 4: Iteration: 1000 / 2000 [ 50%]  (Warmup)
#> Chain 4: Iteration: 1001 / 2000 [ 50%]  (Sampling)
#> Chain 4: Iteration: 1200 / 2000 [ 60%]  (Sampling)
#> Chain 4: Iteration: 1400 / 2000 [ 70%]  (Sampling)
#> Chain 4: Iteration: 1600 / 2000 [ 80%]  (Sampling)
#> Chain 4: Iteration: 1800 / 2000 [ 90%]  (Sampling)
#> Chain 4: Iteration: 2000 / 2000 [100%]  (Sampling)
#> Chain 4: 
#> Chain 4:  Elapsed Time: 0.30643 seconds (Warm-up)
#> Chain 4:                0.300302 seconds (Sampling)
#> Chain 4:                0.606732 seconds (Total)
#> Chain 4:

If rstan is available and correctly working this function should sample the posterior distribution and create a bpc object.

To see a summary of the results we can run the summary function. Here we get three tables:

  1. The parameters of the model
  2. The probabilities of one player beating the other (this probability is based on the predictive posterior distribution)
  3. The rank of the player based on their abilities (this rank is based on the predictive posterior ranks).
summary(m)
#> Estimated baseline parameters with HPD intervals:
#> 
#> 
#> Parameter               Mean   HPD_lower   HPD_higher    n_eff   Rhat
#> --------------------  ------  ----------  -----------  -------  -----
#> lambda[Seles]           0.51       -2.37         3.34   663.50   1.01
#> lambda[Graf]            0.94       -1.85         3.71   622.34   1.01
#> lambda[Sabatini]       -0.33       -3.12         2.50   662.40   1.01
#> lambda[Navratilova]     0.04       -2.87         2.91   637.48   1.01
#> lambda[Sanchez]        -1.12       -4.04         1.64   670.86   1.01
#> NOTES:
#> * A higher lambda indicates a higher team ability
#> 
#> 
#> Posterior probabilities:
#> These probabilities are calculated from the predictive posterior distribution
#> for all player combinations
#> 
#> 
#> i             j              i_beats_j
#> ------------  ------------  ----------
#> Graf          Navratilova         0.70
#> Graf          Sabatini            0.77
#> Graf          Sanchez             0.86
#> Graf          Seles               0.60
#> Navratilova   Sabatini            0.58
#> Navratilova   Sanchez             0.74
#> Navratilova   Seles               0.41
#> Sabatini      Sanchez             0.66
#> Sabatini      Seles               0.32
#> Sanchez       Seles               0.20
#> 
#> 
#> Rank of the players' abilities:
#> The rank is based on the posterior rank distribution of the lambda parameter
#> 
#> 
#> Parameter              MedianRank   MeanRank   StdRank
#> --------------------  -----------  ---------  --------
#> lambda[Graf]                    1       1.42      0.64
#> lambda[Seles]                   2       2.08      0.88
#> lambda[Navratilova]             3       3.04      0.91
#> lambda[Sabatini]                4       3.65      0.83
#> lambda[Sanchez]                 5       4.81      0.48

Features of the bpcs package

  • Bayesian computation of different variations of the Bradley-Terry (including with home advantage, random effects and the generalized model).
  • Bayesian computation of different variations of the Davidson model to handle ties in the contest (including with home advantage, random effects and the generalized model).
  • Accepts a column with the results of the contest or the scores for each player.
  • Customize a normal prior distribution for every parameter.
  • Compute HDP interval for every parameter with the get_hpdi_parameters function
  • Compute rank of the players with the get_rank_of_players function.
  • Compute all the probability combinations for one player beating the other with the get_probabilities function.
  • Convert aggregated tables of results into long format (one contest per row) with the expand_aggregated_data.
  • Obtain the posterior distribution for every parameter of the model with the get_sample_posterior function.
  • Easy predictions using the predict function.
  • We do not reinforce any table or plotting library! Results are returned as data frames for easier plotting and creating tables
  • We reinforce the need to manually specify the model to be used.

Models available

  • Bradley-Terry (bt) (Bradley and Terry 1952)
  • Davidson model (davidson) for handling ties (Davidson 1970)

Options to add to the models:

  • Order effect (-ordereffect). E.g. for home advantage (Davidson and Beaver 1977)
  • Generalized models (-generalized). When we have contestant specific predictors (Springall 1973)
  • Intercept random effects (-U). For example, to compensate clustering or repeated measures (Böckenholt 2001)

E.g.:

  • Simple BT model: bt
  • Davidson model with random effects: davidson-U
  • Generalized BT model with order effect: bt-generalized-ordereffect

Notes:

  • The model type should be first
  • The order of the options do not matter: bt-U-ordereffect is equivalent to bt-ordereffect-U
  • The - is mandatory

Roadmap

Goals for bpcs 1.1.0 (Before June 2021)

  • Add new models for modeling time effects (newer contests have higher impact on the ability than older contests)
  • Add model for the Bayesian ELO-type rating system
  • Add models for cumulative comparisons
  • Improve test coverage
  • Add ties to expand_aggregated_data.

Vignettes

This package provides a series of small and self contained vignettes that exemplify the use of each model. In the vignettes, we also provide examples of code for data transformation, tables and plots.

Below we list all our vignettes with a short description:

  • Getting Started: This vignette shows a basic example on tennis competition data, covering how to run a Bradley-Terry model, MCMC diagnostics, posterior predictive values, ranking, predict new matches

  • Ties and home advantage: This vignette covers a soccer example from the Brazilian soccer league. Here, we first model the results using a Bradley-Terry model and the Davidson model to handle ties. Then, we extend both models to include for order effects, this allows us to investigate the home advantage in and without the presence of ties.

  • Bradley-Terry with random effects: This vignette covers the problem of ranking black-box optimization algorithms based on benchmarks. Since in benchmarking we often run the same optimization algorithm more than once with the same benchmark problem, we need to compensate for the repeated measures effect. We deal with this utilizing a simple Bradley-Terry model with random effects.

Contributing and bugs

If you are interested you are welcome to contribute to the repository through pull requests.

We have a short contributing guide vignette.

If you find bugs, please report it in https://github.com/davidissamattos/bpcs/issues

Icon credits

  • Boxing gloves image by “surang” from “flaticons.com”
  • Hex Sticker created with the hexSticker package

References

Agresti, Alan. 2003. Categorical Data Analysis. Vol. 482. John Wiley & Sons.

Böckenholt, Ulf. 2001. “Hierarchical Modeling of Paired Comparison Data.” Psychological Methods 6 (1): 49.

Bradley, Ralph Allan, and Milton E Terry. 1952. “Rank Analysis of Incomplete Block Designs: I. The Method of Paired Comparisons.” Biometrika 39 (3/4): 324–45.

Davidson, Roger R. 1970. “On Extending the Bradley-Terry Model to Accommodate Ties in Paired Comparison Experiments.” Journal of the American Statistical Association 65 (329): 317–28.

Davidson, Roger R, and Robert J Beaver. 1977. “On Extending the Bradley-Terry Model to Incorporate Within-Pair Order Effects.” Biometrics, 693–702.

Springall, A. 1973. “Response Surface Fitting Using a Generalization of the Bradley-Terry Paired Comparison Model.” Journal of the Royal Statistical Society: Series C (Applied Statistics) 22 (1): 59–68.

Stan Development Team. 2020. “RStan: The R Interface to Stan.” https://mc-stan.org/.

Copy Link

Version

Install

install.packages('bpcs')

Monthly Downloads

8

Version

1.0.0

License

MIT + file LICENSE

Issues

4

Pull Requests

0

Stars

12

Forks

0

Maintainer

David Issa Mattos

Last Published

December 9th, 2020

Functions in bpcs (1.0.0)

HPDI_from_stanfit

Calculate HPDI for all parameters from a stanfit object Here we use the coda package
HPD_higher_from_column

Returns the higher value of the HPD interval for a data frame column
check_if_there_are_ties

Check if a data frame column contains ties
bpc

Bayesian Paired comparison regression models in Stan
check_numeric_predictor_matrix

Check if all values in the predictor matrix are numeric and not NA. Note that TRUE will be cast to 1 and FALSE will be cast to 0
check_predictors_df_contains_all_players

Check if the predictor df contains all players and only those
HPD_lower_from_column

Returns the lower value of the HPD interval for a data frame column
check_if_there_are_na

Check for NA in the specfic columns and returns T or F is there is at least 1 NA in those columns
bpcs-package

bpcs - A package for Bayesian Paired Comparison analysis with Stan
brasil_soccer_league

This is a dataset with the results matches fromo the first league of the Brazilian soccer championship from 2017-2019. It was reduced and translatedfrom the adaduque/Brasileirao_Dataset repository
create_index_with_existing_lookup_table

Create two columns with the indexes for the names Here we use an existing lookup table. Should be used in predicting
get_rank_of_players

Generate a ranking of the ability based on sampling the posterior distribution of the ranks.
create_predictor_matrix_with_player_lookup_table

Receives a predictor dataframe, a string with the column of the player, a vector of strings with the columns for the predictors and a lookup table and returns an ordered matrix for Stan To be used with the predictors data frame
get_sample_posterior

Get the posterior samples for a parameter of the model.
create_index_lookuptable

Create a lookup table of names and indexes Note that the indexes will be created in the order they appear. For string this doesnt make much difference but for numbers the index might be different than the actual number that appears in names
create_bpc_object

Defines the class bpc and creates the bpc object. To create we need to receive some defined parameters (the arguments from the bpc function), a lookup table and a the stanfit object generated from the rstan sampling procedure
sample_stanfit

Return a data frame by resampling the posterior from a stanfit Here we select a parameter, retrieve the all the posterior from the stanfit and then we resample this posterior n times
replace_parameter_index_with_names

Replace the name of the parameter from index to name using a lookup_table Receives a data frame and returns a dataframe.
create_array_of_par_names

Create an array with the parameter name and to what player/cluster it refers to in the order stan presents
create_index_predictors_with_lookup_table

Receives one column with player names and returns a data frame with the relevant index columns based on a given lookup table To be used with the predictors data frame
predict.bpc

Predict results for new data.
compute_scores

Giving a player0 an player1 scores, this functions adds one column to the data frame containing who won (0= player0 1=player1 2=tie) and another if it was a tie. The ties column superseeds the y column. If it was tie the y column does not matter y column: (0= player0 1=player1 2=tie) ties column (0=not tie, 1=tie)
print.bpc

Print method for the bpc object.
check_result_column

Check if a data frame column contains only the values 1 0 and 2. Used to check the format of the results
compute_ties

Giving a result column we create a new column with ties (0 and 1 if it has)
optimization_algorithms

Dataset containing an example of the performance of different optimization algorithms against different benchmark functions. This is a reduced version of the dataset presented at the paper: "Statistical Models for the Analysis of Optimization Algorithms with Benchmark Functions.". For details on how the data was collected we refer to the paper.
get_model_parameters

Return all the name of parameters in a model from a bpc_object. Here we exclude the log_lik and the lp__ since they are not parameters of the model
get_stanfit

Retrieve the stanfit object generated by rstan.
check_z_column

Check if a data frame column contains only the values 1 or 0. For the z column
get_stanfit_summary

Get stanfit summary table of all parameters excluding log_lik.
get_probabilities

Get the empirical win/draw probabilities based on the ability/strength parameters. Instead of calculating from the probability formula given from the model we create a predictive posterior distribution for all pair combinations and calculate the posterior wins/loose/draw The function returns the mean value of win/loose/draw for the player i. To calculate for player j the probability is 1-p_i
create_predictors_lookup_table

Receives a vector with predictors strings (the column names) and returns a predictor_lookup_table
create_index

Create two columns with the indexes for the names of the players Here we create a new lookup table. Should be used when sampling the parameters
create_index_cluster_lookuptable

Create a lookup table of names and indexes Note that the indexes will be created in the order they appear. For string this does not make much difference but for numbers the index might be different than the actual number that appears in names
create_cluster_index

Create two columns with the indexes for the names of the players Here we create a new lookup table. Should be used when sampling the parameters
expand_aggregated_data

Expand aggregated data Several datasets for the Bradley-Terry Model aggregate the number of wins for each player in a different column. The models we provide are intended to be used in a long format. A single result for each contest. This function expands datasets that have aggregated data into this long format.
logit

Logit function
launch_shinystan

Tiny wrapper to launch a shinystan app to investigate the MCMC.
get_hpdi_parameters

Return the mean and the HPDI of the parameters of the model
get_loo

Tiny wrapper for the PSIS-LOO-CV method from the loo package.
get_waic

Tiny wrapper for the WAIC method from the loo package.
inv_logit

Inverse logit function
match_cluster_names_to_cluster_lookup_table

Receives a column with cluster names and returns a data frame with the relevant index column based on a given cluster lookup table
create_cluster_index_with_existing_lookup_table

Create two columns with the indexes for the names Here we use an existing lookup table. Should be used in predicting
match_player_names_to_lookup_table

Receives two columns with player names and returns a data frame with the relevant index columns based on a given lookup table
summary.bpc

Summary of the model bpc model.
tennis_agresti

This is the expansion of the tennis data from Agresti (2003) p.449 This data refers to matches for several women tennis players during 1989 and 1990
%>%

Pipe operator