interpretation of residuals in the linear regression Assume we have a population of $N$ people and each people is denoted by $i\in\{1,2,...,N\}$. Among these $N$ people, there are $M$ male and $N-M$ female. 
We have two conditions denoted by $0$ and $1$ and let $D_i\in\{0,1\}$ be the random variable indicates the condition for individual $i$. For each individual $i$, it has a value denoted by $\epsilon_{0i}$ under condition $0$ and it has a value denoted by $\epsilon_{1i}$ under condition $1$. Hence, we can write
$\epsilon_i = \epsilon_{0i}+D_i(\epsilon_{1i}-\epsilon_{0i})$,
as the value of each individual $i$ in the population. Assume $E[\epsilon_{0i}]=E[\epsilon_{1i}]=0$.
Now, in this specific population we just described, the people's gender is fixed and it has $M$ male and $N-M$ female in the population. If we use $X_i$ to denote each individual's gender in this specific population. Is $X_i$ a random variable? or not?
If we run a linear regression 
$\epsilon_i=cX_i+u_i$, where $u_i$ is defined as $\epsilon_i-cX_i$ and $c$ is the OLS coefficient, then what is the random variable $u_i$? What does $u_i$ mean intuitively? Any physical meaning of $u_i$? And what could be the shape (distribution) of $u_i$ in terms of $\epsilon_{0i}$ and $\epsilon_{1i}$ (or truncated version of them)?
 A: The interpretation of residuals in regression models is a philosophical construct more than anything, even under presumptive correct model specification. It all stems from the notion that random errors arise from millions of unobserved and/or unadjusted causes of variation, all averaging out into an "error term". This "error term" often takes a normal distribution because it is a maximum entropy distribution, but it is not actually necessary statistically (or philosophically) for residuals to assume such a distribution. The asymptotic normality of averages ensure approximately correct finite sample CIs and PIs. As a frequentist, I believe that if all such covariates were somehow included in a model for an outcome, it would be deterministic and yield $R^2=1$ each time under repeated experiments.
The gender of participants entered into any kind of study is a natural blocking factor. Thus, when we think of the generalizability of results, we can view the gender as a "fixed" value which, under replications of the study, would not actually vary (the outcome is what would vary). For that reason, I think it is fair to call $X$ a covariate or blocking factor but not a random variable. We can obtain somewhat different (larger) error estimates in models by taking $X$ to be random and assume either an empirical or known density. It is important to note that most conventional inference assumes the regressors to be fixed quantities either by experimental design or population representativeness. 
You describe more fully the impact of $D$ as factor whose influence is not appropriately "handled" by simple adjustment in a multivariate model. Under the assumptions you give, it does not give rise to identically distributed errors, even though the mean difference is 0 for $D=0$ vs $D=1$. Therefore the errors take a mixture distribution. If we know this, there are statistical tools which could address this for better inference: If we observe $D$, a model which stratifies by $D$ may give more efficient estimates of treatment effects, or an appropriate transformation may achieve homoscedasticity. Whether or not we observe $D$ we can estimate the multivariate normal likelihood using an EM approach to maximize a likelihood for mixture normal errors. If we ignore $D$, the result is surprisingly similar to other omitted covariate findings. Basically, the standard error of the residuals may be over or under estimated and won't generalize to another sample in which the proportion of $D$ differs, or $D$'s relation with $X$ differs.
