Regression on residuals within joint model

I have a regression model (part of a larger hierarchical model) where I wish to try construct a regression using the residuals of a regression.

To simplify, say we have $$j$$ in $$n$$ observations, a $$j \times r$$ design matrix $$X$$, and corresponding coefficient vector $$\beta$$, such that $$\begin{equation*} y = X\beta + \epsilon \\ \epsilon \sim N(0, \sigma^2) \end{equation*}$$

Assuming we have equivalent observations and coefficients for each $$k$$ of $$m$$ subjects, we can write $$$$y_k = X_k\beta_k + \epsilon_k$$$$ The residuals, $$\epsilon_k$$ are of interest to me. I would thus like to try to model differences in the residuals between subjects, using a vector of $$p$$ subject-specific covariates, $$z_k$$. This is where things get confusing for me. Now we would have $$$$\epsilon_k = z_k\psi_k \sim N(0,\sigma^2_k)$$$$ I cannot place a prior on $$z_k\psi_k$$ though, since $$z_k$$ can vary arbitrarily. I then thought that perhaps the solution is to have a nested regression for the residuals, such that $$$$\epsilon_k \sim N(z_k\psi_k,\sigma^2_k) \ \text{or} \\ \epsilon_k = z_k\psi_k + \tau_k, \ \ \tau_k\sim N(0, \sigma^2_k)$$$$ However, this seems non-identifiable to me. If we don't restrict $$\epsilon_k$$, then the values of $$X_k\beta_k$$ and $$\epsilon_k$$ can exchange, and the likelihood would be equivalent. There seems to be no way to disentangle the effect of the two sets of covariates.

It seems to me that if I restrict $$\beta_k$$ to be the same across all subjects I have some chance, ie: $$$$y_k = X_k\beta + \epsilon_k$$$$ since this induces some shrinkage towards a global value, but I am unsure if this would be enough to allow me to infer these effects. Is there an alteration I can make, or an alternative approach that would make this tractable?

• Within a frequentist mixed effects framework, one can have a model for the residuals in addition to the typical model for the means. The model for the residuals can have time-varying or person-level predictors. However, I am wondering if I am somehow misunderstanding what you want to do. Commented Mar 5, 2020 at 19:23
• Thanks Erik, could you perhaps link to an example in a paper or book? Commented Mar 5, 2020 at 19:24
• Hmm, upon fitting this model in stan ($y_k \sim N(X_k\beta_k + z_k\psi_k, \sigma^2_k)$ there doesn't seem to be any issues and I can recover $\psi$ almost exactly. Perhaps I'm overthinking things! Commented Mar 5, 2020 at 19:30
• Sure thing. The paper by Leckie et al. uses a Bayesian MCMC approach, albeit with largely diffuse priors. Commented Mar 5, 2020 at 19:43
• Thanks @alanocallaghan. Added a response! Commented Apr 5, 2020 at 21:56

One option is to model the residual variance as a function of various predictors in your regression model, including time-varying and time-invariant predictors. Examples of this approach and its implementation in nlme can be found here and also here. A similar approach can be utilized in Bayesian setups as well, e.g., this article by Leckie et al.
The residual variance is modeled as a log-linear function of the time-varying and time-invariant predictors. In a frequentist mixed effects modeling paradigm, the residual variance ($$\sigma^2$$) is estimated from the model, but in this approach, its log is modeled directly as a function of predictors:
log($$\sigma^2$$) = $$\alpha_{0j}$$ + $$\alpha_{1j}X_{ij}$$.
As noted in the Leckie et al paper, if $$\alpha_{0j}$$ + $$\alpha_{1j}X_{ij}$$ are both equal to 0, the residual variance is constant across occasions (level 1 units in multilevel modeling parlance), which is the typical assumption of linear regression models.