Here's a reformulation of your model broken down by component assumptions.

1. We have an unknown matrix of win probabilities, $P=(p_{ij})$, indexed from 1 to 5, where $0 \leq p_{ij} \leq 1$ and $p_{ij}=1-p_{ji}$ for all $i$, $j$. For every $(i,j)$ with $i \leq j$, we observe $n_{ij}$ independent draws from each corresponding Bernoulli($p_{ij}$) distribution.
2. There is some ranking $r(1),...,r(5)$ of the indices where $r(1)$ is the index of the "best" character and $r(5)$ is the index of the worst character. More precisely, we have (i) $p_{r(i)r(j)}\leq 1/2$ for all $i \leq j$ and (ii) $p_{r(i)r(j)} \geq p_{r(i)r(k)}$ for all $j \leq k$.
3. We observe "size" covariates $s_i$, and would like to find a relationship of the form $p_{ij}= f(s_i-s_j) + \epsilon$ where $f$ is some monotonic function and $\epsilon$ is a small error term.

I will describe a simple way of doing it using Maximum Likelihood, before outlining a more technically difficult Bayesian approach.

Maximum Likelihood

Background: see definition of logit, EM Algorithm.

Let $\tilde{P}=(p_{r^{-1}(i)r^{-1}(j)})$ be the matrix $P$ permuted so that the indices correspond to the true ranking. Reparameterize $\tilde{P}$ by writing $$ \tilde{p}_{ij}=(1/2)\text{logit}(\sum_{k=2}^i \sum_{l=1}^{i-1} e^{u_{kl}})+1/2 $$ for $i < j$. This is done so that your transformed variables $u_{21},u_{31},u_{32},...,u_{54}$ have a 1-1 relationship with valid $\tilde{P}$ but can take any real values.
For every possible ranking, use the EM algorithm to find the likelihood of the ranking, and keep the ranking which has the highest likelihood. For this ranking, use the maximum likelihood estimates of $u_{21},...,u_{54}$ to convert back to the ML estimate of the matrix $P$, which we call $\hat{P}=(\hat{p}_{ij})$.

Now, to determine the "consistency" of the relationship between $p_{ij}$ and $s_i-s_j$, carry out a logistic regression $\hat{p}_{ij} \sim \text{logit}(\beta_0 + \beta_1 (s_i-s_j))$ for all $i < j$. The residuals give you an idea of the consistency. You can try high-order regression formulae, (eg quadratic) but then you may not preserve monotonicity.

Bayesian Approach

Specify a prior on $P$ by assigning uniform (or optionally, beta-weighted) probability to all 5 by 5 matrices with entries in $[0,1]$ which satisfy condition 2, and zero probability to all such matrices which violate condition 2; this gives you the prior density $p(P)$. Sample from the posterior be using resampling (see any intro text on Markov Chain Monte Carlo). From each posterior draw, compute the posterior of the sum or residuals (or whatever measure of consistency you use) based on priors for the regression coefficients $\beta_0, \beta_1$, weighting by the likelihood $$ \exp[-\sum_{i <j}(\text{logit}(\beta_0 + \beta_1 (s_i-s_j))-p_{ij})^2] $$ to get the joint posterior for $P$ and your measure of "consistency."

Here's a reformulation of your model broken down by component assumptions.

1. We have an unknown matrix of win probabilities, $P=(p_{ij})$, indexed from 1 to 5, where $0 \leq p_{ij} \leq 1$ and $p_{ij}=1-p_{ji}$ for all $i$, $j$. For every $(i,j)$ with $i \leq j$, we observe $n_{ij}$ independent draws from each corresponding Bernoulli($p_{ij}$) distribution.
2. There is some ranking $r(1),...,r(5)$ of the indices where $r(1)$ is the index of the "best" character and $r(5)$ is the index of the worst character. More precisely, we have (i) $p_{r(i)r(j)}\leq 1/2$ for all $i \leq j$ and (ii) $p_{r(i)r(j)} \geq p_{r(i)r(k)}$ for all $j \leq k$.
3. We observe "size" covariates $s_i$, and would like to find a relationship of the form $p_{ij}= f(s_i-s_j) + \epsilon$ where $f$ is some monotonic function and $\epsilon$ is a small error term.

I will describe a simple way of doing it using Maximum Likelihood, before outlining a more technically difficult Bayesian approach.

Maximum Likelihood

Background: see definition of logit, EM Algorithm.

Let $\tilde{P}=(p_{r^{-1}(i)r^{-1}(j)})$ be the matrix $P$ permuted so that the indices correspond to the true ranking. Reparameterize $\tilde{P}$ by writing $$ \tilde{p}_{ij}=(1/2)\text{logit}(\sum_{k=i}^5 \sum_{l=1}^{i-1} e^{u_{kl}})+1/2 $$ for $i < j$. This is done so that your transformed variables $u_{21},u_{31},u_{32},...,u_{54}$ have a 1-1 relationship with valid $\tilde{P}$ but can take any real values.
For every possible ranking, use the EM algorithm to find the likelihood of the ranking, and keep the ranking which has the highest likelihood. For this ranking, use the maximum likelihood estimates of $u_{21},...,u_{54}$ to convert back to the ML estimate of the matrix $P$, which we call $\hat{P}=(\hat{p}_{ij})$.

Now, to determine the "consistency" of the relationship between $p_{ij}$ and $s_i-s_j$, carry out a logistic regression $\hat{p}_{ij} \sim \text{logit}(\beta_0 + \beta_1 (s_i-s_j))$ for all $i < j$. The residuals give you an idea of the consistency. You can try high-order regression formulae, (eg quadratic) but then you may not preserve monotonicity.

