Mixing observed and predicted values within explanatory variables - is it ok?

Suppose that out of a population we measure 100% of $x_i$ and 50% of $y_i$.

Now, we use the 50% measured pairwise $x_i$ and $y_i$ to adjust model 1:

$$Y = \beta_0 + \beta_1X + \epsilon_1$$

Based on this model we predict the other 50% $\hat{y_i}$. Then, we have 100% of measured $x_i$, 50% of measured $y_i$ and 50% of predicted $\hat{y_i}$.

If one wants to use $X$ and $Y$ information to model $W$ (model 2):

$$W = \beta_2 + \beta_3X + \beta_4Y + \epsilon_2$$

Is it correct to mix measured and predicted $Y$ on model 2? Or should we use model 1 to predict 100% of $\hat{y_i}$ and use them on model 2?

I understand that with measured $y_i$ there is only one type of error which is the measurement error, while with $\hat{y_i}$ there is also the prediction error.

My intuition says that keeping the 50% of measured $y_i$ to predict $W$ will minimize the errors but I don't know if mixing up $y_i$ and $\hat{y_i}$ on model 2 is incorrect for some statistical reason.

• I've heard this type of model called a mediation model, but I'm not sure if there's consensus on how to deal with missing data on the mediator variable. – Neal Fultz Oct 6 '14 at 18:59

The OP has not specified an estimation procedure, so I will consider the question for the case of linear least-squares regression. I will also ignore the constant term-let's say variables are already centered on their means.

Should we use model 1 to predict 100% of $\hat{y_i}$ and use them on model 2?

If we do that, then the regressor matrix in the regression for $W$ will exhibit perfect collinearity (since the one regressor will be $X$ and the other $Χ \hat \beta$ throughout). In this case we could instead regress $W$ on $X$ alone and obtain estimates of composite coefficients (see this post).

Is it correct to mix measured and predicted Y on model 2?

It does not violate a priori any mathematical rule, and so it is "correct". Whether the estimator has good properties it is a matter to be investigated.

In the regression for $W$ the postulated relation is that $W$ depends on the actual realizations of $Y$. Denoting $Z$ the theoretical regressor matrix, and writing it in block form (for the two halves of the population) we have:

$$\mathbf Z_{n \times 2} = \left [\begin{matrix} \mathbf x_1 & \mathbf y_1 \\ \mathbf x_2 & \mathbf y_2 \end {matrix}\right]$$

while the actual regressor matrix used will be

$$\mathbf {\hat Z} = \left [\begin{matrix} \mathbf x_1 & \mathbf y_1 \\ \mathbf x_2 & \hat \beta\mathbf x_2 \end {matrix}\right],\;\;\mathbf {\hat Z}' = \left [\begin{matrix} \mathbf x_1' & \mathbf x_2' \\ \mathbf y_1' & \hat \beta\mathbf x_2' \end {matrix}\right]$$

Then the $\gamma$ vector that has the coefficients in the regression for $W$ will be estimated by OLS as

