# Bayesian Inference of a Multivariate Regression Model

Authors: Marick S. Sinay; John S. J. Hsu
Journal: Journal of Probability and Statistics
 (2014)

---

## Abstract

We explore Bayesian inference of a multivariate linear regression model with use of a flexible prior for the covariance structure. The commonly adopted Bayesian setup involves the conjugate prior, multivariate normal distribution for the regression coefficients and inverse Wishart specification for the covariance matrix. Here we depart from this approach and propose a novel Bayesian estimator for the covariance. A multivariate normal prior for the unique elements of the matrix logarithm of the covariance matrix is considered. Such structure allows for a richer class of prior distributions for the covariance, with respect to strength of beliefs in prior location hyperparameters, as well as the added ability, to model potential correlation amongst the covariance structure. The posterior moments of all relevant parameters of interest are calculated based upon numerical results via a Markov chain Monte Carlo procedure. The Metropolis-Hastings-within-Gibbs algorithm is invoked to account for the construction of a proposal density that closely matches the shape of the target posterior distribution. As an application of the proposed technique, we investigate a multiple regression based upon the 1980 High School and Beyond Survey.


---

## Body

## 1. Introduction

The multivariate multiple regression model is a natural extension of the univariate multiple regression model. The key difference, as the name implies, is that the univariate response variable is instead a multivariate response vector. By utilizing the multivariate multiple regression model the covariance of the response vector can be modeled. From an estimation standpoint van der Merwe and Zidek [1] suggest an intrinsic role to be played by the covariance structure, whereas, in the case of separate univariate multiple regression models, the covariance of the distinct response variables cannot be modeled. Although optimal point estimates of any linear combination of the means of the various response variables can still be obtained, an appropriate estimate of the variance of said estimator cannot be obtained without fully incorporating the covariance amongst the multivariate response vector. Under this framework multivariate analysis is required to most appropriately produce an estimate of the standard error.As a particular example, in educational testing data when multiple subject area exams are administered it is common practice to simply report the sum of the individual exam scores as a total score. For instance, the metric most commonly associated with the Scholastic Aptitude Test (SAT) is simply the sum of the student’s verbal score and the mathematical score. Other exams, such as the ACT exam, have even more than two subject areas and report a composite score which is the arithmetic average of the individual subject area scores. In these instances, multivariate analysis, by capturing the covariance amongst the various subject exams, is required to properly estimate the standard error of the final score.Formal Bayesian analysis has long been used for multivariate multiple regression models [2]. The inverse Wishart is widely used in this respect, since it is a conjugate prior distribution for the multivariate normal covariance matrix [3–5]. Also see Dawid [6] for a general discussion of the inverse Wishart and Wishart distributions. However, in contrast to traditional Bayesian methods we will not make use of the standard inverse Wishart conjugate prior for the covariance matrix. The reason is that the inverse Wishart is a rather restrictive distribution in its ability to capture prior information that may be available to the practitioner. See Leonard and Hsu [7], Hsu et al. [8], and Sinay et al. [9] for a more detailed explanation of the disadvantages of the inverse Wishart as a prior distribution for the covariance matrix.Leonard and Hsu [7] presented an alternative approach that remedies the shortcomings of the inverse Wishart and allows for greater flexibility in the prior specification. In a univariate normal model setting, the normal distribution has been used as a prior for the logarithm of the variance parameter. In this same vein, Leonard and Hsu [7] consider the matrix logarithm transformation of the covariance matrix for the multivariate case. Making use of a result from Bellman [10, page 171], it can be demonstrated that the exponential terms of a multivariate normal likelihood function can be expressed in the form of a Volterra linear integral equation. An approximation to a function, that is, proportional to the likelihood, can then be obtained via Bellman’s iterative solution to the Volterra integral equation. The resulting approximation has a multivariate normal form, with respect to the unique elements of the matrix logarithm of the covariance matrix. This allows a normal prior specification to act as a conjugate prior distribution, thereby yielding an approximate normal posterior for the covariance structure.One of the primary benefits of such a technique is the ability to specify varying degrees of confidence in each element of the normal prior hyperparameter mean vector via the variance terms of the prior hyperparameter covariance matrix. Obviously, larger variance terms in the prior hyperparameter covariance matrix indicate a lack of confidence in the corresponding prior location hyperparameter. Another chief advantage of this method is the ability to model beliefs about any possible interdependency between the covariance parameters. This can be accomplished by specifying the covariance terms of the prior hyperparameter covariance matrix. Note that in this way both the interrelationships and the strength of prior beliefs with respect to the covariance parameters can be modeled.Bayesian estimates of all relevant parameters of interest are calculated using Markov chain Monte Carlo (MCMC) techniques. With respect to the covariance structure, since an approximate posterior distribution is used as a proposal density, we employ the Metropolis-Hastings-within-Gibbs (MHWG) algorithm [11, page 291], to correctly estimate the true target posterior density.Having laid out the general outline we now move on to the body of the paper. We begin by introducing and defining the standard multivariate multiple regression model. We follow that by making a distributional assumption about the error matrix of the multivariate regression model. This provides us with a mechanism to state the likelihood function for the model. In turn, we then go through the formal analytical Bayesian derivations. Subsequent to the Bayesian analysis we outline the MCMC procedure and discuss how the posterior means and standard errors are numerically calculated. We conclude with an application to the High School and Beyond Survey [12].


## 2. Multivariate Multiple Regression Model

We consider the standard multivariate multiple regression model:(1)Y11⋯Y1p⋮⋱⋮Yn1⋯Ynpn×p=x11⋯x1k⋮⋱⋮xn1⋯xnkn×kβ11⋯β1p⋮⋱⋮βk1⋯βkpk×p+ϵ11⋯ϵ1p⋮⋱⋮ϵn1⋯ϵnpn×p.
Notationally in matrix form we can succinctly write
(2)Y=Xβ+ϵ,
where Y is the (n×p) matrix of response variables, X is the (n×k) matrix of explanatory variables, β is the (k×p) matrix of unknown regression coefficients, and ϵ is the (n×p) matrix of random errors. In Section 2.1 we introduce the matrix normal distribution and make a general assumption about the distribution of the n×p random error matrix ϵ in (2). The matrix normal representation will greatly facilitate the Baysian analysis that follows. Based upon the matrix normal, we proceed in Section 2.2 to develop the joint likelihood function for β and Σ. Hierarchical prior specifications are discussed in Section 2.3 and the joint posterior distribution is reported in Section 2.4.


### 2.1. Distributional Assumptions and Matrix Normal Distribution

The matrix normal distribution is closely related to, and is a generalization of the multivariate normal. In particular, then×p random matrix M~MNn×pΦ,Σ,Ω, if and only if, the np×1 random vector Vec⁡M~NnpVec⁡Φ,Σ⊗Ω, where MNn×p denotes the n×p dimensional matrix normal distribution, Φ is a n×p location matrix, Σ is a p×pfirst covariance matrix, and Ω is a n×nsecond covariance matrix [13, page 54]. Vec⁡· and ⊗ are the standard vector operator and Kronecker product, respectively.We make the distributional assumption that, conditional on thep×p covariance matrix Σ, the n×p random error matrix ϵ=(ϵ1T,ϵ2T,…,ϵnT)T, in (2), follows a matrix normal with n×p zero mean matrix and covariance matrices given by Σ and In, where In is a n×n identity matrix. Formally, we have ϵ∣Σ~MNn×p0,Σ,In, or equivalently, the error terms ϵ1,ϵ2,…,ϵn are independent and identically distributed normal random vectors each with mean vector 0 and covariance matrix Σ. The probability density function for the error matrix is given by
(3)fϵϵ∣Σ=2π-np/2Σ-n/2exp⁡-12tr⁡ϵTϵΣ-1,
where tr⁡· is the standard trace function.


### 2.2. Likelihood Function forΣ Conditional on β

From the multivariate multiple regression model (2) we can write ϵ=Y-Xβ. Therefore, the joint likelihood function for β and Σ is given by the following:
(4)Lβ,Σ∣Y=2π-np/2Σ-n/2×exp⁡-12tr⁡Y-XβTY-XβΣ-1.For a given value ofβ let S=n-1Y-XβTY-Xβ. Note that S is a p×p symmetric almost surely positive definite matrix. Then the likelihood function (4) for Σ conditional on β can be written as
(5)LΣ∣Y,β=2π-np/2Σ-n/2exp⁡-n2tr⁡SΣ-1.
In Bayesian analysis for a univariate normal model, the logarithm of the variance parameter has been modeled by a univariate normal prior distribution. In a multivariate setting the matrix logarithm of a covariance matrix has also been investigated by Chiu et al. [14]. Along these same lines, we consider the matrix logarithm of Σ and S:
(6)Ap×p=log⁡Σ=Elog⁡DET,Λp×p=log⁡S=E0log⁡D0E0T,
where E is a p×p orthonormal matrix whose columns are normalized eigenvectors and D is a p×p diagonal matrix of the corresponding normalized eigenvalues associated with Σ. E0 and D0 are defined analogously for S. Using the fact that A=log⁡Σ from (6) and noting that Σ=exp⁡tr⁡⁡A, we can express the exact likelihood function (5) in the following equivalent fashion:
(7)LA∣Y,β=2π-np/2exp⁡-n2tr⁡A+Sexp⁡-A.
We now define the following unconventional matrix operator Vec*·. Let aij be the element in the ith row and jth column of the matrix A, and then
(8)αq×1=Vec*A=a11,a22,…,app∣a12,a23,…ap-1,p∣…,a1,p-1,a2p∣a1pT,
where α=α1,…,αqT is a q×1 vector and q=1/2pp+1. We analogously define λ=Vec*Λ, which will appear again in Section 3.4. Moving forward we will use α, which captures the unique elements of the matrix logarithm of the covariance matrix, to model the covariance structure.


### 2.3. Prior Distributional Specfications

