Constructing counterfactual analysis with generated quantities or another model

Modeling
1

Hi all. This might be a basic/stupid question, but it’d be nice to check if it makes sense. I’m working on a hierarchical model where I want to generate counterfactual simulations by modifying fitted parameters. What’s the best practice for this? I’m currently recomputing the model outside Stan playing with the parameters and the draws, but it becomes a bit messy.

My idea now would be to fit two models. First, I’d fit the original model and get the fitted parameters. And in the second one I’d use the parameters as data (manipulating them to construct counterfactuals) and the original data as parameters.

To exemplify this: suppose I have this original model that I fit with my data:

data {
  int<lower=0> N;
  array[N] int<lower=0> y;
}

parameters {
  real<lower=0,upper=1> theta;  // success probability
  real<lower=0> lambda;          // exposure intensity
}

model {
  theta ~ beta(2, 2);
  lambda ~ exponential(0.5);
  
  for (n in 1:N) {
    y[n] ~ poisson(lambda * theta);
  }
}

After fitting this, I would run the following counterfactual model, feeding stan with manipulated parameters as data:

data {
  int<lower=0> N;
  
  // Parameters as data for counterfactuals
  real<lower=0,upper=1> theta_cf;
  real<lower=0> lambda_cf;
}

generated quantities {
  array[N] int y_cf;
  
  // What if theta was 0.8?
  for (n in 1:N) {
    y_cf[n] = poisson_rng(lambda_cf * theta_cf);
  }
}

My questions are:

  1. Is this the recommended approach for counterfactual analysis: creating a separate Stan file where parameters become data?
  2. Or is there a more elegant way using generated quantities, for instance?

Thanks in advance!

2

For using fitted parameters for post-fitted calculations, cmdstanr::generate_quantities() works well. You pass in the fitted draws and data, and it returns the things calculated in generated quantities.

3

I see. But the problem is I’d like to play around with the fitted parameters. So suppose i have parameters alpha, beta, theta. So I’d like to see (post-fitted) what would happen if i forced alpha = beta or beta = theta. Since I’d be doing a bunch of these computations, I thought the way I proposed would be easier to perform these counterfactuals.