Help needed in building a causal prediction model using rstan

Modeling
1

Hi rstan community!

I’m working on a causal prediction challenge where I need to build a prediction model using rstan. I have an outcome variable, Y, which follows a poisson distribution. My goal is to predict Y using a combination of binary, categorical, and continuous predictors. These predictors all have a causal structure, and I want to ensure that the model accounts for the relationships between the predictors and the outcome accurately (see attached DAG).

Rplot01
Rplot01885×758 46.5 KB

I’ve made progress in constructing the base model in stan, thanks to online resources including this forum. However, I’m unsure if this model captures my DAG structure above. Any assistance, advice, or resources you could provide would be greatly appreciated. Thank you in advance!

data {
  int<lower=0> N;
  int Y[N];
  int chl[N];
  vector[N] sam;
  int msl[N];  // Binary variable for msl
  vector[N] rnfl;
  vector[N] m;
  vector[N] p;
  vector[N] dr;
  vector[N] ptime;
  vector[N] wt;
  int aww[N];  // Multinomial variable for aww
  int mlr[N];  // Multinomial variable for mlr
  int aca[N];  // Multinomial variable for aca
}

parameters {
  real alpha_Y;
  real alpha_chl;
  real alpha_sam;
  real alpha_msl;
  real alpha_rnfl;
  real alpha_m;
  real alpha_p;
  real alpha_dr;
  
  real beta_Y;
  real beta_chl;
  real beta_sam;
  real beta_msl;
  real beta_rnfl;
  real beta_m;
  real beta_p;
  real beta_dr;
  
  real alpha_aww[3];  // Coefficients for aww multinomial regression
  real alpha_mlr[3];  // Coefficients for mlr multinomial regression
  real alpha_aca[3];  // Coefficients for aca multinomial regression
}

model {
  sam ~ normal(alpha_sam, 1);
  rnfl ~ normal(alpha_rnfl, 1);
  m ~ normal(alpha_m, 1);
  p ~ normal(alpha_p, 1);
  dr ~ normal(alpha_dr, 1);
  
  for (n in 1:N) {
    chl[n] ~ bernoulli_logit(alpha_chl);
    msl[n] ~ bernoulli_logit(alpha_msl);
    Y[n] ~ poisson_log(alpha_Y + beta_Y * wt[n] + log(ptime[n]));
    
    aww[n] ~ categorical_logit(alpha_aww);  // Multinomial likelihood for aww
    mlr[n] ~ categorical_logit(alpha_mlr);  // Multinomial likelihood for mlr
    aca[n] ~ categorical_logit(alpha_aca);  // Multinomial likelihood for aca
  }
}

generated quantities {
  int predicted_Y[N];
  
  for (n in 1:N) {
    predicted_Y[n] = poisson_log_rng(alpha_Y + beta_Y * wt[n] + log(ptime[n]));
  }
}