Why divergences in reparameterization of ACE model? solution?

Modeling
#1

I have a model that works great and an equivalent model that encounters divergences. Here’s the data simulation script (2.3 KB) and shell script to run the sampler (s04.txt, 740 Bytes). This is the working model gen2c.stan (1.7 KB). After I reparameterize the variance then I get lots of divergences, gen2c.stan (2.0 KB).

Maybe the easiest way to highlight the difference between the models is to look at the diff,

diff --git a/ace/gen2c.stan b/ace/gen2c.stan
index dbe7d82..57ad0d1 100644
--- a/ace/gen2c.stan
+++ b/ace/gen2c.stan
@@ -6,6 +6,8 @@ data {
   real score[N, 2*numPhenotypes];
 }
 transformed data {
+  vector[3] varPropPrior;
+  for (x in 1:3) varPropPrior[x] = 1./3.;
 }
 parameters {
   cholesky_factor_corr[numPhenotypes] aCorCF;
@@ -13,13 +15,20 @@ parameters {
   cholesky_factor_corr[numPhenotypes] eCorCF;
   row_vector[numPhenotypes] iAC[N];
   row_vector[numPhenotypes] iC[N];
-  row_vector<lower=0>[numPhenotypes] aScale;
-  row_vector<lower=0>[numPhenotypes] cScale;
-  row_vector<lower=0>[numPhenotypes] eScale;
+  simplex[3] varProp[numPhenotypes];
+  row_vector<lower=0>[numPhenotypes] totalVar;
+}
+transformed parameters {
+  row_vector[numPhenotypes] aScale;
+  row_vector[numPhenotypes] cScale;
+  row_vector[numPhenotypes] eScale;
+  aScale = sqrt(to_row_vector(varProp[,1]) .* totalVar);
+  cScale = sqrt(to_row_vector(varProp[,2]) .* totalVar);
+  eScale = sqrt(to_row_vector(varProp[,3]) .* totalVar);
 }
 model {
-  row_vector[numPhenotypes] totalVar;
-  
+  for (px in 1:numPhenotypes) varProp[px,] ~ dirichlet(varPropPrior);
+
   for (n in 1:N) {
     row_vector[numPhenotypes] errScale;
     matrix[numPhenotypes,numPhenotypes] errCov;
@@ -39,7 +48,6 @@ model {
   cCorCF ~ lkj_corr_cholesky(3.0);
   eCorCF ~ lkj_corr_cholesky(3.0);
 
-  totalVar = square(aScale) + square(cScale) + square(eScale);
   totalVar ~ lognormal(0,1);
 }
 generated quantities {

Instead of totalVar = aScale^2 + cScale^2 + eScale^2, I want to redefine aScale, cScale, and eScale as proportions of totalVar. This should be equivalent, but sampling doesn’t work very well. The reason I want to reparameterize like this is because sometimes I want to fix totalVar to 1 and only estimate the proportions. Maybe this isn’t feasible in Stan? Suggestions?

#2

bayesplot::mcmc_parcoord and / or pairs with pars = c("varProp", "totalVar") to see where the divergences are coming from.

2 Likes
#3

It probably doesn’t matter but aScale/cScale/eScare are all sampled on the standard deviation scale, whereas totalVar is sampled on the variance scale.

#4

mcmc_parcoord is helpful. Also, in retrospect, maybe I don’t need to discompose the variance components as in the model that was diverging. Thank you.