3.1  Introduction

In this chapter, we review several basic features of the multivariate normal distribution. Section 3.2 considers general properties of the multivariate normal density and Section 3.3 considers the sampling distribution of the maximum likelihood estimators (MLEs) of the mean vector and variance–covariance matrix based on iid (independent, identically distributed) random sampling of a multivariate normal distribution. Section 3.4 reviews the standard conjugate distributions for Bayesian inference with the multivariate normal distribution and Section 3.5 considers various generalizations and robustifications of the normal model. The properties of the multivariate normal distribution are well known and available in many places; our primary sources are the texts by Johnson and Wichern (1998) and Morrison (2005). A classic and comprehensive treatment is given by Anderson’s (2003) text.

In item response theory (IRT), the multivariate normal and its generalizations are most often used as the underlying variables distribution for a data-augmentation version of the normal-ogive model (Bartholomew and Knott, 1999; Fox, 2010), as a population distribution for the proficiency parameters θ _p , for persons p = 1, ..., P, and as a prior distribution for other model parameters (e.g., difficulty parameters b_i , log-discrimination parameters log a_i , etc.) whose domain is the entire real line or Euclidean space. Especially, in the latter two cases, the multivariate normal distribution serves to link standard IRT modeling with hierarchical linear models (HLMs) and other structures that incorporate various group structures and other dependence on covariates, to better model item responses in terms of the contexts in which they are situated (e.g., Fox, 2003, 2005a,b, 2010).

Because of the wide variety of applications of the normal distribution in IRT, we will depart somewhat from the notation used in the rest of the book in this chapter. We use X = (X ₁, X ₂, ..., X_K ) ^T to represent a generic K-dimensional random vector and X to represent a generic random variable. Their observed values will be denoted x and x, respectively. Also, because very many densities will be discussed, we will use the generic notation f () to denote a density. The particular role or nature of f () will be clear from context.

3.2  Multivariate Normal Density

The univariate normal density with mean μ and variance σ² for a random variable X is as follows:

As it is well known, E[X] = μ, and Var(X) = E[(X − μ)²] = σ². We often write X ∼ N (μ, σ²) to convey that the random variable X has this density. The quantity in the exponent of the normal density thus measures the square of the distance between realizations of the random variable, X, and the mean μ, scaled in units of the standard deviation σ.

The multivariate normal density, for X ∈ ℜ ^K , is as follows:

where again E[X] = μ = (μ₁, ..., μ _K ) ^T is the mean vector and Var(X) = E[(X − μ) (X − μ) ^T ] = Σ is the K × K symmetric nonnegative-definite variance–covariance matrix. We often write X ∼ N_K (μ, Σ), or just X ∼ N(μ, Σ) if the dimension K is clear from context, to indicate that X follows the multivariate normal density. The quantity (x − μ) ^T ∑ ⁻¹(x − μ) in the exponent is the squared distance from x to μ, again scaled by the variance–covariance matrix Σ. In other contexts, this is called the Mahalanobis distance between x and μ (Morrison, 2005).

The following theorem gives some properties of the multivariate normal random variables.

3.3  Sampling from a Multivariate Normal Distribution

3.4  Conjugate Families

Recall that if the likelihood for data X is (any function proportional to) the conditional density of X given (possibly multidimensional) parameter η, f(x|η), then a conjugate family of prior distributions for f(x|η) is a parametric family of densities f(η; τ) for η with (possibly multidimensional) hyperparameter τ, such that for any member of the conjugate family, the posterior distribution f(η|x; τ) = f(x|η) f(η; τ)/ ∫_h f(x|h) f(h; τ) d h can be rewritten as f(η|τ ^∗), a member of the same parametric family as f(η|τ), where τ ^∗ = τ ^∗(x, τ) is some function of x and τ. Since only the form of f(η|τ ^∗) as a function of η matters, in verifying that f(η|τ ^∗) and f(η; τ) belong to the same parametric family, it is usual to ignore multiplicative constants that do not depend on η.

For example, if X ₁, ..., X_N are iid N(μ, σ²), then the likelihood is

as a function of μ, as would be expected, since $\overline{x}$

is sufficient for μ.

If we assume σ² is known and place a normal $f(\mu ;{\mu }_0,{\tau }_0^2)=(1/\sqrt{2\pi }{\tau }_0)e^{-(1/2{\tau }_0^2){(\mu -{\mu }_0)}^2}$

prior on μ, then the posterior density for μ will be

after completing the square, collecting terms, and identifying the normalizing constant, where

In this case, the posterior mean is ${\mu }_N={\rho }_N\overline{x}+(1-{\rho }_N){\mu }_0$

