Hello! I am wondering the best way to use a derived quantity from a one model as a response variable in a second model (while maintaining the error associated with each value). Specifically, I have a simple linear regression as my first model and then I take a new dataset and calculate the difference between observed values of y from predicted values of y based on the model parameters. I then want to use these differences as a response variable in a new model: a simple ANOVA-style where I look at delta y ~ treatment. Is it possible to do this in a joint model in STAN? How would I code it? I have my current code below and I added the code for the second model in a second model statement at the bottom. Thank you for any advice!

```
data {
int < lower = 1 > N; // Sample size for control regression
int < lower = 1 > n_treatment; // number of treatments
int < lower = 1 > n_PlateID; // number of plates
vector[N] y; // Response
vector[N] x; // Predictor
vector[N] x2; // Covariate
int < lower = 1 > n_resid; // Sample size for residual model
int<lower = 1, upper = n_treatment> Treatment[n_resid]; // Treatment variable
vector[n_resid] PlateID; // PlateID for random effect of plate
vector[n_resid] AllPR; // All final PR
vector[n_resid] AllInitial; // All initial PR
vector[n_resid] AllpHmin; // All pHmin
}
parameters {
real alpha; // Intercept
real beta; // Slope (regression coefficients for x)
real beta2; // Slope for x2
real < lower = 0 > sigma; // Error SD
real < lower = 0 > sigma_resid; // Error SD for residual model
real beta_resid; // beta for second model
real alpha_plate; // intercept for second model
}
model {
// conditional mean
vector[N] mu;
// linear combination
mu = alpha + x * beta + x2*beta2;
// likelihood function
y ~ normal(mu, sigma);
}
generated quantities {
vector[N] y_rep;
vector[n_resid] PR_expected;
vector[n_resid] PR_Resid;
vector[n_resid] y_resid;
for (i in 1:N) {
// generate predictions
real y_hat = alpha + x[i] * beta + x2[i] *beta2; // Posterior predictions
// generate replication values
y_rep[i] = normal_rng(y_hat, sigma);
} // The posterior predictive distribution
// calculate predicted PR based on pH from plate data
for (i in 1: n_resid){
PR_expected[i] = alpha + AllpHmin[i] * beta + AllInitial[i]* beta2; // calculate the expected value
PR_Resid[i] = AllPR[i] - PR_expected[i]; // calculate the residuals
}
}
// This is the second model statement where PR_Resid from above is the new y value
model {
// conditional mean
real mu_resid;
// generate model for residuals ~ treatment
mu_resid = alpha_plate + beta_resid*Treatment[i];
// likelihood
PR_Resid ~ normal(mu_resid, sigma_resid);
}
```