Using regression weights when $Y$ might be measured with nonzero-mean measurement error

Question

Suppose we observe data $Y, X$ and would like to fit a regression model for $\mathbf{E}[Y \,|\, X]$. Unfortunately, $Y$ is sometimes measured with errors whose mean is nonzero.

Let $Z \in \left\{\text{unbiased}, \text{biased}\right\}$ indicate whether $Y$ is measured with classical zero mean errors or nonzero-mean errors, respectively. We would like to estimate $\mathbf{E}[Y \,|\, X, Z = \text{unbiased}]$. Unfortunately, $Z$ is generally not observed, and $\mathbf{E}[Y \,|\, X, Z = \text{unbiased}] \neq \mathbf{E}[Y \,|\, X]$. If we fit a regression of $Y$ on $X$, we'll get biased predictions.

Suppose we cannot generally observe $Z$, but have access to a model for $\Pr[Z \,|\, X,Y]$ (because we manually learned $Z$ on a small training set and fit a classification model with $Z$ as the target variable). Does fitting a regression of $Y$ on $X$ using $\Pr[Z = \text{unbiased} \,|\, X,Y]$ as regression weights produce an unbiased estimate of $\mathbf{E}[Y \,|\, X, Z = \text{unbiased}]$ (or, failing that, a less biased estimate than we would get without using weights)? Is this method used in practice, and does it have a name?

Clarification: the goal is to fit a model that minimizes mean squared error on unseen data (test data) where $Z = \text{unbiased}$. The optimal predictor for that objective is $\mathbf{E}[Y \,|\, X, Z = \text{unbiased}]$, so that is the function we are trying to estimate. Methods for solving this problem should be ranked in terms of how well they achieve that objective.

---

Small example in R with df$y_is_unbiased playing the role of $Z$ and df$y_observed playing the role of $Y$:

library(ggplot2)
library(randomForest)

set.seed(12345)

get_df <- function(n_obs, constant, beta, sd_epsilon, mismeasurement) {
    df <- data.frame(x1=rnorm(n_obs), x2=rnorm(n_obs), epsilon=rnorm(n_obs, sd=sd_epsilon))

    ## Value of Y if measured correctly
    df$y_unbiased <- constant + as.matrix(df[c("x1", "x2")]) %*% beta + df$epsilon

    ## Value of Y if measured incorrectly
    df$y_biased <- df$y_unbiased + sample(mismeasurement, size=n_obs, replace=TRUE)

    ## Y is equally likely to be measured correctly or incorrectly
    df$y_is_unbiased<- sample(c(TRUE, FALSE), size=n_obs, replace=TRUE)
    df$y_observed <- ifelse(df$y_is_unbiased, df$y_unbiased, df$y_biased)

    return(df)
}

## True coefficients
constant <- 5
beta <- c(1, 5)

df <- get_df(n_obs=2000, constant=constant, beta=beta, sd_epsilon=1.0, mismeasurement=c(-10.0, 5.0))

ggplot(df, aes(x=x1, y=y_observed, color=y_is_unbiased)) + geom_point() + scale_color_manual(values=c("#ff7f00", "#377eb8"))

## For facet_wrap title
df$string_y_is_unbiased <- paste0("y_is_unbiased: ", df$y_is_unbiased)

## Notice that Pr[Y | Z = biased] differs from Pr[Y | Z = unbiased]
ggplot(df, aes(x=y_observed)) + geom_histogram(color="black", fill="grey", binwidth=0.5) + facet_wrap(~ string_y_is_unbiased, ncol=1)

## Recover true constant and beta (plus noise) when using y_unbiased
summary(lm(y_unbiased ~ x1 + x2, data=df))

## Biased estimates when using y_biased (constant is biased downward)
summary(lm(y_biased ~ x1 + x2, data=df))

## Also get biased estimates when using y_observed (constant is biased downward)
summary(lm(y_observed ~ x1 + x2, data=df))

## Now image that we "rate" subset of the data (manually check/research whether y was measured with or without bias)
n_rated <- 1000
df_rated <- df[1:n_rated, ]

## Use a factor so that randomForest does classification instead of regression
df_rated$y_is_unbiased <- factor(df_rated$y_is_unbiased)

model_pr_unbiased <- randomForest(formula=y_is_unbiased ~ y_observed + x1 + x2, data=df_rated, mtry=2)

