Bayesian Reduced Rank Regression in Econometrics

Caleb Jin / 2018-05-09

Foreword

After I read this paper, Bayesian Reduced Rank Regression in Econometrics (Geweke 1996), I write down the nodes of the key idea and R code to realize it.

Introduction

In contemporary sociey, a big amount of data are generated and collected more easily and routinely in various academic and industrial areas, such as engineering, politics, B2C e-ccommerce, genomics, etc. Many problems could be cast into statistical problems under the framework of multivariate linear regression model, which is characterized by that both reponse variables and predictors are high dimensionality.

We assume $n$ independent observations of the response $y_{i} \in R^{q}$ with predictor vector $x_{i} \in R^{p}$ , $i = 1, 2, \dots, n$ . Cosider multivariate linear regression model as follows: $\begin{array}{rcl} (1) & Y = X C + E, \end{array}$ where $Y = (y_{1}, y_{2}, \dots, y_{n})^{⊤} \in R^{n \times q}$ is the response matrix, $X = (x_{1}, x_{2}, \dots, x_{n})^{⊤} \in R^{n \times p}$ is the design matrix, $C \in R^{p \times q}$ is the coefficient matrix, and $E = (e_{1}, e_{2}, \dots, e_{n})^{⊤} \in R^{n \times q}$ is the disturbance matrix with $e_{i}$ ’s $\overset{i i d}{\sim} N_{q} (0, Σ_{e})$ . We assume $Σ_{e} = σ^{2} I_{q}$ . Therefore, we have $Y \sim M N (X C, Σ_{e}, I_{n})$ .

Bayesian Rank Reduced Regression model

We consider to decompose the coefficient matrix into two part as follows: $\begin{array}{rcl} C = A B^{⊤}, \end{array}$ where $A \in R^{p \times r}$ , $B \in R^{r \times q}$ and known $r \leq min (p, q)$ . However, this decomposition is not unique, because with a $r \times r$ nonsingular matrix $Q, C = A B^{⊤} = A Q Q^{- 1} B^{⊤} = \tilde{A} {\tilde{B}}^{⊤} .$ In order to indentify it, we further decompose $A$ . $A^{⊤} = [I_{r}, {A^{*}}^{⊤}]$ . The author assumes that $p (A, B) \propto \exp (- \frac{τ^{2}}{2} t r a c e {A^{⊤} A + B^{⊤} B})$ and $σ^{2} \sim I G (\frac{a}{2}, \frac{b}{2})$ .

Posterior Distribution

Let ${\tilde{a}}_{k}$ and ${\tilde{b}}_{k}$ denote the $k$ th column of $A$ and $B$ , respectively. Let $a_{j}^{⊤}$ and $b_{l}^{⊤}$ denote the $y$ th row of $A$ and $l$ th row of $B$ , respectively.

$\begin{array}{rcl} p (Y | A, B, Σ_{e}) \\ \propto & | Σ_{e} |^{- \frac{q}{2}} \exp (- \frac{1}{2} t r a c e {(Y - X A B^{⊤}) (Y - X A B^{⊤})^{⊤} Σ_{e}^{- 1}}) \\ = & (σ^{2})^{- \frac{n q}{2}} \exp (- \frac{1}{2 σ^{2}} t r a c e {(Y - X A B^{⊤}) (Y - X A B^{⊤})^{⊤}}) \\ = & (σ^{2})^{- \frac{n q}{2}} \exp (- \frac{1}{2 σ^{2}} t r a c e {(Y - X_{(\tilde{j})} A_{(j)} B^{⊤}) (Y - X_{(\tilde{j})} A_{(j)} B^{⊤})^{⊤}}) \\ \times & \exp (- \frac{1}{2 σ^{2}} t r a c e {({\tilde{x}}_{j} a_{j}^{⊤} B^{⊤} B a_{j} {\tilde{x}}_{j}^{⊤})}) \\ \times & \exp (\frac{1}{σ^{2}} t r a c e {({\tilde{x}}_{j} a_{j}^{⊤} B^{⊤} (Y - X_{(\tilde{j})} A_{(j)} B^{⊤})^{⊤})}) \\ = & (σ^{2})^{- \frac{n q}{2}} \exp (- \frac{1}{2 σ^{2}} t r a c e {(Y - X_{(\tilde{j})} A_{(j)} B^{⊤}) (Y - X_{(\tilde{j})} A_{(j)} B^{⊤})^{⊤}}) \\ \times & \exp (- \frac{1}{2 σ^{2}} a_{j}^{⊤} B^{⊤} B a_{j} {\tilde{x}}_{j}^{⊤} {\tilde{x}}_{j}) \times \exp ((- 2) \frac{- 1}{2 σ^{2}} a_{j}^{⊤} B^{⊤} (Y - X_{(\tilde{j})} A_{(j)} B^{⊤})^{⊤} {\tilde{x}}_{j}) . \end{array}$

