A hierarchical model can involve more than one coefficient that varies across groups.
Whenever several coefficients differ by group, a multivariate prior distribution becomes necessary.
1. Basic hierarchical structure with multiple varying coefficients
Assume the existence of $J$ groups.
Within each group $j$, the data vector $y^{(j)}$ is modeled by a likelihood of the form:
$p(y^{(j)} \mid \theta^{(j)})$
The parameter vector $\theta$ contains multiple coefficients that vary from group to group.
A prior distribution is assigned to each group-specific parameter vector:
$p(\theta_j \mid \phi)$
where $\phi$ represents hyperparameters.
This structure supports any model involving group-dependent intercepts, group-dependent slopes, or both.
2. Example: Varying-intercept, varying-slope regression
A common form is a linear regression where each group has its own intercept and its own slope:
$y_{ij} \sim \text{Normal}(\alpha_j + x_{ij}\,\beta_j, \sigma_y^2)$
To describe how intercepts and slopes vary across groups, a bivariate normal prior is used:
$\begin{pmatrix} \alpha_j \\ \beta_j \end{pmatrix} \sim \text{Normal} \left( \begin{pmatrix} \mu_\alpha \\ \mu_\beta \end{pmatrix}, \begin{pmatrix} \sigma_\alpha^2 & \rho\,\sigma_\alpha\sigma_\beta \\ \rho\,\sigma_\alpha\sigma_\beta & \sigma_\beta^2 \end{pmatrix} \right)$
This prior specifies:
- a mean intercept $\mu_\alpha$
- a mean slope $\mu_\beta$
- two variances $\sigma_\alpha^2$, $\sigma_\beta^2$
- a correlation $\rho$ between intercepts and slopes
Hyperprior distributions are required for the parameters
$(\mu_\alpha, \mu_\beta, \sigma_\alpha, \sigma_\beta, \rho)$.
Typical choices start with uniform distributions (with restrictions such as positive scales and a correlation between −1 and 1).
Additional information can be introduced when the number of groups is small.
3. Extending to higher-dimensional varying coefficients
When more than two coefficients vary by group, the model takes the following multivariate form:
$y_{ij} \sim \text{Normal}(X_j \,\beta^{(j)}, \sigma_y^2)$
$\beta^{(j)} \sim \text{Normal}(\mu_\beta, \Sigma_\beta)$ (15.5)
- $X_j$: design matrix for group $j$
- $\beta^{(j)}$: $K$-dimensional coefficient vector for group $j$
- $\mu_\beta$: population-level mean vector
- $\Sigma_\beta$: $K \times K$ group-level covariance matrix
Additional layers can be added—such as group-dependent variances, structured mean vectors, or multiple grouping factors—but the crucial component is the prior on $\Sigma_\beta$.
The mean vector $\mu_\beta$ and the data variance $\sigma_y^2$ also require hyperpriors, but these quantities are often well estimated, making the exact choice of prior less critical.
4. Inverse-Wishart prior for $\Sigma_\beta$
Let $K$ be the number of varying coefficients.
The full matrix $\Sigma_\beta$ is $K \times K$.
A conditionally conjugate choice for $\Sigma_\beta$ is the inverse-Wishart distribution:
$\Sigma_\beta \sim \text{Inv-Wishart}_{K+1}(I)$
Although this prior is convenient algebraically, it places strong restrictions on the diagonal elements (the variances).
Such restrictions make the model less flexible for representing large uncertainty in those variances.
5. Scaled inverse-Wishart approach for more flexibility
A more flexible parameterization rewrites (15.5) using a decomposition:
$\beta^{(j)} = \mu_\beta + \xi \otimes \eta^{(j)}$
where:
- $\xi$: a $K$-dimensional scale vector
- $\eta^{(j)}$: a $K$-dimensional standardized random vector
- $\otimes$: element-wise multiplication
The new model becomes:
$y_{ij} \sim \text{Normal}(X_j(\mu_\beta + \xi \otimes \eta^{(j)}), \sigma_y^2)$
$\eta^{(j)} \sim \text{Normal}(0, \Sigma_\eta)$
From this decomposition:
$\Sigma_\beta = \text{Diag}(\xi)\, \Sigma_\eta\, \text{Diag}(\xi)$
Advantages of this reparameterization include:
- separate prior distributions for $\xi$ and $\Sigma_\eta$
- rich representations of uncertainty
- weak prior constraints on variances
A noninformative version sets:
- $\Sigma_\eta \sim \text{Inv-Wishart}_{K+1}(I)$
- $\xi$ governed by independent uniform priors
- informative priors can later be added to $\xi$ if needed
6. Business-school example (prediction problem with sparsity)
A practical application involves predicting first-year grades in U.S. business schools.
Grade predictions often use a regression based on:
- GMAT verbal score (GMAT-V)
- GMAT quantitative score (GMAT-Q)
- undergraduate GPA (UGPA)
- an intercept term
These predictions are used for admission decisions, so accuracy is crucial.
Concern about subgroup bias
A study examined whether prediction errors differed between black students and non-black students.
Data came from:
- 59 business schools
- approximately 8500 students over two years
- roughly 4% black students
Separate regressions were fit for each school:
$y_{ij} = \alpha_j + \beta_{j,\,V}(\text{GMAT-V})_{ij} + \beta_{j,\,Q}(\text{GMAT-Q})_{ij} + \beta_{j,\,G}(\text{UGPA})_{ij} + \epsilon_{ij}$
Residual analysis indicated systematic overprediction for black students.
Expanded regression to allow subgroup differences
A single binary indicator for race was considered too restrictive.
A more flexible formulation doubled the predictors:
- original four variables
- the same four multiplied by the race indicator
Thus, each school required eight coefficients:
- four for predicting non-black students
- four additional terms representing adjustments for black students
Total parameters per school
For school $j$:
- 888 regression coefficients
- 111 residual standard deviation $\sigma_j$
Total: 9 parameters per school.
Challenge: sparse data for black students
Some schools had no black students or only a handful.
Classical regression with noninformative priors cannot estimate nine parameters per school under such sparsity.
A hierarchical model resolves this issue by pooling information across schools.
Hierarchical solution
A multivariate prior is assigned:
$\beta_j \sim \text{Normal}(\alpha, \Lambda_\beta)$
for each of the 59 schools.
Unknown quantities include:
- all coefficients $\beta_j$
- the population mean vector $\alpha$
- the population covariance $\Lambda_\beta$
- 59 variance parameters $\sigma_1, \ldots, \sigma_{59}$
A uniform prior is placed on $\alpha$, $\Lambda_\beta$, and the log-variances.
Initial values for $\alpha$ and the variances were produced by pooling the data across schools.
An EM algorithm was then used to locate the posterior mode.
Conclusions of the study
The hierarchical multivariate regression produced accurate predictions for both black and non-black students, outperforming the standard non-hierarchical regression.
A systematic difference was discovered:
- Conditioning on the mean test scores for black students, non-black students were predicted to have higher first-year GPAs in approximately 85% of schools.
- More than 60% of these differences exceeded one posterior standard deviation.
- About 20% exceeded two posterior standard deviations.
Such differences cannot be estimated reliably without a hierarchical model because many schools have extremely small numbers of black students.