Contest models: from the multinomial to competitive breeding values

This vignette derives the two estimators BreedingR uses for competitive ability measured through grouped contests: the Laplace estimator of the genetic variance of contest strength, and the composite-likelihood test for a rival-dependent (differential) indirect effect. The object throughout is a group contest: several competitors share a group, a fixed number of discrete outcomes is divided among them, and the division is observed.

1. The contest likelihood

A group $g$ holds competitors $i \in \mathcal{P}_g$. Competitor $i$ enters with a known exposure $s_{i}>0$ (a dose, a number of trials) and carries a latent strength $\lambda_i>0$. The $N_g$ outcomes of the group are divided by a biased lottery: each outcome falls to $i$ with probability

$$p_{ig} \;=\; \frac{s_{i}\,\lambda_i}{\sum_{k\in\mathcal{P}_g} s_{k}\,\lambda_k}, \qquad (y_{1g},\dots,y_{kg}) \;\sim\; \mathrm{Multinomial}\!\left(N_g;\, p_{1g},\dots,p_{kg}\right).$$

Write $\eta_i=\log\lambda_i$ and take the exposure as an offset $o_i=\log s_i$. With two competitors this is the Bradley-Terry model (Bradley and Terry, 1952); with more it is Plackett-Luce (Luce, 1959; Plackett, 1975), equivalently the conditional logit. The strengths are defined only up to a common scale: multiplying every $\lambda$ in a group by the same constant leaves every $p_{ig}$ unchanged, because the constant factors out of numerator and denominator. Only differences of $\eta$ are identifiable, and an anchor is required (a zero-mean random term supplies one).

2. Reduction to a pair, by exact conditioning

Take a pair $\{i,j\}\subseteq\mathcal{P}_g$ and condition on their joint count $m=y_{ig}+y_{jg}$. The claim is that the conditional distribution of $y_{ig}$ is Binomial and depends on no other competitor. Start from the multinomial mass function and hold every count except $y_{ig},y_{jg}$ fixed:

