| Title: | Efficient Bayesian Inference for the Bradley--Terry Model |
|---|---|
| Description: | A suite of functions that allow a full, fast, and efficient Bayesian treatment of the Bradley--Terry model. Prior assumptions about the model parameters can be encoded through a multivariate normal prior distribution. Inference is performed using a latent variable representation of the model. |
| Authors: | Rowland Seymour [aut, cre, cph] (ORCID: <https://orcid.org/0000-0002-8739-3921>), Christine Stawitz [aut] (ORCID: <https://orcid.org/0000-0003-3122-4501>), Heather Turner [aut] (ORCID: <https://orcid.org/0000-0002-1256-3375>) |
| Maintainer: | Rowland Seymour <[email protected]> |
| License: | GPL (>= 3) |
| Version: | 1.0.0.9000 |
| Built: | 2026-09-23 17:34:14 UTC |
| Source: | https://github.com/rowlandseymour/speedybbt |
This function fits the Bradley-Terry model with comparison and player specific effects. Each comparison can be assigned a real value to allow for a specific effect for the comparison, such as bias, ordering or home/away effect. The value of this effect is denoted kappa. The player specific effects are described through a formula and data.frame containing the value. The function places a normal prior distribution on both kappa and the player specific parameters beta.
BBTm( outcome, player1, player2, lambda.initial = NULL, player.prior.var = NULL, beta.initial = NULL, n.iter = 1000, formula = NULL, data = NULL, advantage = NULL, kappa.initial = NULL, kappa.var = NULL, hyperparameter = TRUE, chi = 0.01, psi = 0.01, burn.in = 100, n.thin = 1, verbose = interactive() )BBTm( outcome, player1, player2, lambda.initial = NULL, player.prior.var = NULL, beta.initial = NULL, n.iter = 1000, formula = NULL, data = NULL, advantage = NULL, kappa.initial = NULL, kappa.var = NULL, hyperparameter = TRUE, chi = 0.01, psi = 0.01, burn.in = 100, n.thin = 1, verbose = interactive() )
outcome |
Vector of outcomes. 1 if player2 is the winner, 0 if player1 is the winner. |
player1 |
Vector of first players. |
player2 |
Vector of second players. |
lambda.initial |
(optional) Vector containing the values of the player parameters for the first MCMC iteration. |
player.prior.var |
(optional) Matrix specifying the prior covariance of the player correlation parameters. |
beta.initial |
(optional) Vector containing the values of the player specific parameters for the first MCMC iteration. |
n.iter |
Number of MCMC samples to be drawn. |
formula |
Formula with no left-hand-side specifying the player specific effects. |
data |
Data frame with a row corresponding to each player and a column corresponding to each covariate. |
advantage |
(optional) A vector with the value of the comparisons specific effect for each comparison. |
kappa.initial |
(optional) An initial value for the comparison specific value kappa. |
kappa.var |
(optional) The prior variance of the comparison specific value kappa. |
hyperparameter |
Boolean indicating if inference should be performed for the prior variance hyperparameter. If TRUE the prior variance (main diagonal of the covariance matrix) must be set to 1. |
chi |
Rate parameter for the inverse-gamma prior distribution on the hyperparameter. |
psi |
Shape parameter for the inverse-gamma prior distribution on the hyperparameter. |
burn.in |
The number of iterations to use for a burn.in, default is 100. |
n.thin |
(optional) The number of iterations to thin the MCMC samples by. Default is 1. |
verbose |
(optional) A boolean indicating if progress should be printed to the console.
Default is |
If player.prior.var is omitted, independent and identical
N(0, 5^2) prior distributions are placed on each object quality parameter.
If beta.initialis omitted, it is set to a vector of zeroes.
If kappa.var is omitted, it is set to N(0, 5^2), if kappa.initial is omitted
it is set to 0.5.
A "mcmc" object containing samples from the posterior distribution.
##################### ## Wimbledon 2019 ## #################### # Fit model where the quality of each player depends on their rank # and the number of points they had immediately before the tournament. # Allow an effect for a match being in the first or second week. wimbledonModel <- BBTm( outcome = wimbledon$matches$outcome, player2 = wimbledon$matches$loser, player1 = wimbledon$matches$winner, advantage = wimbledon$matches$secondWeek, formula = ~ rank + points, data = wimbledon$players, n.iter = 1000, verbose = FALSE ) # Plot posterior distributions hist(parameter(wimbledonModel, "kappa"), main = "", xlab = expression(kappa), freq = FALSE )##################### ## Wimbledon 2019 ## #################### # Fit model where the quality of each player depends on their rank # and the number of points they had immediately before the tournament. # Allow an effect for a match being in the first or second week. wimbledonModel <- BBTm( outcome = wimbledon$matches$outcome, player2 = wimbledon$matches$loser, player1 = wimbledon$matches$winner, advantage = wimbledon$matches$secondWeek, formula = ~ rank + points, data = wimbledon$players, n.iter = 1000, verbose = FALSE ) # Plot posterior distributions hist(parameter(wimbledonModel, "kappa"), main = "", xlab = expression(kappa), freq = FALSE )
This function uses MCMC to sample from the posterior distribution of the Bradley–Terry model with ties.A multivariate normal prior distribution on the player quality parameters can be specified. An exponential prior distribution is placed on the tie parameter theta, and a Metropolis- Hasting random walk algorithm is used to update this parameter.
BBTm.ties( n.objects, outcome, player1, player2, player.prior.var = NULL, theta.initial = NULL, lambda.initial = NULL, n.iter = 1000, hyperparameter = TRUE, chi = 0.01, psi = 0.01, rw.sd = 0.1, theta.rate = 0.01, burn.in = 100, n.thin = 1, verbose = interactive() )BBTm.ties( n.objects, outcome, player1, player2, player.prior.var = NULL, theta.initial = NULL, lambda.initial = NULL, n.iter = 1000, hyperparameter = TRUE, chi = 0.01, psi = 0.01, rw.sd = 0.1, theta.rate = 0.01, burn.in = 100, n.thin = 1, verbose = interactive() )
n.objects |
The number of objects in the study. |
outcome |
Vector of outcomes. 0 if player 1 is the winner, 1 if player 2 is the winner, and 2 if it is a tie. |
player1 |
Vector of first players. |
player2 |
Vector of second players. |
player.prior.var |
(optional) Matrix specifying the prior covariance of the player correlation parameters. |
theta.initial |
(optional) Value of the tied parameter there for the first MCMC iteration. |
lambda.initial |
(optional) Vector containing the values of the player parameters for the first MCMC iteration. |
n.iter |
Number of MCMC samples to be drawn. |
hyperparameter |
Boolean indicating if inference should be performed for the prior variance hyperparameter. If TRUE the prior variance (main diagonal of the covariance matrix) must be set to 1. |
chi |
(optional) Rate parameter for the inverse-gamma prior distribution on the hyperparameter. |
psi |
(optional) Shape parameter for the inverse-gamma prior distribution on the hyperparameter. |
rw.sd |
(optional) Number describing the standard deviation of normal distribution proposal distribution for theta. |
theta.rate |
(optional) The rate parameter of the exponential prior distribution placed on theta. |
burn.in |
(optional) The number of iterations to use as a burn-in period. Default is 100. |
n.thin |
(optional) The number of iterations to thin the MCMC samples by. Default is 1. |
verbose |
(optional) A boolean indicating if progress should be printed to the console.
Default is |
If player.prior.var is omitted, independent and identical
N(0, 5^2) prior distributions are placed on each object quality parameter.
If lambda.initial is omitted, it is set to a vector of zeroes.
A "mcmc" object containing samples from the posterior distribution
############################################ ## Deprivation in Dar es Salaam, Tanzania ## ## Seymour et al (2022) ## ############################################ # Construct covariance matrix based on spatial informartion prior.var <- expm::expm(darEsSalaam$adjacencyMatrix) prior.var <- diag(diag(prior.var)^-0.5) %*% prior.var %*% diag(diag(prior.var)^-0.5) # Fit BT model with ties # Using `n.iter = 3` here to reduce model runtime, you will need # a larger number of iterations for valid inference. darTiedModel <- BBTm.ties( n.objects = nrow(darEsSalaam$adjacencyMatrix), outcome = darEsSalaam$comparisons$outcome, player1 = darEsSalaam$comparisons$subward1, player2 = darEsSalaam$comparisons$subward2, player.prior.var = prior.var, hyperparameter = TRUE, rw.sd = 0.005, n.iter = 3, burn.in = 0, verbose = FALSE ) # Get posterior means darTiedModelLambda <- parameter(darTiedModel, "lambda") - colMeans(parameter(darTiedModel, "lambda")) lambda.mean <- rowMeans(darTiedModelLambda) # Generate trace plots plot(lambda.mean) plot(parameter(darTiedModel, "theta"), type = "l")############################################ ## Deprivation in Dar es Salaam, Tanzania ## ## Seymour et al (2022) ## ############################################ # Construct covariance matrix based on spatial informartion prior.var <- expm::expm(darEsSalaam$adjacencyMatrix) prior.var <- diag(diag(prior.var)^-0.5) %*% prior.var %*% diag(diag(prior.var)^-0.5) # Fit BT model with ties # Using `n.iter = 3` here to reduce model runtime, you will need # a larger number of iterations for valid inference. darTiedModel <- BBTm.ties( n.objects = nrow(darEsSalaam$adjacencyMatrix), outcome = darEsSalaam$comparisons$outcome, player1 = darEsSalaam$comparisons$subward1, player2 = darEsSalaam$comparisons$subward2, player.prior.var = prior.var, hyperparameter = TRUE, rw.sd = 0.005, n.iter = 3, burn.in = 0, verbose = FALSE ) # Get posterior means darTiedModelLambda <- parameter(darTiedModel, "lambda") - colMeans(parameter(darTiedModel, "lambda")) lambda.mean <- rowMeans(darTiedModelLambda) # Generate trace plots plot(lambda.mean) plot(parameter(darTiedModel, "theta"), type = "l")
This function constructs a win matrix from a data frame of comparisons. It is needed for the MCMC functions.
comparisons_to_matrix(n.objects, comparisons)comparisons_to_matrix(n.objects, comparisons)
n.objects |
The number of areas in the study. |
comparisons |
An N x 2 data frame, where N is the number of comparisons. Each row should correspond to a judgment. The first column is the winning object, the second column is the more losing object. The areas should be labeled from 1 to n.objects. |
A matrix where the i, j th element is the number of times object i beat object j.
# Generate some sample comparisons comparisons <- data.frame("winner" = c(1, 3, 2, 2), "loser" = c(3, 1, 1, 3)) # Create matrix from comparisons win.matrix <- comparisons_to_matrix(3, comparisons)# Generate some sample comparisons comparisons <- data.frame("winner" = c(1, 3, 2, 2), "loser" = c(3, 1, 1, 3)) # Create matrix from comparisons win.matrix <- comparisons_to_matrix(3, comparisons)
A comparative judgment data set on deprivation in subwards in Dar es Salaam, Tanzania.Citizens were shown pairs of subwards at random and asked which was more deprived.If they said they were equal, one of the pair was chosen at random to be more deprived.The data was collected in August 2018. The sex of each judge is also included.
darEsSalaamdarEsSalaam
A list with three elements.
A dataframe containing the comparisons. Each row corresponds
to a judgement made by a single judge. This differs from the data
in the BBT package as it explicitly includes ties rather than randomly
allocating a winner. The columns are:
The outcome of the comparison. 1 if item 2 is the winner, 0 if item 1 is the winner, and 2 if it was a tie (although there are no instances of this happening).
The first item in the comparison.
The second item in the comparison.
The sex of the judge making the comparison.
A dataframe containing the names and shapefiles of the subwards with columns:
The name of the subward.
A list containing the shapefile of the subward.
An adjacency matrix of the subwards formed from the shapefiles.This considers subwards as nodes and places edges between adjacent subwards. Two additional edges have been manually included to allow for crossings of the Kurasini creek.
This data set was collected by Madeleine Ellis, James Goulding, Bertrand Perrat, Gavin Smith and Gregor Engelmann. We gratefully acknowledge the Rights Lab at the University of Nottingham for supporting funding for the comprehensive ground truth survey. We also acknowledge HumanitarianStreet Mapping Team (HOT) for providing a team of experts in data collection to facilitate the surveys. This work was also supported by the EPSRC Horizon Centre for Doctoral Training - My Life in Data (EP/L015463/1) and EPSRC grant Neodemographics (EP/L021080/1).
A comparative judgment data set for risk of forced marriage at ward level in Nottinghamshire. There are 12 judges and 76 wards.
forcedMarriageforcedMarriage
A list with three elements.
A dataframe containing the comparisons. Each row corresponds to a judgement made by a single judge. The columns are:
The ID of the judge.
The time the comparison was made.
The ID of the ward that was judged to have a higher risk of forced marriage.
The ID of the ward that was judged to have a lower risk of forced marriage.
A dataframe containing information about each of the wards. The columns are:
The name of the ward.
The code of the ward.
The electoral division of the ward.
The name of the file containing the shapefile for the ward.
The final element is an adjacency matrix, where the wards are nodes and edges are placed between adjacent wards.
The data was collected using support from the Engineering and Physical Sciences Research Council (grant reference EP/R513283/1), the Economic and Social Sciences Research Council (ES/V015370/1) and the Research England Policy Support Fund. The data was collected following ethical approval from the University of Nottingham School of Politics and International Relations ethics committee.
A comparative judgment data set for risk of honour based abuse in Oxford and Banbury
oxon.comparisonsoxon.comparisons
A data frame with 1,167 comparisons.
The ID of the comparison.
The ID of the user who made the comparisons.
The ID of the first area involved in the comparison.
The ID of the second area involved in the comparison.
The ID of the selected area. If the
comparison was tied, the selected_item_id is NA.
The state of the outcome. selected indicates that a judgment was made, skipped indicates that the comparison was skipped, and tied indicates that the comparison was tied.
The data was collected following ethical approval the University of Birmingham's Science, Engineering and Maths Ethics Committee.
This function helps the user extract relevant estimates from parameters fit using a Bayesian Bradley-Terry model.
parameter(model_output, parameter_name, indices_to_extract = NULL)parameter(model_output, parameter_name, indices_to_extract = NULL)
model_output |
The mcmc object that contains the draws output from the BBT model. |
parameter_name |
The name of the parameter to extract draws for. |
indices_to_extract |
A single column index or a vector of column indices to extract from the draws matrix. |
A vector containing the draws of the parameter.
prior.var <- expm::expm(forcedMarriage$adjacencyMatrix) prior.var <- diag(diag(prior.var)^-0.5) %*% prior.var %*% diag(diag(prior.var)^-0.5) # Fit model # Using `n.iter = 3, burn.in = 0` to reduce runtime # - you should use more iterations for reliable posterior estimates forcedMarriageModel <- speedyBBTm( outcome = rep(1, length(forcedMarriage$comparisons$win)), player1 = forcedMarriage$comparisons$win, player2 = forcedMarriage$comparisons$lost, player.prior.var = prior.var, n.iter = 3, burn.in = 0 ) parameter(forcedMarriageModel, "lambda", c(10, 20, 30, 40)) parameter(forcedMarriageModel, "alpha.sq")prior.var <- expm::expm(forcedMarriage$adjacencyMatrix) prior.var <- diag(diag(prior.var)^-0.5) %*% prior.var %*% diag(diag(prior.var)^-0.5) # Fit model # Using `n.iter = 3, burn.in = 0` to reduce runtime # - you should use more iterations for reliable posterior estimates forcedMarriageModel <- speedyBBTm( outcome = rep(1, length(forcedMarriage$comparisons$win)), player1 = forcedMarriage$comparisons$win, player2 = forcedMarriage$comparisons$lost, player.prior.var = prior.var, n.iter = 3, burn.in = 0 ) parameter(forcedMarriageModel, "lambda", c(10, 20, 30, 40)) parameter(forcedMarriageModel, "alpha.sq")
Generates a plot of the posterior mean and 95% credible interval of the quality parameter estimates.
plot_qualities( model_output, player_names, flip = FALSE, main = "Posterior Mean and 95% Credible \n Interval of Player Quality", xlab = "Player", ylab = "Posterior Mean Quality", quality_label = NULL, ... )plot_qualities( model_output, player_names, flip = FALSE, main = "Posterior Mean and 95% Credible \n Interval of Player Quality", xlab = "Player", ylab = "Posterior Mean Quality", quality_label = NULL, ... )
model_output |
The "mcmc" model output object containing the draws. |
player_names |
A character vector containing the names of the items (e.g., players or wards). |
flip |
Whether to flip the plot to put the quality estimates on the x-axis instead of the y-axis; defaults to FALSE. |
main |
The text to use for the plot title. |
xlab |
The text to use for the x-axis label. |
ylab |
The text to use for the y-axis label. |
quality_label |
A vector of three custom labels to use instead of numbers on the quality axis. |
... |
other arguments passed to |
players <- wimbledon$players$name wimbledonModel <- BBTm( outcome = wimbledon$matches$outcome, player2 = wimbledon$matches$loser, player1 = wimbledon$matches$winner, advantage = wimbledon$matches$secondWeek, formula = ~ rank + points, data = wimbledon$players, n.iter = 200 ) plot_qualities(model_output = wimbledonModel, player_names = players)players <- wimbledon$players$name wimbledonModel <- BBTm( outcome = wimbledon$matches$outcome, player2 = wimbledon$matches$loser, player1 = wimbledon$matches$winner, advantage = wimbledon$matches$secondWeek, formula = ~ rank + points, data = wimbledon$players, n.iter = 200 ) plot_qualities(model_output = wimbledonModel, player_names = players)
This function uses MCMC to sample from the posterior distribution of the standard Bradley–Terry model. Standard model means that there are no tied comparisons and no player or comparison specific variables. This provides a fast implementation of the standard model. A multivariate normal prior distribution on the player quality parameters can be specified.
speedyBBTm( outcome = NULL, player1 = NULL, player2 = NULL, win.matrix = NULL, player.prior.var = NULL, lambda.initial = NULL, n.iter = 1000, hyperparameter = TRUE, chi = 0.01, psi = 0.01, burn.in = 100, n.thin = 1, verbose = interactive() )speedyBBTm( outcome = NULL, player1 = NULL, player2 = NULL, win.matrix = NULL, player.prior.var = NULL, lambda.initial = NULL, n.iter = 1000, hyperparameter = TRUE, chi = 0.01, psi = 0.01, burn.in = 100, n.thin = 1, verbose = interactive() )
outcome |
Vector of outcomes. 1 if player 2 is the winner, 0 if player 1 is the winner. |
player1 |
Vector of first players. |
player2 |
Vector of second players. |
win.matrix |
(optional) A win-loss matrix where the i,j th element is the number of times object i beat object j. |
player.prior.var |
(optional) A matrix specifying the prior covariance of the player correlation parameters. |
lambda.initial |
(optional) A vector containing the values of the player correlation parameters for the first MCMC iteration. |
n.iter |
The number of MCMC samples to be drawn. |
hyperparameter |
A boolean indicating if inference should be performed for the prior variance hyperparameter. If TRUE the prior variance (main diagonal of the covariance matrix) must be set to 1. |
chi |
(Optional) The rate parameter for the inverse-gamma prior distribution on the hyperparameter. Default is 0.01. |
psi |
(Optional) The shape parameter for the inverse-gamma prior distribution on the hyperparameter. Default is 0.01. |
burn.in |
(optional) The number of iterations to discard as burn-in. Default is 100. |
n.thin |
(optional) The number of iterations to thin the MCMC samples by. Default is 1. |
verbose |
(optional) A boolean indicating if progress should be printed to the console.
Default is |
If player.prior.var is omitted, independent and identical
N(0, 1^2) prior distributions are placed on each object quality parameter.
If lambda.initial is ommitted, it is set to a vector of zeroes.
An "mcmc" object containing samples from the posterior distribution.
######################################## ## Forced Marriage in Nottinghamshire ## ######################################## # Construct covariance matrix based on spatial information prior.var <- expm::expm(forcedMarriage$adjacencyMatrix) prior.var <- diag(diag(prior.var)^-0.5) %*% prior.var %*% diag(diag(prior.var)^-0.5) # Fit model # Using `n.iter = 3` here to reduce runtime, you will need more # iterations for inference. forcedMarriageModel <- speedyBBTm( outcome = rep(1, length(forcedMarriage$comparisons$win)), player1 = forcedMarriage$comparisons$win, player2 = forcedMarriage$comparisons$lost, player.prior.var = prior.var, n.iter = 3, burn.in = 0 ) # Plot results plot(forcedMarriageModel[, "lambda[1]"], xlab = "Iteration", ylab = expression(lambda[i]) )######################################## ## Forced Marriage in Nottinghamshire ## ######################################## # Construct covariance matrix based on spatial information prior.var <- expm::expm(forcedMarriage$adjacencyMatrix) prior.var <- diag(diag(prior.var)^-0.5) %*% prior.var %*% diag(diag(prior.var)^-0.5) # Fit model # Using `n.iter = 3` here to reduce runtime, you will need more # iterations for inference. forcedMarriageModel <- speedyBBTm( outcome = rep(1, length(forcedMarriage$comparisons$win)), player1 = forcedMarriage$comparisons$win, player2 = forcedMarriage$comparisons$lost, player.prior.var = prior.var, n.iter = 3, burn.in = 0 ) # Plot results plot(forcedMarriageModel[, "lambda[1]"], xlab = "Iteration", ylab = expression(lambda[i]) )
A comparative judgment data set for risk of female genital mutilation at ward level in South Yorkshire.
sy.comparisonssy.comparisons
A data frame with 877 comparisons.
The ID of each comparison.
The ID of the user who made the comparisons.
The ID of the first area involved in the comparison.
The ID of the second area involved in the comparison.
The ID of the selected area. If the
comparison was tied, the selected_item_id is NA.
The state of the outcome. selected indicates that a judgment was made, skipped indicates that the comparison was skipped, and tied indicates that the comparison was tied.
The data was collected following ethical approval the University of Birmingham's Science, Engineering and Maths Ethics Committee.
The outcomes of all 127 men's singles matches in the 2019 Wimbledon champtionship.
wimbledonwimbledon
A list containing a dataframe with the outcomes of the matches and a dataframe describing the players.
Each row of the matchs dataframe corresponds to a match. The columns are:
The id of the winner of the match.
The id of the loser of the match.
A boolean indicating if the match was played in the second week of the tournament.
The outcome of the match. 1 if the winner is item 2, 0 if the winner is item 1.
Data frame describing the players. The columns are:
The name of the player.
The rank of the player in the ATP league table.
The number of points received so far in the ATP 2019 tour prior to Wimbledon starting.
The ID of the player.
http://tennis-data.co.uk/alldata.php