class: center, middle, inverse, title-slide # <a href="https://github.com/EllaKaye/BradleyTerryScalable">BradleyTerryScalable</a> ## Ranking items scalably with the Bradley-Terry model ### <div class="line-block">Ella Kaye, with David Firth<br /> University of Warwick</div> ### July 6th, 2017 --- class: inverse, center, middle # The Bradley-Terry model --- background-image: url("football.png") background-size: cover --- background-image: url(http://authorservices.taylorandfrancis.com/custom/uploads/2015/10/Maths-Stats-Journal-Finder-1170x827.jpg) background-size: cover --- - Let there be `\(K\)` items. The Bradley-Terry probability that item `\(i\)` beats item `\(j\)` is `$$p_{ij} = \frac{\pi_i}{\pi_i + \pi_j},$$` where `\(\pi_k\)` is a strength parameter for player `\(k\)`, `\(1 \leq k \leq K\)`. - Equivalently `$$\text{logit}(p_{ij}) = \lambda_i - \lambda_j,$$` where `\(\lambda_k = \log(\pi_k)\)`. - A constraint is needed: - `\(\sum_{i=1}^K{\pi_i} = K\)`, - `\(\text{mean}({\lambda_i}) = 0\)`. - The log-likelihood is given by `$$\ell(\pi) = \sum_{i = 1}^{K}\sum_{j = 1}^{K} [w_{ij}(\log \pi_i - \log(\pi_i + \pi_j))],$$` where `\(w_{ij}\)` is the number of times item `\(i\)` beats item `\(j\)`. --- # Comparison graph/Existence of the MLE <img src="useR_presentation_files/figure-html/unnamed-chunk-1-1.png" style="display: block; margin: auto;" /> The MLE exists and is finite whenever the comparison graph is fully connected (there is a directed path from `\(i\)` to `\(j\)` for every `\(i\)` and every `\(j\)`). --- # Options for fitting the Bradley-Terry model - If the comparison graph is fully connected, can find the MLE (using MM-algorithm of Hunter (2004)) - If the comparison graph is not fully connected: - find MLE on each fully-connected component - find MAP estimate using Bayesian approach of Caron and Doucet (2012) - this exists and is finite, and allows ranking of all items, even when the comparison graph is not fully connected. - N.B. Can also include a prior when the comparison graph is fully connected --- # Fitting the model - Using the MM-algorithm to find the MLE: `$$\pi_i^{(n+1)} = \frac{W_i}{\sum_{j=1}^K \frac{n_{ij}}{\pi_i^{(n)} + \pi_j^{(n)}}}$$` where `\(W_i = \sum_{j=1}^K w_{ij}\)` is the number of wins for item `\(i\)` and `\(n_{ij} = w_{ij} + w_{ji}\)` is the number of comparisons between item `\(i\)` and item `\(j\)`. - Using the EM-algorithm to find the MAP estimate: `$$\pi_i^{(n+1)} = \frac{a - 1 + W_i}{b + \sum_{j=1}^K \frac{n_{ij}}{\pi_i^{(n)} + \pi_j^{(n)}}},$$` where `\(a\)` and `\(b\)` are the shape and rate parameters of a gamma-distributed prior on `\(\pi\)`: `\(p(\pi) = \prod_{i=1}^K \mathcal{G}(\pi_i; a, b)\)`. --- class: inverse, center, middle # `BradleyTerryScalable` --- # Aims - Fit the Bradley-Terry model to large and (potentially) sparse data sets - Has to be fast *enough* ('Scalable' is in the package name!) - Has to be able to handle large data sets - Has to be able to deal with cases when the comparison graph is not fully connected - Do one thing well, within a complete package - Well thought-out user interface - Complement to (NOT a replacement for) `BradleyTerry2` package --- # Installation ```r # CRAN install.packages("BradleyTerryScalable") # GitHub devtools::install_github("EllaKaye/BradleyTerryScalable") ``` ```r library(BradleyTerryScalable) ``` --- # Workflow - `btdata(x)` to create object of class `"btdata"` - may need to call `codes_to_counts()` first - `summary(btdata)` - `select_components(btdata, subset)` - `btfit(btdata, a)` to fit model and create object of class `"btfit"` - S3 methods for `btfit` object: - `summary` - `fitted` - `coef` - `vcov` - `simulate` - `btprob(object)` for Bradley-Terry probabilities `\(p_{ij}\)` --- # `btdata(x)` - `x` can be a three or four column data frame with counts of wins, an `igraph` object or a square matrix or contigency table ```r data(citations) citations ``` ``` ## citing ## cited Biometrika Comm Statist JASA JRSS-B ## Biometrika 714 730 498 221 ## Comm Statist 33 425 68 17 ## JASA 320 813 1072 142 ## JRSS-B 284 276 325 188 ``` ```r citations_data <- btdata(citations) summary(citations_data) ``` ``` ## Number of items: 4 ## Density of wins matrix: 1 ## Fully-connected: TRUE ``` --- ```r data(toy_data) head(toy_data, 4) ``` ``` ## # A tibble: 4 x 3 ## player1 player2 outcome ## <chr> <chr> <chr> ## 1 Cyd Amy W1 ## 2 Amy Ben D ## 3 Ben Eve W2 ## 4 Cyd Dan W2 ``` ```r toy_data_4c <- codes_to_counts(toy_data, c("W1", "W2", "D")) toy_btdata <- btdata(toy_data_4c) summary(toy_btdata) ``` ``` ## Number of items: 8 ## Density of wins matrix: 0.25 ## Fully-connected: FALSE ## Number of fully-connected components: 3 ## Summary of fully-connected components: ## Component size Freq ## 1 1 1 ## 2 3 1 ## 3 4 1 ``` --- # `btfit(btdata, a)` ```r btfit(btdata, a) ``` - If `a = 1`, finds MLE on each fully connected component, using the MM-algorithm - If `a > 1`, finds the MAP estimate of `\(\pi\)` - Returns a `btfit` S3 object ```r btfit(btdata, a, MAP_by_component = FALSE, subset = NULL, maxit = 10000, epsilon = 0.001) ``` --- ```r toy_fit_MAP <- btfit(toy_btdata, a = 1.1) summary(toy_fit_MAP) ``` ``` ## Warning: `data_frame()` is deprecated as of tibble 1.1.0. ## Please use `tibble()` instead. ## This warning is displayed once every 8 hours. ## Call `lifecycle::last_warnings()` to see where this warning was generated. ``` ``` ## $call ## btfit(btdata = toy_btdata, a = 1.1) ## ## $item_summary ## # A tibble: 8 x 3 ## component item estimate ## <chr> <chr> <dbl> ## 1 full_dataset Eve 1.90 ## 2 full_dataset Cyd 0.472 ## 3 full_dataset Han 0.245 ## 4 full_dataset Amy -0.0766 ## 5 full_dataset Gal -0.102 ## 6 full_dataset Ben -0.423 ## 7 full_dataset Dan -0.536 ## 8 full_dataset Fin -1.48 ## ## $component_summary ## # A tibble: 1 x 4 ## component num_items iters converged ## <chr> <int> <int> <lgl> ## 1 full_dataset 8 101 TRUE ``` --- ```r toy_fit_MLE <- btfit(toy_btdata, a = 1) summary(toy_fit_MLE, SE = TRUE) ``` ``` ## $call ## btfit(btdata = toy_btdata, a = 1) ## ## $item_summary ## # A tibble: 7 x 4 ## component item estimate SE ## <chr> <chr> <dbl> <dbl> ## 1 2 Han 0.696 0.911 ## 2 2 Gal 0.413 0.768 ## 3 2 Fin -1.11 1.05 ## 4 3 Cyd 0.592 0.991 ## 5 3 Amy 0.0325 0.699 ## 6 3 Ben -0.243 0.944 ## 7 3 Dan -0.382 0.712 ## ## $component_summary ## # A tibble: 2 x 4 ## component num_items iters converged ## <chr> <int> <int> <lgl> ## 1 2 3 6 TRUE ## 2 3 4 10 TRUE ``` --- ```r toy_data %>% codes_to_counts(c("W1", "W2", "D")) %>% btdata() %>% btfit(1) %>% summary(ref = "Amy", subset = function(x) "Amy" %in% names(x)) ``` ``` ## $call ## btfit(btdata = ., a = 1) ## ## $item_summary ## # A tibble: 4 x 3 ## component item estimate ## <chr> <chr> <dbl> ## 1 3 Cyd 0.560 ## 2 3 Amy 0 ## 3 3 Ben -0.276 ## 4 3 Dan -0.414 ## ## $component_summary ## # A tibble: 1 x 4 ## component num_items iters converged ## <chr> <int> <int> <lgl> ## 1 3 4 10 TRUE ``` --- # `btprob(object)` Gives the Bradley-Terry probabilities `\(\frac{\pi_i}{\pi_i + \pi_j}\)` ```r citations %>% btdata() %>% btfit(1, epsilon = 1e-6) %>% btprob() ``` ``` ## 4 x 4 sparse Matrix of class "dgCMatrix" ## citing ## cited JRSS-B Biometrika JASA Comm Statist ## JRSS-B . 0.56683657 0.67885745 0.9615071 ## Biometrika 0.43316343 . 0.61764636 0.9502196 ## JASA 0.32114255 0.38235364 . 0.9219760 ## Comm Statist 0.03849288 0.04978036 0.07802401 . ``` --- ```r citations %>% btdata() %>% btfit(1, epsilon = 1e-6) %>% btprob(as_df = TRUE) ``` ``` ## # A tibble: 6 x 5 ## component cited citing prob1wins prob2wins ## <chr> <chr> <chr> <dbl> <dbl> ## 1 full_dataset JRSS-B Biometrika 0.567 0.433 ## 2 full_dataset JRSS-B JASA 0.679 0.321 ## 3 full_dataset Biometrika JASA 0.618 0.382 ## 4 full_dataset JRSS-B Comm Statist 0.962 0.0385 ## 5 full_dataset Biometrika Comm Statist 0.950 0.0498 ## 6 full_dataset JASA Comm Statist 0.922 0.0780 ``` --- # A proper example ```r comp_cites ``` ``` ## # A tibble: 415,496 x 3 ## cited_pdpass citing_pdpass count ## <chr> <chr> <dbl> ## 1 10032645 10030580 2 ## 2 10083166 10030580 4 ## 3 10170329 10030580 2 ## 4 10191810 10030580 2 ## 5 10223803 10030580 2 ## 6 10254226 10030580 1 ## 7 11490383 10030580 2 ## 8 11720482 10030580 2 ## 9 11811663 10030580 1 ## 10 11862620 10030580 4 ## # … with 415,486 more rows ``` --- ```r system.time(comp_cites_data <- btdata(comp_cites)) ``` ``` ## user system elapsed ## 2.794 0.147 2.947 ``` ```r summary(comp_cites_data) ``` ``` ## Number of items: 27137 ## Density of wins matrix: 0.0005642131 ## Fully-connected: FALSE ## Number of fully-connected components: 11297 ## Summary of fully-connected components: ## Component size Freq ## 1 1 11285 ## 2 2 10 ## 3 3 1 ## 4 15829 1 ``` --- ```r comp_cites_data <- btdata(comp_cites) system.time(comp_cites_fit <- btfit(comp_cites_data, 1)) ``` ``` ## user system elapsed ## 7.212 0.644 7.998 ``` ```r system.time(comp_cites_fit <- btfit(comp_cites_data, 1.1)) ``` ``` ## user system elapsed ## 21.500 2.328 24.159 ``` --- <img src="useR_presentation_files/figure-html/unnamed-chunk-19-1.png" width="80%" height="80%" style="display: block; margin: auto;" /> --- class: center, middle # Further information CRAN [https://CRAN.R-project.org/package=BradleyTerryScalable](https://CRAN.R-project.org/package=BradleyTerryScalable) GitHub [https://github.com/EllaKaye/BradleyTerryScalable]([https://github.com/EllaKaye/BradleyTerryScalable]) Documentation (`pkgdown` site) [https://ellakaye.github.io/BradleyTerryScalable]([https://ellakaye.github.io/BradleyTerryScalable]) --- class: inverse, center, middle # Thank you! ## Any questions? ## I'd love to hear from you! [E.Kaye.1@warwick.ac.uk](mailto:E.Kaye.1@warwick.ac.uk) [@ellamkaye](https://twitter.com/ellamkaye) [ellakaye.rbind.io](https://ellakaye.rbind.io) [github.com/EllaKaye](https://github.com/ellakaye)