As part of my graduate thesis (area: psychology) I have gathered preference data. The data includes approximately 50000 heads-up comparison between elementX and elementX. I have a total of 15 elements. The participant would be shown 36 comparisons in random order, of which he/she would then choose the preferred element. I need to be able to calculate a meaningful "interval" scale-value for each element in order to see if the actual preferences of each element can predict theoretically generated scale values using regression analysis.

The Data: 14 elements, forming a 14x14 matrix table 91 unique possible comparisons ((14*14-14)/2) = 50000 comparison (NOT evenly distributed across all comparisons)

The problem: Since showing the participant all 91 possible comparison would be too time consuming I've just presented each participant 36 comparisons. Moreover, each participant was not eligible for each element, the eligible elements was determined if the participant owned the element.

For example: ParticipantX owns element; 1, 5, 6, 8 and 10. And would therefore only be shown comparisons that includes those elements.

Furthermore, as owning some elements were considered more rare, I made a priority list for each element, so that element1 > element2 > .... > element14

For example: If the participant would be eligible for 50 possible comparisons (but restricted by the 36 comparison-cap) the comparisons shown to the participant would be based on the rarity order.

I'm afraid I destroyed the data by making this rarity order, since it resulted in some some elements being compared unequally to some elements; For example: of 1000 comparisons that element 1 was included in, 30% was with element3, 10% with element4 etc ...

My questions: What method would you suggest would be the best for calculating "scale" values for each element? Currently I've calculated the basic probability that a given element has been chosen over another element. For example; element1 being chosen 95% of times to be preferred over any other element. However, I would prefer a more refined method to calculate the element values.

  • 1
    $\begingroup$ Bradley-Terry model is supposed to be the way to do this, but I don't understand it. Did you have any luck since you posted the question? $\endgroup$
    – endolith
    Commented Dec 6, 2014 at 20:06

1 Answer 1


Hopefully this isn't too late, but your data is still valuable using the Bradley-Terry model. This paper is a nice reference MM Algorithms for generalized Bradley-Terry models. Equation 4 is the one you want.

Brief explanation as I guess you are not a mathematition, and I'm being nice as I'd really like a copy of your data to test my implementation of this algorithm!

Explanation of Bradley-Terry Model

Basically Bradley-Terry assigns a value to each object, which is the answer you want - e.g. if you are asking "Which of these cars looks sexier", the value would measure sexiness. Then it says if you have two cars with sexinesses $s_i$ and $s_j$, the probability of someone choosing car $i$ over car $j$ is:

$$P(o_i\succ o_j)=\frac{e^{s_i}}{e^{s_i}-e^{s_j}}$$

That actually only depends on $s_j-s_i$ as follows:

bt function

So if they have the same sexiness the chances are 50-50 (random guessing), and if one has a higher sexiness it is more likely to be chosen.

This model is mostly plucked out of the air for its plausibility and nice mathematical properties. I don't know how realistic it really is. It is widely use though.

Maximum Likelihood

BT then uses something called Maximum Likelihood to estimate the values of $s_i$. Put in words, it finds the sexinesses of all the cars that would make the choices that you measured ("Dave says car A is sexier than car B" etc.) as likely as possible, according to our choice model above. So if people choose car A over car B a lot, then giving $s_A$ a lower sexiness than $s_B$ would mean that their choices (A > B) are very unlikely. So those sexiness values are not a maximum likelihood estimate.

To actually find all the $s_i$'s you can use a generic optimisation method (e.g. fminsearch in matlab). Or there is an iterative algorithm that they claim is guaranteed to converge on a "unique maximum" (which I assume means global?).

That algorithm is equation 4 in the linked paper (plus normalisation).


Basically you count the number of times each object is beaten by every other. Store that in $w_{ij}$. The total number of comparisons between any two objects is $N_{ij}=w_{ij}+w_{ji}$. Define new numbers $\gamma_i=e^{s_i}$. Then start with all $\gamma_i=0$ and iteratively update, looping repeatedly through the $i$'s.

$$\gamma_i\leftarrow W_{i}\left[\sum_{j\neq i}\frac{N_{ij}}{\gamma_{i}+\gamma_{j}}\right]^{-1}$$

You don't need to complete a whole loop before updating each $i$. By which I mean you don't need to store old values of $\gamma_i$. As soon as a new value is calculated you can use it in subsequent updates.

$W_i$ is the total number of wins by $i$: $\sum_{j}w_{ij}=W_{i}$

You also need to normalise $\gamma$ after each loop over it:

$$\gamma \leftarrow \frac{\gamma}{\mathrm{mean}(\gamma)}$$


You probably won't need this since your data set is so large and dense. In some cases you need to regularise by adding "phantom" comparisons of every object with an extra virtual one. Apparently the standard method is to count $\lambda$ wins for the real objects vs virtual, and $\lambda$ losses. I've found $\lambda$=0.1 to work fine.

Please can you send me your data and expected results if you still have it? :-)

  • $\begingroup$ Imagine a comparison between $A$ and $B$ as sampling an observed score for $a$ and $b$ from two different underlying distributions, then comparing the scores. Then you get the same sort of sigmoidal curve for $P(A > B)$. However, the implied distribution here is somewhat longer-tailed than a normal distribution. More in next comment as I ran out of space. $\endgroup$
    – TLW
    Commented Apr 16, 2023 at 19:39
  • $\begingroup$ Imagine $x_1 \sim \mathcal{N}(\mu_1, 1)$ and $x_2 \sim \mathcal{N}(\mu_2, 1)$ and look at the asymptotics when $\mu_2 >>> \mu_1$. I get $P(x_1 > x_2) \approx \frac{1}{\sqrt{2\pi}}\frac{e^{\frac{-(\mu_2-\mu_1)^2}{2}}}{\mu_2-\mu_1}$. Compare the tail behavior of BT, where $P(x_1 > x_2) \approx e^{-(\mu_2-\mu_1)}$. Note the linear exponent for BT and quadratic for Gaussian. This means that BT is modelling a heavier-tailed distribution. BT also assumes that all values have the same $\sigma$, which may not be a good approximation. $\endgroup$
    – TLW
    Commented Apr 16, 2023 at 19:42
  • 1
    $\begingroup$ This does make me wonder if there's a simple extension to include $(\mu, \sigma)$ for each variable not just $\mu$. $\endgroup$
    – TLW
    Commented Apr 16, 2023 at 19:44
  • $\begingroup$ I am glossing over the normalization step and the log/exp here I realize. $\endgroup$
    – TLW
    Commented Apr 16, 2023 at 19:44

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.