1. Multivariate outcomes: meta-analysis with binomial data

1.1 Goal of the example

  • There are 22 clinical trials, each comparing two treatments (treatment 0 vs 1) with outcomes death vs survival.
  • For each trial $j$:
    • Treatment $i = 0,1$
    • Data: number of deaths $y_{ij}$ out of $n_{ij}$ patients
  • Interest: the distribution of treatment effects on the probability of death across studies.
  • Earlier analysis (in a previous section) used:
    • A normal approximation to the log-odds ratio in each study (univariate model).
    • That approximation is okay with large samples, but:
      • It ignores the average death rate level in each study.
      • That can matter if the baseline risk correlates with the treatment effect.
  • Now they build a more exact binomial–hierarchical model that treats both treatment arms jointly as a bivariate outcome per study.

1.2 Data model

  • For each treatment $i$ and study $j$:
    • $y_{ij} \mid n_{ij}, p_{ij} \sim \text{Bin}(n_{ij}, p_{ij})$
    • $p_{ij}$ is the probability of death for treatment $i$ in study $j$.
  • There are 44 parameters $p_{ij}$ (2 per study × 22 studies).
  • To model these, they move to the logit scale so that:
    • $\text{logit}(p_{ij})$ lies in $(-\infty, \infty)$
    • Normal distributions become reasonable as population models.

1.3 Reparameterization: average and effect per study

Instead of modeling $\text{logit}(p_{0j})$ and $\text{logit}(p_{1j})$ directly, they define for each study $j$:

  • Average logit level:
    • $ \beta_{1j} = \big(\text{logit}(p_{0j}) + \text{logit}(p_{1j})\big) / 2 $
  • Treatment effect (difference in logits):
    • $ \beta_{2j} = \text{logit}(p_{1j}) – \text{logit}(p_{0j}) $

Interpretation:

  • $\beta_{2j}$ is the log-odds ratio (treatment vs control) in study $j$.
  • $\beta_{1j}$ represents the average risk level (baseline severity) in that study.

So each study $j$ has a pair $(\beta_{1j}, \beta_{2j})$.

1.4 Hierarchical model for study-level parameters

  • The 22 study-specific pairs $(\beta_{1j}, \beta_{2j})$ are treated as exchangeable and modeled as draws from a bivariate normal distribution:
    • $(\beta_{1j}, \beta_{2j})^\top \sim \mathcal{N}\big((\alpha_1, \alpha_2)^\top,\ \Lambda\big)$
  • Here:
    • $\alpha_1$ = population mean of average logit (baseline death level)
    • $\alpha_2$ = population mean of treatment effect
    • $\Lambda$ = $2 \times 2$ covariance matrix with entries:
      • $\Lambda_{11}$: variance of $\beta_{1j}$ across studies
      • $\Lambda_{22}$: variance of $\beta_{2j}$ across studies
      • $\Lambda_{12}$: covariance between baseline level and treatment effect
    • Correlation: $\rho_{12} = \Lambda_{12} / (\Lambda_{11} \Lambda_{22})^{1/2}$

They note: this is equivalent to a normal model in $(\text{logit}(p_{0j}), \text{logit}(p_{1j}))$, just linearly reparameterized. The chosen transformation is convenient because it tends to make $\beta_{1}$ and $\beta_{2}$ closer to independent in their population distribution.

1.5 Hyperpriors

  • For the means:
    • $\alpha_1, \alpha_2$ get uniform (flat) priors.
  • For the covariance matrix $\Lambda$:
    • No universally “noninformative” prior exists.
    • Here they use:
      • Independent uniform priors on the variances $\Lambda_{11}$ and $\Lambda_{22}$.
      • Independent uniform prior on the correlation $\rho_{12}$.
  • This choice leads to a proper posterior (they mention an exercise that proves this).

1.6 Posterior computation