$\begin{array}{rcl} p (a_{j} | A_{(j)}, B, Σ_{e}, Y) \\ \propto & p (A | B, Σ_{e}, Y) \propto p (Y | A, B, Σ_{e}) p (A, B) \\ \propto & p (Y | A, B, Σ_{e}) \exp (- \frac{τ^{2}}{2} t r a c e {A^{⊤} A + B^{⊤} B}) \\ \propto & p (Y | A, B, Σ_{e}) \exp (- \frac{τ^{2}}{2} \sum_{j = 1}^{r} a_{j}^{⊤} a_{j}) \\ \propto & \exp (- \frac{1}{2 σ^{2}} a_{j}^{⊤} B^{⊤} B a_{j} {\tilde{x}}_{j}^{⊤} {\tilde{x}}_{j}) \times \exp ((- 2) \frac{- 1}{2 σ^{2}} a_{j}^{⊤} B^{⊤} (Y - X_{(\tilde{j})} A_{(j)} B^{⊤})^{⊤} {\tilde{x}}_{j}) \\ \times & \exp (- \frac{τ^{2}}{2} a_{j}^{⊤} a_{j}) \\ = & \exp {- \frac{1}{2} (a_{j}^{⊤} (σ^{- 2} B^{⊤} B {\tilde{x}}_{j}^{⊤} {\tilde{x}}_{j} + I_{r} τ^{2}) a_{j} - 2 a_{j}^{⊤} B^{⊤} (Y - X_{(\tilde{j})} A_{(j)} B^{⊤})^{⊤} {\tilde{x}}_{j} σ^{- 2})} \\ = & \exp {- \frac{1}{2} (a_{j}^{⊤} {Σ_{j}^{A}}^{- 1} a_{j} - 2 a_{j}^{⊤} {Σ_{j}^{A}}^{- 1} Σ_{j}^{A} B^{⊤} (Y - X_{(\tilde{j})} A_{(j)} B^{⊤})^{⊤} {\tilde{x}}_{j} σ^{- 2})} \\ = & \exp {- \frac{1}{2} (a_{j}^{⊤} {Σ_{j}^{A}}^{- 1} a_{j} - 2 a_{j}^{⊤} {Σ_{j}^{A}}^{- 1} μ_{j}^{A})} \\ \propto & \exp {- \frac{1}{2} ((a_{j} - μ_{j}^{A})^{⊤} {Σ_{j}^{A}}^{- 1} (a_{j} - μ_{j}^{A}))}, \end{array}$ where $Σ_{j}^{A} = {(B^{⊤} B {\tilde{x}}_{j}^{⊤} {\tilde{x}}_{j} σ^{- 2} + I_{r} τ^{2})}^{- 1}$ and $μ_{j}^{A} = Σ_{j}^{A} B^{⊤} (Y - X_{(\tilde{j})} A_{(j)} B^{⊤})^{⊤} {\tilde{x}}_{j} σ^{- 2}$ .

