8
$\begingroup$

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 a systematic bias (i.e. errors whose mean is nonzero).

Let $Z \in \left\{\text{unbiased}, \text{biased}\right\}$ indicate whether $Y$ is measured with bias or not. We would actually 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_}')
$\endgroup$
  • 1
    $\begingroup$ This almost sounds like "Instrumental Variables", where you observe some variables that are correlated with the error in your regression. I'm afraid that doesn't help too much, though. $\endgroup$ – shabbychef Feb 15 at 5:59
  • $\begingroup$ @shabbychef True, but there are no instruments available in this setup. $\endgroup$ – Adrian Feb 15 at 21:19
  • 1
    $\begingroup$ With your updated problem, the bias is now a function of the $x_i$, and we should expect the regression coefficients to change. That is, the 'bias' term is $-0.25p$ where $p=\operatorname{logistic}(x_1 + 2x_2)$. I would expand $p$ with a Taylor expansion to show that there is extra linear dependence of $y$ on $x_1$ and $x_2$. Going back to your original question, the bias term you see, and the $z$ variable you observe, do not change variance in any way, but do change the expected value. So they should be hacked into the linear specification, I would think, and not into the weights. $\endgroup$ – shabbychef Feb 16 at 1:54
2
$\begingroup$

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.).

$\endgroup$
  • 1
    $\begingroup$ Thank you for this detailed answer (+1), I'll read it carefully and get back to you. I agree that a more correct description of the problem might be "Y is sometimes measured with nonzero-mean error", rather than "Y is biased." $\endgroup$ – Adrian Feb 18 at 5:02
  • $\begingroup$ Interesting! One downside / tricky detail here is that this requires estimating the full distribution of $Y \,|\, X$, rather than just the conditional mean. Normality is a common assumption, but that may not hold in my application. $\endgroup$ – Adrian Feb 20 at 1:39
  • $\begingroup$ Yes, but that should not really be tricky. Whatever regression model you use will have a clear model form that determines this conditional distribution. $\endgroup$ – Ben Feb 20 at 2:54
  • 1
    $\begingroup$ Ben, are you sure that's true in a machine learning / prediction setting? Fitting a linear regression to minimize mean squared error on test data does not require making any assumptions about the normality of the residuals (as long as you aren't doing finite sample hypothesis tests, aren't interested in prediction intervals, don't care whether your estimator is a maximum likelihood estimator, etc). I'm interested in estimating a conditional mean (of $Y$ given $X$) and haven't made assumptions about the distribution of $Y$. $\endgroup$ – Adrian Feb 20 at 3:07
  • 1
    $\begingroup$ I don't think Adrian is necessarily interested in the full distribution of (Y|X,Z=0), merely in its mean. If one wants to know the full distribution of Y|X,Z=0 then of course the exact distributions of Y and X are very relevant, but if one only wants to estimate E[Y|X,Z=0] then this is not necessary: the law of large numbers works for arbitrary distributions. $\endgroup$ – user1111929 Feb 20 at 13:17
1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Interesting idea (+1)! I updated my post with a second, slightly more complex example in which $\Pr[Z = \text{unbiased} \,|\, X]$ varies with $X$ instead of being constant. In this case, both the constant and the coefficient on x2 are biased when regressing $Y$ on $X$, and I don't think your method works as well (since it only removes bias from the constant). Have a look and let me know what you think! $\endgroup$ – Adrian Feb 15 at 21:21
1
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ I agree that the weighted regression generally won't produce unbiased estimates unless the classifier is 100% accurate. Is it guaranteed to reduce bias? Is your cutoff approach necessarily better, or does that depend on sample sizes and the classifier's accuracy? $\endgroup$ – Adrian Feb 17 at 5:28
  • $\begingroup$ Not really the classifier's accuracy, mostly the resulting amount of rows matters, as for very small datasets one can sometimes not afford to discard a large number of rows. But given enough data I don't see a reason for doubt: your approach weighs rows with 50% chance of bias as half as relevant as rows with 0% chance of bias. Personally, I'd never do that, I'd highly prefer 1000 guaranteed clean rows over 2000 rows that each have 50% chance of bias. Note that you can also combine both approaches, to apply a more gradual system than a plain cutoff, I've updated my answer to elaborate on this. $\endgroup$ – user1111929 Feb 18 at 16:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.