Modelling hyper priors for sigma with partially exponentiated components

Modeling
1

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!

2

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

3

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