Hence, $a_{j} | A_{(j)}, B, Σ_{e}, Y \sim N_{r} (μ_{j}^{A}, Σ_{j}^{A})$ . $\begin{array}{rcl} p (Y | A, B, Σ_{e}) \\ \propto & | Σ_{e} |^{- \frac{q}{2}} \exp (- \frac{1}{2} t r a c e {(Y - X A B^{⊤})^{⊤} Σ_{e}^{- 1} (Y - X A B^{⊤})}) \\ = & (σ^{2})^{- \frac{n q}{2}} \exp (- \frac{1}{2 σ^{2}} t r a c e {(Y - X A B^{⊤})^{⊤} (Y - X A B^{⊤})}) \\ = & (σ^{2})^{- \frac{n q}{2}} \exp (- \frac{1}{2 σ^{2}} t r a c e {(Y_{(\tilde{l})} - X A (B^{⊤})_{(\tilde{l})})^{⊤} (Y_{(\tilde{l})} - X A (B^{⊤})_{(\tilde{l})})}) \\ \times & \exp (- \frac{1}{2 σ^{2}} [b_{l}^{⊤} (X A)^{⊤} X A b_{l} - 2 b_{l}^{⊤} (X A)^{⊤} {\tilde{y}}_{l} + {\tilde{y}}_{l}^{⊤} {\tilde{y}}_{l}]) . \end{array}$

$\begin{array}{rcl} p (b_{l} | A, (B^{⊤})_{(\tilde{l})}, Y, Σ_{e}) \\ \propto & p (B | A, Y, Σ_{e}) \propto p (Y | A, B, Σ_{e}) p (A, B) \\ \propto & p (Y | A, B, Σ_{e}) \exp (- \frac{τ^{2}}{2} t r a c e {A^{⊤} A + B^{⊤} B}) \\ \propto & p (Y | A, B, Σ_{e}) \exp (- \frac{τ^{2}}{2} \sum_{l = 1}^{r} b_{l}^{⊤} b_{l}) \\ \propto & \exp (- \frac{1}{2 σ^{2}} [b_{l}^{⊤} (X A)^{⊤} X A b_{l} - 2 b_{l}^{⊤} (X A)^{⊤} {\tilde{y}}_{l} + {\tilde{y}}_{l}^{⊤} {\tilde{y}}_{l}]) \exp (- \frac{τ^{2}}{2} b_{l}^{⊤} b_{l}) \\ = & \exp (- \frac{1}{2} [b_{l}^{⊤} (X A)^{⊤} X A b_{l} σ^{- 2} - 2 b_{l}^{⊤} (X A)^{⊤} {\tilde{y}}_{l} σ^{- 2} + {\tilde{y}}_{l}^{⊤} {\tilde{y}}_{l} σ^{- 2}]) \exp (- \frac{τ^{2}}{2} b_{l}^{⊤} b_{l}) \\ = & \exp {- \frac{1}{2} (b_{l}^{⊤} ((X A)^{⊤} X A σ^{- 2} + I_{r} τ^{2}) b_{l} - 2 b_{l}^{⊤} (X A)^{⊤} {\tilde{y}}_{l} σ^{- 2})} \\ = & \exp {- \frac{1}{2} (b_{l}^{⊤} {Σ_{j}^{B}}^{- 1} b_{l} - 2 b_{l}^{⊤} {Σ_{j}^{B}}^{- 1} Σ_{j}^{B} (X A)^{⊤} {\tilde{y}}_{l} σ^{- 2})} \\ = & \exp {- \frac{1}{2} (b_{l}^{⊤} {Σ_{j}^{B}}^{- 1} b_{l} - 2 b_{l}^{⊤} {Σ_{j}^{B}}^{- 1} μ_{j}^{B})} \\ \propto & \exp {- \frac{1}{2} ((b_{l} - μ_{j}^{B})^{⊤} {Σ_{j}^{B}}^{- 1} (b_{l} - μ_{j}^{B}))}, \end{array}$ where $Σ_{j}^{B} = ((X A)^{⊤} X A σ^{- 2} + I_{r} τ^{2})^{- 1}$ and $μ_{j}^{B} = Σ_{j}^{B} (X A)^{⊤} {\tilde{y}}_{l} σ^{- 2}$ .

