Menggunakan bobot regresi ketika mungkin diukur dengan kesalahan pengukuran bukan-rata-rata

8

Misalkan kita mengamati data dan ingin mencocokkan model regresi untuk . Sayangnya, kadang-kadang diukur dengan kesalahan yang rata-rata bukan nol.Y,XE[Y|X]Y

Biarkan menunjukkan apakah diukur dengan kesalahan rata-rata nol klasik atau kesalahan bukan-rata. Kami ingin memperkirakan . Sayangnya, Z umumnya tidak diamati, dan \ mathbf {E} [Y \, | \, X, Z = \ text {bias}] \ neq \ mathbf {E} [Y \, | \, X] . Jika kita cocok dengan regresi Y pada X , kita akan mendapatkan prediksi yang bias.Z{unbiased,biased}YE[Y|X,Z=unbiased]ZE[Y|X,Z=unbiased]E[Y|X]YX

Misalkan kita secara umum tidak dapat mengamati Z , tetapi memiliki akses ke model untuk Pr[Z|X,Y] (karena kita secara manual mempelajari Z pada set pelatihan kecil dan menyesuaikan model klasifikasi dengan Z sebagai variabel target) . Apakah sesuai dengan regresi Y pada X menggunakan Pr[Z=unbiased|X,Y] karena bobot regresi menghasilkan estimasi yang tidak bias dari E[Y|X,Z=unbiased] (atau, gagal itu, estimasi yang kurang bias dari yang akan kita dapatkan tanpa menggunakan bobot)? Apakah metode ini digunakan dalam praktik, dan apakah ada nama?

Klarifikasi: tujuannya adalah untuk mencocokkan model yang meminimalkan kesalahan kuadrat rata-rata pada data yang tidak terlihat (data uji) di mana Z=tidak bias . Prediktor optimal untuk tujuan itu adalah E[Y|X,Z=unbiased] , jadi itulah fungsi yang kami coba perkirakan. Metode untuk memecahkan masalah ini harus diberi peringkat dalam hal seberapa baik mereka mencapai tujuan itu.


Contoh kecil dalam R dengan df$y_is_unbiasedmemainkan peran Z dan df$y_observedmemainkan peran 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))

Dalam contoh ini, model adalah hutan acak dengan . Jika model ini sangat akurat, itu akan menghasilkan bobot 1,0 di mana tidak bias, 0,0 di mana bias, dan regresi tertimbang jelas akan tidak bias. Apa yang terjadi ketika model untuk memiliki presisi pengujian dan penarikan yang tidak sempurna (akurasi <100%)? Apakah regresi berbobot dijamin kurang bias dibandingkan dengan regresi tertimbang pada ?Pr[Z=unbiased|X,Y]formula=y_is_unbiased ~ y_observed + x1 + x2YYPr[Z=unbiased|X,Y]YX


Contoh yang sedikit lebih rumit di mana bervariasi dengan (sebagai lawan dari contoh sederhana yang saya posting di atas, di mana ):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))

Dalam contoh ini, regresi tertimbang pada terlihat kurang bias dibandingkan dengan regresi tanpa bobot. Apakah itu benar secara umum? Saya juga mencoba saran shabbychef (lihat jawaban di bawah) pada contoh ini, dan tampaknya lebih buruk daripada regresi berbobot.YX


Bagi mereka yang lebih suka Python ke R, inilah simulasi kedua dengan 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
sumber
1
Ini hampir terdengar seperti "Variabel Instrumental", di mana Anda mengamati beberapa variabel yang berkorelasi dengan kesalahan dalam regresi Anda. Saya khawatir itu tidak banyak membantu.
shabbychef
@shabbychef Benar, tetapi tidak ada instrumen yang tersedia di pengaturan ini.
Adrian
1
Dengan masalah Anda yang diperbarui, bias sekarang merupakan fungsi dari , dan kami harus mengharapkan koefisien regresi berubah. Artinya, istilah 'bias' adalah mana . Saya akan memperluas dengan ekspansi Taylor untuk menunjukkan bahwa ada ketergantungan linear ekstra pada dan . Kembali ke pertanyaan awal Anda, istilah bias yang Anda lihat, dan variabel Anda amati, jangan ubah varians dengan cara apa pun, tetapi ubah nilai yang diharapkan. Jadi mereka harus diretas ke dalam spesifikasi linier, saya pikir, dan tidak ke bobot. xsaya-0,25halhal=logistik(x1+2x2)halyx1x2z
shabbychef

