Finding priors for multilevel time-series model (response surface on L2)

Modeling
1

Hey everyone,

I’m currently working on finding weakly informative priors for a multilevel time-series model that includes a response surface analysis on L2. I expect the scaled and centered values to mostly fall between –2 and 2, but they’re often out of bounds and show an asymmetric tendency toward positive values instead of being roughly centered around zero.

Here are the current quantiles:

  • q05: –43.6
  • q25: –3.25
  • q75: 5.72
  • q95: 49.4

I suspect the main issue lies in the polynomial terms. One way I managed to bring the values into a more reasonable range was by scaling the polynomial coefficients of mu and lambda by 0.5, as well as scaling the entire exponential term of sigma. However, this feels more like a hack than a sound modeling practice.

I’d really appreciate any advice on how to specify priors that set more reasonable bounds and ideally reduce the asymmetry.

data {
  int<lower=1> N;              
  int<lower=1> Nobs;          
  array[Nobs] int<lower=1, upper=N> subj;
  vector[Nobs] lag_y;
  vector[N] S;
  vector[N] O;
}

parameters {
  vector[6] beta_mu;
  vector[6] beta_lambda;
  vector[6] beta_e;
  array[N] vector[2] z_u;
  vector<lower=0>[2] tau;
}

transformed parameters {
  array[N] vector[2] u;
  for (i in 1:N) {
    u[i,1] = tau[1] * z_u[i,1];
    u[i,2] = tau[2] * z_u[i,2];
  }
}

model {
  beta_mu     ~ normal(0, 1); 
  beta_lambda ~ normal(0, 1);      
  beta_e      ~ normal(0, 0.5);   

  tau[1]      ~ normal(0, 0.5);  
  tau[2]      ~ normal(0, 0.05);

  for (i in 1:N) 
  z_u[i]      ~ normal(0, 1);
}


generated quantities {
  // Simulate random effects
  array[N] vector[2] z_u_rng;
  array[N] vector[2] u_rng;

  for (i in 1:N) {
    z_u_rng[i,1] = normal_rng(0, 1);
    z_u_rng[i,2] = normal_rng(0, 1);
    u_rng[i,1] = tau[1] * z_u_rng[i,1];
    u_rng[i,2] = tau[2] * z_u_rng[i,2];
  }

  // Squared and interaction terms
  vector[N] S2 = S .* S;
  vector[N] O2 = O .* O;
  vector[N] SO = S .* O;

  vector[Nobs] mu_i;
  vector[Nobs] lambda_i;
  vector[Nobs] sigma_i;
  vector[Nobs] y_sim;

  for (n in 1:Nobs) {
    int i = subj[n];

mu_i[n] = beta_mu[1] 
         + beta_mu[2]*S[i] 
         + beta_mu[3]*O[i] 
         + beta_mu[4]*S2[i]   
         + beta_mu[5]*SO[i] 
         + beta_mu[6]*O2[i] 
         + u_rng[i,1];

lambda_i[n] = beta_lambda[1] 
            + beta_lambda[2]*S[i] 
            + beta_lambda[3]*O[i] 
            + beta_lambda[4]*S2[i] 
            + beta_lambda[5]*SO[i]
            + beta_lambda[6]*O2[i] 
            + u_rng[i,2];

    
sigma_i[n] = exp(beta_e[1] 
                       + beta_e[2]*S[i] 
                       + beta_e[3]*O[i] 
                       + beta_e[4]*S2[i] 
                       + beta_e[5]*SO[i] 
                       + beta_e[6]*O2[i]);


    
    y_sim[n] = normal_rng(mu_i[n] + lambda_i[n] * lag_y[n], sigma_i[n]);
  }
}
2

Welcome to the Stan forums. It’s late and I might be missing something, but I see you’re transforming with the u array, but I don’t see it or anything else in the model block, it looks like you’d just be sampling from the prior.