Hence, $b_{l} | A, (B^{⊤})_{(\tilde{l})}, Y, Σ_{e} \sim N_{r} (μ_{j}^{B}, Σ_{j}^{B})$ .

We know that the element in $k$ th row and $k$ th column of $Σ_{e}$ is $σ^{2}, k = 1, 2, \dots, q$ . $\begin{array}{rcl} p (σ^{2} | A, B, Y, (Σ_{e})_{(k) (\tilde{k})}) \\ \propto & p (Σ_{e} | A, B, Y) \propto p (Y | A, B, Σ_{e}) p (Σ_{e}) \propto p (Y | A, B, Σ_{e}) p (σ^{2}) \\ = & (σ^{2})^{- \frac{n q}{2}} \exp (- \frac{1}{2 σ^{2}} t r a c e {(Y - X A B^{⊤}) (Y - X A B^{⊤})^{⊤}}) \\ \times & (σ^{2})^{- \frac{a}{2} - 1} \exp (- \frac{b}{2 σ^{2}}) \\ = & (σ^{2})^{- \frac{n q + a}{2} - 1} \exp (- \frac{1}{2 σ^{2}} (t r a c e {(Y - X A B^{⊤}) (Y - X A B^{⊤})^{⊤}} + b)) . \end{array}$ Hence, $σ^{2} | A, B, Y, (Σ_{e})_{(k) (\tilde{k})} \sim I G (\frac{n q + a}{2}, \frac{1}{2} (t r a c e {(Y - X A B^{⊤}) (Y - X A B^{⊤})^{⊤}} + b))$ .

Gibbs Sampling

The algorithm is easy to construct. Begin with initial values $A^{(0)}, Σ^{(0)}$ , then for t = 1,2,…

Step (1) Given $A_{(j)}^{(t - 1)}, B^{⊤ (t - 1)}, {σ^{2}}^{(t - 1)}$ , draw $a_{j}^{(t)}$ from $N_{r} (μ_{j}^{A}, Σ_{j}^{A}), j = r + 1, r + 2, \dots, p$ ;
Step (2) Given $A^{(t)}, B_{(\tilde{l})}^{⊤ (t - 1)}, {σ^{2}}^{(t - 1)}$ , draw $b_{l}^{(t)}$ from $N_{r} (μ_{j}^{B}, Σ_{j}^{B}), l = 1, 2, \dots, q$ ;
Step (3) Given $B^{(t)}, A^{(t)}$ , draw ${σ^{2}}^{(t)}$ from $I G (a^{*}, b^{*})$ .
Step (4) Repeat step (1) to step (3) until convergence.

Note that, the step (1) samples the $a_{j}$ from $j = r + 1$ , because first part of $A$ is $I_{r}$ , which is known. Hence, step (1) samples the $A^{*}$ and then insert $I_{r}$ together to construct $A$ .

Model Selection

To this point we have proceeded as if r were known. In most applications this will not be true and so analysis to this point is conditional.When r is unknown, the analysis may be carried out for several alternative values of r. Under the Bayesian framework, model selection can be implemented by using the model posterior probability conditioning the data, $p (M_{r} | y) = \frac{p (y | M_{r}) p (M_{r})}{\sum_{r \in M} p (y | M_{r}) p (M_{r})} = \frac{p (y | M_{r})}{\sum_{r \in M} p (y | M_{r})},$ where $M = 1, 2, \dots, m i n (p, q)$ , and $p (M_{r}) = \frac{1}{c a r d (M)}$ .

Hence, as the prior of rank(model) is flat, $p (M_{r} | y)$ is only determined by marginal likelihood ( $p (y | M_{r})$ ). The model selection problem here is then converted to find out the rank which maximizes the marginal likelihood ( $p (y | M_{r})$ ). In order to calculate the $p (y | M_{r})$ , I used Laplace and Gelfand and Dey (GD) methods, and then compare with DIC.

Laplace

