# Linear multivariate model

**URL:** <https://discourse.mc-stan.org/t/linear-multivariate-model/5407>\
**Category:** Modeling\
**Tags:** fitting-issues\
**Created:** [August 28, 2018, 11:19am UTC](https://discourse.mc-stan.org/t/linear-multivariate-model/5407 "2018-08-28T11:19:10Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![jroon](https://avatars.discourse-cdn.com/v4/letter/j/fbc32d/32.png) [@jroon](https://discourse.mc-stan.org/u/jroon)\
**Post date:** [August 28, 2018, 11:19am UTC](https://discourse.mc-stan.org/t/linear-multivariate-model/5407/1 "2018-08-28T11:19:10Z")

</div>

Hi all, I am trying to specify a linear multivariate model and I’m having some difficulties: 1) I’m not recovering parameters correctly from a simulated data set and 2) I need to speed it up for my real data.

First here is a simulated data script: [gen\_test\_data\_y2\_x3.R](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/9/9f86dd3676277442cd493a258cf13075a013f0cb.R) (530 Bytes)  
And here is the script to call stan from r: [run\_test\_model.R](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/7/72132bfcc621e7c1474e3b482779670eabaccc66.R) (311 Bytes)

My model is based on section 9.15 of the Stan manual (save to linear\_multivariate2.stan):

```data
    int<lower=1> NvarsY;
    int<lower=1> NvarsX;
    int<lower=0> N;
    vector[NvarsX] x[N];
    vector[NvarsY] y[N];
}

parameters {
    matrix[NvarsY, NvarsX] beta;
    cholesky_factor_corr[NvarsY] L_Omega;
    vector<lower=0>[NvarsY] L_sigma;
}

model {
    vector[NvarsY] mu;
    matrix[NvarsY, NvarsY] L_Sigma;
    
    //Priors
    L_Sigma = diag_pre_multiply(L_sigma, L_Omega);
    to_vector(beta) ~ normal(0, 5);
    L_Omega ~ lkj_corr_cholesky(4);
    L_sigma ~ cauchy(0, 2.5);

    //Model
    for (n in 1:N){
        mu = beta * x[n];
    }

    y ~ multi_normal_cholesky(mu, L_Sigma);
}

```

My recovered parameters from running this were:

> summary(params$beta[,1,])  
> V1 V2 V3  
> Min. :-3.953 Min. :-18.506 Min. :-13.2165  
> 1st Qu.: 3.760 1st Qu.: -4.397 1st Qu.: -0.1258  
> Median : 5.483 Median : -1.051 Median : 2.8229  
> Mean : 5.465 Mean : -1.053 Mean : 2.8344  
> 3rd Qu.: 7.095 3rd Qu.: 2.287 3rd Qu.: 5.7141  
> Max. :14.481 Max. : 18.401 Max. : 20.0149  
> summary(params$beta[,2,])  
> V1 V2 V3  
> Min. : 0.7665 Min. :-19.278 Min. :-11.595  
> 1st Qu.: 9.0719 1st Qu.: -5.228 1st Qu.: 2.816  
> Median :10.7551 Median : -2.049 Median : 5.885  
> Mean :10.7435 Mean : -2.046 Mean : 5.835  
> 3rd Qu.:12.4121 3rd Qu.: 1.115 3rd Qu.: 8.754  
> Max. :20.2864 Max. : 15.609 Max. : 24.415

… which don’t match the correlations I specified in the generated data.

My second problem is that for my real dataset (~2000 rows, 23 Y variables, 40 X variables) it takes about 24 hours to run. Perhaps I’m simply expecting too much giving the dimensionality of the data and not huge number of rows. I can reduce the number of X variables conceivably based on prior knowledge, but there would still be perhaps 20 X variables, and the number of Y variables are fixed.

---

<div class="post-metadata">

**Author:** ![martinmodrak](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/martinmodrak/32/133_2.png) [@martinmodrak](https://discourse.mc-stan.org/u/martinmodrak)\
**Post date:** [September 6, 2018, 12:38pm UTC](https://discourse.mc-stan.org/t/linear-multivariate-model/5407/2 "2018-09-06T12:38:42Z")

</div>

This might be because you are not simulating directly from your model: in general you should simulate by first drawing parameters from priors (e.g. `beta`, `L_Omega` and `L_sigma`) and follow exactly your model (e.g. calculate `mu` from `beta` and `x`). It is well possible that `sigma` or `mu` you hardcoded in your simulation script actually have very low prior mass and thus cannot be recovered. Also simulating from the prior helps you understand the implications of your prior and check that it makes sense (as in the [Visualisation paper](https://arxiv.org/abs/1709.01449))

---

<div class="post-metadata">

**Author:** ![jroon](https://avatars.discourse-cdn.com/v4/letter/j/fbc32d/32.png) [@jroon](https://discourse.mc-stan.org/u/jroon)\
**Post date:** [September 6, 2018, 1:11pm UTC](https://discourse.mc-stan.org/t/linear-multivariate-model/5407/3 "2018-09-06T13:11:48Z")

</div>

Hi @martinmodrak - thanks for your answer. I actually figured this out yesterday by using the model from this thread and removing the subject specific terms: [Divergent transitions in a multivariate multilevel regression model](https://discourse.mc-stan.org/t/divergent-transitions-in-a-multivariate-multilevel-regression-model/5312/18)

Also, in the simulation r script I set `mvrnorm` command I set `empircal=TRUE`.

The new model I’m using is below:

```
data{
    int<lower=0> NvarsX; // num independent variables
    int<lower=0> NvarsY; // num dependent variables
    int<lower=0> N; //Total number of rows
    matrix[N, NvarsX] x;
    vector[NvarsY] y [N];
}

parameters{
    cholesky_factor_corr[NvarsY] L_Omega; //For the correlations between outcome variables
    vector<lower=0>[NvarsY] L_sigma; //Residual SD of outcome variables

    vector[NvarsY] Beta0; //intercepts
    vector[NvarsX] Beta [NvarsY]; //slopes
}

transformed parameters{
    vector[NvarsY] mu[N];
    matrix[NvarsY, NvarsY] L_Sigma;

    L_Sigma = diag_pre_multiply(L_sigma, L_Omega);

    for (i in 1:N){
        for (dv in 1:NvarsY){
            mu[i, dv] = Beta0[dv] + dot_product(to_vector(Beta[dv, 1:NvarsX]), x[i, 1:NvarsX]);
        }
    }
}

model{
    //Priors
    L_Omega ~ lkj_corr_cholesky(1);
    L_sigma ~ cauchy(0,5);

    for (i in 1:N){
        y[i, 1:NvarsY] ~ multi_normal_cholesky(mu[i, 1:NvarsY], L_Sigma);
    }

}

generated quantities{
    matrix[NvarsY, NvarsY] Omega;
    matrix[NvarsY, NvarsY] Sigma;
    Omega = multiply_lower_tri_self_transpose(L_Omega);
    Sigma = quad_form_diag(L_Omega, L_sigma);
} 

```

This model successfully recovered the correct parameter values.

I would still like to make it faster. I’ll explore reparameterising it when I get some free time next. (In the mean time…back to writing grant applications… :/ )