Jawaban:

2

Saya akan menggunakan 'probabilitas bias' sebagai variabel dummy dalam regresi; mungkin bisa 'menopang' bias yang ada dalam kasus bias. Menggunakan contoh Anda, (tetapi menelepon set.seed(1234)sebelum panggilan ke get_df), saya mencoba

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

dan mendapatkan:

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

Koefisien untuk istilah 1-pr_unbiasedharus menjadi ukuran bias.

shabbychef
sumber
Ide menarik (+1)! Saya memperbarui posting saya dengan contoh kedua yang sedikit lebih kompleks di mana bervariasi dengan alih-alih konstan. Dalam hal ini, baik konstanta dan koefisien pada x2 bias ketika melakukan regresi pada , dan saya tidak berpikir metode Anda bekerja juga (karena hanya menghilangkan bias dari konstanta). Lihatlah dan beri tahu saya apa yang Anda pikirkan! Pr[Z=tidak bias|X]XYX
Adrian
2

Ini adalah masalah variabel yang dihilangkan di mana Anda memiliki variabel indikator yang tidak teramati, tetapi yang memiliki hubungan dengan variabel respons. Karena "bias" adalah properti dari estimator, bukan variabel regresi, saya akan membingkai ulang pertanyaan Anda sebagai salah satu di mana Anda ingin menemukan fungsi regresi yang benar tergantung pada menggunakan data regresi yang tidak termasuk variabel ini, dan satu set data pelatihan regresi terpisah yang digunakan untuk memperkirakan probabilitas .ZZ=0hal0(x,y)P(Z=0|X=x,Y=y)

Misalkan menunjukkan kepadatan bersyarat dari variabel respons dalam masalah regresi dengan variabel respons dan variabel penjelas (tetapi tidak termasuk ). Dari aturan probabilitas kondisional, target distribusi bunga dapat ditulis sebagai:halY|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).

Dengan demikian, kita dapat melihat bahwa cukup untuk dapat memperkirakan fungsi regresi dalam model regresi dengan dihilangkan, dan juga memperkirakan fungsi probabilitas yang Anda miliki sebagai penduga terpisah dari data pelatihan Anda. Yang pertama dapat diestimasi menggunakan estimasi OLS tanpa memaksakan bobot apa pun. "Pembobotan" terjadi setelah estimasi fungsi ini, dengan substitusi ke dalam persamaan di atas.pY|XZp0

Kita bisa melihat bahwa itu tidak perlu (atau diinginkan) untuk menggunakan bobot apapun dalam regresi pada , karena itu sudah cukup untuk memperkirakan bersyarat kepadatan tanpa pertimbangan . Estimasi OLS dari koefisien regresi ini memberikan estimator , dan karena Anda juga memiliki estimator terpisah maka Anda memiliki:YXpY|XZp^Y|Xp^0

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

Setelah Anda mengganti penaksir ini, yang tersisa hanyalah mencoba menentukan tetapan penskalaan yang menghasilkan fungsi kerapatan yang tepat. Ini dapat dilakukan dengan serangkaian metode integrasi numerik (misalnya, aturan Simpson, quadrature, Metropolis-Hastings, dll.).