Let $θ = (A, B, σ^{2})$ $\begin{array}{rcl} p (Y | M_{r}) & = & \int \int \int p (Y | θ, M) p (θ) d A d B d σ^{2} \\ \approx & p (Y | \hat{θ}, M) p (\hat{θ} | M) | (n q)^{- 1} {\hat{Σ}}_{M} |^{1 / 2} (2 π)^{(k_{M} / 2)}, \end{array}$ where $(\hat{A}, \hat{B}, {\hat{σ}}^{2}) = \arg max p (Y | A, B, σ^{2}, M) p (A, B, σ^{2} | M)$ .

Hence, $\begin{array}{rcl} \log p (Y | M_{r}) \\ \approx & \log p (Y | \hat{θ}, M) + \log p (\hat{θ} | M) - \frac{1}{2} k_{M} \log n q + \frac{1}{2} | Σ_{M} | + \frac{k_{M}}{2} \log 2 π \\ = & - \frac{1}{2} (- 2 \log (p (Y | \hat{A}, \hat{B}, \hat{σ^{2}}, M)) + k_{M} \log n q) + C \\ = & - \frac{1}{2} (n q \log (2 π \hat{σ^{2}}) - \frac{1}{2 \hat{σ^{2}}} t r a c e {(Y - X \hat{C})^{⊤} (Y - X \hat{C})} \\ + [r (p + q - r) + 1] \log (n q)) + C \\ = & - \frac{1}{2} B I C + C \\ (2) & = & - \frac{1}{2} B I C, as n \to \infty . \end{array}$

The problem we are facing when we use Laplace here is that we need to find the mode of $p (A, B, σ^{2} | Y)$ , in which $(A, B, σ^{2})$ is high-dimensional. To address this problem, Besag proposed an iterated conditional modes (ICM) algorithm. The ICM obtains the local maximum of the joint posterior by iteratively maximizing the full conditionals as follows:

Begin with initial values $A^{(0)}, Σ^{(0)}$ , then for t = 1,2,…
Given $A_{(j)}^{(t - 1)}, B^{⊤ (t - 1)}, {σ^{2}}^{(t - 1)}$ , $a_{j}^{(t)} \leftarrow μ_{j}^{A}, j = r + 1, r + 2, \dots, p$ ;
Given $A^{(t)}, B_{(\tilde{l})}^{⊤ (t - 1)}, {σ^{2}}^{(t - 1)}$ , $b_{l}^{(t)} \leftarrow μ_{j}^{B}, l = 1, 2, \dots, q$ ;
Given $B^{(t)}, A^{(t)}$ , ${σ^{2}}^{(t)} \leftarrow \frac{b^{*}}{a^{*} - 1}$ .

We obtain $C^{(t)} = A^{(t)} B^{⊤ (t)}$ and estimate ${\hat{C}}_{I C M} = T^{- 1} \sum_{t = 1}^{T} C^{(t)}$ . Put the ${\hat{C}}_{I C M}$ in the Eq.(2) to calculate $- \frac{1}{2} B I C$ . Hence, the Laplace here is actually propotional to BIC.

DIC

Let $θ = (A, B, σ^{2})$ , $\begin{array}{rcl} D (θ) \\ = & - 2 \log p (Y | θ, M) \\ = & n q \log (2 π σ^{2}) + \frac{1}{σ^{2}} t r a c e {(Y - X C)^{⊤} (Y - X C)} \end{array}$

Then DIC can be computed by $D I C \approx 2 T^{- 1} \sum_{t = 1}^{T} D (θ^{(t)}) - D (T^{- 1} \sum_{t = 1}^{T} θ^{(t)}),$ where $θ^{(t)}$ is a MCMC sample generated from $p (θ | Y)$ .

Note that $D I C = D (\bar{θ}) + 2 P_{D}$ , which is analogous to AIC.

Gelfand & Dey (GD)

