# Cholesky decomposition and correlation among random effects

**URL:** <https://discourse.pymc.io/t/cholesky-decomposition-and-correlation-among-random-effects/1210>\
**Category:** Questions\
**Created:** [May 11, 2018, 11:05pm UTC](https://discourse.pymc.io/t/cholesky-decomposition-and-correlation-among-random-effects/1210 "2018-05-11T23:05:44Z")\
**Posts on this page:** 1\
**Page:** 2

<div class="post-metadata">

**Author:** ![chartl](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/chartl/32/1515_2.png) [@chartl](https://discourse.pymc.io/u/chartl)\
**Post date:** [March 11, 2019, 9:05pm UTC](https://discourse.pymc.io/t/cholesky-decomposition-and-correlation-among-random-effects/1210/21 "2019-03-11T21:05:58Z")

</div>

I’m digging this up because it was referenced here: `https://discourse.pymc.io/t/random-intercepts-and-slopes-correlated-convergence-issues/2865/15`. I suspect that the issue has to do with the conflict between:

> [@Jack\_Caster](#):
>
> The covariance matrix is \mathbf{\Sigma} = \left(\begin{array}{cc} 23.8 & 0.82 \\ 0.82 & 5.75\end{array}\right)

and

> [@Jack\_Caster](#):
>
> This value of ρ=0.07, however, is not correct. In fact, the Pearson coefficient of each participant is about -0.5:

This to me sounded like a Simpson’s paradox. Having dug up the [notebook](https://github.com/JackCaster/GLM_with_PyMC3/blob/master/notebooks/LME%20-%20Sleep%20study.ipynb), it seems that this is the case. Plotting the posteriors of model `Random effect on the intercepts and slopes` gives strongly anti-correlated within-subject estimates:

![image](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/c/c72b5881442a4fd8a8bb2adc5b162c2d26cc2c2d.png)

but looking across samples

```auto
tis = trace_intercept_slope
for j in range(3):
    sns.scatterplot(tis['gamma_Z_slope'][:, j], tis['gamma_Z_intercept'][:, j], alpha=0.2)
    plt.figure()
    
colors = ['red', 'black', 'orange', 'yellow', 'green', 'blue', 'purple', 'magenta']
all_vals = [[],[]]
for j in range(tis['gamma_Z_slope'].shape[1]):
    sns.scatterplot(tis['gamma_Z_slope'][:150, j],
                    tis['gamma_Z_intercept'][:150, j],
                    color=colors[j % len(colors)],
                    alpha=0.2)
    all_vals[0].extend(tis['gamma_Z_slope'][:150, j])
    all_vals[1].extend(tis['gamma_Z_intercept'][:150, j])

np.cov(np.array(all_vals))

```

![image](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/8/843fd5565a3bbe20b50612bcc3bf5b0dc3e0c99f.png)

Or using a bivariate normal approximation:

![image](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/4/430a65952d5cd91ff13edb1c9b0386b3591b586e.png)

The overall covariance (using the posterior samples `all_vals`) between slope and intercept is

\mathbf{\Sigma} = \left(\begin{array}{cc} 36.7 & 3.5 \\ 3.5 & 645.2\end{array}\right)

which lines up with

> [@Jack\_Caster](#):
>
> I ran the notebook again. The sampling is way faster, but the correlation is still off (ρ=0.11\rho=0.11). Mmmm, maybe the variance in the intercept is quite high, which may suggest that something is off
> 
> ```auto
> [[596.98838465, 16.46747156],
> [16.46747156, 36.57017179] ]
> 
> ```

[Previous page](https://discourse.pymc.io/t/cholesky-decomposition-and-correlation-among-random-effects/1210.md?page=1)