We will assumea priori that β is independent of Σ. More specifically, with respect to β we make the assumption of a uniform prior distribution:
(9)πβ∝1.
Note that the uniform prior assumption for β can be viewed as a limiting case of the informative multivariate normal prior specification, which this modeling approach could accommodate.We will assumea priori that, given η and Υ, α follows a q dimensional normal distribution with mean location hyperparameter vector η and covariance hyperparameter matrix Υ. The multivariate normal provides a very rich and flexible family of prior distributions for the matrix logarithm of the covariance structure. This adds far greater flexibility than the conventional inverse Wishart prior specification. Since the multivariate normal is fully parameterized by a mean vector and covariance matrix, we have the ability to model more complex prior information. In particular, we can specify different prior mean values for each element of α via the elements of the location hyperparameter η. Moreover, we have the ability to model varying degrees of strength of the prior belief in each of the q elements of η through the q diagonal elements of the covariance hyperparameter matrix Υ. Additionally, with the multivariate normal prior we are able to model potential interdependency among the elements of α because we can specify nontrivial covariance terms in the covariance hyperparameter matrix. That is, the off diagonal elements of Υ can be used to specify any potential correlations amongst the elements of α.We are now able to craft a more complex and accurate prior specification for the covariance structure. A subjective Bayesian may in fact wish to specify allq+1/2qq+1 hyperparameters. In this way, the practitioner can fully take advantage of any relevant prior information through use of the flexible multivariate normal prior specification for the covariance structure. Alternatively, we can opt to model η=ημ and Υ=Υσ, where μ and σ are of smaller order than η and Υ, respectively. That is,a priori we may wish to only model certain subsets of the covariance structure. An obvious choice is to consider the variance components as one subset and the covariance components as another. However, we stress the point that the fully general multivariate normal prior specification can be utilized in its totality.Here we will consider the intraclass matrix form for the prior specification as an example of the fully generalized multivariate normal prior distribution. Specifically, we will consider the firstp elements of α separate from the remaining q-p terms. That is, we wish to model the variance components separately from the covariance components of α. Formally, we assume α∣μ,Δ~NqJμ,Δ for the prior distribution. We have the following prior distributional form:
(10)πα∣μ,Δ∝Δ-1/2exp⁡-12α-JμTΔ-1α-Jμ,
where the 2×1 vector μ=μ1,μ2T and
(11)Jq×2=1⋯10⋯00⋯01⋯1T,Δq×q=σ12Ip00σ22Iq-p.
Note that J is a q×2 matrix, whose first p elements of the first column are equal to one and the remaining q-p terms of the first column are equal to zero. The second column of J consists of the first p elements equal to zero and the remaining q-p elements equal to one. In Δ, Ip and Iq-p are indicator matrices of dimension p and q-p, respectively.Thus,μ1 and σ12 are the location and variance hyperparameters, respectively, for the variance components of α. Analogously, μ2 and σ22 are the location and variance hyperparameters for the covariance components of α. In this way we can specify two different location hyperparameters and two levels of confidence.For the hyperparameters,μ=μ1,μ2T and Δ=hσ12,σ22, we will assume the following vague prior distribution:
(12)πμ,Δ∝1.
Note here that the uniform prior specification can be viewed as a limiting case of a multivariate normal and inverse Wishart prior specification for μ and Δ, respectively. Furthermore, the analysis could in fact accommodate such nontrivial specifications quite easily. Having stated all the prior distributional assumptions we now turn to the posterior Bayesian analysis. We begin this by first examining the exact joint posterior distribution.


### 2.4. Exact Joint Posterior Distribution

The exact joint posterior distribution for all parameters and hyperparameters will be proportional to the product of (4) the exact likelihood function, (9) the prior distribution for β, (10) the prior distribution for α, and (12) the prior distribution for μ and Δ. Note that here we will use α interchangeably with Σ:
(13)πα,β,μ,Δ∣Y∝Δ-1/2Σ-1/2×exp⁡-12tr⁡Y-XβTY-XβΣ-1×exp⁡-12α-JμTΔ-1α-Jμ.
The μ term in (13) can be integrated out, so that
(14)πα,β,Δ∣Y∝Δ-1/2Σ-1/2×exp⁡-12tr⁡Y-XβTY-XβΣ-1×JTΔ-1Jexp⁡-12α-Jμ^TΔ-1α-Jμ^,
where μ^=(XTΔ-1X)-1XTΔ-1α.We clearly see that the exact joint posterior distribution is in fact not tractable. This is the driving motivation behind the implementation of the numerical MHWG sampling techniques. Rather than working with the cumbersome exact joint posterior distribution it is much easier to consider the so-called full conditional distributions for each of the parameters and hyperparameters.


## 2.1. Distributional Assumptions and Matrix Normal Distribution

The matrix normal distribution is closely related to, and is a generalization of the multivariate normal. In particular, then×p random matrix M~MNn×pΦ,Σ,Ω, if and only if, the np×1 random vector Vec⁡M~NnpVec⁡Φ,Σ⊗Ω, where MNn×p denotes the n×p dimensional matrix normal distribution, Φ is a n×p location matrix, Σ is a p×pfirst covariance matrix, and Ω is a n×nsecond covariance matrix [13, page 54]. Vec⁡· and ⊗ are the standard vector operator and Kronecker product, respectively.We make the distributional assumption that, conditional on thep×p covariance matrix Σ, the n×p random error matrix ϵ=(ϵ1T,ϵ2T,…,ϵnT)T, in (2), follows a matrix normal with n×p zero mean matrix and covariance matrices given by Σ and In, where In is a n×n identity matrix. Formally, we have ϵ∣Σ~MNn×p0,Σ,In, or equivalently, the error terms ϵ1,ϵ2,…,ϵn are independent and identically distributed normal random vectors each with mean vector 0 and covariance matrix Σ. The probability density function for the error matrix is given by
(3)fϵϵ∣Σ=2π-np/2Σ-n/2exp⁡-12tr⁡ϵTϵΣ-1,
where tr⁡· is the standard trace function.


## 2.2. Likelihood Function forΣ Conditional on β

From the multivariate multiple regression model (2) we can write ϵ=Y-Xβ. Therefore, the joint likelihood function for β and Σ is given by the following:
(4)Lβ,Σ∣Y=2π-np/2Σ-n/2×exp⁡-12tr⁡Y-XβTY-XβΣ-1.For a given value ofβ let S=n-1Y-XβTY-Xβ. Note that S is a p×p symmetric almost surely positive definite matrix. Then the likelihood function (4) for Σ conditional on β can be written as
(5)LΣ∣Y,β=2π-np/2Σ-n/2exp⁡-n2tr⁡SΣ-1.
In Bayesian analysis for a univariate normal model, the logarithm of the variance parameter has been modeled by a univariate normal prior distribution. In a multivariate setting the matrix logarithm of a covariance matrix has also been investigated by Chiu et al. [14]. Along these same lines, we consider the matrix logarithm of Σ and S:
(6)Ap×p=log⁡Σ=Elog⁡DET,Λp×p=log⁡S=E0log⁡D0E0T,
where E is a p×p orthonormal matrix whose columns are normalized eigenvectors and D is a p×p diagonal matrix of the corresponding normalized eigenvalues associated with Σ. E0 and D0 are defined analogously for S. Using the fact that A=log⁡Σ from (6) and noting that Σ=exp⁡tr⁡⁡A, we can express the exact likelihood function (5) in the following equivalent fashion:
(7)LA∣Y,β=2π-np/2exp⁡-n2tr⁡A+Sexp⁡-A.
We now define the following unconventional matrix operator Vec*·. Let aij be the element in the ith row and jth column of the matrix A, and then
(8)αq×1=Vec*A=a11,a22,…,app∣a12,a23,…ap-1,p∣…,a1,p-1,a2p∣a1pT,
where α=α1,…,αqT is a q×1 vector and q=1/2pp+1. We analogously define λ=Vec*Λ, which will appear again in Section 3.4. Moving forward we will use α, which captures the unique elements of the matrix logarithm of the covariance matrix, to model the covariance structure.


## 2.3. Prior Distributional Specfications

We will assumea priori that β is independent of Σ. More specifically, with respect to β we make the assumption of a uniform prior distribution:
(9)πβ∝1.
Note that the uniform prior assumption for β can be viewed as a limiting case of the informative multivariate normal prior specification, which this modeling approach could accommodate.We will assumea priori that, given η and Υ, α follows a q dimensional normal distribution with mean location hyperparameter vector η and covariance hyperparameter matrix Υ. The multivariate normal provides a very rich and flexible family of prior distributions for the matrix logarithm of the covariance structure. This adds far greater flexibility than the conventional inverse Wishart prior specification. Since the multivariate normal is fully parameterized by a mean vector and covariance matrix, we have the ability to model more complex prior information. In particular, we can specify different prior mean values for each element of α via the elements of the location hyperparameter η. Moreover, we have the ability to model varying degrees of strength of the prior belief in each of the q elements of η through the q diagonal elements of the covariance hyperparameter matrix Υ. Additionally, with the multivariate normal prior we are able to model potential interdependency among the elements of α because we can specify nontrivial covariance terms in the covariance hyperparameter matrix. That is, the off diagonal elements of Υ can be used to specify any potential correlations amongst the elements of α.We are now able to craft a more complex and accurate prior specification for the covariance structure. A subjective Bayesian may in fact wish to specify allq+1/2qq+1 hyperparameters. In this way, the practitioner can fully take advantage of any relevant prior information through use of the flexible multivariate normal prior specification for the covariance structure. Alternatively, we can opt to model η=ημ and Υ=Υσ, where μ and σ are of smaller order than η and Υ, respectively. That is,a priori we may wish to only model certain subsets of the covariance structure. An obvious choice is to consider the variance components as one subset and the covariance components as another. However, we stress the point that the fully general multivariate normal prior specification can be utilized in its totality.Here we will consider the intraclass matrix form for the prior specification as an example of the fully generalized multivariate normal prior distribution. Specifically, we will consider the firstp elements of α separate from the remaining q-p terms. That is, we wish to model the variance components separately from the covariance components of α. Formally, we assume α∣μ,Δ~NqJμ,Δ for the prior distribution. We have the following prior distributional form:
(10)πα∣μ,Δ∝Δ-1/2exp⁡-12α-JμTΔ-1α-Jμ,
where the 2×1 vector μ=μ1,μ2T and
(11)Jq×2=1⋯10⋯00⋯01⋯1T,Δq×q=σ12Ip00σ22Iq-p.
Note that J is a q×2 matrix, whose first p elements of the first column are equal to one and the remaining q-p terms of the first column are equal to zero. The second column of J consists of the first p elements equal to zero and the remaining q-p elements equal to one. In Δ, Ip and Iq-p are indicator matrices of dimension p and q-p, respectively.Thus,μ1 and σ12 are the location and variance hyperparameters, respectively, for the variance components of α. Analogously, μ2 and σ22 are the location and variance hyperparameters for the covariance components of α. In this way we can specify two different location hyperparameters and two levels of confidence.For the hyperparameters,μ=μ1,μ2T and Δ=hσ12,σ22, we will assume the following vague prior distribution:
(12)πμ,Δ∝1.
Note here that the uniform prior specification can be viewed as a limiting case of a multivariate normal and inverse Wishart prior specification for μ and Δ, respectively. Furthermore, the analysis could in fact accommodate such nontrivial specifications quite easily. Having stated all the prior distributional assumptions we now turn to the posterior Bayesian analysis. We begin this by first examining the exact joint posterior distribution.