## Examine OOB confusion matrix (error rate < 5%)
print(model_pr_unbiased)

## Use the model to get Pr[Y is unbiased | X, observed Y] on unrated data
df_unrated <- df[(n_rated+1):nrow(df), ]
df_unrated$pr_unbiased <- as.vector(predict(model_pr_unbiased, newdata=df_unrated, type="prob")[, "TRUE"])

## Train a model on unrated data, using pr_unbiased as regression weights -- is this unbiased?
summary(lm(y_observed ~ x1 + x2, data=df_unrated, weights=df_unrated$pr_unbiased))

In this example, the model $\Pr[Z = \text{unbiased} \,|\, X,Y]$ is a random forest with formula=y_is_unbiased ~ y_observed + x1 + x2. If this model were perfectly accurate, it would generate weights of 1.0 where $Y$ is unbiased, 0.0 where $Y$ is biased, and the weighted regression would clearly be unbiased. What happens when the model for $\Pr[Z = \text{unbiased} \,|\, X,Y]$ has test precision and recalls that aren't perfect (<100% accuracy)? Is the weighted regression guaranteed to be less biased than an unweighted regression of $Y$ on $X$?

---

Slightly more complex example in which $\Pr[Z = \text{unbiased} \,|\, X]$ varies with $X$ (as opposed to the simpler example I posted above, where $\Pr[Z = \text{unbiased} \,|\, X] = \frac{1}{2} \; \forall X$):

library(ggplot2)
library(randomForest)

set.seed(12345)

logistic <- function(x) {
    return(1 / (1 + exp(-x)))
}

pr_y_is_unbiased <- function(x1, x2) {
    ## This function returns Pr[ Z = unbiased | X ]
    return(logistic(x1 + 2*x2))
}

get_df <- function(n_obs, constant, beta, sd_epsilon, mismeasurement) {
    df <- data.frame(x1=rnorm(n_obs), x2=rnorm(n_obs), epsilon=rnorm(n_obs, sd=sd_epsilon))

    ## Value of Y if measured correctly
    df$y_unbiased <- constant + as.matrix(df[c("x1", "x2")]) %*% beta + df$epsilon

    ## Value of Y if measured incorrectly
    df$y_biased <- df$y_unbiased + sample(mismeasurement, size=n_obs, replace=TRUE)

    ## Note: in this example, Pr[ Z = biased | X ] varies with X
    ## In the first (simpler) example I posted, Pr[ Z = biased | X ] = 1/2 was constant with respect to X
    df$y_is_unbiased <- runif(n_obs) < pr_y_is_unbiased(df$x1, df$x2)

    df$y_observed <- ifelse(df$y_is_unbiased, df$y_unbiased, df$y_biased)

    return(df)
}

## True coefficients
constant <- 5
beta <- c(1, 5)

df <- get_df(n_obs=2000, constant=constant, beta=beta, sd_epsilon=1.0, mismeasurement=c(-10.0, 5.0))

ggplot(df, aes(x=x1, y=y_observed, color=y_is_unbiased)) + geom_point() + scale_color_manual(values=c("#ff7f00", "#377eb8"))

## For facet_wrap title
df$string_y_is_unbiased <- paste0("y_is_unbiased: ", df$y_is_unbiased)

## Notice that Pr[Y | Z = biased] differs from Pr[Y | Z = unbiased]
ggplot(df, aes(x=y_observed)) + geom_histogram(color="black", fill="grey", binwidth=0.5) + facet_wrap(~ string_y_is_unbiased, ncol=1)

## Recover true constant and beta (plus noise) when using y_unbiased
summary(lm(y_unbiased ~ x1 + x2, data=df))

## Biased estimates when using y_biased (constant is biased downward)
summary(lm(y_biased ~ x1 + x2, data=df))

## Also get biased estimates when using y_observed
## Note: the constant is biased downward _and_ the coefficient on x2 is biased upward!
summary(lm(y_observed ~ x1 + x2, data=df))

## Now image that we "rate" subset of the data (manually check/research whether y was measured with or without bias)
n_rated <- 1000
df_rated <- df[1:n_rated, ]

## Use a factor so that randomForest does classification instead of regression
df_rated$y_is_unbiased <- factor(df_rated$y_is_unbiased)

model_pr_unbiased <- randomForest(formula=y_is_unbiased ~ y_observed + x1 + x2, data=df_rated, mtry=2)

