Utilisation de poids de régression lorsque

8

Supposons que nous observions les données et que nous souhaitions ajuster un modèle de régression pour . Malheureusement, est parfois mesuré avec des erreurs dont la moyenne est non nulle.Y,XE[Y|X]Y

Soit indiquer si est mesuré avec des erreurs moyennes nulles classiques ou des erreurs non nulles, respectivement. Nous aimerions estimer . Malheureusement, n'est généralement pas observé, et . Si nous ajustons une régression de sur , nous obtiendrons des prédictions biaisées.Z{unbiased,biased}OuiE[Oui|X,Z=impartial]ZE[Oui|X,Z=impartial]E[Oui|X]OuiX

Supposons que nous ne pouvons généralement pas observer , mais que nous ayons accès à un modèle pour (parce que nous avons appris manuellement Z sur un petit ensemble d'apprentissage et ajusté un modèle de classification avec Z comme variable cible) . Le fait d'ajuster une régression de Y sur X en utilisant \ Pr [Z = \ text {impartial} \, | \, X, Y] comme poids de régression produit une estimation non biaisée de \ mathbf {E} [Y \, | \, X, Z = \ text {impartiale}] (ou, à défaut, une estimation moins biaisée que celle que nous obtiendrions sans utiliser de pondérations)? Cette méthode est-elle utilisée dans la pratique et a-t-elle un nom?ZPr[Z|X,Oui]ZZOuiXPr[Z=impartial|X,Oui]E[Oui|X,Z=impartial]

Clarification: l'objectif est d'adapter un modèle qui minimise l'erreur quadratique moyenne sur des données invisibles (données de test) où Z=impartial . Le prédicteur optimal pour cet objectif est E[Oui|X,Z=impartial] , c'est donc la fonction que nous essayons d'estimer. Les méthodes permettant de résoudre ce problème doivent être classées en fonction de leur efficacité à atteindre cet objectif.


Petit exemple en R avec df$y_is_unbiasedle rôle de Z et df$y_observedle rôle de Oui :

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

Dans cet exemple, le modèle est une forêt aléatoire avec . Si ce modèle était parfaitement précis, il générerait des pondérations de 1,0 où est sans biais, 0,0 où est biaisé, et la régression pondérée serait clairement sans biais. Que se passe-t-il lorsque le modèle de a une précision de test et des rappels qui ne sont pas parfaits (<100% de précision)? La régression pondérée est-elle garantie d'être moins biaisée qu'une régression non pondérée de sur ?Pr[Z=impartial|X,Oui]formula=y_is_unbiased ~ y_observed + x1 + x2OuiOuiPr[Z=impartial|X,Oui]OuiX


Exemple un peu plus complexe dans lequel varie avec (contrairement à l'exemple plus simple que j'ai posté ci-dessus, où ):Pr[Z=impartial|X]XPr[Z=impartial|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))

Dans cet exemple, la régression pondérée de sur semble moins biaisée que la régression non pondérée. Est-ce vrai en général? J'ai également essayé la suggestion de shabbychef (voir la réponse ci-dessous) sur cet exemple, et elle semble faire pire que la régression pondérée.OuiX


Pour ceux qui préfèrent Python à R, voici la deuxième simulation en 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
la source
1
Cela ressemble presque à des "variables instrumentales", où vous observez certaines variables qui sont corrélées à l'erreur de votre régression. J'ai bien peur que cela n'aide pas trop.
shabbychef
@shabbychef True, mais aucun instrument n'est disponible dans cette configuration.
Adrian
1
Avec votre problème mis à jour, le biais est désormais fonction du , et nous devrions nous attendre à ce que les coefficients de régression changent. Autrement dit, le terme de «biais» est où . Je développerais avec une expansion de Taylor pour montrer qu'il y a une dépendance linéaire supplémentaire de sur et . Pour revenir à votre question initiale, le terme de biais que vous voyez et la variable vous observez ne changent en rien la variance, mais changent la valeur attendue. Ils devraient donc être piratés dans la spécification linéaire, je pense, et non dans les poids. Xje-0,25pp=la logistique(X1+2X2)pyX1X2z
shabbychef