## 2.4. Exact Joint Posterior Distribution

The exact joint posterior distribution for all parameters and hyperparameters will be proportional to the product of (4) the exact likelihood function, (9) the prior distribution for β, (10) the prior distribution for α, and (12) the prior distribution for μ and Δ. Note that here we will use α interchangeably with Σ:
(13)πα,β,μ,Δ∣Y∝Δ-1/2Σ-1/2×exp⁡-12tr⁡Y-XβTY-XβΣ-1×exp⁡-12α-JμTΔ-1α-Jμ.
The μ term in (13) can be integrated out, so that
(14)πα,β,Δ∣Y∝Δ-1/2Σ-1/2×exp⁡-12tr⁡Y-XβTY-XβΣ-1×JTΔ-1Jexp⁡-12α-Jμ^TΔ-1α-Jμ^,
where μ^=(XTΔ-1X)-1XTΔ-1α.We clearly see that the exact joint posterior distribution is in fact not tractable. This is the driving motivation behind the implementation of the numerical MHWG sampling techniques. Rather than working with the cumbersome exact joint posterior distribution it is much easier to consider the so-called full conditional distributions for each of the parameters and hyperparameters.


## 3. Markov Chain Monte Carlo Approach

As already noted, the joint posterior distribution in (14) is not analytically tractable with respect to drawing inference for the relevant parameters and hyperparameters. This will give rise to our consideration in the subsequent subsections of the full conditional posterior distributions. The conditional posterior distributions we derive below will provide the framework for the MHWG sampling techniques.


### 3.1. Exact Conditional Posterior Distribution forβ

Recall thatβ^=XTX-1XTY is the maximum likelihood estimator of β. Furthermore, we can rewrite the exact likelihood function (4) as
(15)Lβ,Σ∣Y=2π-np/2Σ-n/2×exp⁡-12tr⁡Y-Xβ^TY-Xβ^Σ-1IIIIIIIIIIIIIIIIIII12IIIIIIIIIIIIIIIIII+β-β^TXTXβ-β^Σ-1.
Therefore, the posterior distribution for β conditional on Σ will be proportional, with respect to only the terms involving β, to the product of (15) the exact likelihood function multiplied by (9) the prior distribution for β:
(16)πβ∣Y,Σ∝exp⁡-12tr⁡β-β^TXTXβ-β^Σ-1.
We recognize that the posterior distribution of β conditional on Σ is of a matrix normal form. Making use of the relationship between the matrix normal and the multivariate normal distributions as stated above in Section 2.1, we have
(17)Vec⁡β∣Y,Σ~NkpVec⁡β^,Σ⊗XTX-1.


### 3.2. Exact Conditional Posterior Distribution forα

The joint prior distribution forα, μ, and Δ is given by the product of (10) the prior distribution for α and (12) the joint prior distribution for μ and Δ. We can derive the joint prior distribution for just α and Δ by integrating over the joint prior distribution for α, μ, and Δ with respect to μ. Upon completion of the integration we have the following joint prior distribution for α and Δ:(18)πα,Δ∝Δ-1/2JTΔ-1J-1/2exp⁡-12αTG*α,
where
(19)G*q×q=Iq-JJTΔ-1J-1JTΔ-1T×Δ-1Iq-JJTΔ-1J-1JTΔ-1
and Iq is a q×q identity matrix. By integrating μ out we will help to facilitate the MCMC procedure both in terms of speed and simplification of the algorithm.The exact posterior distribution will be proportional to the product of (7) the exact likelihood function, multiplied by (18) the joint prior distribution for α and Δ. Note that here we will use α interchangeably with A:
(20)πα∣Y,β,Δ∝exp⁡-n2tr⁡A+Sexp⁡-A-12αTG*α.
Note that the above exact posterior distribution is not of a known form and cannot be directly simulated from in an easy fashion. This will motivate us to use a proposal density that closely matches this target density in a MHWG sampling routine.


### 3.3. Exact Conditional Posterior Distribution forΔ

The exact posterior distribution forΔ=hσ12,σ22 conditional on β and α will be proportional to (18) the joint prior distribution for α and Δ. Note that the exact likelihood function (7) does not depend upon Δ and thus can be omitted entirely:
(21)πΔ∣Y,α,β∝σ1-(p-1)/2σ2-(q-p-1)/2×exp⁡-12σ12∑i=1p‍αi-α¯v2-12σ22∑i=p+1q‍αi-α¯c2,
where α¯v=1/p∑i=1p‍αi and α¯c=1/(q-p)∑i=p+1q‍αi are the arithmetic means of the variance and covariance components of α, respectively. We recognize that the posterior distributions for σ12 and σ22 conditional on α are independent Inverse Gamma random variables:
(22)σ12∣Y,α,β~InverseGammap-32,12∑i=1p‍αi-α¯v2,(23)σ22∣Y,α,β~InverseGammaq-p-32,12∑i=p+1q‍αi-α¯c2.
This result is intuitive and theoretically appealing in that the posterior distribution for σ12, the variance hyperparameter for the variance components of α, depends only on the number of variance terms p and α1,…,αp, whereas the posterior distribution for σ22, the variance hyperparameter for the covariance components of α, depends only on the number of covariance terms q-p and αp+1,…,αq. This draws out the point of modeling the variance components separate from the covariance components. In addition, the Inverse Gamma is highly tractable and lends itself to the numerical procedures in the subsequent section.


### 3.4. Approximate Conditional Posterior Distribution forα

In order to implement the MHWG sampling algorithm, we have derived all the full conditional posterior distributionsπβ∣Y,Σ, πα∣Y,β,Δ, and πΔ∣Y,α,β in (17), (20), (22), and (23), respectively. However, we clearly see that the simulation of α based upon the true conditional posterior distribution in (20), the target distribution, is not tractable. Therefore, the Metropolis-Hastings algorithm is employed and a proposal density needs to be constructed. The algorithm works best if the proposal density closely matches the shape of the target distribution.We can construct an approximation to a function, that is, proportional to the likelihood function by utilizing the linear Volterra integral. The key advantage of this is that the approximation can be written as a multivariate normal with respect to the unique elements of the matrix logarithm of the covariance matrix. This allows for a multivariate normal to act as conjugate prior. Hence, we have a multivariate normal posterior, that is, a good proxy for the true posterior, and the proposal can be easily simulated from. Interested readers should refer to Leonard and Hsu [7] and Hsu et al. [8] for a detailed exposition of how this is performed.Leonard and Hsu [7] show how we can use Bellman’s solution of a Volterra integral equation [10, page 171] to derive the following approximation, that is, proportioanl to the likelihood function of α given Y and β:
(24)L*α∣Y,β∝exp⁡-12α-λTQα-λ.
Recall from Section 2.2 that λ=Vec*Λ, where Λ is as defined in (6). The q×q symmetric almost surely positive definite matrix Q is the likelihood information matrix of α and is a function of the normalized eigenvalues and normalized eigenvectors of S. In particular,
(25)Qq×q=n2∑i=1p‍fiifiiT+n∑pppi<jp∑ppξijfijfijT,
where fijq×1=ei*ej and
(26)ξij=di-dj2didjlog⁡di-log⁡dj2
where dj and ej for j=1,…,p are the jth normalized eigenvalue and eigenvector, respectively, of S. fij denotes the q×1 vector that satisfies the condition αTei*ej=eiTAej. We see that the approximate likelihood function (24) is a multivariate normal form with respect to α. This functional form of the approximate likelihood function in (24) will be the driving mechanism in the Bayesian analysis for α. Again, for details of the derivation of the approximate likelihood function please refer to Leonard and Hsu [7] and Hsu et al. [8].Equation (24) provides an excellent approximation to a function, that is, proportional to the exact likelihood (7), when the sample size n is large. To illustrate the effects of n, the sample size, and p, the dimension of the covariance matrix Σ, we conduct the following exercise. Without loss of generality, we consider a simplified model, when Y1,…,Yn is a random sample from Np0,Σ. In our illustrative example, three dimensional sizes (p=3, 5, and 10), and four sample sizes (n=20, 100, 500, and 5000) were considered for comparison. For a fair comparison, the same sample covariance is used for all four different sample sizes, for each p. The sample covariance matrix S is assumed to consist of elements sij for the ith row and jth column, where sij=1.0-i-j×0.1. For example, when p=3, then
(27)S=1.00.90.80.91.00.90.80.91.0.The histograms, in Figure1, represent the exact likelihood from (7) of the (1, 1)th element, a11, of A, where the histogram is normalized for comparison purposes, and the dotted curves represent the univariate normal densities based on approximation (24). Please note that these histograms are computed according to an importance sampling method using (24) as the importance function. For an overview of importance sampling methods, please see, for example, Rubinstein [15] and Leonard et al. [16]. It can be seen from Figure 1 that the approximation is better when the sample size n is bigger and the dimension size p is smaller. Similar patterns were found for other variance and covariance elements of the covariance matrix A. The approximation is fairly accurate when n is 500 or greater.Figure 1
Comparison of the approximate likelihood (normalized) with the exact likelihood (normalized) for covariance matrices of dimensionsp = 3, 5, and 10, and sample sizes n = 20, 100, 500, and 5000. The histogram is the normalized likelihood and the dashed curve is the approximate normalized likelihood.The approximate joint posterior distribution forα and Δ conditional on β will be proportional to the product of the approximate likelihood function (24) and the joint prior distribution (18):
(28)π*α,Δ∣Y,β∝σ1-(p-1)/2σ2-(q-p-1)/2×exp⁡-12α-λTQα-λ+αTG*α.Completing the square for the two terms in the exponent of (28) yields the following approximate posterior distribution for α conditional on β and Δ:
(29)π*α∣Y,β,Δ∝exp⁡-12α-α*TQ+G*α-α*,
where the q×1 vector α*=Q+G*-1Qλ. Recall that Q and G* are as defined in (25) and (19), respectively. Thus, we have the following approximate posterior distribution for α conditional on β and Δ:
(30)α∣Y,β,Δ~Nqα*,Q+G*-1.
This demonstrates the conjugacy of utilizing the approximate likelihood function. Equation (30) provides an efficient proposal distribution for implementing a MHWG algorithm. In short, we have developed a highly flexible while at the same time tractable Bayesian methodology for the covariance structure.