## Examine OOB confusion matrix (error rate < 5%)
print(model_pr_unbiased)

## Use the model to get Pr[Y is unbiased | X, observed Y] on unrated data
df_unrated <- df[(n_rated+1):nrow(df), ]
df_unrated$pr_unbiased <- as.vector(predict(model_pr_unbiased, newdata=df_unrated, type="prob")[, "TRUE"])

## Train a model on unrated data, using pr_unbiased as regression weights -- is this unbiased? If not, is it _less_ biased than the unweighted model?
summary(lm(y_observed ~ x1 + x2, data=df_unrated, weights=df_unrated$pr_unbiased))

## What happens if we use pr_unbiased as a feature (aka predictor) in the regression, rather than a weight?
## In this case the weighted regression seems to do better, but neither is perfect
## Note: copied from shabbychef's answer
summary(lm(formula = y_observed ~ x1 + x2 + I(1 - pr_unbiased), data = df_unrated))

In this example, the weighted regression of $Y$ on $X$ looks less biased than the unweighted regression. Is that true in general? I also tried shabbychef's suggestion (see answer below) on this example, and it appears to do worse than the weighted regression.

---

For those who prefer Python to R, here's the second simulation in Python:

import numpy as np
import pandas as pd
from sklearn.ensemble import RandomForestClassifier
from sklearn.linear_model import LinearRegression

def logistic(x):
    return 1 / (1 + np.exp(-x))

def pr_y_is_unbiased(x1, x2):
    # This function returns Pr[ Z = unbiased | X ]
    return logistic(x1 + 2*x2)

def get_df(n_obs, constant, beta, sd_epsilon, mismeasurement):
    df = pd.DataFrame({
        'x1': np.random.normal(size=n_obs),
        'x2': np.random.normal(size=n_obs),
        'epsilon': np.random.normal(size=n_obs, scale=sd_epsilon),
    })

    df['y_unbiased'] = constant + np.dot(np.array(df[['x1', 'x2']]), beta) + df['epsilon']

    # Note: df['y_biased'].mean() will differ from df['y_unbiased'].mean() if the mismeasurements have a nonzero mean
    df['y_biased'] = df['y_unbiased'] + np.random.choice(mismeasurement, size=n_obs)

    df['y_is_unbiased'] =  np.random.uniform(size=n_obs) < pr_y_is_unbiased(df['x1'], df['x2'])

    df['y_observed'] = df.apply(lambda row: row['y_unbiased'] if row['y_is_unbiased'] else row['y_biased'], axis=1)

    return df


constant = 5
beta = np.array([1, 5])
print(f'true coefficients:\n constant = {constant}, beta = {beta}')

n_obs = 2000

# Note: the mean of the possible mismeasurements is nonzero (this is the source of the bias)
df = get_df(n_obs=n_obs, constant=constant, beta=beta, sd_epsilon=1.0, mismeasurement=[-10.0, 5.0])

lr = LinearRegression()
lr.fit(X=df[['x1', 'x2']], y=df['y_observed'])

print(f'estimates from unweighted regression of Y on X ({df.shape[0]} obs):\n constant = {lr.intercept_}, beta = {lr.coef_}')

# Note: pretend that we only observe y_is_unbiased on a "rated" subset of the data
n_rated = n_obs // 2
df_rated = df.iloc[:n_rated].copy()
df_unrated = df.iloc[n_rated:].copy()

rf = RandomForestClassifier(n_estimators=500, max_features=2, oob_score=True)
rf_predictors = ['y_observed', 'x1', 'x2']

rf.fit(X=df_rated[rf_predictors], y=df_rated['y_is_unbiased'])

print(f'random forest classifier OOB accuracy (for predicting whether Y is unbiased): {rf.oob_score_}')

df_unrated['pr_y_is_unbiased'] = rf.predict_proba(df_unrated[rf_predictors])[:, 1]

lr.fit(X=df_unrated[['x1', 'x2']], y=df_unrated['y_observed'], sample_weight=df_unrated['pr_y_is_unbiased'])
print(f'estimates from weighted regression of Y on X ({df_unrated.shape[0]} obs):\n constant = {lr.intercept_}, beta = {lr.coef_}')

Answer 1

I would use the 'probability of biased' as a dummy variable in the regression; it can possibly 'sop up' the bias present in the biased case. Using your example, (but calling set.seed(1234) before the call to get_df), I tried

