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 :
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.