Hi all,

I am trying to run a very complicated Beta Fay-Herriot style model and would love some assistance with the necessary Jacobian adjustment. Below I present a highly simplified version of my model.

The general setup is as follows. y \in (0,1), with \text{Var}(y_i) = \psi_i and E(y_i) = \mu_i. Using the mean-precision parametization of the Beta distribution, the model setup is as follows.

To ensure, a_i>0, b_i>0, the formula for \phi_i implies the following bounds on \mu_i.

Now given that \psi_i is a random variable in our model, the bounds on \mu_i are also random. Using a combination of the non-centered parameterization and the necessary transform we can sample a standard normal parameter Z \sim N(0,1) and then transform this parameter to the correct bounds on \eta_i which are implied by taking the \text{logit}^{-1} transform of the bounds on \mu_i. On the bottom of this page, the transformation required for variables \in (L, U) specifies that a Jacobian adjustment must be added if L and U are random, which in this case they are.

The current stan model is below.

```
data {
int<lower=0> n;
vector<lower=0,upper=1>[n] y;
int<lower=0> q;
matrix[n, q] x;
vector[2] boundpsi;
}
parameters {
vector[q] B;
real<lower=0> sigma_v;
vector<lower=boundpsi[1],upper=boundpsi[2]>[n] psi;
vector[n] Z;
}
transformed parameters{
// DECLARACTIONS
vector[n] mu;
vector[n] eta;
vector[n] alpha;
vector[n] etaLower;
vector[n] etaUpper;
vector[n] phi;
// FULL VECTORS
etaLower = logit(0.5 * (1 - sqrt(1 - 4 * psi)));
etaUpper = logit(0.5 * (1 + sqrt(1 - 4 * psi)));
// Non-centered parameterization: https://mc-stan.org/docs/2_28/stan-users-guide/reparameterization.html
eta = Z * sigma_v + x*B;
// adjust the bounds on eta and transform to mu: https://mc-stan.org/docs/2_28/stan-users-guide/vectors-with-varying-bounds.html
mu = inv_logit(etaLower + (etaUpper - etaLower) .* inv_logit(eta));
phi = ((mu .* (1 - mu)) ./ (psi)) - 1;
alpha = mu .* phi;
}
model{
// draw from a standard normal
target += std_normal_lpdf(Z);
// priors for the model parameters
for(j in 1:q){
target += std_normal_lpdf(B[j]);
}
target += cauchy_lpdf(sigma_v | 0,2);
// increment log-probability for Jacobian
for(f in 1:n){
// https://mc-stan.org/docs/2_27/reference-manual/logit-transform-jacobian-section.html
target += log( (etaUpper[f] - etaLower[f]) * inv_logit(eta[f]) * (1 - inv_logit(eta[f]) ) );
}
// Beta model
target += beta_lpdf(y | alpha, phi - alpha);
}
```

Note a few caveats that cannot be edited to ensure this simple example still operates in my more complex version. Firstly, I wish to use the non-centered parameterization for the linear predictor, as I get significant pathological behaviour if I don’t. In this example, `psi`

is a parameter, whilst in my complex model, `psi`

is a transformed parameter. This means that the bounds on `eta`

cannot be defined in the parameters block using the `<lower=., upper=.>`

syntax.

The simple model converges very well, and I find close to no differences when I exclude/include the Jacobian adjustment. However, my complex model is having significant convergence problems so I want to be sure that the fault is not because of this bounded parameter.

Thankyou in advance for any assistance.