where ${\rho }_N={\tau }_0^2/({\tau }_0^2+{\sigma }^2/n)$ is the classical reliability coefficient. Thus, if the prior distribution is μ ∼ N(μ₀, τ₀), then the posterior distribution will be μ|x ₁, ..., x_N ∼ N(μ _N , τ _N ); this shows that the normal distribution is the conjugate prior for a normal mean μ, when the variance σ² is known.

Further calculation (as shown in Gelman et al., 2004) shows that when both the mean μ and the variance σ² are unknown, and the data X ₁, ..., X_N are an iid sample from N(μ, σ²), the joint distribution

is the conjugate prior, where the notation “ ${\sigma }^2\sim \text{Inv}-{\chi }^2({\nu }_0,{\sigma }_0^2)$

” means that ${\sigma }_0^2/{\sigma }^2\sim {\chi }_{\nu 0}^2$ . Here, κ₀ and ν₀ are hyperparmeters that function as “prior sample sizes”—the larger κ₀ and ν₀, the greater the influence of the prior on the posterior. In this case, the joint posterior distribution for μ, σ² is

where

with $S^2=(1/N-1){\mathrm{\Sigma }}_{n=1}^N{(x_n-\overline{x})}^2$

.

Although this is the conjugate family when μ and σ² are both unknown, the forced dependence between μ and σ² in the prior is often awkward in applications. Common alternatives are as follows:

1. To replace σ²/κ₀ with an arbitrary τ² in the conditional prior for μ|σ², forcing independence. In this case, the conditional posterior for μ|σ² is as in the “σ² known” case above, but the marginal posterior for σ² is neither conjugate nor in closed form (it is not difficult to calculate however) or
2. To make the conditional prior for μ|σ² very flat/uninformative by letting κ₀ → 0. In this case, the posterior distributions for μ|σ² and for σ² mimic the sampling distributions of the MLEs

If one wishes for a noninformative prior for σ² (analogous to the flat prior choice for μ in (b) above), a choice that preserves conjugacy would be to take the degrees of freedom ν₀ = 0 in the $\text{Inv}-{\chi }^2({\nu }_0,{\sigma }_0^2)$

prior for σ². This leads to a prior f(σ²) ∝ 1/σ². Another noninformative choice is the Jeffreys prior for σ² (proportional to the square root of the Fisher information for the parameter), which in this case leads to the prior f(σ) ∝ 1/σ. The Jeffreys prior for log σ² (or equivalently log σ)issimply f(σ²) ∝ 1. All of these choices are improper priors and care must be taken that the posterior turns out to be proper. One common way to avoid this issue is to force the prior to be proper, say, by taking σ² ∼ Unif(0, M) for some suitably large number M.

The conjugate prior distribution for a multivariate normal distribution with parameters μ and Σ has a form similar to that of the univariate case, but with multivariate normal and inverse-Wishart densities replacing the univariate normal and inverse-χ-squared densities. In particular, assuming sampling x ₁, ..., x _N iid from N(μ, Σ), the joint prior for μ and Σ is of the following form:

where the notation “ Σ ∼ Inv-Wishart $({\nu }_0,{\mathrm{\Sigma }}_0^{-1})$

” means that ${\mathbf{\Sigma }}_0{\mathbf{\Sigma }}^{-1}\sim {\omega }_{{\nu }_0}\left({\mathbf{\Sigma }}_0\right)$ the usual Wishart distribution with ν₀ degrees of freedom and parameter matrix Σ ₀. Again, κ₀ and ν₀ function as prior sample sizes. Then, the joint posterior distribution will be as follows:

where

with $\mathbf{S}=(1/N-1){\mathrm{\Sigma }}_{n=1}^N({\mathbf{x}}_n-\overline{\mathbf{x}}){({\mathbf{x}}_n-\overline{\mathbf{x}})}^T$

(Gelman et al., 2004). Once again the joint conjugate prior forces some prior dependence between μ and Σ that may be awkward in practice. And once again the common fixes are as follows:

1. Replace Σ/κ₀ with an arbitrary Λ ₀ in the conditional prior for μ|Σ, forcing independence or
2. Make the conditional prior for μ|Σ very flat/uninformative by letting κ₀ → 0

The details, generally analogous to the univariate normal case, are worked out in many places; in particular, see Gelman et al. (2004). Gelman et al. (2004) suggest another way to reduce the prior dependence between μ and ∑ . Sun and Berger (2007) provide an extensive discussion of several “default” objective/noninformative prior choices for the multivariate normal.

3.5  Generalizations of the Multivariate Normal Distribution

