 # Modelling hyper priors for sigma with partially exponentiated components

Hi everyone,

I’m struggling to implement priors on sigma when exponentiating sone components of it for the following model
y_i \sim \mathcal{N}(p_i,\sigma_i^2)
with a mean equation:
logit(p_i) = .....
and the variance equation, where n_i represents the sample size and x_1, x_2, x_3 are observed data:
\sigma_i^2 = p_i*(1-p_i)/n_i +exp(\tau_{group[i]} + \beta_{1group[i]} * x_{1i} + \beta_2*x_{2i} + \beta_{3}*x_{3group[i]})

Previously, there was no exp() function in \sigma, hence the the distribution of the mean and standard deviation was simply set to \mathcal{N}(0,0.05). But exponentiating was necessary to secure \sigma being positive.
Now I want to set priors on \tau,\beta_1, \beta_2,\ beta_3 such that they allow a 5% change, but I don’t know how to specify the priors. The model looks as follows:

data{
...
}

parameters{
...
// Hyper-parameters

real mu_beta1;
real<lower=0> sig_beta1;
real mu_beta2;
real<lower=0> sig_beta2;
real mu_beta3;
real<lower=0> sig_beta3;
real<lower=0> sig_tau;

// Scaled Parameters

vector[group] beta1_sc;
real beta2_sc;
real beta3_sc;
vector<lower=0>[SY] tau_sq_sc;
}

transformed parameters {
...
// Parameters
vector[group]beta1;
realbeta2;
real beta3;
vector<lower=0>[SY] tau_sq;

beta1 = mu_beta1 + sig_beta1 * beta1_sc;
beta2 = mu_beta2 + sig_beta2 * beta2_sc;
beta3 = mu_beta3 + sig_beta3 * beta3_sc;
tau_sq = sig_tau * tau_sq_sc;
}

model {
vector[N] logit_p;
vector[N] p;
...

// Hyper priors
mu_beta1 ~ normal(0,0.05);
sig_beta1 ~ normal(0,0.0.5) T[0,];

mu_beta2 ~ normal(0,0.0.5);
sig_beta2 ~ normal(0,0.0.5) T[0,];

mu_beta3 ~ normal(0,0.0.5);
sig_beta3 ~ normal(0,0.0.5) T[0,];

sig_tau ~ normal(0,0.2) T[0,];

// Priors
beta1_sc ~ normal(0,1);
beta2_sc ~ normal(0,1);
beta3_sc ~ normal(0,1);

for(i in 1:SY){
tau_sq_sc[i] ~ normal(0,1) T[0,];
}

for(i in 1:N){
logit_p[i] = ... ;
p[i] = inv_logit(logit_p[i]);
y[i] ~ normal(p[i], sqrt(p[i] * (1 - p[i]) /n[i] +
exp(tau_sq[group[i]] +
beta1[group[i]] * x1[i] +
beta2 * x2[i] +
beta3 * x3[group[i]])));
}
}


Would be great, if you have suggestions on how to model this and thanks a lot in advance!

Hey there! Sorry it took us some time to respond.

So, I guess you are modelling proportions? Maybe you can take a look at the Beta distribution (proportion parameterization), where you could model \kappa with a log link (using exp).

Also, I don’t really understand what the indexing of x3 is doing? And when \beta_2 is not varying, you don’t need to implement it using a NCP. In fact, you would not learn \sigma_{\beta_2} anyway (because it is not varying). Do you have a working version of your model already? In any case, the usual recommendation is to start with something simple (that works) and then increase complexity to see when/where stuff breaks.

Let me know if the Beta proportion suggestion makes sense to you, or if you have more questions.

Cheers,
Max

Thank you very much. Yes, I do model proportions. For now I sticked with the normal distribution but changed the variance equation to:

log(\sigma_i^2 ) = log(p_i * (1-p_i)/n_i) + \tau_{group[i]} + \beta_{1group[i]}* x_{1[i]} + \beta_{2} * x_{2[i]} + \beta_{3} * x_{3group[i]}

This ensures the total variance to be positive since \sigma_i^2=exp(log(sigma_i^2)). The model is working now. You are right, a NCP for \beta_2 is unnecessary. x_3 is measured at the group level, this is why I used the indexing.
I’ll have a look into the Beta proportion.

Thanks again,
Sina

1 Like