Unexpected Lack of Convergence of Model with Multicollinearity

Hi!

I’m working on Richard Mcelreath’s book Statistical Rethinking for self-study and I’m answering end of chapter exercises using PyMC3.

One question involves creating a model with highly correlated data. The idea is to predict height with the left leg and the right leg. Data was generated as follows:

height = np.random.normal(loc=10, scale=2, size=100)
leg_prop = np.random.uniform(low=0.4, high=0.5, size=100)
leg_left = leg_prop * height + np.random.normal(loc=0, scale=0.02, size=100)
leg_right = leg_prop * height + np.random.normal(loc=0, scale=0.02, size=100)

highly_corr_df = pd.DataFrame({
    'height': height,
    'leg_left': leg_left,
    'leg_right': leg_right
})

The first model is simple and converges:

with pm.Model() as highly_corr_model:
    bL = pm.Normal('bL', 2, 10)
    bR = pm.Normal('bR', 2, 10)
    alpha = pm.Normal('alpha', 10, 100)
    sigma = pm.HalfCauchy('sigma', 1)
    mu = bL * highly_corr_df['leg_left'] + bR * highly_corr_df['leg_right'] + alpha
    height = pm.Normal('height', mu, sigma, observed=highly_corr_df['height'])

    trace_highly_corr_model = pm.sample(4000, tune=1000, start={
        'bL': 0,
        'bR': 0,
        'alpha': 10,
        'sigma': 1
    })

However, the question asks what happens if bR was only positive. Here’s the R equivalent:

m5.8s2 <- map2stan(
    alist(
        height ~ dnorm(mu, sigma),
        mu <- a + bl*leg_left + br*leg_right,
        a ~ dnorm(10,100),
        bl ~ dnorm(2, 10),
        br ~ dnorm(2, 10) & T[0,],
        sigma ~ dcauchy(0, 1),
    ),
    data=d, chains=4,
    start=list(a=10, bl=0, br=0, sigma=1)
)

Here’s my code that uses Bound to ensure that bR can only be positive:

with pm.Model() as highly_corr_model_truncated:
    bL = pm.Normal('bL', 2, 10, testval=0)
    BoundedNormal = pm.Bound(pm.Normal, lower=0.0)
    bR = BoundedNormal('bR', mu=2, sd=10, testval=0.1) 
    alpha = pm.Normal('alpha', 10, 100, testval=10)
    sigma = pm.HalfCauchy('sigma', 1, testval=1)
    mu = bL * highly_corr_df['leg_left'] + bR * highly_corr_df['leg_right'] + alpha
    height = pm.Normal('height', mu, sigma, observed=highly_corr_df['height'])

    trace_highly_corr_model_truncated = pm.sample(4000, tune=1000)

This model has convergence issues: !

Any ideas why this is happening? I don’t know if this is a bug or if I wrote code incorrectly. Any help would be appreciated! Thanks.

The implementation looks correct to me, so the unidentifiablitly might be an issue of the model itself.
Did you try to fit the same model/data in R using rethinking?