Usando pesos de regressão quando pode ser medido com erro de medição médio diferente de zero

8

Suponha que observemos os dados e gostaríamos de ajustar um modelo de regressão para . Infelizmente, às vezes é medido com erros cuja média é diferente de zero.Y,XE[Y|X]Y

Deixe indicando se é medido com erros médios clássicos de zero ou erros médios diferentes de zero, respectivamente. Gostaríamos de estimar . Infelizmente, geralmente não é observado e . Se ajustarmos uma regressão de em , obteremos previsões tendenciosas.Z{unbiased,biased}YE[Y|X,Z=unbiased]ZE[Y|X,Z=unbiased]E[Y|X]YX

Suponha que geralmente não podemos observar , mas temos acesso a um modelo para (porque aprendemos Z manualmente em um pequeno conjunto de treinamento e ajustamos um modelo de classificação com Z como variável de destino) . Ajustar uma regressão de Y em X usando \ Pr [Z = \ text {imparcial} \, | \, X, Y] como pesos de regressão produz uma estimativa imparcial de \ mathbf {E} [Y \, | \, X, Z = \ text {imparcial}] (ou, na sua falta, uma estimativa menos tendenciosa do que obteríamos sem o uso de pesos)? Este método é usado na prática e possui um nome?ZPr[Z|X,Y]ZZYXPr[Z=unbiased|X,Y]E[Y|X,Z=unbiased]

Esclarecimento: o objetivo é ajustar um modelo que minimize o erro quadrático médio em dados não vistos (dados de teste) em que Z=unbiased . O preditor ideal para esse objetivo é E[Y|X,Z=unbiased] , portanto essa é a função que estamos tentando estimar. Os métodos para resolver esse problema devem ser classificados em termos de quão bem eles atingem esse objetivo.


Pequeno exemplo em R com df$y_is_unbiasedo papel de Z e df$y_observedo papel de 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))

Neste exemplo, o modelo é uma floresta aleatória com . Se esse modelo fosse perfeitamente preciso, ele geraria pesos de 1,0, onde é imparcial, 0,0, onde é tendencioso, e a regressão ponderada seria claramente imparcial. O que acontece quando o modelo para possui precisão de teste e recalls que não são perfeitos (<100% de precisão)? A regressão ponderada é garantida como menos tendenciosa do que uma regressão não ponderada de em ?Pr[Z=unbiased|X,Y]Y Y Pr [ Z = imparcialformula=y_is_unbiased ~ y_observed + x1 + x2YYPr[Z=unbiased|X,Y]YX


Exemplo ligeiramente mais complexo no qual varia com (em oposição ao exemplo mais simples que eu postei acima, onde ):Pr[Z=unbiased|X]XPr[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))

Neste exemplo, a regressão ponderada de em parece menos tendenciosa que a regressão não ponderada. Isso é verdade em geral? Eu também tentei a sugestão de shabbychef (veja a resposta abaixo) neste exemplo, e ela parece pior do que a regressão ponderada.YX


Para quem prefere o Python ao R, aqui está a segunda simulação no 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
fonte
1
Isso quase soa como "Variáveis ​​instrumentais", onde você observa algumas variáveis ​​correlacionadas com o erro em sua regressão. Receio que isso não ajude muito, no entanto.
21819 shabbychef
@shabbychef Verdadeiro, mas não há instrumentos disponíveis nesta configuração.
Adrian
1
Com o seu problema atualizado, o viés agora é uma função do e devemos esperar que os coeficientes de regressão sejam alterados. Ou seja, o termo 'viés' é que . Eu expandiria com uma expansão de Taylor para mostrar que existe uma dependência linear extra de em e . Voltando à sua pergunta original, o termo de viés que você vê e a variável você observa, não alteram a variação de forma alguma, mas alteram o valor esperado. Então, eles deveriam ser invadidos pela especificação linear, eu pensaria, e não pelos pesos. - 0,25 pxi0.25pp y x 1 x 2 zp=logistic(x1+2x2)pyx1x2z
shabbychef

Respostas:

2

Eu usaria a 'probabilidade de parcialidade' como uma variável dummy na regressão; possivelmente 'absorve' o viés presente no caso tendencioso. Usando o seu exemplo (mas ligando set.seed(1234)antes da ligação para get_df), tentei

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

e pegou:

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

O coeficiente para o termo 1-pr_unbiaseddeve ser o tamanho do viés.

shabbychef
fonte
Ideia interessante (+1)! Atualizei minha postagem com um segundo exemplo um pouco mais complexo, no qual varia com vez de ser constante. Nesse caso, tanto a constante quanto o coeficiente em x2 são enviesados ​​ao fazer a regressão de em , e não acho que seu método funcione tão bem (já que ele remove o viés da constante). Dê uma olhada e deixe-me saber o que você pensa! XPr[Z=unbiased|X]XXYX
Adrian
2

Esse é um problema de variável omitida em que você tem uma variável indicadora que não é observada, mas que tem um relacionamento com a variável de resposta. Como "viés" é uma propriedade de um estimador, não uma variável de regressão, vou reformular sua pergunta como aquela em que você deseja encontrar a verdadeira função de regressão condicional em usando dados de regressão que não incluem essa variável e um conjunto separado de dados de treinamento de regressão usado para estimar as probabilidades .ZZ=0p0(x,y)P(Z=0|X=x,Y=y)

