Impact of substituting default values for predictors in Regression Model Suppose I have a regression model (it could be any model, but I am mainly interested in multiple linear regression) that is designed to predict Y as a function of a set of X predictors.
In my field (civil engineering), I often encounter situations where there is a complex, generalized prediction equation involving around 3 or more predictors, but when it comes to implementation in practice, it is common to not have any data for 2 to 3 of the predictors. That is, there may be some missing data on individual points, but for some predictors there are NO data available. It is common practice in such a case to assume fixed "default values" for missing predictors.
I am of the opinion that (a) using fixed default values will underestimate true variance in Y and thus underestimate risk; and (b) it would be better to use the available (site specific) data (which is quite sparse), to fit a simplified model. The data is normally highly variable and the original data set (on which the model was developed) is not available. So using Imputation to get missing values is not really viable.
My question: (a) Is my suspicions above generally correct? (b) are there any landmark studies around this issue to clarify risks or best practice around missing predictor variables.
 A: With respect to (a), yes it will underestimate the true variance in $Y$.
Let's assume that we live in a world where we actually have the 'true' multiple linear regression model (that is we know what $\beta_0,\beta_1,\beta_2$ are), and there are only two predictor variables. The multiple linear regression equation looks like this:
$$Y_i=\beta_0+\beta_1X_{i,1} + \beta_2X_{i,2} + \epsilon_t$$
To keep things simple assume $X_{i,1}, X_{i,2}$ are independent for all $i$, and that $$\begin{align}
X_{i,1} &\sim \mathcal{N}(\mu_1,\sigma^2_1)\\
X_{i,2} &\sim \mathcal{N}(\mu_2,\sigma^2_2)\\
\epsilon_i &\sim \mathcal{N}(0,\sigma^2_\epsilon)
\end{align}$$
We can look at two cases, where we have the value of $X_{i,2}$ and when we do not. Remember that the predictions your regression model is making is the expected value of $Y$ given that we know what all the covariates are. That is
$$
\begin{align}
E[Y_i|X_{i,1},X_{i,2}] &= E[\beta_0+\beta_1X_{i,1} + \beta_2X_{i,2} + \epsilon_t|X_{i,1},X_{i,2}]\\
&= \beta_0+\beta_1X_{i,1} + \beta_2X_{i,2} \\
&= Y_i\\\\
E[Y_i|X_{i,1}] &= E[\beta_0+\beta_1X_{i,1} + \beta_2X_{i,2} + \epsilon_t|X_{i,1}]\\
&= \beta_0+\beta_1X_{i,1} + \beta_2 \mu_2 \\
\end{align}
$$
So we first see that even if we have the perfect model our point estimate is less accurate (unless $X_{i,2} = \mu_2$).
With respect to the conditional variance we see something similar (remember the variance of a constant is 0)
$$
\begin{align}
\text{Var}[Y_i| X_{i,1},X_{i,2} ] &= \text{Var}[\beta_0+\beta_1X_{i,1} + \beta_2X_{i,2} + \epsilon_t|X_{i,1},X_{i,2}]\\
&= \sigma^2_\epsilon\\\\
\text{Var}[Y_i| X_{i,1}] &=\text{Var}[\beta_0+\beta_1X_{i,1} + \beta_2X_{i,2} + \epsilon_t|X_{i,1}]\\
&= \beta_2^2 \sigma^2_2 + \sigma^2_\epsilon 
\end{align}
$$
So yes we see that if we are conditioning on fewer predictor variables (i.e. we don't know some of them), we become less certain in what we are trying to predict (due to the higher variance).
A: I gather that you are interested in $\text{Var}(Y|X_1=x_1)$ when you have $\text{P}(Y|X_1=x_1,X_2)$
The determination of the variance of random variable from its conditional dependence on an additional random variable is described by the LoTV, the Law of Total Variance, also known as Eve's Law due to the way the Es and Vs line up in the expression. 
In your case, where $X_1$ is specified but $X_2$ is not, Eve's Law gives you 
$$
\text{Var}(Y|X_1=x_1)=\text{E}[\text{Var}[Y|X_1=x_1,X_2]] + \text{Var}[\text{E}[Y|X_1=x_1,X_2]).
$$
As it is quite common to be confused over what is being averaged over, so I'll rewrite the same law to try to make it explicit
$$
\text{Var}_{Y|X1}(Y|X_1=x_1)=\text{E}_{X_2}[\text{Var}_{Y|x_1,X_2}[Y|X_1=x_1,X_2]] + \text{Var}_{X_2}[\text{E}_{Y|x_1,X_2}[Y|X_1=x_1,X_2]).
$$
where references to $\text{P}(Y|X_1,X_2)$ is really just $\text{P}(\epsilon)\sim N(0,\sigma^2)$ since once you've specified  $X_1$ and $X_2$, the only moving part left is the irreducible noise. Also, I've simplified $\text{P}(X_2|X_1=x_1)=P(X_2)$ as I assume that the predictors are independent. (Is that the case?)
The first term is computed in two phases. First, you need to find the variance of the fully conditional model. That is, you'll find the variance of the underlying distribution of $Y$ when $X_1$ and $X_2$ are specified. This variance should just be your irreducible noise $\sigma^2$, which is presumably not a function of the predictors $X_1$ and $X_2$. (If you don't satisfy this Gauss-Markov assumption, then the linear model itself must change.) When you finish the first phase,  you need to calculate the expectation of that variance over a credible distribution of $X_2$ values. If, as we posit with the Gauss-Markov assumptions,  you benefit from homoskedasticity w.r.t. $X_2$, then $\text{E}[\sigma^2|X_1,X_2]=\sigma^2(X_1)=\sigma^2_{\epsilon}$. (David uses this result and all associated assumptions in his answer.)
The second term is also evaluated in two phases. First you find how $\text{E}[Y|X_1, X_2]$ depends on $X_2$, then you find the variance of that function respect to the underyling density $\text{P}(X_2),$ which is very likely a significant contributor to variance without violating assumptions. This is what has you worried! And rightly!
The crux of the problem is finding some trustworthy $\text{P}(X_2).$ (Perhaps $\text{P}(X_2|X_1)?)$ You don't have $X_2$ and you are loath to impute it, but I hope you are comfortable selecting (and getting general consensus from  your colleagues for) a prior $P(X_2)=\pi(X_2)$ which can be used to finish your calculation. This is a classic problem in Bayesian inference: to find $\text{Var}(Y|X_1=x_1)$ when you have $\text{P}(Y|X_1=x_1,X_2)$ and $X_2$ remains in play. This problem is ripe for numerical evaluation, especially if $X_1$ and $X_2$ covary. That is, you can calculate the variance by sampling from the posterior of $\text{P}(Y|X_1=x_1)$. You can sample from the posterior by sampling from the prior $\text{P}(X_2=x_2)$ and computing $\text{P}(Y|X_1=x_1,X_2)P(X_2=x_2).$ 
While you may find many papers that use conditional variance, I don't think you'll find a landmark paper about conditional variance. It is sufficient to find a section of a textbook about that.
As an aside, the popularity of the name Eve's Law has led to the renaming of its mate LoTE Law of Total Expectation as Adam's Law.
A: Are you familiar with metrology or measurement science, and the concept of uncertainty? Yes, if you take a default value to be a constant, without variance or uncertainty, then the use of the default value in place of your incomplete data will underestimate the uncertainty in your outcome. However, if you properly account for the uncertainty in the default value, then this uncertainty will become part of the overall uncertainty or quantified doubt about the value of your outcome measurand. Even so, whether you have greater uncertainty from using the default value or from using the incomplete predictor data is a matter of your particular situation.
A: As per our conversation in the comments section of your original post, it is extremely important to note, that the training data has been very carefully curated at great expense, and thus your unlabelled test examples, which contain missing values for some features, can be assumed to have been "sampled from the same distribution as the training set" 
For this reason, I would suggest the following if computationally feasible. When making a prediction for any particular data point for which you have access to a subset of the features, you train a model using your training data, but only giving yourself access at train time to the exact set of features which are provided in the data point you wish to label. 
This does, provided there are enough features and your data can be missing in enough different permutations, mean you can't feasibly pre-train all of the models you might need, so you'll have to train a model every time you want to make a prediction. 
If this isn't computationally feasible, I'd suggest training one model on all features in your training set, and then filling the missing values in the data point you need to make a prediction for. The simplest way to fill these values is with their mean/median (probably best to use mean if you're using an analytic method like linear regression and median if you're using tree-based methods).  That said, you will likely get better results if you fill your missing values based on conditional means (i.e. based on the features that you do have, what are the likely values of those that you don't have). There are many ways of doing this based on deep learning and it's an active field of research, but a more simple way of doing this would be to use KNN. A sketch of how this works:
For your data point with some missing features, which you wish to make a prediction for:
1: Calculate its Euclidean distance from all points in your training data, based on the features you do have access to
2: Sort your training data by that feature and select the K closest values
3: For each of the missing features, fill with the mean of these K training instances 
This is clearly not super computationally efficient, having to calculate the distance between a point, and every training example, before being able to make a prediction with it, but it's likely faster than training a model every time you wish to make a prediction as per my first suggestion, so it's a decent half-way house. 