Let $θ = (A, B, σ^{2})$ ,then the GD estimator is $p (Y | M) \approx {[T^{- 1} \sum_{t = 1}^{T} \frac{g (θ^{(t)})}{p (Y | θ^{(t)}) p (θ^{(t)})}]}^{- 1},$ where $θ^{(t)}$ is a MCMC sample generated from $p (θ | Y)$ . Define $g (θ) = N (θ | \tilde{θ}, \tilde{Σ})$ , where $\tilde{θ}$ and $\tilde{Σ}$ are MCMC sample mean and variance, respectively. I use formular as follows: $\begin{array}{rcl} \log p (Y | M) & \approx & \log {[\frac{1}{T} \sum_{t = 1}^{T} \frac{g^{(t)}}{f^{(t)}}]}^{- 1} \\ = & \log [\frac{T}{\sum_{t = 1}^{T} \frac{g^{(t)}}{f^{(t)}}}] \\ = & \log T - \log (\sum_{t = 1}^{T} \frac{g^{(t)}}{f^{(t)}}) \\ (3) & = & \log T - \log (\sum_{t = 1}^{T} \exp (\log g^{(t)} - \log f^{(t)})) \\ = & \log T - \log (\sum_{t = 1}^{T} \exp c_{t}) \\ = & \log T - \log (e^{c_{1}} e^{\sum_{t = 1}^{T} (c_{t} - c_{1})}) \\ = & \log T - [c_{1} + \log (\sum_{t = 2}^{T} (c_{t} - c_{1}))], \end{array}$ where $g^{(t)} = g (θ), f^{(t)} = p (Y | θ^{(t)}) p (θ^{(t)}), c_{t} = \log g^{(t)} - \log f^{(t)}$ .

Note that when calculating $\log p (Y | M)$ , there is a computation problem in formula (3), $\exp (\log g^{(t)} - \log f^{(t)})$ goes to infinity, due to a very large magnitude of $\log g^{(t)} - \log f^{(t)}$ . To solve this problem, I used $\log T - [c_{1} + \log (\sum_{t = 2}^{T} (c_{t} - c_{1}))]$ to calculate $\log p (Y | M)$ .

Simulation Study

Data Generation

In the simulation study, my goal is to find out the true model or true rank of coefficient matrix( $C$ ) based on Laplace, DIC and GD method. I set $n = 100, q = 12, p = 7$ and $σ^{2} = 2$ . The coefficient matrix is as follows:

$C = [\begin{matrix} 1 & 0 & 0 & 2 & - 1 & 0 & 0 & 0 & 0 & 0 & 1 & - 1 \\ 0 & 1 & 0 & 0 & 0 & - 3 & 2 & 0 & 0 & 0 & - 1 & 3 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 3 & - 3 & 4 & 2 & 2 \\ 1 & 0 & 0 & 2 & - 1 & 0 & 0 & 0 & 0 & 0 & 1 & - 1 \\ 0 & 1 & 0 & 0 & 0 & - 3 & 2 & 0 & 0 & 0 & - 1 & 3 \\ 0 & 1 & 0 & 0 & 0 & - 3 & 2 & 0 & 0 & 0 & - 1 & 3 \\ 0 & 0 & 1 & 0 & 0 & 0 & 0 & 3 & - 3 & 4 & 2 & 2 \end{matrix}]$ Hence, the true rank of $C$ is 3.

Generate data based on the model from Eq.(1). For the prior, $p (A, B) \propto \exp (- \frac{τ^{2}}{2} t r a c e {A^{⊤} A + B^{⊤} B})$ and $σ^{2} \sim I G (\frac{a}{2}, \frac{b}{2})$ , where $τ = 10^{- 3}$ and $a = b = 1$ .

After running the MCMC simulation, the result based on 1 replication is shown in Table 1.

Table 1: Model Selection among Laplace, DIC and GD Based on 1 Replication.
Method( $\log$ )	1	2	3	4	5	6	7
Laplace	-3007	-2535	-2223	-2257	-2285	-2310	-2329
DIC	84300	17779	7302	7317	7325	7363	7405
GD	-3062	-2631	-2347	-2397	-2436	-2470	-2495
MSE	1.032	0.186	0.0120	0.014	0.0170	0.0192	0.0217

