brat.RdFits a Bradley Terry model (intercept-only model) by maximum likelihood estimation.
brat(refgp = "last", refvalue = 1, ialpha = 1)Integer whose value must be from the set {1,...,\(M+1\)}, where there are \(M+1\) competitors. The default value indicates the last competitor is used—but don't input a character string, in general.
Numeric. A positive value for the reference group.
Initial values for the \(\alpha\)s. These are recycled to the appropriate length.
The Bradley Terry model involves \(M+1\) competitors
who either win or lose against each other (no draws/ties
allowed in this implementation–see bratt
if there are ties). The probability that Competitor
\(i\) beats Competitor \(j\) is \(\alpha_i /
(\alpha_i+\alpha_j)\),
where all the \(\alpha\)s are positive.
Loosely, the \(\alpha\)s can be thought of as
the competitors' `abilities'. 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+1\)
by \(M+1\) 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 \(M\) linear predictors. The log link ensures that the \(\alpha\)s are positive.
The Bradley Terry model can be fitted by logistic regression, but this approach is not taken here. 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.
An object of class "vglmff"
(see vglmff-class).
The object is used by modelling functions
such as vglm.
Agresti, A. (2013). Categorical Data Analysis, 3rd ed. Hoboken, NJ, USA: Wiley.
Stigler, S. (1994). Citation patterns in the journals of statistics and probability. Statistical Science, 9, 94–108.
The BradleyTerry2 package has more comprehensive capabilities than this function.
The function Brat is useful for coercing
a \(M+1\) by \(M+1\) matrix of counts into a one-row
matrix suitable for brat. 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
1-row matrix with \(M(M+1)\) columns.
Only an intercept model is recommended with brat.
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.
Presently, the residuals are wrong, and the prior weights
are not handled correctly. Ideally, the total number of
counts should be the prior weights, after the response has
been converted to proportions. This would make it similar
to family functions such as multinomial
and binomialff.
# 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)
fit <- vglm(Brat(mat) ~ 1, brat(refgp = 1), trace = TRUE)
#> Iteration 1: loglikelihood = -20597.93
#> Iteration 2: loglikelihood = -20525.99
#> Iteration 3: loglikelihood = -20520.45
#> Iteration 4: loglikelihood = -20520.38
#> Iteration 5: loglikelihood = -20520.38
fit <- vglm(Brat(mat) ~ 1, brat(refgp = 1), trace = TRUE, crit = "coef")
#> Iteration 1: coefficients =
#> -1.88777996, -0.36440031, 0.22664468
#> Iteration 2: coefficients =
#> -2.63602167, -0.46252561, 0.26673847
#> Iteration 3: coefficients =
#> -2.91276757, -0.47880148, 0.26897129
#> Iteration 4: coefficients =
#> -2.94851677, -0.47956414, 0.26895435
#> Iteration 5: coefficients =
#> -2.94907236, -0.47956977, 0.26895406
#> Iteration 6: coefficients =
#> -2.94907250, -0.47956977, 0.26895406
summary(fit)
#>
#> Call:
#> vglm(formula = Brat(mat) ~ 1, family = brat(refgp = 1), trace = TRUE,
#> crit = "coef")
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept):1 -2.94907 0.10255 -28.759 < 2e-16 ***
#> (Intercept):2 -0.47957 0.06059 -7.915 2.47e-15 ***
#> (Intercept):3 0.26895 0.07083 3.797 0.000146 ***
#> ---
#> Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#>
#> Names of linear predictors: loglink(alpha2), loglink(alpha3), loglink(alpha4)
#>
#> Log-likelihood: -20520.38 on 0 degrees of freedom
#>
#> Number of Fisher scoring iterations: 6
#>
#> Warning: Hauck-Donner effect detected in the following estimate(s):
#> '(Intercept):1'
#>
c(0, coef(fit)) # Log-abilities (in order of "journal")
#> (Intercept):1 (Intercept):2 (Intercept):3
#> 0.0000000 -2.9490725 -0.4795698 0.2689541
c(1, Coef(fit)) # Abilities (in order of "journal")
#> alpha2 alpha3 alpha4
#> 1.00000000 0.05238827 0.61904967 1.30859502
fitted(fit) # Probabilities of winning in awkward form
#> Comm.Statist>Biometrika JASA>Biometrika JRSS-B>Biometrika
#> 1 0.04978037 0.3823537 0.5668361
#> Biometrika>Comm.Statist JASA>Comm.Statist JRSS-B>Comm.Statist Biometrika>JASA
#> 1 0.9502196 0.921976 0.961507 0.6176463
#> Comm.Statist>JASA JRSS-B>JASA Biometrika>JRSS-B Comm.Statist>JRSS-B
#> 1 0.078024 0.678857 0.4331639 0.03849296
#> JASA>JRSS-B
#> 1 0.321143
(check <- InverseBrat(fitted(fit))) # Probabilities of winning
#> Biometrika Comm.Statist JASA JRSS-B
#> Biometrika NA 0.9502196 0.6176463 0.43316389
#> Comm.Statist 0.04978037 NA 0.0780240 0.03849296
#> JASA 0.38235372 0.9219760 NA 0.32114304
#> JRSS-B 0.56683611 0.9615070 0.6788570 NA
check + t(check) # Should be 1's 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