Models for multivariate and multinomial responses

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.