### 3.5. Metropolis-Hastings within Gibbs Sampling Procedure

Based upon the theoretical results derived above we outline the following procedure for implementing the MHWG algorithm. Specifically, from (17), (30), (20), (22), and (23) we have a formal setup for implementing a MCMC procedure with a Metropolis-Hastings step. Below we outline the specific steps involved in implementing the MHWG algorithm.(1)
Simulateσ12t and σ22t from (22) and (23), respectively. Initial starting values for α may be set equal to λ^=Vec*log⁡S^, where
(31)S^=n-1Y-Xβ^TY-Xβ^andβ^ is the maximum likelihood estimator for β. Subsequent simulations of σ12t and σ22t will be based upon simulated values of αt.(2)
Simulateβt from (17). Initial starting values for Σ may be based upon S^. Subsequent simulations of βt will be based upon simulated values of αt.(3)
Simulate a candidate valueα~ from (30) based upon σ12t and σ22t and βt from steps (1) and (2), respectively. Then let
(32)αt+1=α~,withprobabilitymin⁡ρ,1,αt,otherwise,whereρ is given by following expression:
(33)ρ=πα~∣Y,βt,Δtπ*α~∣Y,βt,Δtπ*αt∣Y,βt,Δtπαt∣Y,βt,Δtandπ*·∣· and π·∣· are as defined in (29) and (20), respectively.The last procedure in step(3) is the Metropolis-Hastings algorithm. We employ this procedure since we are utilizing an approximation to the exact posterior distribution [11, page 291].Posterior moments for any parameter of interest can be calculated easily upon the MHWG results. It is usually the case in MHWG estimation procedures that the firstl simulations are not included in the estimates, due to the fact that the Markov chain has not yet reached a steady state. In the parlance of numerical procedures, l is usually referred to as the number ofburn in iterations.In addition to the burn in value, we also explored the potential need to perform a so-calledthinning procedure. A thinning step entails only retaining every rth simulated value of the MHWG sampling, where r is chosen large enough so that any autocorrelation is removed. By examining the sample autocorrelation function (ACF) plots, practitioners can decide if thinning is necessary. We investigated several plots for numerous parameters. For illustrative purposes, in Figure 2, we present the sample ACF plots for β11, β12, α1, and σ12. We can see that autocorrelation is not a significant concern in our particular application. Other sample ACF plots looked quite similar to numerous other parameters. Practitioners should explore the need to use a thinning procedure in their specific applications.Figure 2
The sample ACF plots demonstrate that in our application a thinning procedure was not necessary. Other sample ACF plots for various parameters looked quite similar. Practitioners should check for this in their applications.


## 3.1. Exact Conditional Posterior Distribution forβ

Recall thatβ^=XTX-1XTY is the maximum likelihood estimator of β. Furthermore, we can rewrite the exact likelihood function (4) as
(15)Lβ,Σ∣Y=2π-np/2Σ-n/2×exp⁡-12tr⁡Y-Xβ^TY-Xβ^Σ-1IIIIIIIIIIIIIIIIIII12IIIIIIIIIIIIIIIIII+β-β^TXTXβ-β^Σ-1.
Therefore, the posterior distribution for β conditional on Σ will be proportional, with respect to only the terms involving β, to the product of (15) the exact likelihood function multiplied by (9) the prior distribution for β:
(16)πβ∣Y,Σ∝exp⁡-12tr⁡β-β^TXTXβ-β^Σ-1.
We recognize that the posterior distribution of β conditional on Σ is of a matrix normal form. Making use of the relationship between the matrix normal and the multivariate normal distributions as stated above in Section 2.1, we have
(17)Vec⁡β∣Y,Σ~NkpVec⁡β^,Σ⊗XTX-1.


## 3.2. Exact Conditional Posterior Distribution forα

The joint prior distribution forα, μ, and Δ is given by the product of (10) the prior distribution for α and (12) the joint prior distribution for μ and Δ. We can derive the joint prior distribution for just α and Δ by integrating over the joint prior distribution for α, μ, and Δ with respect to μ. Upon completion of the integration we have the following joint prior distribution for α and Δ:(18)πα,Δ∝Δ-1/2JTΔ-1J-1/2exp⁡-12αTG*α,
where
(19)G*q×q=Iq-JJTΔ-1J-1JTΔ-1T×Δ-1Iq-JJTΔ-1J-1JTΔ-1
and Iq is a q×q identity matrix. By integrating μ out we will help to facilitate the MCMC procedure both in terms of speed and simplification of the algorithm.The exact posterior distribution will be proportional to the product of (7) the exact likelihood function, multiplied by (18) the joint prior distribution for α and Δ. Note that here we will use α interchangeably with A:
(20)πα∣Y,β,Δ∝exp⁡-n2tr⁡A+Sexp⁡-A-12αTG*α.
Note that the above exact posterior distribution is not of a known form and cannot be directly simulated from in an easy fashion. This will motivate us to use a proposal density that closely matches this target density in a MHWG sampling routine.


## 3.3. Exact Conditional Posterior Distribution forΔ

The exact posterior distribution forΔ=hσ12,σ22 conditional on β and α will be proportional to (18) the joint prior distribution for α and Δ. Note that the exact likelihood function (7) does not depend upon Δ and thus can be omitted entirely:
(21)πΔ∣Y,α,β∝σ1-(p-1)/2σ2-(q-p-1)/2×exp⁡-12σ12∑i=1p‍αi-α¯v2-12σ22∑i=p+1q‍αi-α¯c2,
where α¯v=1/p∑i=1p‍αi and α¯c=1/(q-p)∑i=p+1q‍αi are the arithmetic means of the variance and covariance components of α, respectively. We recognize that the posterior distributions for σ12 and σ22 conditional on α are independent Inverse Gamma random variables:
(22)σ12∣Y,α,β~InverseGammap-32,12∑i=1p‍αi-α¯v2,(23)σ22∣Y,α,β~InverseGammaq-p-32,12∑i=p+1q‍αi-α¯c2.
This result is intuitive and theoretically appealing in that the posterior distribution for σ12, the variance hyperparameter for the variance components of α, depends only on the number of variance terms p and α1,…,αp, whereas the posterior distribution for σ22, the variance hyperparameter for the covariance components of α, depends only on the number of covariance terms q-p and αp+1,…,αq. This draws out the point of modeling the variance components separate from the covariance components. In addition, the Inverse Gamma is highly tractable and lends itself to the numerical procedures in the subsequent section.


## 3.4. Approximate Conditional Posterior Distribution forα

In order to implement the MHWG sampling algorithm, we have derived all the full conditional posterior distributionsπβ∣Y,Σ, πα∣Y,β,Δ, and πΔ∣Y,α,β in (17), (20), (22), and (23), respectively. However, we clearly see that the simulation of α based upon the true conditional posterior distribution in (20), the target distribution, is not tractable. Therefore, the Metropolis-Hastings algorithm is employed and a proposal density needs to be constructed. The algorithm works best if the proposal density closely matches the shape of the target distribution.We can construct an approximation to a function, that is, proportional to the likelihood function by utilizing the linear Volterra integral. The key advantage of this is that the approximation can be written as a multivariate normal with respect to the unique elements of the matrix logarithm of the covariance matrix. This allows for a multivariate normal to act as conjugate prior. Hence, we have a multivariate normal posterior, that is, a good proxy for the true posterior, and the proposal can be easily simulated from. Interested readers should refer to Leonard and Hsu [7] and Hsu et al. [8] for a detailed exposition of how this is performed.Leonard and Hsu [7] show how we can use Bellman’s solution of a Volterra integral equation [10, page 171] to derive the following approximation, that is, proportioanl to the likelihood function of α given Y and β:
(24)L*α∣Y,β∝exp⁡-12α-λTQα-λ.
Recall from Section 2.2 that λ=Vec*Λ, where Λ is as defined in (6). The q×q symmetric almost surely positive definite matrix Q is the likelihood information matrix of α and is a function of the normalized eigenvalues and normalized eigenvectors of S. In particular,
(25)Qq×q=n2∑i=1p‍fiifiiT+n∑pppi<jp∑ppξijfijfijT,
where fijq×1=ei*ej and
(26)ξij=di-dj2didjlog⁡di-log⁡dj2
where dj and ej for j=1,…,p are the jth normalized eigenvalue and eigenvector, respectively, of S. fij denotes the q×1 vector that satisfies the condition αTei*ej=eiTAej. We see that the approximate likelihood function (24) is a multivariate normal form with respect to α. This functional form of the approximate likelihood function in (24) will be the driving mechanism in the Bayesian analysis for α. Again, for details of the derivation of the approximate likelihood function please refer to Leonard and Hsu [7] and Hsu et al. [8].Equation (24) provides an excellent approximation to a function, that is, proportional to the exact likelihood (7), when the sample size n is large. To illustrate the effects of n, the sample size, and p, the dimension of the covariance matrix Σ, we conduct the following exercise. Without loss of generality, we consider a simplified model, when Y1,…,Yn is a random sample from Np0,Σ. In our illustrative example, three dimensional sizes (p=3, 5, and 10), and four sample sizes (n=20, 100, 500, and 5000) were considered for comparison. For a fair comparison, the same sample covariance is used for all four different sample sizes, for each p. The sample covariance matrix S is assumed to consist of elements sij for the ith row and jth column, where sij=1.0-i-j×0.1. For example, when p=3, then
(27)S=1.00.90.80.91.00.90.80.91.0.The histograms, in Figure1, represent the exact likelihood from (7) of the (1, 1)th element, a11, of A, where the histogram is normalized for comparison purposes, and the dotted curves represent the univariate normal densities based on approximation (24). Please note that these histograms are computed according to an importance sampling method using (24) as the importance function. For an overview of importance sampling methods, please see, for example, Rubinstein [15] and Leonard et al. [16]. It can be seen from Figure 1 that the approximation is better when the sample size n is bigger and the dimension size p is smaller. Similar patterns were found for other variance and covariance elements of the covariance matrix A. The approximation is fairly accurate when n is 500 or greater.Figure 1
Comparison of the approximate likelihood (normalized) with the exact likelihood (normalized) for covariance matrices of dimensionsp = 3, 5, and 10, and sample sizes n = 20, 100, 500, and 5000. The histogram is the normalized likelihood and the dashed curve is the approximate normalized likelihood.The approximate joint posterior distribution forα and Δ conditional on β will be proportional to the product of the approximate likelihood function (24) and the joint prior distribution (18):
(28)π*α,Δ∣Y,β∝σ1-(p-1)/2σ2-(q-p-1)/2×exp⁡-12α-λTQα-λ+αTG*α.Completing the square for the two terms in the exponent of (28) yields the following approximate posterior distribution for α conditional on β and Δ:
(29)π*α∣Y,β,Δ∝exp⁡-12α-α*TQ+G*α-α*,
where the q×1 vector α*=Q+G*-1Qλ. Recall that Q and G* are as defined in (25) and (19), respectively. Thus, we have the following approximate posterior distribution for α conditional on β and Δ:
(30)α∣Y,β,Δ~Nqα*,Q+G*-1.
This demonstrates the conjugacy of utilizing the approximate likelihood function. Equation (30) provides an efficient proposal distribution for implementing a MHWG algorithm. In short, we have developed a highly flexible while at the same time tractable Bayesian methodology for the covariance structure.