The MSE in the table is defined as follows: $M S E = \frac{t r a c e {(\hat{C} - C)^{⊤} (\hat{C} - C)}}{p q} .$ MSE indicates the goodness of fit in the reduced rank regression model. We can observe that the MSE is minimized when the rank is 3. This demonstrates that parameter estimation is good enough. For the Laplace, DIC and GD method, in this case, all of them select the true rank. Given rank $(C) = 3$ , the estimated $C$ is as follows:

1.07	0.17	0.01	2.34	-0.89	-0.12	-0.11	-0.08	-0.07	0.03	1.26	-1.10
-0.02	1.00	0.04	-0.14	-0.17	-2.83	2.10	0.10	0.03	0.03	-1.04	3.17
0.13	-0.03	0.99	0.18	-0.02	0.08	0.00	3.06	-2.85	4.03	2.10	2.06
0.94	0.14	0.03	2.05	-0.78	-0.08	-0.11	-0.03	-0.10	0.07	1.13	-0.95
0.04	0.99	0.04	-0.00	-0.21	-2.79	2.05	0.10	0.01	0.04	-0.95	3.05
0.12	1.05	0.00	0.18	-0.29	-2.93	2.14	-0.02	0.12	-0.12	-0.98	3.03
0.08	-0.01	0.98	0.08	0.01	0.03	0.04	3.03	-2.81	3.98	2.01	2.14

In addition, when selected rank is larger than 3, the MSE is not worse than I expect. In other words, the model based on larger model (larger rank in $C$ ), the overfiteed models still have a good estimation for the $C$ . The values of DIC, Laplace and GD are not far away from the value at true rank. This is very reasonable for overfitted model. But when selected rank smaller than true rank, the MSE increase significantly, and the corresponding values of DIC, Laplace and GD are far away from the value at true rank. This means the model miss some important variables. The analysis is based on one replication.

In Table 2, the result is based on 1000 replication. I want to see the perfermance of methods for model selection. In Table 2, Laplace and GD always select true rank, however, the DIC tends to select larger model, because the average of selected rank in DIC is about 3.5, whereas, that of Laplace and GD are 3. Besides, for the successful probability of selecting true rank, the Laplace and GD, of course, is 100%, but is 61.7% for DIC. This makes sense because Laplace and GD are approximated and exact calculation for marginal likehihood, respectively. They work well in model fitting. The DIC is analogous to AIC, which tends to select larger model and better for short-term prediction.

Table 2: Comparison among Laplace, DIC, GD Based on 1000 Replication.
Method	Mean of Selected rank	Selection Probabiliy
Laplace	3	1
DIC	3.525	0.617
GD	3	1

R Code for Bayeian reduced rank regression with DIC used in the simulation study