Bayesian Approach

Specify a prior on $P$ by assigning uniform (or optionally, beta-weighted) probability to all 5 by 5 matrices with entries in $[0,1]$ which satisfy condition 2, and zero probability to all such matrices which violate condition 2; this gives you the prior density $p(P)$. Sample from the posterior be using resampling (see any intro text on Markov Chain Monte Carlo). From each posterior draw, compute the posterior of the sum or residuals (or whatever measure of consistency you use) based on priors for the regression coefficients $\beta_0, \beta_1$, weighting by the likelihood $$ \exp[-\sum_{i <j}(\text{logit}(\beta_0 + \beta_1 (s_i-s_j))-p_{ij})^2] $$ to get the joint posterior for $P$ and your measure of "consistency."

Here's a reformulation of your model broken down by component assumptions.

1. We have an unknown matrix of win probabilities, $P=(p_{ij})$, indexed from 1 to 5, where $0 \leq p_{ij} \leq 1$ and $p_{ij}=1-p_{ji}$ for all $i$, $j$. For every $(i,j)$ with $i \leq j$, we observe $n_{ij}$ independent draws from each corresponding Bernoulli($p_{ij}$) distribution.
2. There is some ranking $r(1),...,r(5)$ of the indices where $r(1)$ is the index of the "best" character and $r(5)$ is the index of the worst character. More precisely, we have (i) $p_{r(i)r(j)}\leq 1/2$ for all $i \leq j$ and (ii) $p_{r(i)r(j)} \geq p_{r(i)r(k)}$ for all $j \leq k$.
3. There exists a "size" function $q$ with $q(r(i)) \geq q(r(j))$ for all $i \leq j$, such that $p_{ij} \geq p_{kl}$ if and only if $q(i)-q(j) \geq q(k)-q(l)$.

With regards to the model, all size functions are equivalent with respect to the ranks of the differences $d_{ij}=q(r(i))-q(r(j))$. Letting $o_{ij}$ be the ordinal rank of $d_{ij}$ (that is, $o_{ij} \leq o_{kl}$ if and only if $d_{ij} \geq d_{kl}$), the theoretical difficulty lies in understanding what matrices of rankings $O=(o_{ij})$ are possible. (to be contd.)

Here's a reformulation of your model broken down by component assumptions.

1. We have an unknown matrix of win probabilities, $P=(p_{ij})$, indexed from 1 to 5, where $0 \leq p_{ij} \leq 1$ and $p_{ij}=1-p_{ji}$ for all $i$, $j$. For every $(i,j)$ with $i \leq j$, we observe $n_{ij}$ independent draws from each corresponding Bernoulli($p_{ij}$) distribution.
2. There is some ranking $r(1),...,r(5)$ of the indices where $r(1)$ is the index of the "best" character and $r(5)$ is the index of the worst character. More precisely, we have (i) $p_{r(i)r(j)}\leq 1/2$ for all $i \leq j$ and (ii) $p_{r(i)r(j)} \geq p_{r(i)r(k)}$ for all $j \leq k$.
3. We observe "size" covariates $s_i$, and would like to find a relationship of the form $p_{ij}= f(s_i-s_j) + \epsilon$ where $f$ is some monotonic function and $\epsilon$ is a small error term.

I will describe a simple way of doing it using Maximum Likelihood, before outlining a more technically difficult Bayesian approach.

Maximum Likelihood

Background: see definition of logit, EM Algorithm.

Let $\tilde{P}=(p_{r^{-1}(i)r^{-1}(j)})$ be the matrix $P$ permuted so that the indices correspond to the true ranking. Reparameterize $\tilde{P}$ by writing $$ \tilde{p}_{ij}=(1/2)\text{logit}(\sum_{k=2}^i \sum_{l=1}^{i-1} e^{u_{kl}})+1/2 $$ for $i < j$. This is done so that your transformed variables $u_{21},u_{31},u_{32},...,u_{54}$ have a 1-1 relationship with valid $\tilde{P}$ but can take any real values.
For every possible ranking, use the EM algorithm to find the likelihood of the ranking, and keep the ranking which has the highest likelihood. For this ranking, use the maximum likelihood estimates of $u_{21},...,u_{54}$ to convert back to the ML estimate of the matrix $P$, which we call $\hat{P}=(\hat{p}_{ij})$.

Now, to determine the "consistency" of the relationship between $p_{ij}$ and $s_i-s_j$, carry out a logistic regression $\hat{p}_{ij} \sim \text{logit}(\beta_0 + \beta_1 (s_i-s_j))$ for all $i < j$. The residuals give you an idea of the consistency. You can try high-order regression formulae, (eg quadratic) but then you may not preserve monotonicity.

Bayesian Approach

Specify a prior on $P$ by assigning uniform (or optionally, beta-weighted) probability to all 5 by 5 matrices with entries in $[0,1]$ which satisfy condition 2, and zero probability to all such matrices which violate condition 2; this gives you the prior density $p(P)$. Sample from the posterior be using resampling (see any intro text on Markov Chain Monte Carlo). From each posterior draw, compute the posterior of the sum or residuals (or whatever measure of consistency you use) based on priors for the regression coefficients $\beta_0, \beta_1$, weighting by the likelihood $$ \exp[-\sum_{i <j}(\text{logit}(\beta_0 + \beta_1 (s_i-s_j))-p_{ij})^2] $$ to get the joint posterior for $P$ and your measure of "consistency."