## 3.5. Metropolis-Hastings within Gibbs Sampling Procedure

Based upon the theoretical results derived above we outline the following procedure for implementing the MHWG algorithm. Specifically, from (17), (30), (20), (22), and (23) we have a formal setup for implementing a MCMC procedure with a Metropolis-Hastings step. Below we outline the specific steps involved in implementing the MHWG algorithm.(1)
Simulateσ12t and σ22t from (22) and (23), respectively. Initial starting values for α may be set equal to λ^=Vec*log⁡S^, where
(31)S^=n-1Y-Xβ^TY-Xβ^andβ^ is the maximum likelihood estimator for β. Subsequent simulations of σ12t and σ22t will be based upon simulated values of αt.(2)
Simulateβt from (17). Initial starting values for Σ may be based upon S^. Subsequent simulations of βt will be based upon simulated values of αt.(3)
Simulate a candidate valueα~ from (30) based upon σ12t and σ22t and βt from steps (1) and (2), respectively. Then let
(32)αt+1=α~,withprobabilitymin⁡ρ,1,αt,otherwise,whereρ is given by following expression:
(33)ρ=πα~∣Y,βt,Δtπ*α~∣Y,βt,Δtπ*αt∣Y,βt,Δtπαt∣Y,βt,Δtandπ*·∣· and π·∣· are as defined in (29) and (20), respectively.The last procedure in step(3) is the Metropolis-Hastings algorithm. We employ this procedure since we are utilizing an approximation to the exact posterior distribution [11, page 291].Posterior moments for any parameter of interest can be calculated easily upon the MHWG results. It is usually the case in MHWG estimation procedures that the firstl simulations are not included in the estimates, due to the fact that the Markov chain has not yet reached a steady state. In the parlance of numerical procedures, l is usually referred to as the number ofburn in iterations.In addition to the burn in value, we also explored the potential need to perform a so-calledthinning procedure. A thinning step entails only retaining every rth simulated value of the MHWG sampling, where r is chosen large enough so that any autocorrelation is removed. By examining the sample autocorrelation function (ACF) plots, practitioners can decide if thinning is necessary. We investigated several plots for numerous parameters. For illustrative purposes, in Figure 2, we present the sample ACF plots for β11, β12, α1, and σ12. We can see that autocorrelation is not a significant concern in our particular application. Other sample ACF plots looked quite similar to numerous other parameters. Practitioners should explore the need to use a thinning procedure in their specific applications.Figure 2
The sample ACF plots demonstrate that in our application a thinning procedure was not necessary. Other sample ACF plots for various parameters looked quite similar. Practitioners should check for this in their applications.


## 4. Application: High School and Beyond Survey

In 1980 the National Education Longitudinal Studies program of the National Center for Education Statistics administered the High School and Beyond (HSB) Survey [12]. The HSB study contains both a 1980 senior class cohort and a 1980 sophomore class cohort. Within the sophomore class cohort we have a total sample size of n = 14,667 students. The HSB study has been analyzed extensively by many, for example, Astone and McLanahan [17], Grogger [18], St. John [19], and Zwick and Sklar [20].


### 4.1. Description of Data

The HSB study contains a myriad of data and variables. In particular, for the sophomore class cohort a total ofp=7 exams were administered in the areas of vocabulary, reading, two exams in mathematics, science, writing, and civics. As is often the case with educational testing data the test scores were standardized to have a mean of fifty and a standard deviation equal to ten. We will use these standardized test scores for the sophomore class cohort as our multivariate response variables.Grade point average data was also obtained from official transcripts that were included in the survey. In addition, a number of demographic variables were collected on both the school and individual student levels. These variables included school type, relative urbanization of school environment, geographic region of the country, gender, and race. All of these variables will serve as categorical or qualitative explanatory variables in the multivariate multiple regression. Table1 provides a description of the categorical explanatory variables as well as the associated number of students per category. Note that the original HSB study did not include a school type four.Table 1
Categorical explanatory variable descriptions.
Categorical variables
Levels
Level description
Students per level
School type
School type 0
Regular public
9534
School type 1
Alternative
437
School type 2
Cuban/Hispanic public
156
School type 3
Other Hispanic public
1483
School type 5
Regular catholic
1341
School type 6
Black catholic
758
School type 7
Cuban/Hispanic catholic
216
School type 8
Private elite
294
School type 9
Private nonelite
448
Urbanization
Urban 1
Urban
3451
Urban 2
Suburban
7325
Urban 3
Rural
3891
Geographic region
Region 1
New England
736
Region 2
Middle Atlantic
2518
Region 3
South Atlantic
2178
Region 4
East South Central
713
Region 5
West South Central
1750
Region 6
East North Central
2931
Region 7
West North Central
1128
Region 8
Mountain
696
Region 9
Pacific
2017
Race
Race 1
Hispanic/Spanish
3122
Race 2
American Indian/Alaskan native
224
Race 3
Asian/Pacific Islander
345
Race 4
Black
2050
Race 5
White
8880
Race 6
Other
46
Gender
Gender 1
Male
7265
Gender 2
Female
7402As can be expected with educational data there was some moderate degree of missing data. We employed a data augmentation technique for missing data imputation. The specific procedure was invoked using themice library [21] in R. This particular data augmentation procedure fits nicely within the context of our research here since it employs a multivariate imputation by chained equations technique. That is, it uses a multivariate Gibbs sampler procedure to augment the missing data set. In particular, incomplete columns of data, that is, variables with missing data, are augmented by generating appropriate values of data given the values of the other columns of variables.


### 4.2. Treatment Contrast for Categorical Variables

Characterization of the explanatory data or design matrixX to account for categorical explanatory variables is not unique [22, page 173]. There are actually severalcontrast methods. In our analysis, school type zero, which corresponds to a regular public school, will serve as the base level and will not in fact have its own regressor. In the same fashion, urban will serve as the base level for the relative urbanization categorical explanatory variable. Also, the New England region will act as the base level for the geographic variable. Hispanic or Spanish is the base level for the race identifier variable. Finally, we will designate male as the base gender level. Thus, after properly accounting for the base levels in our particular application k=26.


### 4.3. Classical Multivariate Multiple Regression Results

We consider the standard multivatiate regression model described in (2) for the HSB survey data within the Sophomore class cohort. A total of n = 14,667 students in that cohort participated in the study. We estimated the model for the seven standardized exams VOCAB, READ, MATH 1, MATH 2, SCI, WRITE, and CIVIC regressed on GPA and the other explanatory variables described in the previous section. Most of the regression coefficient estimates are highly significant. In particular, GPA is highly significant for all seven exams. Although, certain particular levels for various categorical variables are not significant for some of the response variables. All of the associated F statistics are highly significant for all response variables. Moreover, based upon the analysis of variance that was performed we can conclude that nearly all the explanatory variables are highly significant for all seven subject area exams.


### 4.4. Posterior Estimates of the Model Parameters

