Using the R2D2M2 prior on a hierarchical ordered logistic model using brms

After working on it for a few days, I’m getting closer but I’m still stuck trying to generalize the code generated by brms. Looking at the code from rstanarm, I have a general idea of where the solution might be.

Basically, rstanarm states that, since we set the latent residual standard deviation to \sigma_\epsilon = 1, then the standard deviation of the latent variable y^* is \sigma_{y^*} = 1/\sqrt{1- R^2}.

From this relation, rstanarm computes the intercepts (cutpoints) and the effect coefficients (beta) like this :

parameters {
  simplex[J] pi;
  array[K > 1] unit_vector[K > 1 ? K : 2] u;
  real<lower=(K > 1 ? 0 : -1), upper=1> R2;
}
transformed parameters {
  vector[K] beta;
  vector[J - 1] cutpoints;
  {
    real Delta_y;
    if (K > 1) {
      Delta_y = inv_sqrt(1 - R2);
      beta = u[1] * sqrt(R2) * Delta_y * sqrt_Nm1;
    } else {
      Delta_y = inv_sqrt(1 - square(R2));
      beta[1] = R2 * Delta_y * sqrt_Nm1;
    }
    cutpoints = make_cutpoints(pi, Delta_y, link);
  }
}

For reference, the function make_cutpoints() is :

  vector make_cutpoints(vector probabilities, real scale, int link) {
    int C = rows(probabilities) - 1;
    vector[C] cutpoints;
    real running_sum = 0;
    if (link == 1)
      for (c in 1 : C) {
        running_sum += probabilities[c];
        cutpoints[c] = logit(running_sum);
      }
    // other links ...
    else
      reject("invalid link");
    return scale * cutpoints;
  }

Notice that the cutpoints are multiplied by \sigma_{y^*}. As far as I understand, it means
that in order to implement the induced dirichlet prior on the intercepts, you must generate cutpoints from the vector of probabilities \pi and then multiply them by 1/\sqrt{1-R^2}.

All the above, I think I understand. The part where I’m stuck is computing the coefficients. According to this vignette, rstanarm defines the effect coefficients as :

\mathbf{\beta} = \mathbf{R^{-1}} * \mathbf{u}*\sigma_y * \sqrt{R^2 * (N - 1)}

This corresponds to beta = u[1] * sqrt(R2) * Delta_y * sqrt_Nm1 in the code above.

First off, I don’t know where the \mathbf{R^{-1}} went. More importantly, I’m not sure how to plug this back into the code generated by brms. I think the unit vector u in the rstanarm implementation serves the same purpose as the simplex R2D2_phi in the brms code but that’s pretty much all I understand.

I’m sorry I don’t have a more specific question to ask but this is where I’ve been stuck for a while now and some input would be a huge help.