The model is moderately high-dimensional (parameters for all $\beta_{1j}$, $\beta_{2j}$, plus $\alpha$, plus $\Lambda$), so they use MCMC:

  • Parameterization used for computation:
    • $\beta$ (all study-specific parameters)
    • $\alpha_1, \alpha_2$
    • $\log(\Lambda_{11}), \log(\Lambda_{22})$
    • Fisher’s $z$ transform of the correlation:
      • $z = \tfrac{1}{2}\log\big((1 + \rho_{12}) / (1 – \rho_{12})\big)$
        to map $\rho_{12} \in (-1,1)$ into $z \in \mathbb{R}$.
  • Strategy:
    1. Construct an approximate posterior based on conditional modes.
    2. Draw initial points from this approximation and perform importance resampling to get starting values for ten parallel Metropolis chains.
    3. Use a multivariate normal jumping kernel:
      • Mean: posterior mode.
      • Covariance: proportional to the inverse Hessian (“curvature matrix”) at the mode.
      • Scaled by $2.4 / \sqrt{49}$ because the jump is in 49-dimensional space.
    4. Run chains for 40,000 iterations.
    5. Diagnose convergence using the scale reduction statistic $\hat{R}$:
      • All parameters had $\hat{R} < 1.2$; most $< 1.1$.
    6. Keep the second half of the simulations, resulting in 200,000 draws used for posterior summaries.

1.7 Posterior results

Table reports posterior quantiles for:

  • Individual study effects $\beta_{2j}$ (treatment effects),
  • One example of a baseline parameter $\beta_{1,1}$,
  • Hyperparameters:
    • $\alpha_1$ (mean of average logits),
    • $\sqrt{\Lambda_{11}}$ (sd of $\beta_{1j}$),
    • $\alpha_2$ (mean treatment effect),
    • $\sqrt{\Lambda_{22}}$ (sd of treatment effects),
    • $\rho_{12}$ (correlation between baseline and effect).

Key points:

  • The posterior of $\rho_{12}$ is centered around 0.21, but with substantial uncertainty.
  • This means: there is some evidence that baseline risk and treatment effect are positively associated, but not strongly.
  • Compared to the earlier univariate meta-analysis (which modeled only log-odds ratios), the multivariate model:
    • Fixes the approximation issue (no normal approximation to log-odds is needed).
    • But does not substantially change the posterior for:
      • Study-specific treatment effects,
      • Overall mean effect,
      • Between-study variance of effects.
  • Conclusion: the multivariate binomial analysis is more exact, but in this case, the simpler univariate normal-approximation model was already quite adequate for the main inferential targets.

2. Multinomial responses: extensions of logistic models

2.1 Extending logistic regression to multinomial outcomes

For a unit $i$ with $k$ possible outcome categories:

  • Data: $y_i = (y_{i1}, \dots, y_{ik})$ with $\sum_{j=1}^k y_{ij} = n_i$.
  • Model: $y_i \sim \text{Multin}(n_i; \alpha_{i1}, \dots, \alpha_{ik})$, with $\sum_{j=1}^k \alpha_{ij} = 1$.
  • A common parameterization:
    • Choose a baseline category, say $j = 1$.
    • Define for $j = 2, \dots, k$:
      • $\log(\alpha_{ij} / \alpha_{i1}) = \eta_{ij} = X_i \beta_{(j)}$
    • Here:
      • $\beta_{(j)}$ is a coefficient vector for category $j$.
      • Set $\beta_{(1)} = 0$ so that $\eta_{i1} = 0$.
  • Likelihood:
    • $p(y \mid \beta) \propto \prod_{i=1}^n \prod_{j=1}^k \left( \frac{e^{\eta_{ij}}}{\sum_{\ell=1}^k e^{\eta_{i\ell}}} \right)^{y_{ij}}$.
  • Interpretation:
    • $\beta_{(j)}$ describes how covariates $X_i$ shift the log odds of category $j$ relative to category 1.