In Tables2, 3, 4, 5, and 6 we present the Bayesian posterior means and the associated standard errors for the matrices β and Σ. In practice, for the MCMC procedure we found that a total iteration size of T = 500,000 was quite sufficient to establish convergence and we used l = 2,000 as the burn in value.Table 2
Posterior mean forβ.
VOCAB
READ
MATH 1
MATH 2
SCI
WRITE
CIVIC
(Intercept)
37.0582
35.6555
33.6377
37.7119
37.3335
31.2303
37.9283
GPA
4.7810
5.2711
5.9424
4.8139
4.4674
5.7163
4.3733
School type 1
−1.1971
−0.3382
−1.0521
−0.5196
−1.1388
−0.9455
−0.5326
School type 2
−0.8747
−1.9121
−0.6682
−0.4325
−2.4736
−1.4395
−2.4039
School type 3
−1.0713
−0.9146
−0.3638
−0.3816
−1.1779
−0.9105
−0.7742
School type 5
2.7359
2.1862
2.3550
1.5275
0.7704
2.5247
2.1433
School type 6
1.8095
1.8877
0.9502
0.2953
0.2920
2.0570
2.1941
School type 7
3.0532
1.6308
1.5283
0.0003
0.5107
2.1525
1.7716
School type 8
10.3236
8.7564
8.8743
8.7145
5.8165
7.0452
5.5020
School type 9
2.9475
1.9583
2.9381
2.9743
2.1760
1.9975
1.7927
Urban 2
0.5558
0.0821
0.6002
0.9824
0.7407
0.5115
0.2248
Urban 3
−0.8006
−0.5359
−0.3611
−0.2156
0.4344
0.0087
−0.3363
Region 2
−1.0396
−0.7905
−0.3661
−1.0051
−0.4256
−0.6938
−0.9750
Region 3
−2.2086
−1.1256
−1.6328
−1.7664
−1.4603
−1.1176
−1.7760
Region 4
−3.9350
−1.8827
−2.4900
−2.8335
−2.0910
−1.3902
−2.3088
Region 5
−2.6954
−1.5560
−1.5315
−2.1349
−1.1494
−0.8252
−1.6827
Region 6
−2.1648
−1.5579
−0.9690
−1.7886
−0.6775
−0.5598
−1.5524
Region 7
−1.9171
−0.4627
0.1163
−0.3902
0.2802
0.2396
−0.3970
Region 8
−1.6480
−0.9278
−1.6167
−2.4756
−0.6947
−0.4478
−1.4824
Region 9
−0.4866
−0.7586
−0.5338
−1.0931
−0.6208
−0.0413
−1.5015
Gender 2
−1.5605
−1.2414
−1.8966
−1.5490
−3.0722
3.0841
0.2048
Race 2
−0.9813
0.5513
−0.8029
−0.2627
0.5303
0.8640
−0.9432
Race 3
1.7286
1.7006
4.9363
3.9146
2.4132
3.5702
1.3489
Race 4
−0.8534
0.0337
−0.6765
−0.5625
−1.3975
−0.3582
0.4884
Race 5
4.5403
3.7077
3.8777
2.3674
4.5873
3.9272
2.7858
Race 6
−0.5276
0.3656
−0.0162
−0.1581
0.4371
0.3194
0.9498Table 3
Posterior standard deviations forβ.
VOCAB
READ
MATH 1
MATH 2
SCI
WRITE
CIVIC
(Intercept)
0.4421
0.4482
0.4195
0.4585
0.4383
0.4201
0.4678
GPA
0.0920
0.0930
0.0872
0.0951
0.0913
0.0874
0.0972
School type 1
0.4327
0.4379
0.4102
0.4483
0.4297
0.4110
0.4568
School type 2
0.7111
0.7199
0.6736
0.7352
0.7047
0.6761
0.7518
School type 3
0.2797
0.2823
0.2651
0.2902
0.2773
0.2662
0.2954
School type 5
0.2561
0.2589
0.2428
0.2659
0.2540
0.2438
0.2707
School type 6
0.3358
0.3392
0.3185
0.3480
0.3327
0.3193
0.3553
School type 7
0.6116
0.6187
0.5786
0.6331
0.6054
0.5806
0.6452
School type 8
0.5164
0.5229
0.4901
0.5359
0.5125
0.4912
0.5457
School type 9
0.4171
0.4222
0.3963
0.4325
0.4138
0.3972
0.4414
Urban 2
0.1864
0.1882
0.1764
0.1927
0.1844
0.1769
0.1964
Urban 3
0.2135
0.2158
0.2020
0.2206
0.2114
0.2027
0.2253
Region 2
0.3623
0.3672
0.3441
0.3759
0.3593
0.3438
0.3825
Region 3
0.3706
0.3759
0.3520
0.3848
0.3681
0.3524
0.3915
Region 4
0.4512
0.4574
0.4286
0.4678
0.4480
0.4290
0.4774
Region 5
0.3877
0.3928
0.3681
0.4020
0.3850
0.3686
0.4097
Region 6
0.3538
0.3587
0.3358
0.3669
0.3513
0.3365
0.3743
Region 7
0.4056
0.4103
0.3846
0.4208
0.4022
0.3858
0.4286
Region 8
0.4647
0.4709
0.4411
0.4819
0.4614
0.4426
0.4917
Region 9
0.3783
0.3830
0.3586
0.3925
0.3753
0.3598
0.3998
Gender 2
0.1421
0.1436
0.1347
0.1470
0.1411
0.1351
0.1500
Race 2
0.5944
0.6018
0.5630
0.6159
0.5897
0.5652
0.6273
Race 3
0.4938
0.4996
0.4679
0.5114
0.4903
0.4692
0.5216
Race 4
0.2599
0.2639
0.2467
0.2697
0.2579
0.2470
0.2752
Race 5
0.2061
0.2085
0.1952
0.2136
0.2042
0.1959
0.2175
Race 6
1.2666
1.2806
1.1995
1.3091
1.2548
1.2040
1.3359Table 4
Posterior mean forΣ.
VOCAB
READ
MATH 1
MATH 2
SCI
WRITE
CIVIC
VOCAB
72.3406
45.3069
34.3079
25.0339
39.1280
35.7151
32.6781
READ
45.3069
74.0616
37.4238
29.2657
40.8084
36.8560
34.3226
MATH 1
34.3079
37.4238
65.0307
39.9597
36.9083
34.2467
27.4186
MATH 2
25.0339
29.2657
39.9597
77.6078
28.6053
25.5057
19.7419
SCI
39.1280
40.8084
36.9083
28.6053
71.1835
34.7257
31.7008
WRITE
35.7151
36.8560
34.2467
25.5057
34.7257
65.3904
31.9590
CIVIC
32.6781
34.3226
27.4186
19.7419
31.7008
31.9590
80.8041Table 5
Classical frequentist estimate forΣ.
VOCAB
READ
MATH 1
MATH 2
SCI
WRITE
CIVIC
VOCAB
72.3308
45.3567
34.2973
24.9684
39.1487
35.7217
32.6627
READ
45.3567
74.0529
37.4264
29.2202
40.8335
36.8651
34.3098
MATH 1
34.2973
37.4264
64.9922
39.9886
36.9119
34.2393
27.3645
MATH 2
24.9684
29.2202
39.9886
77.6275
28.5646
25.4523
19.6396
SCI
39.1487
40.8335
36.9119
28.5646
71.1592
34.7254
31.6765
WRITE
35.7217
36.8651
34.2393
25.4523
34.7254
65.3498
31.9434
CIVIC
32.6627
34.3098
27.3645
19.6396
31.6765
31.9434
80.8225Table 6
Posterior standard deviations forΣ.
VOCAB
READ
MATH 1
MATH 2
SCI
WRITE
CIVIC
VOCAB
0.8456
0.7119
0.6329
0.6507
0.6759
0.6399
0.6872
READ
0.7119
0.8657
0.6514
0.6705
0.6886
0.6513
0.6990
MATH1
0.6329
0.6514
0.7604
0.6731
0.6389
0.6082
0.6393
MATH2
0.6507
0.6705
0.6731
0.9069
0.6567
0.6251
0.6742
SCI
0.6759
0.6886
0.6389
0.6567
0.8316
0.6326
0.6779
WRITE
0.6399
0.6513
0.6082
0.6251
0.6326
0.7647
0.6556
CIVIC
0.6872
0.6990
0.6393
0.6742
0.6779
0.6556
0.9425From Table4 we observe that the posterior means for the variance and covariance components are each slightlypulled towards central quantities, respectively. To better illustrate this shrinkage property of the posterior mean we present the classical frequentist estimate of Σ in Table 5. If we make the element-wise comparison between the posterior mean in Table 4 and classical frequentist estimate in Table 5 we see that, among all diagonal elements (variances), relatively smaller elements of the classical frequentist estimate are pulled up and relatively larger elements of the classical frequentist estimate are pulled down. For example, the estimated variance for MATH 1 moved up from 64.9922 in Table 5 to 65.0307 in Table 4 and the estimated variance for CIVIC moved down from 80.8225 in Table 5 to 80.8041 in Table 4. Similar phenomenon appeared for the off-diagonal elements (covariances). This is due to the fact that we have assumed the intraclass matrix form for the prior specification of α.To further investigate the shrinkage property we considered an informative prior distribution forσ12 and σ22 instead of the vague prior specification of (12). In particular, we assumeda priori that σ12,σ22~iidInverseGamma5000,1. Table 7 presents the posterior mean for Σ with such informative conjugate prior specification. Under this informative prior specification the variance and covariance components are each pulled more towards central quantities, respectively, in comparison to the elements of Table 4. For example, the estimated variance for MATH 1 moved further up from 65.0307 in Table 4 to 65.7643 in Table 7 and the estimated variance for CIVIC moved further down from 80.8041 in Table 4 to 79.0518 in Table 7. Under the informative prior specification for σ12 and σ22 we observe that the shrinkage property is more pronounced.Table 7
Posterior mean forΣ with an informative prior.
VOCAB
READ
MATH 1
MATH 2
SCI
WRITE
CIVIC
VOCAB
72.3394
44.8244
34.4191
25.5068
38.8961
35.6451
32.5834
READ
44.8244
74.1193
37.4280
29.5260
40.5447
36.7774
34.1974
MATH 1
34.4191
37.4280
65.7643
39.5284
36.9010
34.3695
27.7634
MATH 2
25.5068
29.5260
39.5284
76.2078
28.8305
25.8843
20.4709
SCI
38.8961
40.5447
36.9010
28.8305
71.2920
34.7172
31.6945
WRITE
35.6451
36.7774
34.3695
25.8843
34.7172
65.8718
31.8856
CIVIC
32.5834
34.1974
27.7634
20.4709
31.6945
31.8856
79.0518Specifically, observe that the posterior estimates for the variance of the MATH 2 and CIVIC exams are less than their associated classical frequentist estimates, whereas, all others on diagonal elements of the posterior mean estimate forΣ are greater than their respective classical frequentist estimates. An analogous statement concerning the covariance terms can also be made. This draws our attention to the notion of the Bayesian posterior mean as a compromise between the prior information and the information contained in the data.


### 4.5. Posterior Estimates of the Sum Total

