Biased estimator for regression achieving better results than unbiased one in Error In Variables Model

Question

I am working on some syntatic data for Error In Variable model for some research. Currently I have a single independent variable, and I am assuming I know the variance for the true value of the dependent variable.

So, with this information, I can achieve an unbiased estimator for the coefficient of the dependent variable.

The model:

$\tilde{x} = x + e_1$
$y = 0.5x -10 + e_2$
Where:
$e_1\text{~}N(0,\sigma^2)$ for some $\sigma$
$e_2\text{~}N(0,1)$

Where the values of $y,\tilde{x}$ are known for each sample only, and also the standard deviation of the real value of $x$ for the sample is known: $\sigma_x$.

I get the biased ($\hat{\beta}$) coefficient using OLS, and then making adjustments using:

$\beta' = \hat{\beta} * \frac{\hat{\sigma}_\tilde{x}^2}{\sigma_x^2}$

I see that my new, unbiased estimator for the coefficient is much better (closer to the real value) with this model, but the MSE is getting worse than using the biased estimator.

What is happening? I expected an ubiased estimator to yield better results than the biased one.

Matlab code:

reg_mse_agg = [];
fixed_mse_agg = [];
varMult = 1;
numTests = 60;
for dataNumber=1:8
    reg_mses = [];
    fixed_mses = [];

    X = rand(1000,1);
    X(:,1) = X(:,1) * 10;
    X(:,1) = X(:,1) + 5;

    varX = var(X);
    y = 0.5 * X(:,1) -10;

    y = y + normrnd(0,1,size(y));    
    origX = X;
    X = X + normrnd(0,dataNumber * varMult ,size(X));
    train_size = floor(0.5 * length(y));
    for t=1:numTests,
        idx = randperm(length(y));
        train_idx = idx(1:train_size);
        test_idx = idx(train_size+1:end);
        Xtrain = X(train_idx,:);
        ytrain = y(train_idx);        
        Xtest = X(test_idx,:);
        ytest = y(test_idx);

        b = OLS_solver(Xtrain, ytrain);
        %first arg of evaluate returns MSE, working correctly.
        [ reg_mse, ~ ] = evaluate( b,Xtest,ytest);
        reg_mses = [reg_mses ; reg_mse];

        varInd = var(Xtrain);
        varNoise = varInd - varX;
        bFixed = [0 0]';
        bFixed(1) = b(1) * varInd / varX;
        bFixed(2) = mean(ytrain - bFixed(1)*Xtrain);
        [fixed_mse,~ ] = evaluate( bFixed,Xtest,ytest);
        fixed_mses = [fixed_mses ; fixed_mse];

        dataNumber * varMult
        b
        bFixed

    end
    reg_mse_agg = [reg_mse_agg , reg_mses];
    fixed_mse_agg = [fixed_mse_agg , fixed_mses];

end

mean(reg_mse_agg)
mean(fixed_mse_agg)

Results:

biased estimator's MSE:

ans =

  Columns 1 through 7

    1.2171    1.6513    1.9989    2.3914    2.5766    2.6712    2.5997

  Column 8

    2.8346

Unbiased estimator's MSE:

ans =

  Columns 1 through 7

    1.2308    2.0001    2.9555    4.9727    7.6757   11.3106   14.4283

  Column 8

   11.5653

In addition, printing the values of b and bFixed - I see that bFixed is indeed closer to the real values of 0.5,-10 than the biased estimator (as expected).

P.S. The results of the unbiased being worse than the biased estimator are statistical significant - the test for it is omitted from the code, since it is a simplification of the "full version" code.

