# Modeling hierarchical VAR correlated random effects

**URL:** <https://discourse.mc-stan.org/t/modeling-hierarchical-var-correlated-random-effects/14679>\
**Category:** Modeling\
**Created:** [April 28, 2020, 6:58pm UTC](https://discourse.mc-stan.org/t/modeling-hierarchical-var-correlated-random-effects/14679 "2020-04-28T18:58:00Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![Reg\_John](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/reg_john/32/6732_2.png) [@Reg\_John](https://discourse.mc-stan.org/u/Reg_John)\
**Post date:** [April 28, 2020, 6:58pm UTC](https://discourse.mc-stan.org/t/modeling-hierarchical-var-correlated-random-effects/14679/1 "2020-04-28T18:58:00Z")

</div>

I have a VAR model defined as:  
Y\_{1,it} = \alpha\_{1} + \beta\_{1i}\*Y\_{1,i(t-1)} + \beta\_{2i}\*Y\_{2,i(t-1)} +\epsilon\_{1,it}  
Y\_{2,it} = \alpha\_{2} + \beta\_{3i}\*Y\_{1,i(t-1)} + \beta\_{4i}\*Y\_{2,i(t-1)} +\epsilon\_{2,it}.

\beta\_{ji} is defined as \beta\_{ji} = \beta\_{j} + U\_{\beta ji}, for j = 1,2,3,4 where the random effects are defined as  
\epsilon\_{it} \sim MVN(0, \Sigma)  
U\_{\beta j} \sim MVN(0, \Phi).

I will also like to point out that \epsilon\_{1,it} and \epsilon\_{2,it} are correlated. Similarly, all the level-2 random effects U\_{\beta j} are correlated.

In order for the model to be stationary, I simulated U\_{\beta j} from a multivariate normal distribution with very small variances and no covariances. Therefore, my covariance matrix \Phi is a diagonal matrix, and my correlation matrix is almost an identity matrix.

Everything seems to work but I am having trouble recovering the parameter values for the correlation matrix for the level-2 random effects. I am not sure if I should use lkj\_corr\_cholesky(10) as prior for the random effects correlation matrix since some of my values are way off

* * *

My stan code is

```stan
data {
  int N_Y; 
  int N_ID;
  int N_T;
  vector[2] Y[N_ID, N_T];
}

parameters {
  vector[4] UB_I[N_ID]; //random effects
  vector<lower = 0>[4] B_tau; //random effects SD
  cholesky_factor_corr[4] B_ch; //cholesky correlated random effects
  vector[4] Beta; //overall auto and cross lag regression coefficients
  vector[2] alpha; //fixed intercepts
  
  vector<lower = 0>[N_Y] Y_sd; //Level 1 residual sigma
  cholesky_factor_corr[N_Y] Y_ch; //Level 1 residual cholesky corr
}

model {
  vector[N_Y] mu_Y[N_ID,N_T]; //Expected Y
  vector[4] B_I[N_ID]; //Level 2 auto and crossed lag slopes
  
  for(j in 1: 4){
    for(i in 1:N_ID){
      B_I[i,j] = Beta[j] + UB_I[i,j];
    }
  }
  
  for (i in 1:N_ID){
    
    //time = 1
    mu_Y[i,1,1] = alpha[1];
    mu_Y[i,1,2] = alpha[2];
    
    //time > 2
    for(t in 2:N_T){
      mu_Y[i,t,1] = alpha[1] + B_I[i,1]*Y[i,t-1,1] + B_I[i,2]*Y[i,t-1,2];
      mu_Y[i,t,2] = alpha[2] + B_I[i,3]*Y[i,t-1,1] + B_I[i,4]*Y[i,t-1,2];
    }
    
    
    UB_I[i,] ~ multi_normal_cholesky(rep_vector(0,4), diag_pre_multiply(B_tau, B_ch));
    
    for (t in 1:N_T){
      Y[i,t] ~ multi_normal_cholesky(mu_Y[i,t], diag_pre_multiply(Y_sd, Y_ch));
    }
    
  }
  
  //priors
  Y_sd ~ normal(0,3);
  Y_ch ~ lkj_corr_cholesky(1);
  B_ch ~ lkj_corr_cholesky(10);
  B_tau ~ std_normal();
  Beta ~ std_normal();
}

generated quantities {
  corr_matrix[4] B_corr;
  cov_matrix[4] B_cov;
  cov_matrix[N_Y] Y_cov;
  corr_matrix[N_Y] Y_corr;
  
  B_corr = multiply_lower_tri_self_transpose(B_ch);
  B_cov = quad_form_diag(B_corr, B_tau);
  
  Y_corr = multiply_lower_tri_self_transpose(Y_ch);
  Y_cov = quad_form_diag(Y_corr, Y_sd);
}

```

My code for simulating the data is

```stan
N <- 200 #number of individuals
Nt <- 20 #Number of time occasions

eps <- vector("list", length = Nt) #Level-1 residual structure
Y1 <- vector("list", length = Nt)

#Residuals for Structural Model
sigma <- matrix(c(1, .2,.2,1),2)
for (j in 1:Nt) {
  eps[[j]] <- mvtnorm::rmvnorm(N, mean = c(0,0), sigma = sigma)
}

#time 1
  
Y1 <- eps

cov_beta <- matrix(c(.01, 0, 0, 0, 
                     0, .01, 0, 0,
                     0, 0, .02, 0,
                     0, 0, 0, .05), nrow = 4)

ov_beta <- c(.4, .2, .4,.2) #overall beta_lag

#autoregression slopes
beta_lag <- mvtnorm::rmvnorm(N, mean = ov_beta, sigma = cov_beta)

  #time 2 & beyond
for (j in 2:Nt) {
  Y1[[j]][,1] <- Y1[[j-1]][,1]*beta_lag[,1] + Y1[[j-1]][,2]*beta_lag[,2] + eps[[j]][,1]
  Y1[[j]][,2] <- Y1[[j-1]][,1]*beta_lag[,3] + Y1[[j-1]][,2]*beta_lag[,4] + eps[[j]][,2]
}

Y <- array(0, dim = c(N,Nt,2))
for (j in 1:2) {
  for (t in 1:Nt) {
    for(i in 1:N) {
      Y[i,t,j] = Y1[[t]][i,j]
    }    
  }  
}

```

* * *

---

<div class="post-metadata">

**Author:** ![emiruz](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/emiruz/32/5518_2.png) [@emiruz](https://discourse.mc-stan.org/u/emiruz)\
**Post date:** [May 1, 2020, 7:33am UTC](https://discourse.mc-stan.org/t/modeling-hierarchical-var-correlated-random-effects/14679/2 "2020-05-01T07:33:57Z")

</div>

Hey would you mind publishing some of the outputs from the model? What are the Rhats like? Does the model converge? What do the separate chains look like? What happens when you run it a few times? I.e. same problems or different ones?

---

<div class="post-metadata">

**Author:** ![Reg\_John](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/reg_john/32/6732_2.png) [@Reg\_John](https://discourse.mc-stan.org/u/Reg_John)\
**Post date:** [May 2, 2020, 4:30am UTC](https://discourse.mc-stan.org/t/modeling-hierarchical-var-correlated-random-effects/14679/3 "2020-05-02T04:30:57Z")

</div>

Hi emiruz, sorry for the late reply.  
So my stan results are:

```stan
Inference for Stan model: VAR-202005012253-1-a063aa.
3 chains, each with iter=1000; warmup=500; thin=1; 
post-warmup draws per chain=500, total post-warmup draws=1500.

             mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
Beta[1] 0.40 0.00 0.02 0.37 0.39 0.40 0.41 0.44 1110 1.00
Beta[2] 0.17 0.00 0.02 0.14 0.16 0.17 0.18 0.20 1789 1.00
Beta[3] 0.39 0.00 0.02 0.35 0.38 0.39 0.40 0.42 1404 1.00
Beta[4] 0.17 0.00 0.02 0.13 0.15 0.17 0.18 0.21 1384 1.00
B_corr[1,1] 1.00 NaN 0.00 1.00 1.00 1.00 1.00 1.00 NaN NaN
B_corr[1,2] -0.04 0.01 0.14 -0.30 -0.14 -0.04 0.06 0.23 213 1.01
B_corr[1,3] -0.13 0.01 0.14 -0.39 -0.22 -0.13 -0.03 0.14 312 1.00
B_corr[1,4] -0.05 0.01 0.13 -0.30 -0.13 -0.05 0.03 0.20 400 1.01
B_corr[2,1] -0.04 0.01 0.14 -0.30 -0.14 -0.04 0.06 0.23 213 1.01
B_corr[2,2] 1.00 NaN 0.00 1.00 1.00 1.00 1.00 1.00 NaN NaN
B_corr[2,3] 0.03 0.01 0.14 -0.23 -0.06 0.03 0.12 0.30 416 1.00
B_corr[2,4] -0.09 0.01 0.13 -0.34 -0.18 -0.09 0.01 0.17 375 1.00
B_corr[3,1] -0.13 0.01 0.14 -0.39 -0.22 -0.13 -0.03 0.14 312 1.00
B_corr[3,2] 0.03 0.01 0.14 -0.23 -0.06 0.03 0.12 0.30 416 1.00
B_corr[3,3] 1.00 NaN 0.00 1.00 1.00 1.00 1.00 1.00 NaN NaN
B_corr[3,4] 0.04 0.01 0.12 -0.19 -0.04 0.04 0.12 0.28 272 1.01
B_corr[4,1] -0.05 0.01 0.13 -0.30 -0.13 -0.05 0.03 0.20 400 1.01
B_corr[4,2] -0.09 0.01 0.13 -0.34 -0.18 -0.09 0.01 0.17 375 1.00
B_corr[4,3] 0.04 0.01 0.12 -0.19 -0.04 0.04 0.12 0.28 272 1.01
B_corr[4,4] 1.00 NaN 0.00 1.00 1.00 1.00 1.00 1.00 NaN NaN
B_cov[1,1] 0.01 0.00 0.00 0.00 0.01 0.01 0.01 0.02 43 1.02
B_cov[1,2] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 140 1.02
B_cov[1,3] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 381 1.00
B_cov[1,4] 0.00 0.00 0.00 -0.01 0.00 0.00 0.00 0.00 381 1.01
B_cov[2,1] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 140 1.02
B_cov[2,2] 0.01 0.00 0.01 0.00 0.01 0.01 0.01 0.02 34 1.09
B_cov[2,3] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 523 1.00
B_cov[2,4] 0.00 0.00 0.00 -0.01 0.00 0.00 0.00 0.00 360 1.01
B_cov[3,1] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 381 1.00
B_cov[3,2] 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 523 1.00
B_cov[3,3] 0.02 0.00 0.01 0.01 0.01 0.02 0.02 0.03 59 1.07
B_cov[3,4] 0.00 0.00 0.00 -0.01 0.00 0.00 0.00 0.01 365 1.01
B_cov[4,1] 0.00 0.00 0.00 -0.01 0.00 0.00 0.00 0.00 381 1.01
B_cov[4,2] 0.00 0.00 0.00 -0.01 0.00 0.00 0.00 0.00 360 1.01
B_cov[4,3] 0.00 0.00 0.00 -0.01 0.00 0.00 0.00 0.01 365 1.01
B_cov[4,4] 0.04 0.00 0.01 0.03 0.03 0.04 0.05 0.06 492 1.01
Y_corr[1,1] 1.00 NaN 0.00 1.00 1.00 1.00 1.00 1.00 NaN NaN
Y_corr[1,2] 0.22 0.00 0.02 0.19 0.21 0.22 0.23 0.25 2622 1.00
Y_corr[2,1] 0.22 0.00 0.02 0.19 0.21 0.22 0.23 0.25 2622 1.00
Y_corr[2,2] 1.00 NaN 0.00 1.00 1.00 1.00 1.00 1.00 NaN NaN
Y_cov[1,1] 1.00 0.00 0.02 0.96 0.99 1.01 1.02 1.05 834 1.00
Y_cov[1,2] 0.22 0.00 0.02 0.19 0.21 0.22 0.23 0.25 2665 1.00
Y_cov[2,1] 0.22 0.00 0.02 0.19 0.21 0.22 0.23 0.25 2665 1.00
Y_cov[2,2] 0.99 0.00 0.02 0.94 0.97 0.99 1.00 1.04 2054 1.00

```

And my initial parameters for the simulation especially for the covariance and correlation matrices of the beta\_lag effects are:

```
> cor(beta_lag)
            [,1] [,2] [,3] [,4]
[1,] 1.00000000 -0.06905013 -0.10550052 -0.09441254
[2,] -0.06905013 1.00000000 -0.03730771 -0.06898384
[3,] -0.10550052 -0.03730771 1.00000000 0.14342457
[4,] -0.09441254 -0.06898384 0.14342457 1.00000000

> cov(beta_lag)
             [,1] [,2] [,3] [,4]
[1,] 0.010591585 -0.000730314 -0.001462815 -0.002078992
[2,] -0.000730314 0.010561592 -0.000516556 -0.001516892
[3,] -0.001462815 -0.000516556 0.018151328 0.004134475
[4,] -0.002078992 -0.001516892 0.004134475 0.045781003

```

Also, my traceplot for the correlation matrix of the beta\_lag is

 ![Rplot](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/0/06009a770a4bac8a04d85a2c65df976de2647a78.jpeg)  
and the traceplot for the beta\_lag covariance matrix is

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/d/d852b9bce039a161b5dd5d6151b7065900a2cfe6.png)  
As you can see, the estimated parameters for the correlation matrix of the beta\_lag is an issue. I’m not sure now if it’s because the values of the covariance between the lagged coefficients are too small thereby presenting an issue in the estimation process or the simulation process is flawed.

---

<div class="post-metadata">

**Author:** ![Reg\_John](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/reg_john/32/6732_2.png) [@Reg\_John](https://discourse.mc-stan.org/u/Reg_John)\
**Post date:** [May 2, 2020, 4:53am UTC](https://discourse.mc-stan.org/t/modeling-hierarchical-var-correlated-random-effects/14679/4 "2020-05-02T04:53:37Z")

</div>

Just adding the density plots of the correlation plot of the beta\_lag parameters

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/3/3af557ff5a221044e5489c62d043f5c2bf3fd548.png)  
and covariance plot

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/0/0ceb93ebbc23b3105ee850f479ba57e2b91c05da.png)

---

<div class="post-metadata">

**Author:** ![Charles\_Driver](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/charles_driver/32/6524_2.png) [@Charles\_Driver](https://discourse.mc-stan.org/u/Charles_Driver)\
**Post date:** [May 2, 2020, 1:18pm UTC](https://discourse.mc-stan.org/t/modeling-hierarchical-var-correlated-random-effects/14679/5 "2020-05-02T13:18:43Z")

</div>

what are the problems, exactly? I can’t see what’s wrong…  
In case it helps, I wrote up a model like this in: [https://www.researchgate.net/profile/Charles\_Driver/publication/310747801\_Hierarchical\_Bayesian\_Continuous\_Time\_Dynamic\_Modeling/links/5c4a0a4f458515a4c73da33a/Hierarchical-Bayesian-Continuous-Time-Dynamic-Modeling.pdf](https://www.researchgate.net/profile/Charles_Driver/publication/310747801_Hierarchical_Bayesian_Continuous_Time_Dynamic_Modeling/links/5c4a0a4f458515a4c73da33a/Hierarchical-Bayesian-Continuous-Time-Dynamic-Modeling.pdf)

---

<div class="post-metadata">

**Author:** ![Reg\_John](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/reg_john/32/6732_2.png) [@Reg\_John](https://discourse.mc-stan.org/u/Reg_John)\
**Post date:** [May 2, 2020, 6:07pm UTC](https://discourse.mc-stan.org/t/modeling-hierarchical-var-correlated-random-effects/14679/6 "2020-05-02T18:07:55Z")

</div>

I was thinking the correlation estimated values will match closely to the original values in all cases. I guess I’m being too hard on my model :).  
Thanks, very much for sharing your article. It will definitely help me with my model.