Regression on residuals within joint model

Question

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 \begin{equation} y_k = X_k\beta_k + \epsilon_k \end{equation} 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 \begin{equation} \epsilon_k = z_k\psi_k \sim N(0,\sigma^2_k) \end{equation} 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 \begin{equation} \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) \end{equation} 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: \begin{equation} y_k = X_k\beta + \epsilon_k \end{equation} 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?

Answer 1

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.