2.2 Ordered vs nominal categories

  • Nominal categories: no natural order (e.g., disease types).
  • Ordinal categories: ordered (e.g., grades A, B, C, D).

For ordinal outcomes:

  • Often use cumulative probabilities:
    • $\pi_{ij} = \sum_{\ell \le j} \alpha_{i\ell}$ (probability of being in category ≤ j).
    • Model: $\log(\pi_{ij} / (1 – \pi_{ij})) = \delta_j + X_i \beta_{(j)}$.
  • Frequently assume common slopes:
    • $\beta_{(j)} = \beta$ for all $j$.
  • Another option: multinomial probit, which also can be formulated with latent variables.

These models can often be represented in latent-variable form, e.g., with a continuous latent utility or propensity and cutpoints defining category boundaries.

3. Using Poisson models for multinomial data

3.1 Poisson–multinomial relationship

  • If $y_j \sim \text{Poisson}(\lambda_j)$ independently for $j = 1,\dots,k$, and $n = \sum_j y_j$:
    • The conditional distribution of $y$ given $n$ is:
      • $y \mid n, \alpha \sim \text{Multin}(n; \alpha_1,\dots,\alpha_k)$
      • where $\alpha_j = \lambda_j / \sum_i \lambda_i$.
  • This allows one to fit multinomial models using Poisson GLMs, by:
    • Modeling $\lambda_j$ via a log-link,
    • Conditioning on the total.

To ensure the probabilities sum to 1, one can add extra parameters (e.g., $A_{ij}$) and constraints.

4. Example: World Cup chess (paired comparisons with ties)

4.1 Data description

  • 1988–1989 World Cup of chess:
    • Six round-robin tournaments.
    • 29 players total.
    • Each tournament: each player plays each other exactly once.
    • Total of 789 games.
  • For each game:
    • Players $i$ and $j$,
    • Outcome: win / loss / draw,
    • Who had the first move (white), which is known to give an advantage.
  • Table shows a subset for 8 players (Karpov, Kasparov, Korchnoi, Ljubojevic, Seirawan, Short, Spassky, Tal), aggregated over tournaments.
    • Entries like “1-0-1” mean: 1 win, 0 losses, 1 draw for the row player (with white) vs the column player.

4.2 Basic Bradley–Terry model (without draws)

  • Standard paired-comparisons model:
    • Ability parameters: $\alpha_i$ for each player $i$.
    • Probability that $i$ beats $j$:
      • $p_{ij} = \dfrac{\exp(\alpha_i – \alpha_j)}{1 + \exp(\alpha_i – \alpha_j)}$
  • This model:
    • Ignores draws,
    • Ignores first-move advantage.

4.3 Extended model with draws and first-move advantage

They extend the model for games where player $i$ has the first move, using parameters:

  • $\alpha_i$: ability of player $i$,
  • $\alpha_j$: ability of player $j$,
  • $\gamma$: effect of moving first,
  • $\delta$: propensity for draws,
  • $\theta$ will denote the overall parameter vector (but each equation is written in terms of components).

Define probabilities:

  • $p_{ij1} = \Pr(i \text{ defeats } j \mid \theta)$
  • $p_{ij2} = \Pr(j \text{ defeats } i \mid \theta)$
  • $p_{ij3} = \Pr(i \text{ draws with } j \mid \theta)$

and model them as:

  • $p_{ij1} = \dfrac{e^{\alpha_i}}{e^{\alpha_i} + e^{\alpha_j + \gamma} + e^{\delta + \tfrac{1}{2}(\alpha_i + \alpha_j + \gamma)}}$
  • $p_{ij2} = \dfrac{e^{\alpha_j + \gamma}}{e^{\alpha_i} + e^{\alpha_j + \gamma} + e^{\delta + \tfrac{1}{2}(\alpha_i + \alpha_j + \gamma)}}$
  • $p_{ij3} = \dfrac{e^{\delta + \tfrac{1}{2}(\alpha_i + \alpha_j + \gamma)}}{e^{\alpha_i} + e^{\alpha_j + \gamma} + e^{\delta + \tfrac{1}{2}(\alpha_i + \alpha_j + \gamma)}}$