UPDTAE: I added a test that checks $\sum_{\text{for each test}}{(\hat{\beta}-\beta)^2}$ and $\sum_{\text{for each test}}{(\beta'-\beta)^2}$, and the biased estimator is indeed significantly worse (larger value) than the unbiased one according to this metric, even though the MSE of the biased estimator (on test-set) is significantly better.
Where $\beta=0.5$ is the real coefficient of the dependent variable, $\hat{\beta}$ is the biased estimator for $\beta$, and $\beta'$ is the unbiased estimator for $\beta$.

This I believe shows that the reason for the results is NOT the higher variance of the unbiased estimator, as it is still closer to the real value.

Credit: Using Lecture notes of Steve Pischke as resource

Answer 1

Summary: the corrected parameters are for predicting as a function of the true predictor $x$. If $\tilde{x}$ is used in prediction, the original parameters perform better.

Note that there are two different linear prediction models lurking around. First, $y$ given $x$, \begin{equation} \hat{y}_x = \beta\,x + \alpha, \end{equation} second, $y$ given $\tilde{x}$, \begin{equation} \hat{y}_{\tilde{x}} = \tilde{\beta}\,\tilde{x} + \tilde{\alpha}. \end{equation}

Even if we had access to the true parameters, the optimal linear prediction as a functions of $x$ would be different than the optimal linear prediction as a function of $\tilde{x}$. The code in the question does the following

1. Estimate the parameters $\hat{\tilde{\beta}},\hat{\tilde{\alpha}}$
2. Compute estimates $\hat{\beta},\hat{\alpha}$
3. Compare performance of $\hat{y}_1 = \hat{\beta}\,\tilde{x} + \hat{\alpha}$ and $\hat{y}_2 = \hat{\tilde{\beta}}\,\tilde{x} + \hat{\tilde{\alpha}}$

Since in step 3 we are predicting as a function of $\tilde{x}$, not as a function of $x$, using (estimated) coefficients of the second model works better.

Indeed, if we had access to $\alpha$, $\beta$ and $\tilde{x}$ but not $x$, we might substitute a linear estimator of $x$ in the first model, \begin{equation} \hat{\hat{y_x}} = \beta\,\hat{x}(\tilde{x}) + \alpha = \beta\, (\mu_x + (\hat{x}-\mu_x)\,\frac{\sigma^2_{x}}{\sigma^2_{\tilde{x}}}) + \alpha = \frac{\sigma_x^2}{\sigma^2_{\tilde{x}}}\beta + \alpha - \beta(1-\frac{\sigma_x^2}{\sigma^2_{\tilde{x}}})\mu_x. \end{equation} If we first perform the transformation form $\tilde{\alpha},\tilde{\beta}$ to $\alpha,\beta$ and then do the computation in the latest equation, we get back the coefficients $\tilde{\alpha},\tilde{\beta}$. So, if the goal is to do linear prediction given the noisy version of the predictor, we should just fit a linear model to the noisy data. The corrected coefficients $\alpha,\beta$ are of interest if we are interested in the true phenomenon for other reasons than prediction.

Testing

I edited the code in OP to also evaluate MSEs for predictions using the non-noisy version of the prediction (code in the end of the answer). The results are

Reg parameters, noisy predictor
1.3387    1.6696    2.1265    2.4806    2.5679    2.5062    2.5160    2.8684

Fixed parameters, noisy predictor
1.3981    2.0626    3.2971    5.0220    7.6490   10.2568   14.1139   20.7604

Reg parameters, true predictor
1.3354    1.6657    2.1329    2.4885    2.5688    2.5198    2.5085    2.8676

Fixed parameters, true predictor
1.1139    1.0078    1.0499    1.0212    1.0492    0.9925    1.0217    1.2528

That is, when using $x$ instead of $\tilde{x}$, the corrected parameters indeed beat the uncorrected parameters, as expected. Furthermore, the prediction with ($\alpha,\beta,x$), that is, fixed parameters and true predictor, is better than ($\tilde{\alpha},\tilde{\beta},\tilde{x}$), that is, reg parameters and noisy predictor, since obviously the noise harms the prediction accuract somewhat. The other two cases correspond to using the parameters of a wrong model and thus produce weaker results.

Caveat about nonlinearity

Actually, even if the relationship between $y,x$ is linear, the relationship between $y$ and $\tilde{x}$ might not be. This depends on the distribution of $x$. For example, in the present code, $x$ is drawn from the uniform distribution, thus no matter how high $\tilde{x}$ is, we know an upper bound for $x$ and thus the predicted $y$ as a function of $\tilde{x}$ should saturate. A possible Bayesian-style solution would be to posit a probability distribution for $x$ and then plug in $\mathbb{E}(x \mid \tilde{x})$ when deriving $\hat{\hat{y}}_x$ - instead of the linear prediction I used previously. However, if one is willing to posit a probability distribution for $x$, I suppose one should go for a full Bayesian solution instead of an approach based on correcting OLS estimates in the first place.

MATLAB code for replicating the test result

Note that this also contains my own implementations for evaluate and OLS_solver since they were not given in the question.

rng(1)

OLS_solver = @(X,Y) [X ones(size(X))]'*[X ones(size(X))] \ ([X ones(size(X))]' * Y);
evaluate = @(b,x,y)  mean(([x ones(size(x))]*b - y).^2);

reg_mse_agg = [];
fixed_mse_agg = [];
reg_mse_orig_agg = [];
fixed_mse_orig_agg = [];

varMult = 1;
numTests = 60;
for dataNumber=1:8
    reg_mses = [];
    fixed_mses = [];
    reg_mses_orig = [];
    fixed_mses_orig = [];

    X = rand(1000,1);
    X(:,1) = X(:,1) * 10;
    X(:,1) = X(:,1) + 5;

    varX = var(X);
    y = 0.5 * X(:,1) -10;

    y = y + normrnd(0,1,size(y));    
    origX = X;
    X = X + normrnd(0,dataNumber * varMult ,size(X));
    train_size = floor(0.5 * length(y));
    for t=1:numTests,
        idx = randperm(length(y));
        train_idx = idx(1:train_size);
        test_idx = idx(train_size+1:end);
        Xtrain = X(train_idx,:);
        ytrain = y(train_idx);        
        Xtest = X(test_idx,:);
        origXtest = origX(test_idx,:);
        ytest = y(test_idx);

        b = OLS_solver(Xtrain, ytrain);
        %first arg of evaluate returns MSE, working correctly.
        reg_mse = evaluate( b,Xtest,ytest);
        reg_mses = [reg_mses ; reg_mse];

        varInd = var(Xtrain);
        varNoise = varInd - varX;
        bFixed = [0 0]';
        bFixed(1) = b(1) * varInd / varX;
        bFixed(2) = mean(ytrain - bFixed(1)*Xtrain);
        fixed_mse = evaluate( bFixed,Xtest,ytest);
        fixed_mses = [fixed_mses ; fixed_mse];

        reg_mse_orig = evaluate(b, origXtest, ytest);
        reg_mses_orig = [reg_mses; reg_mses_orig];

        fixed_mse_orig = evaluate(bFixed, origXtest, ytest);
        fixed_mses_orig = [fixed_mses_orig; fixed_mse_orig];

    end
    reg_mse_agg = [reg_mse_agg , reg_mses];
    fixed_mse_agg = [fixed_mse_agg , fixed_mses];

    reg_mse_orig_agg = [reg_mse_orig_agg , reg_mses_orig];
    fixed_mse_orig_agg = [fixed_mse_orig_agg , fixed_mses_orig]; 
end

disp('Reg parameters, noisy predictor')
disp(mean(reg_mse_agg))
disp('Fixed parameters, noisy predictor')
disp(mean(fixed_mse_agg))
disp('Reg parameters, true predictor')
disp(mean(reg_mse_orig_agg))
disp('Fixed parameters, true predictor')
disp(mean(fixed_mse_orig_agg))