Assessing multivariate normality in higher dimensions is challenging. A well-known theorem (Johnson and Wichern, 1998; see also Anderson, 2003) states that X is a multivariate normal vector if and only if each linear combination a ^T X of its components is univariate normal, but this is seldom feasible to show in practice. Instead one often only checks the one- and two-dimensional margins of X; for example, examining a Q—Q plot for each of the K components of X, and in addition, examining bivariate scatterplots of each possible pair of variables to determine whether the data points yield an elliptical appearance. (See Johnson and Wichern (1998) for more information on techniques for assessing multivariate normality.)

There is no doubt that true multivariate normality is a rare property for a multivariate dataset. Although the latent ability variable in IRT is often assumed to be normally distributed, the data may not conform well to this assumption at all (Casabianca, 2011; Casabianca et al., 2010; Moran and Dresher, 2007; Woods and Lin, 2009; Woods and Thissen, 2006). In these and other cases, robust alternatives to the multivariate normal distribution can be considered.

An important class of generalizations of the multivariate normal distribution is the family of multivariate elliptical distributions, so named because each of its level sets defines an ellipsoid (Fang and Zhang, 1990; Branco and Dey, 2002). The K-dimensional random variable X has an elliptical distribution if and only if its characteristic function (Billingsley, 1995; Lukacs, 1970) is of the form

where as usual μ is a K-dimensional vector and Σ is a symmetric nonnegative definite K × K matrix, and t ∈ ℜ ^K . When a density f(x) exists for X, it has the form

where g() is itself a univariate density; g() is called the generator density for f(). The density defines a location/scale family with location parameter μ and scale parameter Σ. The parameter μ is the median of f() in all cases, and if E[X] exists, E[X] = μ. If Var(X) exists, Var(X) = −(∂ψ(0)/∂t) · Σ. A little calculation shows that if W = AX + b is elliptically distributed with location μ and scale Σ, then W = AX + b is again elliptically distributed, with location Aμ + b and scale AΣA ^T .

Elliptical distributions are used as a tool for generalizing normal-theory structural equations modeling (e.g., Shapiro and Browne, 1987; Schumacker and Cheevatanarak, 2000). The special case of the K-dimensional multivariate-t distribution on ν degrees of freedom with density

has been used as the error distribution in robust Bayesian and non-Bayesian linear regression modeling since at least Zellner (1976); more recently, robust regression with general elliptical error distributions and the particular case of scale mixtures of normals (of which the univariate-t and multivariate-t are also examples) has also been studied (e.g., Fernandez and Steel, 2000).

A further generalization is the family of skew-elliptical distributions (e.g., Branco and Dey, 2001, 2002). When it exists, the density of a skew-elliptical distribution is of the form

where f _g1() is the density of a multivariate elliptical distribution, and F _g2() is the cumulative distribution function of a (possibly different) univariate elliptical distribution with location parameter 0 and scale parameter 1. The vector parameter λ is a skewness parameter; when λ = 0, the skew-elliptical density reduces to a symmetric elliptical density.

When the generator densities g ₁(x) = g ₂(x) = ϕ(x), the standard normal density, we obtain the special case of the skew-normal distributions. It has been observed (Moran and Dresher, 2007) that the empirical distribution of the latent proficiency variable in IRT models applied to large-scale educational surveys sometimes exhibits some nontrivial skewing, which if unmodeled can cause bias in estimating item parameters and features of the proficiency distribution. Skew-normal and related distributions have been proposed (Xu and Jia, 2011) as a way of accounting for this skewing in the modeling of such data, with as few extra parameters as possible.

Two additional classes of transformed normal distributions used in IRT are the lognormal and logit-normal distributions. A random variable X has a lognormal distribution when its logarithm is normally distributed. The density of the lognormal distribution is of the following form:

where μ and σ² are the mean and variance of the variable’s natural log. This distribution is an alternative to Gamma and Weibull distributions for nonnegative continuous random variables and is used in psychometrics to model examinee response times (e.g., van der Linden, 2006) and latent processes in decision making (e.g., Rouder et al., 2014). It is also used as a prior distribution in Bayesian estimation of nonnegative parameters, for example the discrimination parameter in two- and three-parameter logistic IRT models (van der Linden, 2006).

Similarly, a random variable X has a logit-normal distribution when its logit is normally distributed. The density of the logit-normal distribution is of the following form:

Here, μ and σ² are the mean and variance of the variable’s logit, and x is a proportion, bounded by 0 and 1. A multivariate generalization of the logit-normal distribution (Aitchison, 1985) has been used in latent Dirichlet allocation models for text classification ( Blei and Lafferty, 2007) and in mixed membership models for strategy choice in cognitive diagnosis (Galyardt, 2012).