Verwenden von Regressionsgewichten, wenn

Angenommen, wir beobachten die Daten $Y, X$ und möchten ein Regressionsmodell für $\mathbf{E}[Y \,|\, X]$ . Leider wird$Y$ manchmal mit Fehlern gemessen, deren Mittelwert ungleich Null ist.

Es sei $Z \in \left\{\text{unbiased}, \text{biased}\right\}$ anzugeben, ob $Y$ mit klassischen Fehlern mit einem Mittelwert von Null bzw. mit einem Fehler ungleich Null gemessen wird. Wir möchten E [ Y schätzen$\mathbf{E}[Y \,|\, X, Z = \text{unbiased}]$ . Leider wird$Z$ im Allgemeinen nicht beobachtet und$\mathbf{E}[Y \,|\, X, Z = \text{unbiased}] \neq \mathbf{E}[Y \,|\, X]$ . Wenn wir eine Regression von$Y$ auf$X$ anpassen, erhalten wir voreingenommene Vorhersagen.

Angenommen, wir können $Z$ nicht allgemein beobachten , haben aber Zugriff auf ein Modell für $\Pr[Z \,|\, X,Y]$ (weil wir$Z$ manuellin einem kleinen Trainingssatzgelerntund ein Klassifizierungsmodell mit$Z$ als Zielvariable angepasst haben). Passt eine Regression von$Y$ auf$X$ mit$\Pr[Z = \text{unbiased} \,|\, X,Y]$ als Regressionsgewichte ergeben eine unvoreingenommene Schätzung von$\mathbf{E}[Y \,|\, X, Z = \text{unbiased}]$ (oder, falls dies nicht der Fall ist, eine weniger voreingenommene Schätzung, als wir ohne die Verwendung von Gewichten erhalten würden)? Wird diese Methode in der Praxis angewendet und hat sie einen Namen?

Klarstellung: Ziel ist es, ein Modell anzupassen, das den mittleren quadratischen Fehler bei unsichtbaren Daten (Testdaten) minimiert, wobei $Z = \text{unbiased}$ . Der optimale Prädiktor für dieses Ziel ist $\mathbf{E}[Y \,|\, X, Z = \text{unbiased}]$ , das ist also die Funktion, die wir zu schätzen versuchen. Methoden zur Lösung dieses Problems sollten dahingehend eingestuft werden, wie gut sie dieses Ziel erreichen.

Kleines Beispiel in R mit df$y_is_unbiasedder Rolle von $Z$ und df$y_observedder Rolle von $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 diesem Beispiel ist das Modell $\Pr[Z = \text{unbiased} \,|\, X,Y]$ ist ein zufälliger Wald mitformula=y_is_unbiased ~ y_observed + x1 + x2. Wenn dieses Modell vollkommen genau wäre, würde es Gewichte von 1,0 erzeugen, wenn$Y$ unverzerrt ist, 0,0, wenn$Y$ voreingenommen ist, und die gewichtete Regression wäre eindeutig unverzerrt. Was passiert, wenn das Modell für Pr [ Z = unverzerrt ist ?$\Pr[Z = \text{unbiased} \,|\, X,Y]$ hat Testgenauigkeit und Rückrufe, die nicht perfekt sind (<100% Genauigkeit)? Ist die gewichtete Regression garantiert weniger voreingenommen als eine ungewichtete Regression von$Y$ auf$X$ ?

Etwas komplexeres Beispiel, in dem $\Pr[Z = \text{unbiased} \,|\, X]$ variiert mit$X$ (im Gegensatz zu dem einfacheren Beispiel, das ich oben gepostet habe, wobei$\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 diesem Beispiel sieht die gewichtete Regression von $Y$ auf $X$ weniger voreingenommen aus als die ungewichtete Regression. Ist das im Allgemeinen wahr? Ich habe auch Shabbychefs Vorschlag (siehe Antwort unten) für dieses Beispiel ausprobiert, und es scheint schlechter zu sein als die gewichtete Regression.

Für diejenigen, die Python gegenüber R bevorzugen, ist hier die zweite 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_}')

Ich würde die 'Wahrscheinlichkeit einer Verzerrung' als Dummy-Variable in der Regression verwenden. es kann möglicherweise die im vorgespannten Fall vorhandene Vorspannung "aufheben". Anhand Ihres Beispiels (aber set.seed(1234)vor dem Anruf bei get_df) habe ich es versucht

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

und bekam:

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

Der Koeffizient für den Term 1-pr_unbiasedsollte die Größe der Vorspannung sein.

$Z$$Z=0$$p_0(x,y) \equiv \mathbb{P}(Z=0|X=x,Y=y)$

$p_{Y|X}$$Y$$X$$Z$

$p_{Y|X}$$Z$$p_0$

$Y$$X$$p_{Y|X}$$Z$$\hat{p}_{Y|X}$$\hat{p}_0$

Nachdem Sie diese Schätzer ersetzt haben, müssen Sie nur noch versuchen, die Skalierungskonstante zu bestimmen, die eine ordnungsgemäße Dichtefunktion ergibt. Dies kann durch eine Reihe numerischer Integrationsmethoden erfolgen (z. B. Simpsons Regel, Quadratur, Metropolis-Hastings usw.).

$\alpha<1$$\alpha$

$P(Z=biased|X,Y)$$P(Z=biased|X,Y)<0.01$$N$$M$$n$$m$$f=\frac{n(N+M)}{N(n+m)}$$\frac nN$$\frac mM$$f$$1$

$p=P(Z=biased|X,Y)>\beta$$\beta$$\beta=0.01$$(1-\frac{p}{\beta})^2$