Verwenden von Regressionsgewichten, wenn

8

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

Es sei Z{unbiased,biased}} 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ätzenE[Y.|X.,Z.=unvoreingenommen]] . Leider wirdZ. im Allgemeinen nicht beobachtet undE[Y|X,Z=unbiased]E[Y|X] . Wenn wir eine Regression vonY aufX 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 wirZ manuellin einem kleinen Trainingssatzgelerntund ein Klassifizierungsmodell mitZ als Zielvariable angepasst haben). Passt eine Regression vonY aufX mitPr[Z=unbiased|X,Y] als Regressionsgewichte ergeben eine unvoreingenommene Schätzung vonE[Y|X,Z=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=unbiased . Der optimale Prädiktor für dieses Ziel ist E[Y|X,Z=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=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, wennY unverzerrt ist, 0,0, wennY voreingenommen ist, und die gewichtete Regression wäre eindeutig unverzerrt. Was passiert, wenn das Modell für Pr [ Z = unverzerrt ist ?Pr[Z=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 vonY aufX ?


Etwas komplexeres Beispiel, in dem Pr[Z=unbiased|X] variiert mitX (im Gegensatz zu dem einfacheren Beispiel, das ich oben gepostet habe, wobeiPr[Z=unbiased|X]=12X ):

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_}')
Adrian
quelle
1
Dies klingt fast wie "Instrumentalvariablen", bei denen Sie einige Variablen beobachten, die mit dem Fehler in Ihrer Regression korrelieren. Ich fürchte, das hilft nicht viel.
Shabbychef
@shabbychef Richtig, aber in diesem Setup sind keine Instrumente verfügbar.
Adrian
1
Bei Ihrem aktualisierten Problem ist die Verzerrung nun eine Funktion von , und wir sollten erwarten, dass sich die Regressionskoeffizienten ändern. Das heißt, der 'Bias'-Term ist - 0,25 p, wobei p = logistisch ist ( x 1 + 2 x 2 ) . Ich würde p mit einer Taylor-Erweiterung erweitern, um zu zeigen, dass es eine zusätzliche lineare Abhängigkeit von y von x 1 und x 2 gibt . Gehen Sie zurück zu Ihrer ursprünglichen Frage, dem Bias-Term, den Sie sehen, und dem zxi0.25pp=logistic(x1+2x2)pyx1x2zVariable, die Sie beobachten, ändern Sie die Varianz in keiner Weise, sondern den erwarteten Wert. Ich denke, sie sollten in die lineare Spezifikation gehackt werden und nicht in die Gewichte.
Shabbychef

Antworten:

2

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.

shabbychef
quelle
Pr[Z=unbiased|X]]X.Y.X.
2

ZZ=0p0(x,y)P(Z=0|X=x,Y=y)

pY|XYXZ

p(Y=y|X=x,Z=0)=p(Y=y,Z=0|X=x)p(Z=0|X=x)=p0(x,y)pY|X(y|x)Rp0(x,y)pY|X(y|x) dyyp0(x,y)pY|X(y|x).

pY|XZp0

YXpY|XZp^Y|Xp^0

p^(Y=y|X=x,Z=0)p^0(x,y)p^Y|X(y|x).

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

Ben - Monica wieder einsetzen
quelle
1
Vielen Dank für diese ausführliche Antwort (+1). Ich werde sie sorgfältig lesen und mich bei Ihnen melden. Ich stimme zu, dass eine korrektere Beschreibung des Problems "Y wird manchmal mit einem Fehler ungleich Null gemessen" und nicht "Y ist voreingenommen" sein könnte.
Adrian
Interessant! Ein Nachteil / heikles Detail hier ist, dass dies die Schätzung der vollständigen Verteilung von erfordertY|X
Ja, aber das sollte nicht wirklich schwierig sein. Welches Regressionsmodell Sie auch verwenden, es gibt eine klare Modellform, die diese bedingte Verteilung bestimmt.
Ben - Reinstate Monica
1
YXY
1
Ich denke nicht, dass Adrian notwendigerweise an der vollständigen Verteilung von (Y | X, Z = 0) interessiert ist, sondern nur an seinem Mittelwert. Wenn man die vollständige Verteilung von Y | X wissen will, Z = 0, dann sind natürlich die genauen Verteilungen von Y und X sehr relevant, aber wenn man nur E [Y | X, Z = 0] schätzen will, dann ist dies nicht notwendig: das Gesetz der großen Zahlen gilt für beliebige Verteilungen.
user1111929
1

α<1α

P(Z=biased|X,Y)P(Z=biased|X,Y)<0.01NMnmf=n(N+M)N(n+m)nNmMf1

p=P(Z=biased|X,Y)>βββ=0.01(1pβ)2

user1111929
quelle
Ich bin damit einverstanden, dass die gewichtete Regression im Allgemeinen keine unvoreingenommenen Schätzungen liefert, es sei denn, der Klassifikator ist 100% genau. Wird garantiert, dass die Voreingenommenheit verringert wird? Ist Ihr Cutoff-Ansatz notwendigerweise besser oder hängt dies von der Stichprobengröße und der Genauigkeit des Klassifikators ab?
Adrian
Nicht wirklich die Genauigkeit des Klassifikators, hauptsächlich die resultierende Anzahl von Zeilen ist wichtig, da es sich bei sehr kleinen Datensätzen manchmal nicht leisten kann, eine große Anzahl von Zeilen zu verwerfen. Bei genügend Daten sehe ich jedoch keinen Grund für Zweifel: Ihr Ansatz wiegt Zeilen mit einer 50% igen Verzerrungswahrscheinlichkeit halb so relevant wie Zeilen mit einer 0% igen Verzerrungswahrscheinlichkeit. Persönlich würde ich das niemals tun. Ich würde 1000 garantierte saubere Zeilen gegenüber 2000 Zeilen mit jeweils 50% Verzerrung vorziehen. Beachten Sie, dass Sie auch beide Ansätze kombinieren können, um ein allmählicheres System als einen einfachen Cutoff anzuwenden. Ich habe meine Antwort aktualisiert, um dies näher zu erläutern.
user1111929