$$\hat \gamma = \left(\mathbf {\hat Z}'\mathbf {\hat Z}\right)^{-1}\mathbf {\hat Z}'\mathbf w = \left(\mathbf {\hat Z}'\mathbf {\hat Z}\right)^{-1}\mathbf {\hat Z}'\left(\mathbf Z\gamma + \mathbf u\right)$$

$$=\left(\mathbf {\hat Z}'\mathbf {\hat Z}\right)^{-1}\mathbf {\hat Z}'\mathbf Z\gamma +\left(\mathbf {\hat Z}'\mathbf {\hat Z}\right)^{-1}\mathbf {\hat Z}'\mathbf u \tag{1}$$

Writing $\mathbf y_2 = \hat \beta\mathbf x_2 + \mathbf e_{y_2}$ where $\mathbf e_{y_2}$ is the prediction error We have

$$\mathbf {\hat Z}'\mathbf Z = \mathbf {\hat Z}'\left [\begin{matrix} \mathbf x_1 & \mathbf y_2 \\ \mathbf x_1 & \hat \beta\mathbf x_2 + \mathbf e_{y_2} \end {matrix}\right]=$$

$$=\mathbf {\hat Z}'\left( \left [\begin{matrix} \mathbf x_1 & \mathbf y_2 \\ \mathbf x_2 & \hat \beta\mathbf x_2 \end {matrix}\right] + \left [\begin{matrix} \mathbf 0 & \mathbf 0 \\ \mathbf 0 & \mathbf e_{y_2} \end {matrix}\right]\right) = \mathbf {\hat Z}'\mathbf {\hat Z} + \left [\begin{matrix} 0 & \mathbf x_2'\mathbf e_{y_2} \\ 0 & \hat \beta\mathbf x_2'\mathbf e_{y_2} \end {matrix}\right] \tag{2}$$

Inserting $(2)$ into $(1)$ we get

$$\hat \gamma = \gamma + \left(\mathbf {\hat Z}'\mathbf {\hat Z}\right)^{-1} \left [\begin{matrix} 0 & \mathbf x_2'\mathbf e_{y_2} \\ 0 & \hat \beta\mathbf x_2'\mathbf e_{y_2} \end {matrix}\right]\gamma +\left(\mathbf {\hat Z}'\mathbf {\hat Z}\right)^{-1}\mathbf {\hat Z}'\mathbf u \tag{3}$$

If in the regression of $W$ we have assumed strictly exogenous regressors, then taking the expected value in $(3)$ conditional on the observable random variables, the 3d term (as usual) is zero, but also the 2nd term is zero, since the prediction error is conditionally orthogonal to the predictor. So $\hat \gamma$ is unbiased. It is also consistent. In these aspects it is equivalent to the OLS estimator we would obtain if we could observe all $Y$. But the variance will be different.

• The estimation procedure I'm using is indeed OLS. Tks for your help. – Andre Silva Oct 8 '14 at 1:13

Using half of the predicted $y$s would mean that you strongly underestimate the variance of $\beta_4$ in model 2 as the prediction would not be accounted for!

What you suggest would be some form of single imputation which in general leads to underestimation of variance.

• Welcome to Cross Validated Stack Exchange. Could the underestimation of the variance of $\beta_4$ lead to overfitting? – Andre Silva Oct 7 '14 at 20:23
• @matthew - Can you quantify what you mean by strongly underestimate? Let’s say, in the limiting case where $\epsilon_1=0$, or, more realistically, where $\epsilon_1$ is small. Surely the increase in sample size by taking 50% more samples will more than make up for small ϵ1. The effect will of course depend on the sample size, on $\epsilon_1$ and on $\epsilon_2$. – martino Oct 9 '14 at 11:02

You can think of this in terms of learning a feature of an ensemble in machine learning.

Normally an ensemble would say have 5 models all trained to predict W and then the ensemble model would combine these models to make a better prediction on W. Normally by something simple like a vote. But there is no reason why these 5 predictions could not be combined in a more complex way using any ML algorithm.

Now you want to learn a model to predict Y and then use this as a feature to predict W. Nothing wrong in that.

You are using one model -Y as an input feature to model W.

I think this is also related to semi supervised learning, where you can use proxy measures- To predict the measure you want.

For example in an image classification you train one model to predict sea(Your labelled Y for example), another to predict sand (another model say Z, where you have labels) and then combine these models to predict beach. (Your W where you have limited labelled data). Where W would be a model built out of Y and Z.

• I am not much familiar with machine learning. Is the message from your answer that it would be ok to use the predicted $\hat{y_i}$ together with the measured ones on model 2? Is that it? – Andre Silva Oct 7 '14 at 14:48
• It is okay, but you should make sure you have not 'cheated' in your validation of the performance on model w. So you need to make sure that any data used to learn model Y and then W is separate from your way of measuring the performance of model w as per usual. i.e if you were performing 10 fold cross validation, to measure the performance of model w you would learn a new model y in each fold. (So you might have cross validation inside cross validation) – user27815 Oct 7 '14 at 15:00

Nice question. I can’t say I have come across this before but the following paper on mediation, compounding and suppression might be useful. I can’t see how mixing up $y_i$ and $\hat{y_i}$ in model 2 is incorrect for any statistical reason. Whether or not it leads to better prediction would depend on the relationships between the different variables. To begin with, you can determine an expected error distribution where the predicted $y_i$'s are included and excluded and see where that leads you.

If I was faced with this, I would start by testing the models on a generated dataset to get a better understanding of how the inclusion of the known and predicted $y_i$ values compare to including only the known $y_i$. You can pick values for your model parameters and randomly generate a set of $x_i$ and $y_i$ - discard half of the $y_i$'s. Given this dataset you can calculate the parameters and error distributions for model 1 and model 2. You can repeat this for many generated datasets to help understand the characteristics of the models when the predicted $y_i$’s are included or excluded. The insight gained should help inform your decision on the best approach.

• These are good ideas (+1). I will take a look into that paper. Thank you. – Andre Silva Oct 7 '14 at 14:40
• It’s always worth generating data to gain an understanding of how models work. I expect in your case the answer as to whether it is a good idea to include the predicted $y_i$’s will come down to the relative size of $\epsilon_1$ and $\epsilon_2$ and just how much of the variance is explained in model 2 by $X$ and $Y$. A different problem with different models may come to different conclusions. – martino Oct 7 '14 at 14:51

The problem that you pose reminds me a bit of the classic problem of fitting multivariate Gaussian mixture models (GMM), which is the archetypal application for the expectation maximization (EM) algorithm.

At a high level, it's instructive to think a little bit about how the EM algorithm solves the GMM fitting problem, because some of the basic EM concepts (if not the actual EM algorithm directly itself) are quite applicable and can easily be borrowed and modified to address your situation.

By way of comparison with your own problem, in a GMM fitting context, the basic problem statement is this: you are given a data set consisting of bunch of random vector variables $\vec{x}_{i}$ which are drawn from a probability distribution which is a sum of several multivariate Gaussian distributions. Each "mode" of the distribution has several defining parameters, which you are trying to estimate based upon the data: the mean and covariance of course, plus a relative amplitude parameter which gives the relative size of each mode, compared to the others. These defining parameters are directly analogous to the $\beta_{k}$ in the above problem. The GMM problem also contains an entire second set of what are usually termed latent variables $z_{i}$, which essentially are a postulated set of additional variables that would go a long way toward explaining the observations $\vec{x}_{i}$, if only we could observe them. Unfortunately however, they are essentially "missing information", which is why we refer to them as latent.

In the GMM problem, the information that one imagines residing within the latent variables would theoretically be able to tell you (provided you could actually observe them) which mode within the mixture was specifically responsible for giving rise to each of the observations $\vec{x}_{i}$. If the nodes are numbered $1, 2, 3,..., n$, then each of the $z_{i}$ will be assigned a corresponding value of $1,...,n$ as well, indicating which node the $i$th data point came from. (Variations upon the basic EM scheme also include alternate definitions of the $z_{i}$ which attempt to assign "degree or likelihood of belonging", so that we do not make a hard assignment of each data point to a single mode, but that's starting to get beyond the scope of the discussion here.)

Anyway, in the context of this particular stackexchange question, the $y_{i}$ data in the original question are a sort of loosely analogous to "half-latent" variables in the parlance of the EM algorithm: i.e., some of the values are observed, and some are not, so they are neither fully observed nor fully latent.

For the case with fully latent (i.e, 100% non-observable) variables, the way that the EM algorithm actually solves the GMM problem is to start by randomly assigning a value to each latent variable $z_{i}$, then estimate the best fit parameters $\beta_{k}$ (i.e., the mean, covariance and relative amplitude in GMM), assuming that all of the initial $z_{i}$ assignments were correct (of course they're not correct at all, really, since we just assigned the starting values at random, but don't worry about that yet). Next, using the recently estimated $\hat{\beta}_{k}$, the algorithm estimates the most likely values (i.e., it sort of re-predicts) for the missing/latent variables $\hat{z}_{i}$. Essentially, the algorithm goes back and forth between two complementary steps: first, it tries to estimate the most likely values $\hat{\beta}_{k}$ for the true $\beta_{k}$, based upon its most recent estimates of the $\hat{z}_{i}$, and then it tries to update its estimate of the $\hat{z}_{i}$ based upon the newest estimate of the $\hat{\beta}_{k}$. The algorithm continues to go back and forth in this fashion, until some type of convergence criterion is met for both quantities.

So, in analogy with the EM algorithm, here's what I'd recommend for you: either path that you suggest (either mixing 50% observed $y_{i}$ with 50% predicted $\hat{y}_{i}$, or alternately, using 100% predicted $\hat{y}_{i}$) is equally acceptable, because actually, in either case, it's only an initial step. What you should do next, after you've obtained initial estimates for both the $\hat{\beta}_{k}$ and the $\hat{y}_{i}$, is follow the example of the EM algorithm: go through several alternating stages of refining your estimates of both $\hat{\beta}_{k}$ and $\hat{y}_{i}$, by deriving each new estimate for one set of quantities based upon the previous estimate of the others, continuing back and forth until you reach convergence on both. Of course, in iterating back and forth, when predicting the next update for $\hat{\beta}_{k}$ from the previous $\hat{y}_{i}$, you should substitute the real observed $y_{i}$ wherever they are available, so in that sense, I suppose that my recommendation is indeed to mix together 50% predicted $\hat{y}_{i}$ with 50% observed $y_{i}$, but since the point of the algorithm is to iteratively pursue solution convergence, it doesn't necessarily matter all that much what type of $y$ values you use at the beginning in order to start the ball rolling.

Doing it this way has an additional benefit as well: if anyone ever asks you to justify your approach on theoretical grounds, you can say it's essentially just an extension of the EM algorithm, which relies upon the concept of local convergence in order to reach a stable answer. Since both the EM algorithm as well as local convergence are established principles of iterative/recursive parameter estimation, this methodology rests on solid ground, theoretically speaking.

• It seems an interesting approach and I appreciate the detailed answer, which I'll study. – Andre Silva Oct 7 '14 at 20:10