Bradley Terry Model With Ties

Fits a Bradley Terry model with ties (intercept-only model) by maximum likelihood estimation.

bratt(refgp = "last", refvalue = 1, ialpha = 1, i0 = 0.01)

Arguments

refgp

Integer whose value must be from the set {1,...,\(M\)}, where there are \(M\) competitors. The default value indicates the last competitor is used—but don't input a character string, in general.

refvalue

Numeric. A positive value for the reference group.

ialpha

Initial values for the \(\alpha\)s. These are recycled to the appropriate length.

i0

Initial value for \(\alpha_0\). If convergence fails, try another positive value.

Details

There are several models that extend the ordinary Bradley Terry model to handle ties. This family function implements one of these models. It involves \(M\) competitors who either win or lose or tie against each other. (If there are no draws/ties then use brat). The probability that Competitor \(i\) beats Competitor \(j\) is \(\alpha_i / (\alpha_i+\alpha_j+\alpha_0)\), where all the \(\alpha\)s are positive. The probability that Competitor \(i\) ties with Competitor \(j\) is \(\alpha_0 / (\alpha_i+\alpha_j+\alpha_0)\). Loosely, the \(\alpha\)s can be thought of as the competitors' `abilities', and \(\alpha_0\) is an added parameter to model ties. For identifiability, one of the \(\alpha_i\) is set to a known value refvalue, e.g., 1. By default, this function chooses the last competitor to have this reference value. The data can be represented in the form of a \(M\) by \(M\) matrix of counts, where winners are the rows and losers are the columns. However, this is not the way the data should be inputted (see below).

Excluding the reference value/group, this function chooses \(\log(\alpha_j)\) as the first \(M-1\) linear predictors. The log link ensures that the \(\alpha\)s are positive. The last linear predictor is \(\log(\alpha_0)\).

The Bradley Terry model can be fitted with covariates, e.g., a home advantage variable, but unfortunately, this lies outside the VGLM theoretical framework and therefore cannot be handled with this code.

Value

An object of class "vglmff" (see vglmff-class). The object is used by modelling functions such as vglm.

References

Torsney, B. (2004). Fitting Bradley Terry models using a multiplicative algorithm. In: Antoch, J. (ed.) Proceedings in Computational Statistics COMPSTAT 2004, Physica-Verlag: Heidelberg. Pages 513–526.

Author

T. W. Yee

Note

The function Brat is useful for coercing a \(M\) by \(M\) matrix of counts into a one-row matrix suitable for bratt. Diagonal elements are skipped, and the usual S order of c(a.matrix) of elements is used. There should be no missing values apart from the diagonal elements of the square matrix. The matrix should have winners as the rows, and losers as the columns. In general, the response should be a matrix with \(M(M-1)\) columns.

Also, a symmetric matrix of ties should be passed into Brat. The diagonal of this matrix should be all NAs.

Only an intercept model is recommended with bratt. It doesn't make sense really to include covariates because of the limited VGLM framework.

Notationally, note that the VGAM family function brat has \(M+1\) contestants, while bratt has \(M\) contestants.

Examples

# citation statistics: being cited is a 'win'; citing is a 'loss'
journal <- c("Biometrika", "Comm.Statist", "JASA", "JRSS-B")
mat <- matrix(c( NA, 33, 320, 284,
                730, NA, 813, 276,
                498, 68,  NA, 325,
                221, 17, 142,  NA), 4, 4)
dimnames(mat) <- list(winner = journal, loser = journal)

# Add some ties. This is fictitional data.
ties <- 5 + 0 * mat
ties[2, 1] <- ties[1,2] <- 9

# Now fit the model
fit <- vglm(Brat(mat, ties) ~ 1, bratt(refgp = 1), trace = TRUE,
            crit = "coef")
#> Iteration 1: coefficients = 
#> -1.87315508, -0.35603686,  0.21647979, -4.33749983
#> Iteration 2: coefficients = 
#> -2.62102078, -0.44597533,  0.25173061, -4.53044084
#> Iteration 3: coefficients = 
#> -2.89870009, -0.46059635,  0.25373904, -4.54107239
#> Iteration 4: coefficients = 
#> -2.93502254, -0.46126674,  0.25372840, -4.54144520
#> Iteration 5: coefficients = 
#> -2.93559920, -0.46127181,  0.25372820, -4.54144817
#> Iteration 6: coefficients = 
#> -2.93559935, -0.46127181,  0.25372820, -4.54144817

