What is the correct way of adding bias terms in the residuals of the linear regression model?

First, I fit a linear model:

$y=beta_0 + beta_1x_1+beta_2x_2 + epsilon$

Now I want to visualize $y$ after the effects of $x_1$ and $x_2$ have been removed or adjusted. I can visualize the $y$ vs. $x_1$ or $x_2$ relationship by using only the residuals $epsilon$. The problem is I want to add the bias term in the residuals. Say, I want to plot the adjusted $y$ as a boxplot with respect to another independent variable (e.g. diagnostic group).

For now, I am adjusting the effects of $x_1$, $x_2$ on $y$ as below:

$y_a=y – beta_1(x_1-bar{x_1}) – beta_2(x_2-bar{x_2}) qquad (1)$

Here, for one data point, I am defining the effect of $x_1$ as the change in $y$ caused by the difference of $x_1$ from the mean of $x_1$ i.e. $bar{x_1}$.

After few algebraic manipulations:

$y_a= (beta_0 + beta_1 bar{x_1} + beta_2 bar{x_2}) + epsilon
= bias + residuals qquad (2) $

First, I am not 100% convinced with myself with this technique. However, this article https://surfer.nmr.mgh.harvard.edu/ftp/articles/buckner2004.pdf also uses this technique for covariate adjustment (see Equation 1 on Page 728).

Question1: Is this technique correct? and why if yes/no? Or asking the same question based on equation 2: Is the bias term $(beta_0 + beta_1 bar{x_1} + beta_2 bar{x_2})$ added to the residuals is correct?

Let’s assume the above adjustment technique is correct.
Let’s say $x_1$ is a categorical variable with more than two levels. How to calculate the mean of $x_1$ ($bar{x_1}$)?

Question2:
How to calculate the mean of a categorical variable? To be strict, it doesn’t even makes a sense to calculate the mean or any summary statistic off a categorical variable. Is there any workaround for this?