library(mvtnorm)
library(invgamma)
rm(list=ls())
## Data & Model##
row1 <- c(1,0,0,2,-1,0,0,0,0,0,1,-1)
row2 <- c(0,1,0,0,0,-3,2,0,0,0,-1,3)
row3 <- c(0,0,1,0,0,0,0,3,-3,4,2,2)
C <- rbind(row1,row2,row3,row1,row2,row2,row3)
true.rank <- 3
p <- dim(C)[1]
q <- dim(C)[2]
true.sig2 <- 2
true.SIGe <- diag(true.sig2,q)
n=100
a=b=1
tau2=1e-3
MC.size=5000
burn_in=3000
reptn <- 125
error <- array(NA,dim = c(reptn,p))
DIC <- array(NA,dim = c(reptn,p))
TPR.DIC <- rep(NA,reptn)
TPR.error <- rep(NA,reptn)
hat.rank <- rep(NA,reptn)
v=0
for (i in 1:reptn) {
  t1=Sys.time()
  set.seed(2144+i+125*v)
  X <- matrix(rnorm(n*p,0,1),n,p)
  E <- rmvnorm(n, rep(0,q), true.SIGe, method="chol")
  Y=X%*%C+E
  ## Reduced Rank Regression model ##
  for (r in 1:p) {
    Hat.A <- array(NA,dim = c(p,r,MC.size))
    Hat.B <- array(NA,dim = c(q,r,MC.size))
    Hat.C <- array(NA,dim = c(p,q,MC.size))
    Hat.sig2 <- rep(NA,MC.size)
    # initial values
    hat.C <- coef(lm(Y~X-1))
    if (r==1){
      hat.B <- as.matrix(hat.C[1:r,])
    } else {
      hat.B <- t(hat.C[1:r,])
    }
    if (r==p) {
      hat.A <- diag(1,r)
    } else {
      hat.A <- rbind(diag(1,r),hat.C[(r+1):p,]%*%hat.B%*%solve(crossprod(hat.B)))
    }

    hat.sig2 <- mean(diag(tcrossprod(Y-X%*%hat.A%*%t(hat.B))))  #true.sig2
    ## MCMC ##
    for (jin in 1:MC.size) {
      # Sampling A
      if (r==p){
        Hat.A[,,jin] <- hat.A
      } else {
        for(j in (r+1):p) {
          Sig.A_j <- solve(crossprod(hat.B)*as.numeric(crossprod(X[,j]))/hat.sig2+diag(tau2,r))
          mu.A_j <- Sig.A_j%*%t(hat.B)%*%t(Y-X[,-j]%*%hat.A[-j,]%*%t(hat.B))%*%X[,j]/hat.sig2
          hat.a_j <- rmvnorm(n=1,mean = mu.A_j,sigma = Sig.A_j)
          hat.A[j,] <- hat.a_j
        }
        Hat.A[,,jin] <- hat.A
      }
      # Sampling B
      Sig.B_l <- solve(crossprod(X%*%hat.A)/hat.sig2+diag(tau2,r))
      for(l in 1:q) {
        mu.B_l <- Sig.B_l%*%t(X%*%hat.A)%*%Y[,l]/hat.sig2
        hat.b_l <- t(rmvnorm(n=1,mean = mu.B_l,sigma = Sig.B_l))
        hat.B[l,] <- t(hat.b_l)
      }
      Hat.B[,,jin] <- hat.B
      Hat.C[,,jin] <- hat.A%*%t(hat.B)
      #Sampling Sig2
      hat.sig2 <- rinvgamma(n = 1,shape = (q*n+a)/2,rate = 0.5*(b+sum(diag(tcrossprod(Y-X%*%hat.A%*%t(hat.B))))))
      Hat.sig2[jin] <- hat.sig2
      #print(jin)
    }
    hat.C <- apply(Hat.C[,,-(1:burn_in)],c(1,2),mean)
    hat.sig2 <- mean(Hat.sig2[-(1:burn_in)])
    # print(C)
    D.theta <- rep(NA,MC.size-burn_in)
    j=0
    for (jin in (burn_in+1):MC.size) {
      j=j+1
      D.theta[j] <- n*q*log(2*pi*Hat.sig2[jin])+Hat.sig2[jin]*sum(diag(crossprod(Y-X%*%Hat.C[,,jin])))
    }
    D.hat.theta <- n*q*log(2*pi*hat.sig2)+hat.sig2*sum(diag(crossprod(Y-X%*%hat.C)))
    DIC[i,r] <- 2*mean(D.theta)- D.hat.theta
    error[i,r] <- sum((hat.C-C)^2)/(p*q)
  }
  hat.rank[i] <- which.min(DIC[i,])
  TPR.DIC[i] <- hat.rank[i]==true.rank
  TPR.error[i] <- which.min(error[i,])==true.rank
  print(c(i,hat.rank[i]))
  print(Sys.time()-t1)
}

data.frame(hat.rank=mean(hat.rank),TPR.DIC=mean(TPR.DIC),TPR.ERR=mean(TPR.error))
save(list = ls(), file = paste0('DIC-',v,'.RData'))

Reference

Geweke, John F. 1996. “Bayesian Reduced Rank Regression in Econometrics.” Journal of Econometrics 75 (November). https://doi.org/10.1016/0304-4076(95)01773-9.