summary(fit)
#> 
#> Call:
#> vglm(formula = Brat(mat, ties) ~ 1, family = bratt(refgp = 1), 
#>     trace = TRUE, crit = "coef")
#> 
#> Coefficients: 
#>               Estimate Std. Error z value Pr(>|z|)    
#> (Intercept):1 -2.93560    0.10221 -28.721  < 2e-16 ***
#> (Intercept):2 -0.46127    0.05990  -7.701 1.35e-14 ***
#> (Intercept):3  0.25373    0.07003   3.623 0.000291 ***
#> (Intercept):4 -4.54145    0.17595 -25.810  < 2e-16 ***
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> Names of linear predictors: loglink(alpha2), loglink(alpha3), loglink(alpha4), 
#> loglink(alpha0)
#> 
#> Log-likelihood: -1821.474 on 0 degrees of freedom
#> 
#> Number of Fisher scoring iterations: 6 
#> 
#> Warning: Hauck-Donner effect detected in the following estimate(s):
#> '(Intercept):1', '(Intercept):4'
#> 
c(0, coef(fit))  # Log-abilities (last is log(alpha0))
#>               (Intercept):1 (Intercept):2 (Intercept):3 (Intercept):4 
#>     0.0000000    -2.9355993    -0.4612718     0.2537282    -4.5414482 
c(1, Coef(fit))  #     Abilities (last is alpha0)
#>                alpha2     alpha3     alpha4     alpha0 
#> 1.00000000 0.05309889 0.63048128 1.28882145 0.01065796 

fit@misc$alpha   # alpha_1,...,alpha_M
#> [1] 1.00000000 0.05309889 0.63048128 1.28882145
fit@misc$alpha0  # alpha_0
#> [1] 0.01065796

fitted(fit)  # Probabilities of winning and tying, in awkward form
#>   Comm.Statist>Biometrika JASA>Biometrika JRSS-B>Biometrika
#> 1              0.04991637       0.3841729          0.560484
#>   Biometrika>Comm.Statist JASA>Comm.Statist JRSS-B>Comm.Statist Biometrika>JASA
#> 1               0.9400645         0.9081629           0.9528627       0.6093328
#>   Comm.Statist>JASA JRSS-B>JASA Biometrika>JRSS-B Comm.Statist>JRSS-B
#> 1        0.07648512   0.6677967          0.434881          0.03925753
#>   JASA>JRSS-B
#> 1   0.3266809
#> attr(,"probtie")
#>   Comm.Statist==Biometrika JASA==Biometrika JRSS-B==Biometrika
#> 1               0.01001917      0.006494245        0.004634945
#>   Biometrika==Comm.Statist JASA==Comm.Statist JRSS-B==Comm.Statist
#> 1               0.01001917         0.01535202          0.007879737
#>   Biometrika==JASA Comm.Statist==JASA JRSS-B==JASA Biometrika==JRSS-B
#> 1      0.006494245         0.01535202  0.005522372        0.004634945
#>   Comm.Statist==JRSS-B JASA==JRSS-B
#> 1          0.007879737  0.005522372
predict(fit)
#>      loglink(alpha2) loglink(alpha3) loglink(alpha4) loglink(alpha0)
#> [1,]       -2.935599      -0.4612718       0.2537282       -4.541448
(check <- InverseBrat(fitted(fit)))    # Probabilities of winning
#>              Biometrika Comm.Statist       JASA     JRSS-B
#> Biometrika           NA    0.9400645 0.60933282 0.43488104
#> Comm.Statist 0.04991637           NA 0.07648512 0.03925753
#> JASA         0.38417294    0.9081629         NA 0.32668089
#> JRSS-B       0.56048401    0.9528627 0.66779674         NA
qprob <- attr(fitted(fit), "probtie")  # Probabilities of a tie
qprobmat <- InverseBrat(c(qprob), NCo = nrow(ties))  # Pr(tie)
check + t(check) + qprobmat  # Should be 1s in the off-diagonals
#>              Biometrika Comm.Statist JASA JRSS-B
#> Biometrika           NA            1    1      1
#> Comm.Statist          1           NA    1      1
#> JASA                  1            1   NA      1
#> JRSS-B                1            1    1     NA