# Multiple regression caveats

**URL:** <https://discourse.pymc.io/t/multiple-regression-caveats/1557>\
**Category:** Questions\
**Created:** [July 23, 2018, 3:12pm UTC](https://discourse.pymc.io/t/multiple-regression-caveats/1557 "2018-07-23T15:12:10Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![phildias](https://avatars.discourse-cdn.com/v4/letter/p/a6a055/32.png) [@phildias](https://discourse.pymc.io/u/phildias)\
**Post date:** [July 23, 2018, 3:12pm UTC](https://discourse.pymc.io/t/multiple-regression-caveats/1557/1 "2018-07-23T15:12:10Z")

</div>

Hi there!

I’m trying to get a multiple regression off the ground, but I’m running into some problems.

The general structure of the model is like so:

y\_{i,1} = \mathbf{x}\_{i,1} \cdot \mathbf{\beta}\_{1}' + \epsilon\_{i,1}  
y\_{i,2} = \mathbf{x}\_{i,2} \cdot \mathbf{\beta}\_{2}' + \epsilon\_{i,2}

dim(\mathbf{x}\_{i,1})=dim(\mathbf{\beta}\_{i,1})=(k\_1,1)   
dim(\mathbf{x}\_{i,2})=dim(\mathbf{\beta}\_{i,2})=(k\_2,1)  
dim(y\_{i,1})=dim(y\_{i,2})=(1,1)

\bf{Likelihood}  
 \mathbf{\epsilon}\_i = (\epsilon\_{i,1},\epsilon\_{i,2}) \stackrel{iid}{\sim} MVN(mean=\mathbf{0}\_2,cov=\mathbf{\Sigma\_\epsilon})

\bf{Priors}  
 (\beta\_{1,1},...,\beta\_{{k\_1},1}) \stackrel{iid}{\sim} Normal(mean=0,stdev=10)  
 (\beta\_{1,2},...,\beta\_{{k\_2},2}) \stackrel{iid}{\sim} Normal(mean=0,stdev=10)  
 \mathbf{sd}=(sd\_1,sd\_2) \stackrel{iid}{\sim} HalfCauchy(beta=2.5)  
 Chol(\mathbf{\Sigma}\_\epsilon) \sim LKJCholeskyCov(eta=2, n=2, sd=\mathbf{sd})

In the example I’m giving above, I have 2 regression lines, \mathbf{y}\_1 which is explained by variables \mathbf{x}\_1 and \mathbf{y}\_2 which is explained by variables \mathbf{x}\_2. Notice that this form is general enough that \mathbf{x}\_1 and \mathbf{x}\_2 can be constituted by completely different columns.

What is making my life a bit difficult is translating this into PyMC3. Let me show you my example using some [fake data](https://www.dropbox.com/s/3vl01iz7wygvd27/fake_multiple_multivariate_regression_data_covmtx.csv?dl=0)

```
import pymc3 as pm
import numpy as np
import pandas as pd
import theano
import theano.tensor as tt

# Reading in dataset
dataset = pd.read_csv("fake_multiple_multivariate_regression_data_covmtx.csv")

num_rows = len(dataset)

# Naming the dependent categorical variables
dep_vars = ["y1_cont","y2_cont"]

num_dep_vars = len(dep_vars)

# Exogenous variables impacting each depvar. 
indep_vars = [["uno","x1","x2","x3","x4"],
              ["uno","x3","x4","x5","x6","x7"]]

indep_vars_y1 = len(indep_vars[0])
indep_vars_y2 = len(indep_vars[1])

x_data_y1 = dataset[indep_vars[0]].values
x_data_y2 = dataset[indep_vars[1]].values

y_data = dataset[dep_vars].values

with pm.Model() as model: # model specifications in PyMC3 are wrapped in a with-statement
    # Priors for betas
    betas_y1 = pm.Normal('Betas_Y1', 0, sd=20, shape=indep_vars_y1)
    betas_y2 = pm.Normal('Betas_Y2', 0, sd=20, shape=indep_vars_y2)

    # Priors for error covariance matrix
    sd_dist = pm.HalfCauchy.dist(beta=2.5, shape=num_dep_vars)
    packed_chol = pm.LKJCholeskyCov('chol_cov', eta=2, n=num_dep_vars, sd_dist=sd_dist)
    chol = pm.expand_packed_triangular(num_dep_vars, packed_chol, lower=True)
    
    y1_hat = tt.dot(x_data_y1, betas_y1)
    y2_hat = tt.dot(x_data_y2, betas_y2)

    # Here is my main problem!!!
    mu = pm.Deterministic('mu', tt.stack([y1_hat,y2_hat])) 
    
    # Define a new MvNormal with the given covariance
    vals = pm.MvNormal('vals', mu=mu, chol=chol, observed=y_data)

```

In the example above, I’m having a hard time re-packing `y1_hat` and `y2_hat` in a single array-type element so I can use the `MvNormal`.

Does anyone have any insight on how to do so?

Thank you!  
.

---

<div class="post-metadata">

**Author:** ![bwengals](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/bwengals/32/1237_2.png) [@bwengals](https://discourse.pymc.io/u/bwengals)\
**Post date:** [July 23, 2018, 3:50pm UTC](https://discourse.pymc.io/t/multiple-regression-caveats/1557/2 "2018-07-23T15:50:46Z")

</div>

All you need is a small fix:

```
mu = pm.Deterministic('mu', tt.stack([y1_hat, y2_hat], axis=1)

```

`tt.stack` by default uses axis=0, which takes your two arrays horizontally, so the output shape is (2, 5000) instead of (5000,2).

---

<div class="post-metadata">

**Author:** ![phildias](https://avatars.discourse-cdn.com/v4/letter/p/a6a055/32.png) [@phildias](https://discourse.pymc.io/u/phildias)\
**Post date:** [July 23, 2018, 4:13pm UTC](https://discourse.pymc.io/t/multiple-regression-caveats/1557/3 "2018-07-23T16:13:56Z")

</div>

Thanks, that worked like a charm =)

For some reason, though, another issue popped up: when I try to sample from this model, pymc3 just hangs on me.

So when I call this:

```
with model:
    trace = pm.sample()

```

Python then hangs on this screen after calculating the first quarter of the job:

```
Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [chol_cov, Betas_Y2, Betas_Y1]
Sampling 4 chains: 25%|██▌ | 1000/4000 [03:00<01:02, 47.98draws/s] 

```

Any thoughts?

Thanks again for the help!!!

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [July 23, 2018, 4:31pm UTC](https://discourse.pymc.io/t/multiple-regression-caveats/1557/4 "2018-07-23T16:31:05Z")

</div>

> [@phildias](#):
>
> sd\_dist = pm.HalfCauchy.dist(beta=2.5, shape=num\_dep\_vars)

Try using a more informative prior with lighter tail might help:

```python
sd_dist = pm.HalfNormal.dist(sd=1, shape=num_dep_vars)

```

---

<div class="post-metadata">

**Author:** ![phildias](https://avatars.discourse-cdn.com/v4/letter/p/a6a055/32.png) [@phildias](https://discourse.pymc.io/u/phildias)\
**Post date:** [July 23, 2018, 4:54pm UTC](https://discourse.pymc.io/t/multiple-regression-caveats/1557/5 "2018-07-23T16:54:25Z")

</div>

Thanks for the tip, @junpenglao!

Unfortunately, though, it still hangs… =/

---

<div class="post-metadata">

**Author:** ![bwengals](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/bwengals/32/1237_2.png) [@bwengals](https://discourse.pymc.io/u/bwengals)\
**Post date:** [July 23, 2018, 5:31pm UTC](https://discourse.pymc.io/t/multiple-regression-caveats/1557/6 "2018-07-23T17:31:27Z")

</div>

Might be a problem with the multiprocess sampling. This is likely fixed in the newest version of PyMC3. Before you update, try setting the `cores=1`.

```
with model:
    trace = pm.sample(1000, chains=4, cores=1)
```

---

<div class="post-metadata">

**Author:** ![phildias](https://avatars.discourse-cdn.com/v4/letter/p/a6a055/32.png) [@phildias](https://discourse.pymc.io/u/phildias)\
**Post date:** [July 23, 2018, 8:21pm UTC](https://discourse.pymc.io/t/multiple-regression-caveats/1557/7 "2018-07-23T20:21:26Z")

</div>

Yeah, @bwengals, it really does seem to be a problem with the multiprocess sampling. When I run the snippet you suggested, things work out fine: all four chains get sampled sequentially and no problems are encountered.

I hope the multiprocess issue gets solved in the next version!

Thanks =)