Réponses:

2

J'utiliserais la «probabilité de biais» comme variable muette dans la régression; cela peut éventuellement «éliminer» le biais présent dans le cas biaisé. En utilisant votre exemple, (mais en appelant set.seed(1234)avant l'appel à get_df), j'ai essayé

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

et a obtenu:

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

Le coefficient du terme 1-pr_unbiaseddoit être la taille du biais.

shabbychef
la source
Idée intéressante (+1)! J'ai mis à jour mon article avec un deuxième exemple un peu plus complexe dans lequelPr[Z=impartial|X] varie avec Xau lieu d'être constant. Dans ce cas, la constante et le coefficient sur x2 sont biaisés lors de la régressionOui sur X, et je ne pense pas que votre méthode fonctionne aussi bien (car elle supprime uniquement le biais de la constante). Jetez un œil et faites-moi savoir ce que vous en pensez!
Adrian
2

Il s'agit d'un problème de variable omise où vous avez une variable indicatriceZqui n'est pas observé, mais qui a une relation avec la variable de réponse. Étant donné que le «biais» est une propriété d'un estimateur, et non une variable de régression, je vais reformuler votre question comme celle où vous voulez trouver la vraie fonction de régression conditionnelle àZ=0 en utilisant des données de régression qui n'incluent pas cette variable, et un ensemble distinct de données d'entraînement de régression qui est utilisé pour estimer les probabilités p0(x,y)P(Z=0|X=x,Y=y).