Seja denotar a densidade condicional da variável de resposta no problema de regressão com a variável de resposta e a variável explicativa (mas excluindo ). A partir das regras de probabilidade condicional, a distribuição alvo de interesse pode ser escrita como: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).

Assim, podemos ver que é suficiente estimar a função de regressão no modelo de regressão com omitido e também estimar a função de probabilidade que você possui como estimador separado dos dados de treinamento. O primeiro pode ser estimado usando a estimativa OLS sem impor pesos. A "ponderação" ocorre após a estimativa desta função, por substituição na equação acima.pY|XZp0

Podemos ver que não é necessário (ou desejado) para utilizar quaisquer pesos na regressão de em , uma vez que é suficiente para estimar a densidade condicional sem consideração de . A estimativa do OLS dos coeficientes dessa regressão fornece um estimador , e como você também possui um estimador separado, você tem:YXpY|XZp^Y|Xp^0

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

Depois de substituir esses estimadores, resta apenas tentar determinar a constante de escala que produz uma função de densidade adequada. Isso pode ser feito por vários métodos de integração numérica (por exemplo, regra de Simpson, quadratura, Metropolis-Hastings, etc.).

Ben - Restabelecer Monica
fonte
1
Obrigado por esta resposta detalhada (+1), vou ler com atenção e retornaremos para você. Concordo que uma descrição mais correta do problema pode ser "Y às vezes é medido com erro de média diferente de zero", em vez de "Y é tendencioso".
Adrian
Interessante! Uma desvantagem / detalhe complicado aqui é que isso requer estimar a distribuição completa de , em vez de apenas a média condicional. Normalidade é uma suposição comum, mas isso pode não se aplicar ao meu aplicativo. Y|X
Adrian
Sim, mas isso não deve ser realmente complicado. Qualquer que seja o modelo de regressão usado, terá uma forma clara de modelo que determina essa distribuição condicional.
Ben - Restabelece Monica
1
Ben, você tem certeza de que isso é verdade em um cenário de aprendizado / previsão de máquina? Ajustar uma regressão linear para minimizar o erro quadrático médio nos dados de teste não exige nenhuma suposição sobre a normalidade dos resíduos (desde que você não esteja realizando testes de hipóteses de amostra finita, não esteja interessado em intervalos de previsão, não se importe se seu estimador é um estimador de probabilidade máxima etc.). Estou interessado em estimar uma média condicional (de dado ) e não fizeram suposições sobre a distribuição de . X YYXY
Adrian
1
Não acho que Adrian esteja necessariamente interessado na distribuição completa de (Y | X, Z = 0), apenas em sua média. Se alguém quiser conhecer a distribuição completa de Y | X, Z = 0, é claro que as distribuições exatas de Y e X são muito relevantes, mas se você quiser apenas estimar E [Y | X, Z = 0], então isso é não é necessário: a lei do grande número funciona para distribuições arbitrárias.
user1111929
1

Sua ideia não fornecerá uma estimativa imparcial, a menos que você sempre tenha 100% de certeza se é parcial ou não. Assim que um exemplo tendencioso puder fazer parte do seu conjunto de treinamento com probabilidade diferente de zero, haverá viés, pois você não terá nada para cancelar esse viés. Na prática, seu viés será simplesmente multiplicado por um fator , onde é a probabilidade de um exemplo tendencioso ser detectado como tal.α<1α

Supondo que você tenha dados suficientes, uma abordagem melhor é calcular para cada amostra e remover todas as amostras do conjunto de treinamento em que essa probabilidade exceda um determinado limite. Por exemplo, se for possível treinar seu conjunto de dados apenas nas amostras para as quais , e seu conjunto de dados diminuir de tendencioso e imparcial para tendencioso e imparcial, e o viés será multiplicado pelo fator . Como normalmente será muito menor que , será muito menor queP(Z=biased|X,Y)P(Z=biased|X,Y)<0.01NMnmf=n(N+M)N(n+m)nNmMf1, resultando em uma melhoria significativa.

Observe que as duas técnicas podem ser combinadas: linhas com saem (para alguma opção de , acima de eu usei ) e as linhas que permanecem em obtenha um peso , que deve realmente oferecer o melhor dos dois mundos.p=P(Z=biased|X,Y)>βββ=0.01(1pβ)2

user1111929
fonte
Concordo que a regressão ponderada geralmente não produzirá estimativas imparciais, a menos que o classificador seja 100% exato. É garantido reduzir o viés? Sua abordagem de corte é necessariamente melhor ou depende do tamanho da amostra e da precisão do classificador?
Adrian
Não é exatamente a precisão do classificador, principalmente a quantidade resultante de linhas, pois, para conjuntos de dados muito pequenos, às vezes não é possível descartar um grande número de linhas. Mas, com dados suficientes, não vejo motivo para dúvidas: sua abordagem pesa linhas com 50% de chance de viés e metade da relevância das linhas com 0% de chance de viés. Pessoalmente, eu nunca faria isso, eu preferiria 1000 linhas limpas garantidas em vez de 2000, com 50% de chance de viés. Observe que você também pode combinar as duas abordagens; para aplicar um sistema mais gradual do que um ponto de corte simples, atualizei minha resposta para elaborar isso.
user1111929