Of particular interest to educational testing data is some estimate of the overall summary or composite score of the individual subject area exams for a given set of explanatory variables. An obvious choice is to estimate the sum total of the individual exam scores. In a Bayesian framework, this can easily be accomplished by simply summing the individual posterior estimates of the seven subject area exams. However, the standard error associated with this estimate of the sum total cannot simply be calculated as the square root of the sum of the variance of the individual estimates. This would fail to incorporate the obvious covariance terms that exist amongst the subject area exams. Additionally, we would be overlooking any potential parameter uncertainty.Suppose for an individual we have thep×1 vector of response variables Yh=Yh1,…,YhpT and the associated k×1 vector of explanatory variables xh=xh1,…,xhkT. Note that the linear model for a single observation can be expressed as Yh=βTxh+ϵh. Furthermore, the sum total score is W=1TYh=∑j=1p‍Yhj, where 1=1,…,1T is a p×1 vector of ones. Then, given β and Σ, the total score W follows a univariate normal distribution with mean 1TβTxh and variance 1TΣ1.In order to more fully account for the added variability due to parameter uncertainty we consider the unconditional mean and variance of the sum total for an individual:(34)EW=Eβ,Σ∣YEW∣β,Σ=Eβ,Σ∣Y1TβTxi,Var⁡W=Eβ,Σ∣YVar⁡W∣β,Σ+Varβ,Σ∣YEW∣β,Σ=Eβ,Σ∣Y1TΣ1+Eβ,Σ∣Y1TβTxh2-Eβ,Σ∣Y1TβTxh2.
Here it is understood and the notation implies that the resulting expectations in (34) are in fact taken with respect to the posterior distributions of β and Σ, given Y. Equation (34) can be calculated based upon the results from the MCMC procedures in the following manner:
(35)E^h=1N∑t=lT‍1TβtTxh,V^h=1N∑t=lT‍1TΣt1+1N∑t=lT‍1TβtTxh2-E^h2,
where T is the total number of iterations of the MCMC algorithm, l is the number of burn-in iterations defined in Section 3, and N=T-l.Notice that by calculating the variance of the sum total score in this fashion we fully incorporate the obvious covariance structure amongst the response variables. Furthermore, in a Bayesian sense we capture the added variation due to parameter uncertainty.As a particular example we investigated a hypothetical student. The student is a female of Hispanic or Spanish decent. Her GPA is 3.0 and she attends a private nonelite school in a rural community of the Pacific region of the United States. Thus, for this particular example we have the following characterization of thek×1 vector of explanatory variables xh. xh1=xh10=xh12=xh20=xh21=xh26=1, xh2=3, and all the other remaining elements of xh are equal to zero.We obtained the following estimates of the Bayesian posterior means for the individual area exams, 50.9732, 51.2566, 51.5953, 52.1122, 50.0902, 53.7476, and 52.1576 for the VOCAB, READ, MATH 1, MATH 2, SCI, WRITE, and CIVIC exams, respectively. This results in an estimated sum total of 361.9328 with a standard error of 44.2357.


## 4.1. Description of Data

The HSB study contains a myriad of data and variables. In particular, for the sophomore class cohort a total ofp=7 exams were administered in the areas of vocabulary, reading, two exams in mathematics, science, writing, and civics. As is often the case with educational testing data the test scores were standardized to have a mean of fifty and a standard deviation equal to ten. We will use these standardized test scores for the sophomore class cohort as our multivariate response variables.Grade point average data was also obtained from official transcripts that were included in the survey. In addition, a number of demographic variables were collected on both the school and individual student levels. These variables included school type, relative urbanization of school environment, geographic region of the country, gender, and race. All of these variables will serve as categorical or qualitative explanatory variables in the multivariate multiple regression. Table1 provides a description of the categorical explanatory variables as well as the associated number of students per category. Note that the original HSB study did not include a school type four.Table 1
Categorical explanatory variable descriptions.
Categorical variables
Levels
Level description
Students per level
School type
School type 0
Regular public
9534
School type 1
Alternative
437
School type 2
Cuban/Hispanic public
156
School type 3
Other Hispanic public
1483
School type 5
Regular catholic
1341
School type 6
Black catholic
758
School type 7
Cuban/Hispanic catholic
216
School type 8
Private elite
294
School type 9
Private nonelite
448
Urbanization
Urban 1
Urban
3451
Urban 2
Suburban
7325
Urban 3
Rural
3891
Geographic region
Region 1
New England
736
Region 2
Middle Atlantic
2518
Region 3
South Atlantic
2178
Region 4
East South Central
713
Region 5
West South Central
1750
Region 6
East North Central
2931
Region 7
West North Central
1128
Region 8
Mountain
696
Region 9
Pacific
2017
Race
Race 1
Hispanic/Spanish
3122
Race 2
American Indian/Alaskan native
224
Race 3
Asian/Pacific Islander
345
Race 4
Black
2050
Race 5
White
8880
Race 6
Other
46
Gender
Gender 1
Male
7265
Gender 2
Female
7402As can be expected with educational data there was some moderate degree of missing data. We employed a data augmentation technique for missing data imputation. The specific procedure was invoked using themice library [21] in R. This particular data augmentation procedure fits nicely within the context of our research here since it employs a multivariate imputation by chained equations technique. That is, it uses a multivariate Gibbs sampler procedure to augment the missing data set. In particular, incomplete columns of data, that is, variables with missing data, are augmented by generating appropriate values of data given the values of the other columns of variables.


## 4.2. Treatment Contrast for Categorical Variables

Characterization of the explanatory data or design matrixX to account for categorical explanatory variables is not unique [22, page 173]. There are actually severalcontrast methods. In our analysis, school type zero, which corresponds to a regular public school, will serve as the base level and will not in fact have its own regressor. In the same fashion, urban will serve as the base level for the relative urbanization categorical explanatory variable. Also, the New England region will act as the base level for the geographic variable. Hispanic or Spanish is the base level for the race identifier variable. Finally, we will designate male as the base gender level. Thus, after properly accounting for the base levels in our particular application k=26.


## 4.3. Classical Multivariate Multiple Regression Results

We consider the standard multivatiate regression model described in (2) for the HSB survey data within the Sophomore class cohort. A total of n = 14,667 students in that cohort participated in the study. We estimated the model for the seven standardized exams VOCAB, READ, MATH 1, MATH 2, SCI, WRITE, and CIVIC regressed on GPA and the other explanatory variables described in the previous section. Most of the regression coefficient estimates are highly significant. In particular, GPA is highly significant for all seven exams. Although, certain particular levels for various categorical variables are not significant for some of the response variables. All of the associated F statistics are highly significant for all response variables. Moreover, based upon the analysis of variance that was performed we can conclude that nearly all the explanatory variables are highly significant for all seven subject area exams.


## 4.4. Posterior Estimates of the Model Parameters

