CODEWITHSUNDEEP
rbradleyterry2
Will Cornwell
Will Cornwell

Reputation: 251

Calculate BradleyTerry Ability rankings

I'm trying to output the estimates of ability from a Bradley Terry model using the BradleyTerry2 package in R. I keep getting a very cryptic error. Some of the examples from their documentation work and others return the same error I get using my data. This code uses one of the example analysis from the documentation. So if you load the library the "chameleon" data should already be there 

install.packages("BradleyTerry2")
library (BradleyTerry2)
summary(chameleon.model <- BTm(player1 = winner, player2 = loser,formula = ~ prev.wins.2 + ch.res[ID] + prop.main[ID] + (1|ID), id = "ID",data = chameleons))
BTabilities(chameleon.model)

And the error I get is 

Error in X[, est, drop = FALSE] : (subscript) logical subscript too long

Anyone know how to do this?

Upvotes: 2

Views: 405

Answers (1)

Heather Turner
Heather Turner

Reputation: 3314

I maintain BradleyTerry2. This is a bug that occurs when you have a contest-specific predictor in the ability formula. It should work as documented in ?BTabilities:

... the abilities are computed from the terms of the fitted model that involve player covariates only (those indexed by ‘model$id’ in the model formula). Thus parameters in any other terms are assumed to be zero.

We weren't aware that this wasn't working, so thanks for the bug report. Until it is fixed, you can work out the abilities and standard errors directly:

## get names of relevant coefficients
> nm <- grep("[ID]", names(coef(chameleon.model)),
>            fixed = TRUE, value = TRUE)
> nm
[1] "ch.res[ID]"    "prop.main[ID]"

## get names of corresponding predictors
> IDvar <- gsub("[ID]", "", nm, fixed = TRUE)
> IDvar
[1] "ch.res"    "prop.main"

## create coefficient vector and design matrix
> cf <- coef(chameleon.model)[nm]
> X <- as.matrix(chameleons$predictors[, IDvar])

## compute abilities
> abilities <- X %*% cf
> colnames(abilities) <- "abilities"

## compute standard errors
> V <- vcov(chameleon.model)[nm, nm]
> res <- cbind(abilities = abilities,
>               se = sqrt(diag(X %*% V %*% t(X))))
> head(res)
    abilities       se
C01  3.757229 1.655205
C02  2.404778 1.017782
C03  2.319346 1.133959
C04  1.892671 1.399391
C05  2.253472 1.101628
C06  2.015840 1.075806

This will give NA for individuals with missing values in the predictors, however their ability is modelled separately and returned in the model summary. The above code also assumes continuous covariates, if some of your predictors were factors, you would need to create a model matrix more formally, e.g.

X <- model.matrix(reformulate(IDvar, intercept = FALSE),
                  c(chameleons$predictors))

Upvotes: 2

Related Questions