Laisser pY|X dénoter la densité conditionnelle de la variable de réponse dans le problème de régression avec la variable de réponse Y et variable explicative X (mais à l'exclusion Z). A partir des règles de probabilité conditionnelle, la distribution cible des intérêts peut s'écrire:

p(Oui=y|X=X,Z=0)=p(Oui=y,Z=0|X=X)p(Z=0|X=X)=p0(X,y)pOui|X(y|X)Rp0(X,y)pOui|X(y|X) yyp0(X,y)pOui|X(y|X).

On voit donc qu'il suffit de pouvoir estimer la fonction de régression pOui|X dans le modèle de régression avec Z omis, et également estimer la fonction de probabilité p0que vous avez comme estimateur distinct de vos données d'entraînement. Les premiers peuvent être estimés en utilisant l'estimation OLS sans imposer de pondération. La "pondération" intervient après estimation de cette fonction, par substitution dans l'équation ci-dessus.

Nous pouvons voir qu’il n’est pas nécessaire (ou souhaitable) d’utiliser des poids dans la régression de Oui sur X, car il suffit d'estimer la densité conditionnelle pOui|X sans considération de Z. L'estimation OLS des coefficients de cette régression donne un estimateurp^Oui|X, et puisque vous avez également un estimateur distinct p^0 vous avez alors:

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

Une fois que vous avez substitué ces estimateurs, il ne reste plus qu'à essayer de déterminer la constante de mise à l'échelle qui donne une fonction de densité appropriée. Cela peut être fait par une gamme de méthodes d'intégration numérique (par exemple, la règle de Simpson, la quadrature, Metropolis-Hastings, etc.).

Ben - Réintègre Monica
la source
1
Merci pour cette réponse détaillée (+1), je vais la lire attentivement et vous recontacter. Je conviens qu'une description plus correcte du problème pourrait être «Y est parfois mesuré avec une erreur non nulle», plutôt que «Y est biaisé».
Adrian
Intéressant! Un détail négatif / délicat ici est que cela nécessite d’estimer la distribution complète desOui|X, plutôt que simplement la moyenne conditionnelle. La normalité est une hypothèse courante, mais cela peut ne pas être valable dans ma candidature.
Adrian
Oui, mais cela ne devrait pas vraiment être délicat. Quel que soit le modèle de régression que vous utilisez, le formulaire de modèle clair détermine cette distribution conditionnelle.
Ben - Rétablir Monica le
1
Ben, êtes-vous sûr que c'est vrai dans un contexte d'apprentissage automatique / de prédiction? Ajuster une régression linéaire pour minimiser l'erreur quadratique moyenne sur les données de test ne nécessite pas de faire d'hypothèses sur la normalité des résidus (tant que vous ne faites pas de tests d'hypothèse d'échantillons finis, ne vous intéressez pas aux intervalles de prédiction, ne vous souciez pas si votre estimateur est un estimateur du maximum de vraisemblance, etc.). Je souhaite estimer une moyenne conditionnelle (deOui donné X) et n'ont pas émis d'hypothèses sur la répartition des Oui.
Adrian
1
Je ne pense pas qu'Adrian soit nécessairement intéressé par la distribution complète de (Y | X, Z = 0), simplement par sa moyenne. Si l'on veut connaître la distribution complète de Y | X, Z = 0 alors bien sûr les distributions exactes de Y et X sont très pertinentes, mais si l'on veut seulement estimer E [Y | X, Z = 0] alors c'est pas nécessaire: la loi des grands nombres fonctionne pour les distributions arbitraires.
user1111929
1

Votre idée ne donnera pas d'estimation impartiale, sauf si vous pouvez toujours être sûr à 100% si elle est biaisée ou non. Dès qu'un exemple biaisé peut faire partie de votre ensemble d'entraînement avec une probabilité non nulle, il y aura un biais, car vous n'avez rien pour annuler ce biais. En pratique, votre biais sera simplement multiplié par un facteurα<1, où α est la probabilité qu'un exemple biaisé soit détecté comme tel.

En supposant que vous disposez de suffisamment de données, une meilleure approche consiste à calculer P(Z=bjeunese|X,Oui)pour chaque échantillon, puis supprimez tous les échantillons de l'ensemble d'apprentissage lorsque cette probabilité dépasse un certain seuil. Par exemple, s'il vous est possible de former votre ensemble de données sur les échantillons pour lesquelsP(Z=bjeunese|X,Oui)<0,01et votre jeu de données diminue de N biaisé et M impartial à n biaisé et m des exemples impartiaux, et le biais se multipliera par un facteur F=n(N+M)N(n+m). Depuis généralementnN sera bien inférieur à mM, F sera beaucoup plus petit que 1, résultant en une amélioration significative.

Notez que les deux techniques peuvent être combinées: lignes avec p=P(Z=bjeunese|X,Oui)>β sortir (pour un choix de β, au-dessus j'ai utilisé β=0,01), et les rangées qui restent obtiennent un poids (1-pβ)2, cela devrait vous donner le meilleur des deux mondes.

user1111929
la source
Je suis d'accord que la régression pondérée ne produira généralement pas d'estimations non biaisées à moins que le classifieur ne soit précis à 100%. Est-il garanti de réduire les biais? Votre approche de coupure est-elle nécessairement meilleure ou cela dépend-il de la taille des échantillons et de la précision du classificateur?
Adrian
Ce n'est pas vraiment la précision du classificateur, principalement la quantité de lignes qui en résulte, car pour de très petits ensembles de données, on ne peut parfois pas se permettre de supprimer un grand nombre de lignes. Mais étant donné suffisamment de données, je ne vois pas de raison de douter: votre approche pèse les lignes avec 50% de risque de biais aussi moitié que les lignes avec 0% de risque de biais. Personnellement, je ne ferais jamais ça, je préférerais fortement 1000 lignes propres garanties plutôt que 2000 lignes qui ont chacune 50% de chances de biais. Notez que vous pouvez également combiner les deux approches, pour appliquer un système plus progressif qu'une simple coupure, j'ai mis à jour ma réponse pour développer cela.
user1111929