Priors for covariance matrix for multivariate normal

Questions
1

I’m re-implementing some SAS code (provided by a biostatistician) in PyMC3. The SAS code is as follows:

      prior theta0 ~ normal(2.4347, sd=0.07);
      prior theta1 ~ normal(0.0993, sd=0.04);
      prior theta2 ~ normal(-0.1198, sd=0.02);
      array theta{3} theta0 theta1 theta2;
 
      array tau{3,3} ( 1.0   -0.344   0.365
                      -0.344  1.0     0.055
                       0.365  0.0550  1.0);
      array cov{3,3};
      prior cov ~ iwish(8, tau);

      array theta_i{3} Xi0 Xi1 Xi2;
      random theta_i ~ mvn(theta, cov) subject=subjectID;

I can figure out how to formulate most of it in PyMC3.

    theta = pm.Normal('theta', [2.4347, 0.0993, -0.1198], [0.07, 0.04, 0.02], shape=3)
    cov = np.array([[ 1.0   -0.344   0.365
                     -0.344  1.0     0.055
                      0.365  0.0550  1.0]])
    theta_i = pm.MvNormal('theta_i', t, cov=cov, shape=(n_subjects, 3))

It’s not exactly the same though because in the SAS code, cov is a random variable but in the PyMC3 adaptation, cov is hard-coded.

It’s best that the covariance matrix be a random variable. I know I can use LKJCholeskyCov to set up the priors; however, I’m not sure how I would transform the cov matrix (above) so that LKJCholeskyCov starts with priors based on the predetermined covariance matrix above (we are deriving these covariances from published data and then applying them to analysis of a new dataset). How might I do this? Is the secret in how I handle sd_dist?

    theta = pm.Normal('theta', [2.4347, 0.0993, -0.1198], [0.07, 0.04, 0.02], shape=3)
    sd_dist = pm.HalfCauchy('sd_dist', 1)
    chol, corr, stds = pm.LKJCholeskyCov('L', n=3, eta=5, sd_dist=pm.Exponential.dist(sd_dist), compute_corr=True)
    cov = pm.Deterministic('cov', chol.dot(chol.T))

    theta_i = pm.MvNormal('theta_i', t, chol=chol, shape=(n_subjects, 3))
2

Unfortunately pymc3 does not have a good implementation of InverseWishart, so I don’t think you can port that code exactly without implementing it on your own (PR welcome :wink: ).
I think a reasonable workaround would be to use the LKJCholeskyCov prior you already defined, but add a second observed NvNormal to include the data from the paper.
If that isn’t available you could also randomly generate some that has the right properties. (and verify that the choice of seed does not influence your results much).

3 Likes
3

Also, note that InverseWishart is a pretty bad choice of prior for covariance and since been discouraged by most (you should avoid conjugate priors in general when sampling with HMC) - +1 to use a LKJCholeskyCov and play with eta to get the prior you are happy with

3 Likes
4

That seems like a good solution, thanks!