$$\Pr\!\left(y_{ig}=a,\; y_{jg}=m-a \mid \{y_{kg}\}_{k\neq i,j}\right) \;=\; \frac{\dfrac{N_g!}{a!\,(m-a)!\,\prod_{k\neq i,j} y_{kg}!}\; p_{ig}^{\,a}\,p_{jg}^{\,m-a}\,\prod_{k\neq i,j} p_{kg}^{\,y_{kg}}} {\displaystyle\sum_{a'=0}^{m}\dfrac{N_g!}{a'!\,(m-a')!\,\prod_{k\neq i,j} y_{kg}!}\; p_{ig}^{\,a'}\,p_{jg}^{\,m-a'}\,\prod_{k\neq i,j} p_{kg}^{\,y_{kg}}}.$$

Every factor that does not involve $a$ (the $N_g!$, the $\prod_{k\neq i,j}$ terms, the fixed other counts) is identical in the numerator and in every summand of the denominator, so it cancels:

$$\Pr\!\left(y_{ig}=a \mid y_{ig}+y_{jg}=m\right) \;=\; \frac{\dbinom{m}{a}\,p_{ig}^{\,a}\,p_{jg}^{\,m-a}} {\displaystyle\sum_{a'=0}^{m}\dbinom{m}{a'}\,p_{ig}^{\,a'}\,p_{jg}^{\,m-a'}}.$$

The denominator is a binomial expansion, $\sum_{a'=0}^{m}\binom{m}{a'}p_{ig}^{\,a'}p_{jg}^{\,m-a'}=(p_{ig}+p_{jg})^{m}$, so

$$y_{ig}\mid(y_{ig}+y_{jg}=m)\;\sim\;\mathrm{Binomial}(m,\pi_{ij}), \qquad \pi_{ij}=\frac{p_{ig}}{p_{ig}+p_{jg}}=\frac{s_i\lambda_i}{s_i\lambda_i+s_j\lambda_j}.$$

The result is exact, not an approximation: conditioning on the pair total removes the rest of the group from the problem, because the other competitors’ probabilities do not survive the cancellation. Taking the log-odds,

$$\mathrm{logit}\,\pi_{ij} =\log\frac{p_{ig}}{p_{jg}} =(o_i-o_j)+(\eta_i-\eta_j). \tag{2.1}$$

Section 3 shows that the pair is also the only place where a rival-dependent effect is identifiable.

3. What is estimable: the aliasing of a uniform indirect effect

Give each competitor a direct log-strength $a_i$ and an indirect effect it exerts on its groupmates. Write the realized log-strength of $i$ in group $g$ as its own ability plus the effects the others exert on it:

$$\eta_{ig}=a_i+\sum_{j\in\mathcal{P}_g,\,j\neq i} c_{j\to i},$$

where $c_{j\to i}$ is the effect of $j$ on $i$. Two cases.

Uniform indirect effect. $j$ acts on every groupmate alike, $c_{j\to i}=c_j$. Then

$$\eta_{ig}=a_i+\sum_{j\neq i}c_j =a_i+\Big(\underbrace{\textstyle\sum_{j\in\mathcal{P}_g}c_j}_{\;=\,C_g,\ \text{constant in }i}\Big)-c_i =(a_i-c_i)+C_g.$$

Substituting into the choice probability, the group constant $C_g$ cancels as the common scale did in Section 1:

$$p_{ig}=\frac{s_i\,e^{\,(a_i-c_i)+C_g}}{\sum_{k} s_k\,e^{\,(a_k-c_k)+C_g}} =\frac{s_i\,e^{\,a_i-c_i}}{\sum_{k} s_k\,e^{\,a_k-c_k}},$$

because $e^{C_g}$ multiplies numerator and denominator. Only the combination $b_i=a_i-c_i$ enters. Hence the reparametrization

$$(a_i,\,c_i)\;\longmapsto\;(a_i+t,\;c_i+t)\quad\text{for any constant }t$$

leaves $b_i$, and therefore every probability and the whole likelihood, exactly unchanged. Direct and uniform-indirect effects are perfectly aliased: no estimator separates them, at any sample size. (On the share scale, with outcomes as proportions summing to one, the same collinearity appears as a direct-indirect correlation pinned at $-1$, where the gradient does not vanish and a delta-method standard error is meaningless.)

Differential (rival-dependent) indirect effect. $j$ acts differently depending on who the rival is, for instance by the rival’s lineage. Model the departure from uniform as

$$c_{j\to i}=c_j+\delta_j\,\mathbb{1}[\,\text{class}(i)\neq\text{class}(j)\,].$$

Now form the pair contrast (2.1). The own abilities and the uniform parts enter $\eta_i-\eta_j$; the group constant $C_g$ and any term common to $i$ and $j$ cancel in the difference, but the class-specific term does not, because it is asymmetric in $i$ and $j$:

$$\mathrm{logit}\,\pi_{ij} =(o_i-o_j)+(\eta_i-\eta_j) =(o_i-o_j)+(a_i-a_j) +\delta_i\,\mathbb{1}[\text{class}(j)\neq\text{class}(i)] -\delta_j\,\mathbb{1}[\text{class}(i)\neq\text{class}(j)]. \tag{3.1}$$

The pair is where the differential effect becomes visible: a uniform shift adds the same amount to $\eta_i$ and $\eta_j$ and cancels in the difference, while a rival-dependent shift survives. The estimable quantity $\delta_i$ is the change in $i$’s log-odds of winning an outcome from $j$ when $j$ is of a different class, over and above the two competitors’ abilities.

4. Composite likelihood

Break every group into its pairs and multiply the pair binomials as if independent:

$$c\ell(\boldsymbol{\eta},\boldsymbol{\delta}) =\sum_{\text{pairs }(i,j)}\log\mathrm{Binom}\!\left(y_{ig}\,;\,m_{ij},\,\pi_{ij}\right). \tag{4.1}$$

The pairs are not independent: a competitor appears in many of them, and within a group the counts are constrained to sum to $N_g$. So (4.1) is a composite likelihood (Lindsay, 1988; Varin, Reid and Firth, 2011), not a likelihood. Two properties of a likelihood have to be examined separately: the consistency of the estimator, and the distribution of the ratio statistic.

Consistency holds. Each factor is a genuine conditional probability (Section 2), so each per-pair score has mean zero at the true parameter: $\mathbb{E}\!\left[\nabla\log\mathrm{Binom}_p\right]=0$. The composite score is a sum of these,

$$\nabla c\ell=\sum_p \nabla\log\mathrm{Binom}_p, \qquad \mathbb{E}\!\left[\nabla c\ell\right]=\sum_p \mathbb{E}\!\left[\nabla\log\mathrm{Binom}_p\right]=0,$$

so the estimating equation $\nabla c\ell=0$ is unbiased and its root is consistent under the usual regularity conditions. Unbiasedness of the score requires the marginals to be correct, not the pairs to be independent; this is what justifies the composite construction.

The second Bartlett identity fails. For a true likelihood, the expected negative Hessian equals the variance of the score, $H\equiv-\mathbb{E}[\nabla^2\ell]=\operatorname{Var}(\nabla\ell)\equiv J$. For the composite likelihood the two objects differ:

$$H=-\mathbb{E}\!\left[\nabla^2 c\ell\right]=\sum_p -\mathbb{E}\!\left[\nabla^2\log\mathrm{Binom}_p\right],$$

a sum of per-pair informations, whereas

$$J=\operatorname{Var}\!\left(\nabla c\ell\right) =\sum_p\sum_q \operatorname{Cov}\!\left(\nabla\log\mathrm{Binom}_p,\;\nabla\log\mathrm{Binom}_q\right).$$

The off-diagonal $p\neq q$ terms are nonzero because pairs sharing a competitor are dependent, so $J\neq H$. Two consequences follow.

1. The curvature is not the precision. The variance of $\hat{\boldsymbol{\delta}}$ is the Godambe (sandwich) matrix $G^{-1}=H^{-1}JH^{-1}$ (Godambe, 1960), not $H^{-1}$. A standard error read off the composite Hessian is invalid.

2. The likelihood-ratio statistic is not $\chi^2$. For a true likelihood $2[\ell(\hat\theta)-\ell(\theta_0)]\to\chi^2_q$. For a composite likelihood the same Taylor expansion gives (Kent, 1982)

   $$2\!\left[c\ell(\hat\theta)-c\ell(\theta_0)\right]\;\xrightarrow{d}\;\sum_{k=1}^{q}\lambda_k Z_k^2, \qquad Z_k\stackrel{iid}{\sim}\mathcal N(0,1),$$

   where $\lambda_1,\dots,\lambda_q$ are the eigenvalues of $H^{-1}J$. When $H=J$ (a true likelihood) every $\lambda_k=1$ and the mixture collapses to $\chi^2_q$. Here the $\lambda_k$ are unknown, so the reference distribution is unknown and a $\chi^2$ table does not apply.

5. Calibration by permutation

Because the $\lambda_k$ are unknown, and because overdispersion (a batch or handling effect shared within a group) inflates $J$ beyond what the multinomial predicts, the null distribution of the statistic is obtained empirically. The null must reproduce everything that shapes the statistic (the dependence among pairs, the overdispersion, the group sizes, the exposures) and break only the association under test, between the rival’s class and the outcome. A permutation that relabels which groups are class-mixed, among groups of the same size, with the counts and exposures untouched, does this: it severs class from outcome while preserving the dependence structure and the dispersion. The empirical $p$-value is

$$p=\frac{1+\#\{\,\text{permutation statistic}\ \geq\ \text{observed}\,\}}{1+B},$$

whose floor is $1/(1+B)$ for $B$ permutations. A parametric bootstrap that samples from the fitted multinomial is not a valid null here: it takes the multinomial as the truth and so omits the overdispersion, making the null too tight. It is kept as a diagnostic: the gap between the permutation null and the parametric null measures the overdispersion the multinomial cannot hold.

6. The variance components, by Laplace

The differential test estimates $\boldsymbol{\delta}$ as location parameters. The companion question, how much of contest strength is heritable, puts a pedigree-structured random effect on the strengths and asks for its variance. With $\eta_i=a_i+q_i$, $\mathbf{a}\sim\mathcal N(\mathbf 0,\sigma^2_a\mathbf A)$ (additive, tied by the pedigree) and $\mathbf{q}\sim\mathcal N(\mathbf 0,\sigma^2_q\mathbf I)$ (a permanent competitor effect), the variance is a marginal likelihood over the random effects,

$$L(\sigma^2_a,\sigma^2_q)=\int \Big[\textstyle\prod_g \mathrm{Multinomial}\big(\mathbf y_g\mid \mathbf p(\mathbf u)\big)\Big]\, \mathcal N(\mathbf a\mid\mathbf 0,\sigma^2_a\mathbf A)\,\mathcal N(\mathbf q\mid\mathbf 0,\sigma^2_q\mathbf I)\;d\mathbf u, \qquad \mathbf u=(\mathbf a,\mathbf q).$$

There is no closed form: the integral is in the dimension of the pedigree. In the Gaussian linear model the integral of a normal density against a normal prior is itself normal and closes, which is what makes REML available there; with a multinomial likelihood it does not close and has to be approximated.

Laplace’s method (Tierney and Kadane, 1986). Write the integrand as $e^{-g(\mathbf u)}$ with

$$g(\mathbf u)=-\ell(\mathbf u)+\tfrac12\,\mathbf a'\mathbf A^{-1}\mathbf a/\sigma^2_a+\tfrac12\,\mathbf q'\mathbf q/\sigma^2_q,$$

where $\ell$ is the multinomial log-likelihood. Expand $g$ to second order about its minimizer $\hat{\mathbf u}$ (the joint mode). The gradient vanishes there by definition of the mode, $\nabla g(\hat{\mathbf u})=\mathbf 0$, so the first-order term of the Taylor expansion is zero and the expansion is

$$g(\mathbf u)\approx g(\hat{\mathbf u})+\tfrac12(\mathbf u-\hat{\mathbf u})'\mathbf H(\mathbf u-\hat{\mathbf u}), \qquad \mathbf H=\nabla^2 g(\hat{\mathbf u})=\mathbf W+\begin{pmatrix}\mathbf A^{-1}/\sigma^2_a & \\ & \mathbf I/\sigma^2_q\end{pmatrix},$$

with $\mathbf W$ the multinomial information (below). The remaining integral is Gaussian and closes, $\int e^{-\frac12(\mathbf u-\hat{\mathbf u})'\mathbf H(\mathbf u-\hat{\mathbf u})}d\mathbf u=(2\pi)^{n/2}|\mathbf H|^{-1/2}$, giving

$$-2\log L \;\approx\; \underbrace{-2\ell(\hat{\mathbf u})}_{\text{fit at the mode}} +\underbrace{\mathbf a'\mathbf A^{-1}\mathbf a/\sigma^2_a+\mathbf q'\mathbf q/\sigma^2_q}_{\text{pedigree and permanent penalties}} +\underbrace{n_a\log\sigma^2_a-\log|\mathbf A^{-1}|+n_q\log\sigma^2_q}_{\text{prior normalizers}} +\underbrace{\log|\mathbf H|}_{\text{cost of estimating }\hat{\mathbf u}}. \tag{6.1}$$

The multinomial information $\mathbf W$. From $\ell$, with $\eta$ on the log scale,

$$\frac{\partial \ell}{\partial\eta_i}=\sum_g\big(y_{ig}-N_g p_{ig}\big), \qquad -\frac{\partial^2 \ell}{\partial\eta_i\,\partial\eta_j}=\sum_g N_g\big(\delta_{ij}\,p_{ig}-p_{ig}p_{jg}\big),$$

so each group contributes the block $N_g\big(\mathrm{diag}(\mathbf p_g)-\mathbf p_g\mathbf p_g'\big)$. Two facts follow. First, $\mathbf W$ is sparse: entry $(i,j)$ is nonzero only if $i$ and $j$ shared a group. The sparsity comes from the design, and it is what makes the method affordable: $\mathbf H=\mathbf W+(\text{block-diagonal precision})$ is sparse, and a single sparse Cholesky factorization delivers both the Newton step for the mode and the $\log|\mathbf H|$ that (6.1) needs. Second, $\mathbf W$ is a weighted graph Laplacian: its rows sum to zero ($\sum_j N_g(\delta_{ij}p_{ig}-p_{ig}p_{jg})=N_g(p_{ig}-p_{ig})=0$), so $\mathbf W$ is singular on its own; the precision term $\mathbf A^{-1}/\sigma^2_a$ is what makes $\mathbf H$ positive definite.

The role of $\log|\mathbf H|$. It is the price of having estimated the mode rather than known it; without it $\sigma^2$ comes out systematically underestimated. It is not the REML penalty. Writing the Gaussian restricted likelihood in the standard form,

$$-2\log L_{\text{REML}}=(\mathbf y-\mathbf X\mathbf b)'\mathbf V^{-1}(\mathbf y-\mathbf X\mathbf b)+\log|\mathbf V|+\log|\mathbf X'\mathbf V^{-1}\mathbf X|,$$

the Laplace term corresponds to $\log|\mathbf V|$ (through the identity $\log|\mathbf V|=\log|\mathbf R|+\log|\mathbf G|+\log|\mathbf H|$), which is the maximum-likelihood term. The restricted term proper is $\log|\mathbf X'\mathbf V^{-1}\mathbf X|$, which penalizes the fixed effects for having been estimated, and (6.1) does not carry it. The estimator is therefore marginal (ML) over the random effects, not REML; REML would additionally integrate the fixed effects, which the same expansion could do. When the design has an offset and no free fixed effect, as the contest model does with exposure entering as an offset, the distinction does not change the numbers, but it does as soon as a fixed effect enters.

The variance components are then estimated by minimizing (6.1) over $(\sigma^2_a,\sigma^2_q)$, re-finding the mode at each trial value. The surface of an additive/permanent model has a ridge: trading $\sigma^2_a$ for $\sigma^2_q$ at fixed sum barely moves the fit. A point estimate is therefore not enough, and the interval on the heritable fraction $\sigma^2_a/(\sigma^2_a+\sigma^2_q)$ is read from a profile of (6.1) along that fraction. The profile is valid here because the multinomial is a true likelihood and the profile-likelihood ratio keeps its asymptotic calibration; this is the point at which the composite likelihood of Section 4 and the likelihood of this section part ways, since there $H\neq J$ removes the $\chi^2$ calibration and here it holds.

7. Heritability of a contest-derived effect

Suppose the differential effects $\hat{\boldsymbol{\delta}}$ of Sections 3 to 5 are shown to exist and the question becomes whether they are heritable. The estimates are a first stage, so this is a two-stage problem, and the way the first stage is carried into the second decides the answer.

The measurement-error model. Write

$$\hat{\boldsymbol{\delta}}=\boldsymbol{\delta}+\mathbf e,\qquad \boldsymbol{\delta}\sim\mathcal N(\mathbf 0,\sigma^2_\delta\mathbf A+\sigma^2_{\text{env}}\mathbf I),\qquad \operatorname{Var}(\mathbf e)=\mathbf V_e,$$

where $\mathbf V_e$ is the Godambe sandwich of Section 4, clustered by group, which is what makes it valid under both the pair dependence and any overdispersion. The marginal covariance of the observed vector is then

$$\operatorname{Var}(\hat{\boldsymbol{\delta}})=\sigma^2_\delta\mathbf A+\sigma^2_{\text{env}}\mathbf I+\mathbf V_e. \tag{7.1}$$

$\mathbf V_e$ enters with its coefficient fixed, and that is the whole point. It is tempting to fit $\sigma^2_\delta\mathbf A+\sigma^2_{\text{env}}\mathbf I+\phi\,\mathbf V_e$ with $\phi$ free, on the argument that a free $\phi$ near one would confirm the sandwich. It does not work, and the reason is structural rather than numerical: $\phi\mathbf V_e$ with $\phi$ free is the multiplicative model for heterogeneity and $\sigma^2_{\text{env}}\mathbf I+\mathbf V_e$ is the additive one (Thompson and Sharp, 1999), and the two components $\phi$ and $\sigma^2_{\text{env}}$ are not simultaneously identifiable whenever the $v_i$ are not strongly heterogeneous, because $\mathbf V_e$ is then close to proportional to $\mathbf I$ and the two columns of the design in variance space collapse. Fitted with $\phi$ free on simulated contests, $\phi$ settles between 2.4 and 4.8, $\sigma^2_{\text{env}}$ goes to the zero boundary, and $\hat h^2_\delta$ returns essentially 1.0 against a true 0.6. Fixing the coefficient at one is what Piepho, Möhring, Schulz-Streeck and Ogutu (2012) call the fully efficient two-stage analysis, the only version that reproduces a single-stage fit; carrying a diagonal approximation instead is the standard shortcut of the two-stage literature in breeding (Garrick, Taylor and Fernando, 2009, state plainly that they ignore the off-diagonal of the error covariance).

With the coefficient fixed, (7.1) is an ordinary restricted likelihood in two free components,

$$-2\log L_{\text{REML}}=\log|\mathbf V|+\log|\mathbf X'\mathbf V^{-1}\mathbf X|+\mathbf y'\mathbf P\mathbf y, \qquad \mathbf V=\sigma^2_\delta\mathbf A+\sigma^2_{\text{env}}\mathbf I+\mathbf V_e,$$

with $\mathbf X$ the intercept alone, since $\boldsymbol{\delta}$ is identified only up to an additive constant and the restricted likelihood projects that direction out.

What it delivers depends on how badly the second Bartlett identity fails, so it has to be measured against that severity rather than reported as one number. Taking the median of $\operatorname{diag}(\mathbf H^{-1}\mathbf J\mathbf H^{-1})/\operatorname{diag}(\mathbf H^{-1})$ as the severity index, ten replicates in each of four regimes on simulated contests with a half-sib pedigree give:

regime	severity	$\hat\sigma^2_\delta$ (true 0.30)	$\hat h^2_\delta$ (true)	coverage
no environmental variance	1.4	0.24	0.85 (1.00, a boundary)	0.93
interior	1.4	0.31	0.63 (0.60)	0.93
interior, overdispersed	2.2	0.30	0.69 (0.60)	0.92
interior, strongly overdispersed	4.8	0.18	0.38 (0.40)	0.90

The two ends of that table read in opposite directions, and neither alone would have been safe to generalise. Where the identity fails mildly, $\sigma^2_\delta$ is recovered and $h^2_\delta$ runs optimistic by 0.03 to 0.09, because $\sigma^2_{\text{env}}$ comes back low and $h^2$ cannot exceed one, so replicates that reach the zero boundary are truncated and pull the mean up. Where it fails at the severity a real contest design can reach, $h^2_\delta$ is instead almost unbiased and it is $\sigma^2_\delta$ that deflates, by about forty percent. So on a real design the heritability is the reading to trust and the additive variance is a floor, and on a well-behaved one it is the other way round. Interval coverage holds throughout, and the gap against the unweighted analysis widens exactly where it should: its coverage falls to 0.55 in the last row.

The same severity governs the identifiability discussed above. The free-scale fit breaks (returns $h^2=1$) in the first three regimes and stops breaking in the fourth, because heavy overdispersion makes the $v_i$ genuinely heterogeneous and $\mathbf V_e$ stops being nearly proportional to $\mathbf I$. That is the Thompson and Sharp mechanism seen from the outside: the non-identifiability is a property of the design, not of the fitter.

Three readings that look like diagnostics and are not. Each was measured on simulated contests with a known truth, and each says less than it appears to.

• A weighted animal model with $w_i=1/v_i$, reading the residual variance as a calibration check. The residual it returns is not the scale of $\mathbf V_e$: it satisfies, to within a few percent, the identity $\hat\sigma^2_e\approx\bigl(\operatorname{Var}(\hat{\boldsymbol{\delta}})-\hat\sigma^2_\delta\bigr)/\bar v$. Since $\bar v$ is small next to the dispersion of $\hat{\boldsymbol{\delta}}$, everything the pedigree fails to explain, the environmental variance of $\boldsymbol{\delta}$ included, is divided by a small number and reappears magnified. Replacing $\hat{\boldsymbol{\delta}}$ by a synthetic phenotype whose error is exactly independent $\mathcal N(0,v_i)$ reproduces the same inflated residual, which settles that it never measured calibration.
• The share of the sandwich mass sitting off the diagonal. It can exceed the diagonal mass several times over and still be harmless: the off-diagonal is dominated by one global mode, the non-identifiable constant, and the restricted likelihood projects that direction out. The whitened quadratic forms built from the full matrix and from its diagonal agree to within about seven percent.
• Shrinkage from the ridge on $\boldsymbol{\delta}$. By the Ghosh-Louis result (Louis, 1984), a shrunken set of estimates understates the dispersion of the parameters, so the ridge biases the genetic fraction downward, never upward. Lowering it by two orders of magnitude restores the regression of $\hat{\boldsymbol{\delta}}$ on $\boldsymbol{\delta}$ to a slope of one and leaves the variance components no better.

The coupled alternative, and why it is not simply better. The obvious objection to any two-stage scheme is that a single stage would avoid the transport question altogether: put $\boldsymbol{\delta}$ inside the composite likelihood with the pedigree prior and integrate by the Laplace machinery of Section 6. That construction is available, and on the surface it looks strictly preferable. It is not, and the reason is the same broken identity as in Section 4, arriving in a place where it is much harder to see.

Expand the composite log-likelihood about its maximum, $c\ell(\boldsymbol{\delta})\approx c-\tfrac12(\boldsymbol{\delta}-\mathbf d)'\mathbf H(\boldsymbol{\delta}-\mathbf d)$ with $\mathbf d=\hat{\boldsymbol{\delta}}$. Integrating against the prior $\mathcal N(\mathbf 0,\mathbf G)$, $\mathbf G=\sigma^2_\delta\mathbf A+\sigma^2_{\text{env}}\mathbf I$, gives a closed form,

$$\int e^{\,c\ell(\boldsymbol{\delta})}\,\mathcal N(\boldsymbol{\delta}\mid\mathbf 0,\mathbf G)\,d\boldsymbol{\delta} =e^{c}(2\pi)^{p/2}|\mathbf H|^{-1/2}\, \mathcal N\!\left(\mathbf d\mid\mathbf 0,\;\mathbf G+\mathbf H^{-1}\right),$$

so in this regime the coupled estimator is exactly maximum likelihood on $\hat{\boldsymbol{\delta}}$ with covariance $\mathbf G+\mathbf H^{-1}$. But the sampling truth, by Section 4, is $\mathbf G_{\text{true}}+\mathbf H^{-1}\mathbf J\mathbf H^{-1}$. Equating the two, the coupled estimator solves

$$\hat{\mathbf G}-\mathbf G_{\text{true}} =\mathbf E:=\mathbf H^{-1}(\mathbf J-\mathbf H)\mathbf H^{-1}. \tag{7.2}$$

It replaces the measurement variance by the curvature of the composite likelihood, which is precisely the object the second Bartlett identity would have licensed and which, in a composite likelihood, is wrong by exactly $\mathbf E$. When $\mathbf J=\mathbf H$ the bias vanishes identically; when they differ it is of the order of the quantity being estimated. Which component absorbs it depends on the pattern of $\mathbf E$ and not its size: white $\mathbf E$ loads entirely onto $\sigma^2_{\text{env}}$, $\mathbf E$ carrying the pedigree pattern loads entirely onto $\sigma^2_\delta$, and a numerical projection of (7.2) with the same mean diagonal sends a true $h^2_\delta=0.60$ to $0.375$ in the first case and $0.750$ in the second. The bias does not vanish under the null either: with $\sigma^2_\delta=0$ it is unchanged, so a non-heritable effect can be returned as heritable.

This is why the two-stage form above is preferred here. It injects the sandwich, which is valid under the pair dependence by construction; the coupled form manufactures the same quantity internally out of $\mathbf H$ and gets it wrong. The two-stage error is also visible, whereas the coupled one hides inside a well-behaved number. What the coupled form does genuinely repair is the ridge: with a penalty on $\boldsymbol{\delta}$ the first stage returns $\hat{\boldsymbol{\delta}}=\mathbf B\boldsymbol{\delta}+\boldsymbol{\epsilon}$, and $\mathbf B\mathbf G\mathbf B'$ cannot be written in the form $\sigma^2_\delta\mathbf A+\sigma^2_{\text{env}}\mathbf I$ at all, which is the mechanism behind $\sigma^2_{\text{env}}$ reaching the zero boundary. Coupling generates that shrinkage from the pedigree prior instead, where $\log|\mathbf H|$ already accounts for it.

The construction that would carry both advantages is the coupled form with a curvature adjustment applied to $c\ell$ before the Laplace step, in the manner of Chandler and Bate (2007) or Pace, Salvan and Sartori (2011): rescaling so the effective curvature becomes $\mathbf H\mathbf J^{-1}\mathbf H$ makes the variance the estimator replaces equal to the sandwich, which zeroes $\mathbf E$ by construction rather than calibrating around it. The obstacle is not cost but the rank of $\hat{\mathbf J}=\sum_g\mathbf S_g\mathbf S_g'$, at most the number of groups and exactly singular along the flat directions, so the trace-based adjustments, which never invert $\mathbf J$, are the safe form.

The general theory for what a two-step estimator costs when the first step is itself estimated is Murphy and Topel (1985); the standing warning against using first-stage predictions as a second-stage phenotype without carrying their uncertainty is Hadfield, Wilson, Garant, Sheldon and Kruuk (2010), and carrying $\mathbf V_e$ as above is what answers it. Where only the inference on $\sigma^2_\delta$ is wanted and not the point estimate, a cluster-robust variance in the second stage (Hedges, Tipton and Johnson, 2010), clustered by group as the first-stage sandwich already is, is the cheaper route.

Doing this to your own data

The data you need is one row per competitor per group: who competed, in which group, with how much exposure, and how many outcomes they took.

library(BreedingR)
set.seed(7)

# One row per competitor per group. Substitute your own columns.
n_grupos <- 200
contests <- do.call(rbind, lapply(seq_len(n_grupos), function(g) {
  quem <- sample(1:40, 4)
  dose <- runif(4, 0.5, 2)
  forca <- exp(c(0.9, 0.3, -0.2, -0.6, rep(0, 36))[quem])
  p <- dose * forca / sum(dose * forca)
  data.frame(group = as.character(g), competitor = as.character(quem),
             dose = dose, wins = as.vector(rmultinom(1, 12, p)))
}))
head(contests, 4)

##   group competitor      dose wins
## 1     1         19 1.0100935    1
## 2     1         31 1.9580938    4
## 3     1         28 0.7487832    3
## 4     1          2 1.1886555    4

competition_strength() estimates the strength of each competitor per unit of exposure. It is the fit of Section 1, and it reports who it could not measure.

bt <- competition_strength(contests$wins, contests$group,
                           contests$competitor, exposure = contests$dose)
head(round(sort(bt$log_strength, decreasing = TRUE), 3), 5)

##     1    20     2     9    25 
## 0.994 0.308 0.298 0.281 0.279

Competitors 1 and 2 were planted at 0.9 and 0.3 on this scale and come back near those values, which is the check worth running on any design of your own before trusting a ranking from it: simulate your own group sizes and exposures with strengths you chose, and see whether they come back.

Read identifiable before anything else. A competitor outside the strongly connected component of the contest graph has a strength that came from the prior and not from the data (Ford, 1957). It has no standard error, and it does not belong in a published ranking.

sum(!bt$identifiable)          # how many are prior shrinkage rather than measurement

## [1] 0

The strength then becomes a phenotype like any other, and it carries its own precision, so it goes in as a weight rather than as a bare number:

d2 <- data.frame(id = names(bt$strength), y = bt$log_strength, w = 1 / bt$se^2)
d2 <- d2[bt$identifiable, ]                 # measured competitors only
model(y ~ animal(id), d2, pedigree, weights = "w")

That two-stage route is the one Section 7 dissects, and the warning there applies here: weighting by 1/se^2 alone uses the composite curvature, which overstates the precision. Where the groups are large and the same competitors meet repeatedly, use the cluster-robust variance instead, and read the additive variance rather than the heritability as the finding.

What this design cannot answer, however much data you collect. If you are asking whether a competitor harms its rivals uniformly, Section 3 settles it: that effect is exactly aliased with its own ability, and no estimator separates the two. Only a rival-dependent effect is estimable, and only through the pair contrast.

competition_strength() fits the strengths by the MM algorithm of Hunter (2004), which is the EM-MAP of Caron and Doucet (2012) with a Gamma prior and contains Zermelo’s (1929) iteration as the flat-prior case. It flags competitors outside the strongly connected component of the contest graph, where the estimate is prior shrinkage rather than measurement (Ford, 1957). The variance-component estimator of Section 6 and the differential test of Sections 4 and 5 are analysis-level constructions built on the sparse pieces the package exposes (a_inverse(), sparse_chol(), sparse_solve(), selected_inverse()).

References

• Bradley, R.A. & Terry, M.E. (1952). Rank analysis of incomplete block designs I. Biometrika 39:324-345.
• Caron, F. & Doucet, A. (2012). Efficient Bayesian inference for generalized Bradley-Terry models. Journal of Computational and Graphical Statistics 21:174-196.
• Chandler, R.E. & Bate, S. (2007). Inference for clustered data using the independence loglikelihood. Biometrika 94:167-183.
• Ford, L.R. (1957). Solution of a ranking problem from binary comparisons. American Mathematical Monthly 64:28-33.
• Garrick, D.J., Taylor, J.F. & Fernando, R.L. (2009). Deregressing estimated breeding values and weighting information for genomic regression analyses. Genetics Selection Evolution 41:55.
• Godambe, V.P. (1960). An optimum property of regular maximum likelihood estimation. Annals of Mathematical Statistics 31:1208-1211.
• Hadfield, J.D., Wilson, A.J., Garant, D., Sheldon, B.C. & Kruuk, L.E.B. (2010). The misuse of BLUP in ecology and evolution. The American Naturalist 175:116-125.
• Hedges, L.V., Tipton, E. & Johnson, M.C. (2010). Robust variance estimation in meta-regression with dependent effect size estimates. Research Synthesis Methods 1:39-65.
• Hunter, D.R. (2004). MM algorithms for generalized Bradley-Terry models. Annals of Statistics 32:384-406.
• Kent, J.T. (1982). Robust properties of likelihood ratio tests. Biometrika 69:19-27.
• Lindsay, B.G. (1988). Composite likelihood methods. Contemporary Mathematics 80:221-239.
• Louis, T.A. (1984). Estimating a population of parameter values using Bayes and empirical Bayes methods. Journal of the American Statistical Association 79:393-398.
• Luce, R.D. (1959). Individual Choice Behavior. Wiley, New York.
• Murphy, K.M. & Topel, R.H. (1985). Estimation and inference in two-step econometric models. Journal of Business and Economic Statistics 3:370-379.
• Pace, L., Salvan, A. & Sartori, N. (2011). Adjusting composite likelihood ratio statistics. Statistica Sinica 21:129-148.
• Piepho, H.-P., Möhring, J., Schulz-Streeck, T. & Ogutu, J.O. (2012). A stage-wise approach for the analysis of multi-environment trials. Biometrical Journal 54:844-860.
• Plackett, R.L. (1975). The analysis of permutations. Applied Statistics 24:193-202.
• Thompson, S.G. & Sharp, S.J. (1999). Explaining heterogeneity in meta-analysis: a comparison of methods. Statistics in Medicine 18:2693-2708.
• Tierney, L. & Kadane, J.B. (1986). Accurate approximations for posterior moments and marginal densities. Journal of the American Statistical Association 81:82-86.
• Varin, C., Reid, N. & Firth, D. (2011). An overview of composite likelihood methods. Statistica Sinica 21:5-42.
• Zermelo, E. (1929). Die Berechnung der Turnier-Ergebnisse als ein Maximumproblem der Wahrscheinlichkeitsrechnung. Mathematische Zeitschrift 29:436-460.