Best practices for contraining specific parameters to a constant value

Modeling
1

I’m comparing a series of models that are nested within a maximally complex model. For simplicity, imagine I want to consider two models, one in which the parameter mu is constrained to a single value across levels of stratum, and a “random intercept model” where observations in each stratum have mean mu_stratum, which is normally distributed about mu. Ideally I’d have something like this:


data {
  int<lower=0,upper=1> INCLUDE_RANDOM_INTERCEPT;
...
}

parameters {
  ...
  vector[N] mu_stratum;
  if (INCLUDE_RANDOM_INTERCEPT==1) {
    real<lower=0> sd_mu_stratum;
  }
...
}

transformed parameters {
  if (INCLUDE_RANDOM_INTERCEPT==0) {
    for (i in 1:N) {
      mu_stratum[i] = mu;
    }
  }
  ...
}

model {  
  ...
  if (INCLUDE_RANDOM_INTERCEPT==1) {
     mu_stratum ~ normal(mu, sd_mu_stratum);
  }
   ...
}

However, I can’t have conditional parameter definitions or set transformed parameters to a single global value. The only way I can get something like this to work is to force sd_mu_stratum to be very small when the flag is off, but that feels hacky. Is there a nice way to do this?

2

Stan allows size 0 arrays, so would something like this approach work?

parameters {
  real mu;
  array[N * INCLUDE_RANDOM_INTERCEPT] real mu_stratum;
  array[INCLUDE_RANDOM_INTERCEPT] real<lower=0> sd_mu_stratum ;
}
transformed parameters {
  array[N * (1 - INCLUDE_RANDOM_INTERCEPT)] real mu_equal;
  if (INCLUDE_RANDOM_INTERCEPT == 0) {
    for (n in 1:N) mu_equal[n] = mu;
  }
}
model {
  if (INCLUDE_RANDOM_INTERCEPT) {
    mu_stratum ~ normal(mu, sd_mu_stratum[1]);
  }
}
3

Actually maybe you don’t even need the transformed parameter mu_equal from my example. You could just use mu in the model block if INCLUDE_RANDOM_INTERCEPT == 0.

4

Oh wow, I didn’t know about the length zero array business, that’s really nifty. A very clean way of ensuring you create a bunch of parameters that you have to place priors on but don’t contribute to the likelihood otherwise.

And you’re right, the mu_equal is unnecessary

5

I’m afraid there’s not a nice way to do exactly what you’re asking, which is to pin variables declared as parameters to specific variables. All you can do is the kind of thing @Jonah’s suggesting. There’s a slight hitch—@Jonah gave you the centered parameterization, which will only work if you have tight priors and/or a lot of data. Otherwise, you need to do something like the following, which will use a non-centered parameterization. It takes in arguments for the fixed priors—those will be ignored if vary == 1. Then mu_mu and sd_mu will be of size 0 if vary = 0.

data {
  int<lower=0, upper=1> vary;
  int<lower=0> N;  // number of vars
  real mu_mu_fixed;  // prior location if non-varying
  real sigma_mu_fixed;  // prior scale if non-varying
}
parameters {
  vector[vary ? N : 1] mu_raw;
  vector[vary] mu_mu;  // prior location if varying, if not size = 0
  vector<lower=0>[vary] sd_mu;  // prior sd if varying, if not size = 0
}
transformed parameters {
  vector[N] mu;
  if (vary) {
    mu = mu_mu[1] + sd_mu[1] * mu_raw;
  } else {
    mu = rep_vector(mu_mu_fixed + sd_mu_fixed * mu_raw, N);
  }
}
model {
  mu_raw ~ normal(0, 1);
  ...
}

Then you write the rest of the model using the transformed parameter mu of size N for convenience.

6

I see, this works just fine. You’re right, it is a bit cumbersome relative to some sort of macro for pinning the parameters, but this is much more efficicent than my hack, which was also slow because you end up sampling a bunch of latent variables you don’t actually use

7

@Bob_Carpenter good point about the parameterization. It’s definitely a bit uglier with the non-centered parameterization, but that’s usually the one we want.