Cholesky decomposition and correlation among random effects

Questions
16

I found a brms implementation already made in this blog post. The results indicate that the correlation is, in fact, \rho=0.06 [-0.51, 0.69]. The value I found with my Cholesky implementation is \rho=0.07 [-0.45, 0.59] which seems pretty close (I will run further test myself later).
I found also a rstanarm implementation: here \rho=0.07. In the same post he plotted the individual correlation for each participant

image
image749×569 83.6 KB

which looks like the ones I plotted too (with a correlation of \rho=-0.5). So I would conclude, as you said, that we cannot estimate the correlation between the random effect by looking at the correlation of their posterior samples

I obtained values that seem wrong. The following represents the distribution of \rho computed from each MCMC sample as

rho = []

for i in gamma_Z:
    rho.append(np.corrcoef(i.T)[0, 1])

where gamma_Z has shape (nsubject, 2). np.corrcoef assumes that each row of represents a variable, and each column a single observation of all those variables, hence the transpose T command.

I do not know how to interpret this result. Does it mean that modelling random effects without correlation matrix is wrong? As suggested in the post you linked Covariance of multiple varying coefficients in hierarchical models we are posing a strong prior on the random effects to be independent (but the data attempt to say otherwise).

image

Given the above, \rho_a does not agree well with \rho_b. But \rho_b agrees with the brms results.

I then check whether computing the \rho from the MCMC samples, agrees with \rho computed from the cholesky matrix

    covariance_matrix = pm.Deterministic('covariance_matrix', tt.dot(chol, chol.T))
    
    # Extract the standard deviations and rho
    standard_deviations = pm.Deterministic('standard_deviations', tt.sqrt(tt.diag(covariance_matrix)))
    correlation_matrix = pm.Deterministic('correlation_matrix', tt.diag(standard_deviations**-1).dot(covariance_matrix.dot(tt.diag(standard_deviations**-1))))
    rho = pm.Deterministic('rho', correlation_matrix[np.triu_indices(2, k=1)])

image

I would expect them to be almost identical, but they are not. Maybe, in general, it is not a good idea to compute the correlation from MCMC samples then.

1 Like