Ben - Pasang kembali Monica
sumber
1
Terima kasih atas jawaban terinci ini (+1), saya akan membacanya dengan cermat dan kembali kepada Anda. Saya setuju bahwa deskripsi masalah yang lebih tepat mungkin "Y kadang-kadang diukur dengan kesalahan bukan-nol", daripada "Y bias."
Adrian
Menarik! Satu detail downside / rumit di sini adalah bahwa ini memerlukan estimasi distribusi penuhY|X, bukan hanya rata-rata bersyarat. Normalitas adalah asumsi umum, tetapi itu mungkin tidak berlaku pada aplikasi saya.
Adrian
Ya, tapi itu seharusnya tidak terlalu rumit. Apapun model regresi yang Anda gunakan akan memiliki bentuk model yang jelas yang menentukan distribusi bersyarat ini.
Ben - Reinstate Monica
1
Ben, apakah Anda yakin itu benar dalam pengaturan pembelajaran mesin / prediksi? Menyesuaikan regresi linier untuk meminimalkan kesalahan kuadrat rata-rata pada data uji tidak perlu membuat asumsi tentang normalitas residual (selama Anda tidak melakukan tes hipotesis sampel terbatas, tidak tertarik pada interval prediksi, tidak peduli apakah penaksir Anda adalah penaksir kemungkinan maksimum, dll). Saya tertarik memperkirakan rata-rata bersyarat (dariY diberikan X) dan belum membuat asumsi tentang distribusi Y.
Adrian
1
Saya tidak berpikir Adrian perlu tertarik pada distribusi penuh (Y | X, Z = 0), hanya dalam rata-rata. Jika seseorang ingin mengetahui distribusi penuh Y | X, Z = 0 maka tentu saja distribusi Y dan X yang tepat sangat relevan, tetapi jika seseorang hanya ingin memperkirakan E [Y | X, Z = 0] maka ini adalah tidak perlu: hukum angka besar berlaku untuk distribusi sewenang-wenang.
user1111929
1

Gagasan Anda tidak akan memberikan perkiraan yang tidak bias, kecuali Anda selalu dapat 100% yakin apakah itu bias atau tidak. Segera setelah satu contoh bias dapat menjadi bagian dari rangkaian latihan Anda dengan probabilitas nol, akan ada bias, karena Anda tidak memiliki apa pun untuk membatalkan bias itu. Dalam praktiknya, bias Anda hanya akan dikalikan dengan faktorα<1dimana α adalah probabilitas bahwa contoh bias terdeteksi seperti itu.

Dengan asumsi Anda memiliki cukup data, pendekatan yang lebih baik adalah menghitung P(Z=bsayaSebuahsed|X,Y)untuk setiap sampel, dan kemudian menghapus semua sampel dari set pelatihan di mana probabilitas ini melebihi batas tertentu. Misalnya, jika layak bagi Anda untuk melatih dataset Anda hanya pada sampel yang manaP(Z=bsayaSebuahsed|X,Y)<0,01, dan dataset Anda berkurang dari N bias dan M. tidak bias untuk n bias dan m contoh yang tidak bias, dan bias akan dikalikan dengan faktor f=n(N+M.)N(n+m). Karena biasanyanN akan jauh lebih rendah dari mM, f akan jauh lebih kecil dari 1, menghasilkan peningkatan yang signifikan.

Perhatikan bahwa kedua teknik dapat dikombinasikan: baris dengan p=P(Z=biased|X,Y)>β keluar (untuk beberapa pilihan β, di atas saya gunakan β=0.01), dan baris yang tetap di dapatkan bobot (1pβ)2, yang seharusnya memberi Anda yang terbaik dari kedua dunia.

pengguna1111929
sumber
Saya setuju bahwa regresi berbobot umumnya tidak akan menghasilkan perkiraan yang tidak memihak kecuali penggolongnya 100% akurat. Apakah dijamin mengurangi bias? Apakah pendekatan cutoff Anda tentu lebih baik, atau apakah itu tergantung pada ukuran sampel dan akurasi classifier?
Adrian
Tidak terlalu akurat, terutama sebagian besar jumlah baris yang dihasilkan, karena untuk dataset yang sangat kecil, kadang-kadang orang tidak dapat membuang banyak baris. Tetapi dengan data yang cukup saya tidak melihat alasan untuk ragu: pendekatan Anda menimbang baris dengan peluang bias 50% setengah relevan dengan baris dengan peluang bias 0%. Secara pribadi, saya tidak akan pernah melakukan itu, saya lebih suka 1000 baris bersih dijamin lebih dari 2000 baris yang masing-masing memiliki 50% kemungkinan bias. Perhatikan bahwa Anda juga dapat menggabungkan kedua pendekatan, untuk menerapkan sistem yang lebih bertahap daripada cutoff biasa, saya telah memperbarui jawaban saya untuk menguraikan hal ini.
user1111929