In Tables2, 3, 4, 5, and 6 we present the Bayesian posterior means and the associated standard errors for the matrices β and Σ. In practice, for the MCMC procedure we found that a total iteration size of T = 500,000 was quite sufficient to establish convergence and we used l = 2,000 as the burn in value.Table 2
Posterior mean forβ.
VOCAB
READ
MATH 1
MATH 2
SCI
WRITE
CIVIC
(Intercept)
37.0582
35.6555
33.6377
37.7119
37.3335
31.2303
37.9283
GPA
4.7810
5.2711
5.9424
4.8139
4.4674
5.7163
4.3733
School type 1
−1.1971
−0.3382
−1.0521
−0.5196
−1.1388
−0.9455
−0.5326
School type 2
−0.8747
−1.9121
−0.6682
−0.4325
−2.4736
−1.4395
−2.4039
School type 3
−1.0713
−0.9146
−0.3638
−0.3816
−1.1779
−0.9105
−0.7742
School type 5
2.7359
2.1862
2.3550
1.5275
0.7704
2.5247
2.1433
School type 6
1.8095
1.8877
0.9502
0.2953
0.2920
2.0570
2.1941
School type 7
3.0532
1.6308
1.5283
0.0003
0.5107
2.1525
1.7716
School type 8
10.3236
8.7564
8.8743
8.7145
5.8165
7.0452
5.5020
School type 9
2.9475
1.9583
2.9381
2.9743
2.1760
1.9975
1.7927
Urban 2
0.5558
0.0821
0.6002
0.9824
0.7407
0.5115
0.2248
Urban 3
−0.8006
−0.5359
−0.3611
−0.2156
0.4344
0.0087
−0.3363
Region 2
−1.0396
−0.7905
−0.3661
−1.0051
−0.4256
−0.6938
−0.9750
Region 3
−2.2086
−1.1256
−1.6328
−1.7664
−1.4603
−1.1176
−1.7760
Region 4
−3.9350
−1.8827
−2.4900
−2.8335
−2.0910
−1.3902
−2.3088
Region 5
−2.6954
−1.5560
−1.5315
−2.1349
−1.1494
−0.8252
−1.6827
Region 6
−2.1648
−1.5579
−0.9690
−1.7886
−0.6775
−0.5598
−1.5524
Region 7
−1.9171
−0.4627
0.1163
−0.3902
0.2802
0.2396
−0.3970
Region 8
−1.6480
−0.9278
−1.6167
−2.4756
−0.6947
−0.4478
−1.4824
Region 9
−0.4866
−0.7586
−0.5338
−1.0931
−0.6208
−0.0413
−1.5015
Gender 2
−1.5605
−1.2414
−1.8966
−1.5490
−3.0722
3.0841
0.2048
Race 2
−0.9813
0.5513
−0.8029
−0.2627
0.5303
0.8640
−0.9432
Race 3
1.7286
1.7006
4.9363
3.9146
2.4132
3.5702
1.3489
Race 4
−0.8534
0.0337
−0.6765
−0.5625
−1.3975
−0.3582
0.4884
Race 5
4.5403
3.7077
3.8777
2.3674
4.5873
3.9272
2.7858
Race 6
−0.5276
0.3656
−0.0162
−0.1581
0.4371
0.3194
0.9498Table 3
Posterior standard deviations forβ.
VOCAB
READ
MATH 1
MATH 2
SCI
WRITE
CIVIC
(Intercept)
0.4421
0.4482
0.4195
0.4585
0.4383
0.4201
0.4678
GPA
0.0920
0.0930
0.0872
0.0951
0.0913
0.0874
0.0972
School type 1
0.4327
0.4379
0.4102
0.4483
0.4297
0.4110
0.4568
School type 2
0.7111
0.7199
0.6736
0.7352
0.7047
0.6761
0.7518
School type 3
0.2797
0.2823
0.2651
0.2902
0.2773
0.2662
0.2954
School type 5
0.2561
0.2589
0.2428
0.2659
0.2540
0.2438
0.2707
School type 6
0.3358
0.3392
0.3185
0.3480
0.3327
0.3193
0.3553
School type 7
0.6116
0.6187
0.5786
0.6331
0.6054
0.5806
0.6452
School type 8
0.5164
0.5229
0.4901
0.5359
0.5125
0.4912
0.5457
School type 9
0.4171
0.4222
0.3963
0.4325
0.4138
0.3972
0.4414
Urban 2
0.1864
0.1882
0.1764
0.1927
0.1844
0.1769
0.1964
Urban 3
0.2135
0.2158
0.2020
0.2206
0.2114
0.2027
0.2253
Region 2
0.3623
0.3672
0.3441
0.3759
0.3593
0.3438
0.3825
Region 3
0.3706
0.3759
0.3520
0.3848
0.3681
0.3524
0.3915
Region 4
0.4512
0.4574
0.4286
0.4678
0.4480
0.4290
0.4774
Region 5
0.3877
0.3928
0.3681
0.4020
0.3850
0.3686
0.4097
Region 6
0.3538
0.3587
0.3358
0.3669
0.3513
0.3365
0.3743
Region 7
0.4056
0.4103
0.3846
0.4208
0.4022
0.3858
0.4286
Region 8
0.4647
0.4709
0.4411
0.4819
0.4614
0.4426
0.4917
Region 9
0.3783
0.3830
0.3586
0.3925
0.3753
0.3598
0.3998
Gender 2
0.1421
0.1436
0.1347
0.1470
0.1411
0.1351
0.1500
Race 2
0.5944
0.6018
0.5630
0.6159
0.5897
0.5652
0.6273
Race 3
0.4938
0.4996
0.4679
0.5114
0.4903
0.4692
0.5216
Race 4
0.2599
0.2639
0.2467
0.2697
0.2579
0.2470
0.2752
Race 5
0.2061
0.2085
0.1952
0.2136
0.2042
0.1959
0.2175
Race 6
1.2666
1.2806
1.1995
1.3091
1.2548
1.2040
1.3359Table 4
Posterior mean forΣ.
VOCAB
READ
MATH 1
MATH 2
SCI
WRITE
CIVIC
VOCAB
72.3406
45.3069
34.3079
25.0339
39.1280
35.7151
32.6781
READ
45.3069
74.0616
37.4238
29.2657
40.8084
36.8560
34.3226
MATH 1
34.3079
37.4238
65.0307
39.9597
36.9083
34.2467
27.4186
MATH 2
25.0339
29.2657
39.9597
77.6078
28.6053
25.5057
19.7419
SCI
39.1280
40.8084
36.9083
28.6053
71.1835
34.7257
31.7008
WRITE
35.7151
36.8560
34.2467
25.5057
34.7257
65.3904
31.9590
CIVIC
32.6781
34.3226
27.4186
19.7419
31.7008
31.9590
80.8041Table 5
Classical frequentist estimate forΣ.
VOCAB
READ
MATH 1
MATH 2
SCI
WRITE
CIVIC
VOCAB
72.3308
45.3567
34.2973
24.9684
39.1487
35.7217
32.6627
READ
45.3567
74.0529
37.4264
29.2202
40.8335
36.8651
34.3098
MATH 1
34.2973
37.4264
64.9922
39.9886
36.9119
34.2393
27.3645
MATH 2
24.9684
29.2202
39.9886
77.6275
28.5646
25.4523
19.6396
SCI
39.1487
40.8335
36.9119
28.5646
71.1592
34.7254
31.6765
WRITE
35.7217
36.8651
34.2393
25.4523
34.7254
65.3498
31.9434
CIVIC
32.6627
34.3098
27.3645
19.6396
31.6765
31.9434
80.8225Table 6
Posterior standard deviations forΣ.
VOCAB
READ
MATH 1
MATH 2
SCI
WRITE
CIVIC
VOCAB
0.8456
0.7119
0.6329
0.6507
0.6759
0.6399
0.6872
READ
0.7119
0.8657
0.6514
0.6705
0.6886
0.6513
0.6990
MATH1
0.6329
0.6514
0.7604
0.6731
0.6389
0.6082
0.6393
MATH2
0.6507
0.6705
0.6731
0.9069
0.6567
0.6251
0.6742
SCI
0.6759
0.6886
0.6389
0.6567
0.8316
0.6326
0.6779
WRITE
0.6399
0.6513
0.6082
0.6251
0.6326
0.7647
0.6556
CIVIC
0.6872
0.6990
0.6393
0.6742
0.6779
0.6556
0.9425From Table4 we observe that the posterior means for the variance and covariance components are each slightlypulled towards central quantities, respectively. To better illustrate this shrinkage property of the posterior mean we present the classical frequentist estimate of Σ in Table 5. If we make the element-wise comparison between the posterior mean in Table 4 and classical frequentist estimate in Table 5 we see that, among all diagonal elements (variances), relatively smaller elements of the classical frequentist estimate are pulled up and relatively larger elements of the classical frequentist estimate are pulled down. For example, the estimated variance for MATH 1 moved up from 64.9922 in Table 5 to 65.0307 in Table 4 and the estimated variance for CIVIC moved down from 80.8225 in Table 5 to 80.8041 in Table 4. Similar phenomenon appeared for the off-diagonal elements (covariances). This is due to the fact that we have assumed the intraclass matrix form for the prior specification of α.To further investigate the shrinkage property we considered an informative prior distribution forσ12 and σ22 instead of the vague prior specification of (12). In particular, we assumeda priori that σ12,σ22~iidInverseGamma5000,1. Table 7 presents the posterior mean for Σ with such informative conjugate prior specification. Under this informative prior specification the variance and covariance components are each pulled more towards central quantities, respectively, in comparison to the elements of Table 4. For example, the estimated variance for MATH 1 moved further up from 65.0307 in Table 4 to 65.7643 in Table 7 and the estimated variance for CIVIC moved further down from 80.8041 in Table 4 to 79.0518 in Table 7. Under the informative prior specification for σ12 and σ22 we observe that the shrinkage property is more pronounced.Table 7
Posterior mean forΣ with an informative prior.
VOCAB
READ
MATH 1
MATH 2
SCI
WRITE
CIVIC
VOCAB
72.3394
44.8244
34.4191
25.5068
38.8961
35.6451
32.5834
READ
44.8244
74.1193
37.4280
29.5260
40.5447
36.7774
34.1974
MATH 1
34.4191
37.4280
65.7643
39.5284
36.9010
34.3695
27.7634
MATH 2
25.5068
29.5260
39.5284
76.2078
28.8305
25.8843
20.4709
SCI
38.8961
40.5447
36.9010
28.8305
71.2920
34.7172
31.6945
WRITE
35.6451
36.7774
34.3695
25.8843
34.7172
65.8718
31.8856
CIVIC
32.5834
34.1974
27.7634
20.4709
31.6945
31.8856
79.0518Specifically, observe that the posterior estimates for the variance of the MATH 2 and CIVIC exams are less than their associated classical frequentist estimates, whereas, all others on diagonal elements of the posterior mean estimate forΣ are greater than their respective classical frequentist estimates. An analogous statement concerning the covariance terms can also be made. This draws our attention to the notion of the Bayesian posterior mean as a compromise between the prior information and the information contained in the data.


## 4.5. Posterior Estimates of the Sum Total

Of particular interest to educational testing data is some estimate of the overall summary or composite score of the individual subject area exams for a given set of explanatory variables. An obvious choice is to estimate the sum total of the individual exam scores. In a Bayesian framework, this can easily be accomplished by simply summing the individual posterior estimates of the seven subject area exams. However, the standard error associated with this estimate of the sum total cannot simply be calculated as the square root of the sum of the variance of the individual estimates. This would fail to incorporate the obvious covariance terms that exist amongst the subject area exams. Additionally, we would be overlooking any potential parameter uncertainty.Suppose for an individual we have thep×1 vector of response variables Yh=Yh1,…,YhpT and the associated k×1 vector of explanatory variables xh=xh1,…,xhkT. Note that the linear model for a single observation can be expressed as Yh=βTxh+ϵh. Furthermore, the sum total score is W=1TYh=∑j=1p‍Yhj, where 1=1,…,1T is a p×1 vector of ones. Then, given β and Σ, the total score W follows a univariate normal distribution with mean 1TβTxh and variance 1TΣ1.In order to more fully account for the added variability due to parameter uncertainty we consider the unconditional mean and variance of the sum total for an individual:(34)EW=Eβ,Σ∣YEW∣β,Σ=Eβ,Σ∣Y1TβTxi,Var⁡W=Eβ,Σ∣YVar⁡W∣β,Σ+Varβ,Σ∣YEW∣β,Σ=Eβ,Σ∣Y1TΣ1+Eβ,Σ∣Y1TβTxh2-Eβ,Σ∣Y1TβTxh2.
Here it is understood and the notation implies that the resulting expectations in (34) are in fact taken with respect to the posterior distributions of β and Σ, given Y. Equation (34) can be calculated based upon the results from the MCMC procedures in the following manner:
(35)E^h=1N∑t=lT‍1TβtTxh,V^h=1N∑t=lT‍1TΣt1+1N∑t=lT‍1TβtTxh2-E^h2,
where T is the total number of iterations of the MCMC algorithm, l is the number of burn-in iterations defined in Section 3, and N=T-l.Notice that by calculating the variance of the sum total score in this fashion we fully incorporate the obvious covariance structure amongst the response variables. Furthermore, in a Bayesian sense we capture the added variation due to parameter uncertainty.As a particular example we investigated a hypothetical student. The student is a female of Hispanic or Spanish decent. Her GPA is 3.0 and she attends a private nonelite school in a rural community of the Pacific region of the United States. Thus, for this particular example we have the following characterization of thek×1 vector of explanatory variables xh. xh1=xh10=xh12=xh20=xh21=xh26=1, xh2=3, and all the other remaining elements of xh are equal to zero.We obtained the following estimates of the Bayesian posterior means for the individual area exams, 50.9732, 51.2566, 51.5953, 52.1122, 50.0902, 53.7476, and 52.1576 for the VOCAB, READ, MATH 1, MATH 2, SCI, WRITE, and CIVIC exams, respectively. This results in an estimated sum total of 361.9328 with a standard error of 44.2357.


## 5. Conclusion

In conclusion, we have demonstrated how a flexible prior specification for the covariance structure of a multivariate multiple regression model can provide a richer class of distributions than the inverse Wishart family. We discussed how the likelihood function for the covariance structure can be approximated based upon Bellman’s solution of a linear Volterra integral equation. We discussed the shrinkage properties of the posterior mean of the covariance structure. This highlighted the concept of the posterior means as a compromise between the prior information and the information contained in the data.All posterior estimates were calculated based upon the numerical results of a MHWG procedure. The Metropolis-Hastings algorithm was employed to account for sampling from an approximate posterior distribution for the covariance structure.



---

*Source: 673657-2014-11-24.xml*