Interpretation:

  • $\gamma > 0$ suggests a systematic advantage to moving first.
  • $\delta$ controls how likely draws are, relative to wins and losses.
  • The symmetric term $\tfrac{1}{2}(\alpha_i + \alpha_j + \gamma)$ in the draw exponent reflects that a draw is “in between” the two player strengths.

4.4 Poisson GLM formulation

Let:

  • $y_{ijk}$ be the number of times outcome $k$ ($1 = i$ wins, $2 = j$ wins, $3 = draw$) occurs when $i$ plays $j$ with $i$ having the first move.
  • $n_{ij}$ be the total number of such games.

They assume:

  • $y_{ijk}$ are independent Poisson random variables with means:
    • $\mu_{ijk} = n_{ij} p_{ijk}$.

Then they write a log-link linear predictor for each mean:

  • $\log \mu_{ij1} = \log n_{ij} + \alpha_i – A_{ij}$
  • $\log \mu_{ij2} = \log n_{ij} + \alpha_j + \gamma – A_{ij}$
  • $\log \mu_{ij3} = \log n_{ij} + \delta + \tfrac{1}{2}\alpha_i + \tfrac{1}{2}\alpha_j + \tfrac{1}{2}\gamma – A_{ij}$

where $A_{ij}$ is the log of the common denominator of the three probabilities.

The role of $A_{ij}$:

  • Ensures that $\mu_{ij1} + \mu_{ij2} + \mu_{ij3} = n_{ij}$, so the Poisson model corresponds exactly to the multinomial distribution with probabilities $p_{ij1}, p_{ij2}, p_{ij3}$.
  • $A_{ij}$ is treated as a nuisance parameter to enforce normalization.

4.5 Structure of the Poisson regression (X and β)

  • Outcome vector $y$:
    • Length = $3 \times 29 \times 28$ (three outcomes for each ordered pair $(i,j)$ with $i \ne j$).
  • Mean vector $\mu$:
    • Same length, with components $(\mu_{ij1}, \mu_{ij2}, \mu_{ij3})$.
  • Link: $\log \mu = X \beta$, where:
    • Parameter vector $\beta$ includes:
      • 29 ability parameters $\alpha_i$,
      • $\gamma$ (first-move advantage),
      • $\delta$ (draw behaviour),
      • 29 × 28 nuisance parameters $A_{ij}$.
    • Design matrix $X$ is built from the structures in the log-mean equations:
      • Column for $\alpha_1$:
        • Value 1 in rows for games where player 1 wins,
        • Value 0.5 in rows for draws involving player 1,
        • 0 elsewhere.
      • Column for $\delta$:
        • 1 in rows representing draws,
        • 0 elsewhere.
      • Final block of columns for each $A_{ij}$:
        • 1 in the three rows corresponding to the outcomes of games where $i$ has first move vs $j$,
        • 0 elsewhere.

During posterior simulation they do not care about $A_{ij}$ themselves; they are just needed to make the Poisson model exactly equivalent to the multinomial specification.

4.6 Using an offset and identifiability

  • $\log n_{ij}$ is included as an offset (a predictor with fixed coefficient 1) in the Poisson GLM.
  • With a flat prior on all parameters, the model is overparameterized:
    • Adding a constant $c$ to all $\alpha_i$ and adjusting $A_{ij}$ accordingly leaves the probabilities unchanged.
  • To identify the model, they fix:
    • $\alpha_1 = 0$ (choose player 1 as reference),
    • One of the $A_{ij}$ also set to 0.
  • If one later extends the model to allow abilities to vary (hierarchical ability model), the constraint on $\alpha_1$ may no longer be necessary if identifiability is handled differently.