summary(lm(y_observed ~ x1 + x2 + I(1-pr_unbiased), data=df_unrated))

and got:

Call:
lm(formula = y_observed ~ x1 + x2 + I(1 - pr_unbiased), data = df_unrated)

Residuals:
   Min     1Q Median     3Q    Max 
-9.771 -2.722 -0.386  2.474 11.238 

Coefficients:
                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)           5.515      0.250   22.07   <2e-16 ***
x1                    1.108      0.169    6.54    1e-10 ***
x2                    4.917      0.168   29.26   <2e-16 ***
I(1 - pr_unbiased)   -3.727      0.383   -9.72   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.25 on 996 degrees of freedom
Multiple R-squared:  0.514,     Adjusted R-squared:  0.513 
F-statistic:  351 on 3 and 996 DF,  p-value: <2e-16

The coefficient for the term 1-pr_unbiased should be the size of the bias.

Answer 2

This is an omitted-variable problem where you have an indicator variable $Z$ that is unobserved, but which has a relationship with the response variable. Since "bias" is a property of an estimator, not a regression variable, I am going to reframe your question as one where you want to find the true regression function conditional on $Z=0$ using regression data that does not include this variable, and a separate set of regression training data that is used to estimate the probabilities $p_0(x,y) \equiv \mathbb{P}(Z=0|X=x,Y=y)$.

Let $p_{Y|X}$ denote the conditional density of the response variable in the regression problem with response variable $Y$ and explanatory variable $X$ (but excluding $Z$). From the rules of conditional probability, the target distribution of interest can be written as:

$$\begin{equation} \begin{aligned} p(Y=y|X=x,Z=0) &= \frac{p(Y=y,Z=0|X=x)}{p(Z=0|X=x)} \\[6pt] &= \frac{p_0(x,y) \cdot p_{Y|X}(y|x)}{\int_\mathbb{R} p_0(x,y) \cdot p_{Y|X}(y|x) \ dy} \\[6pt] &\overset{y}{\propto} p_0(x,y) \cdot p_{Y|X}(y|x). \\[6pt] \end{aligned} \end{equation}$$

Thus, we can see that it is sufficient to be able to estimate the regression function $p_{Y|X}$ in the regression model with $Z$ omitted, and also estimate the probability function $p_0$ which you have as a separate estimator from your training data. The former can be estimated using OLS estimation without imposing any weights. The "weighting" occurs after estimation of this function, by substitution into the above equation.

We can see that it is not necessary (or desirable) to use any weights in the regression of $Y$ on $X$, since it is sufficient to estimate the conditional density $p_{Y|X}$ without consideration of $Z$. OLS estimation of the coefficients of this regression gives an estimator $\hat{p}_{Y|X}$, and since you also have a separate estimator $\hat{p}_0$ you then have:

$$\hat{p}(Y=y|X=x,Z=0) \propto \hat{p}_0(x,y) \cdot \hat{p}_{Y|X}(y|x). $$

Once you have substituted these estimators, all that remains is to try to determine the scaling constant that yields a proper density function. This can be done by a range of numerical integration methods (e.g., Simpson's rule, quadrature, Metropolis-Hastings, etc.).

Answer 3

Your idea will not give an unbiased estimate, unless you can always be 100% sure whether it is biased or not. As soon as one biased example can be part of your training set with nonzero probability, there will be bias, as you have nothing to cancel out that bias. In practice, your bias will simply be multiplied by a factor $\alpha<1$, where $\alpha$ is the probability that a biased example is detected as such.

Assuming you have enough data, a better approach is to compute $P(Z=biased|X,Y)$ for each sample, and then remove all samples from the training set where this probability exceeds a certain threshold. For example, if it is feasible for you to train your dataset on only the samples for which $P(Z=biased|X,Y)<0.01$, and your dataset decreases from $N$ biased and $M$ unbiased to $n$ biased and $m$ unbiased examples, and the bias will multiply by a factor $f=\frac{n(N+M)}{N(n+m)}$. Since typically $\frac nN$ will be far lower than $\frac mM$, $f$ will be much smaller than $1$, resulting in a significant improvement.

Note that both techniques can be combined: rows with $p=P(Z=biased|X,Y)>\beta$ go out (for some choice of $\beta$, above I used $\beta=0.01$), and the rows that stay in get a weight $(1-\frac{p}{\beta})^2$, that should actually give you the best of both worlds.