Пример реального анализа данных с помощью DRF
В двух предыдущих статьях я объяснял Distributional Random Forests (DRF), случайный лес, способный оценивать условные распределения, а также расширение метода, позволяющее количественно определять неопределенность, например доверительные интервалы и т. д. Здесь Я представляю пример реального приложения для данных о заработной плате из исследования американского сообщества 2018 года, проведенного Бюро переписи населения США. В первом документе DRF мы получили данные примерно об 1 миллионе штатных сотрудников из опроса американского сообщества 2018 года, проведенного Бюро переписи населения США, из которого мы извлекли информацию о заработной плате и все ковариаты, которые могут иметь отношение к заработной плате. Это богатство данных идеально подходит для экспериментов с таким методом, как DRF (на самом деле мы будем использовать только крошечное подмножество для этого анализа).
При изучении необработанных данных о почасовой оплате труда между двумя полами наблюдается постоянный разрыв в том, что мужчины, как правило, зарабатывают больше. Интересен вопрос, является ли наблюдаемый разрыв в почасовой оплате труда (Вт) мужчин (G=1) и женщин (G=0) только из-за пола или может быть объяснено некоторыми другими мешающими переменными X, которые зависят от пола и, в свою очередь, влияют на заработную плату. То есть мы хотим изучить величину эффекта, соответствующую жирной стрелке на следующем графике причинно-следственных связей:
Например, предположим, что X включает только профессию и что женщины склонны выбирать профессии, не требующие большого денежного вознаграждения, такие как врачи, медсестры или учителя. , в то время как мужчины, как правило, имеют профессиональную работу в азартных играх с непристойной почасовой оплатой. Если бы это само по себе определяло разницу в почасовой оплате труда между полами, мы бы все равно увидели разрыв в оплате труда, глядя непосредственно на данные о почасовой оплате труда. Однако если мы затем зафиксируем профессию (X) только для врачей и сравним два распределения заработной платы там, любое статистически значимое различие может быть связано только с полом.
Мы ориентируемся на двухэтапный анализ:
- Мы фиксируем X на определенном значении и сравниваем распределение заработной платы в двух группах для ковариат, фиксированных на X=х. Это интересно с двух точек зрения: во-первых, если X действительно включает все другие факторы, влияющие на заработную плату и связанные с полом, то фиксация X =x и рассмотрение заработной платы для обоих полов означает, что мы действительно наблюдаем влияние пола на заработную плату. Во-вторых, он позволяет прогнозировать распределение заработной платы в целом для человека с заданными характеристиками x.
- Мы используем приведенный выше предполагаемый причинно-следственный график и правила причинно-следственной связи для оценки гипотетического распределения с помощью DRF: Распределение заработной платы женщин, если бы они рассматривались как мужчины при установлении заработной платы. Если Xсодержит все соответствующие ковариаты и гендерный разрыв в оплате труда отсутствует, это распределение должно быть таким же, как распределение заработной платы для мужчин (без учета статистической случайности).
Эта статья является кульминацией работы нескольких человек: код и набор данных были получены из исходного репозитория DRF, а затем объединены с методами, разработанными в нашей новой статье на arXiv, написанной совместно с Коринн Эменеггер.
Прежде чем продолжить, я хочу отметить, что это лишь пример, иллюстрирующий использование DRF. Я не хочу делать здесь никаких серьезных (причинно-следственных) утверждений просто потому, что анализ, безусловно, в некотором отношении ошибочен, а причинно-следственная диаграмма, которую мы предполагаем ниже, безусловно, неверна. Более того, мы используем только небольшое подмножество доступных данных.
Также обратите внимание, что код работает довольно медленно. Это связано с тем, что, хотя сам DRF закодирован на C, повторная подгонка, необходимая для доверительных интервалов, пока реализована на R.
Тем не менее, давайте углубимся. Далее все изображения, если не указано иное, принадлежат автору.
Данные
Данные PUMS (Общественное использование микроданных) из 1-летнего исследования американского сообщества 2018 г. получены из API Бюро переписи населения США. Опрос рассылается ≈ 3,5 миллионам человек ежегодно и направлен на предоставление более актуальных данных, чем официальная перепись, которая проводится раз в десять лет. Набор данных за 2018 год содержит около 3 миллионов анонимных точек данных для 51 штата и округа Колумбия. Для документа DRF, указанного выше, мы получили только подмножество переменных, которые могут иметь отношение к заработной плате, таких как пол человека, возраст, раса, семейное положение человека, уровень образования и уровень знания английского языка.
Предобработанные данные можно найти здесь. Сначала делаем дополнительную очистку:
##Further data cleaning ##
which = rep(TRUE, nrow(wage))
which = which & (wage$age >= 17)
which = which & (wage$weeks_worked > 48)
which = which & (wage$hours_worked > 16)
which = which & (wage$employment_status == 'employed')
which = which & (wage$employer != 'self-employed')
which[is.na(which)] = FALSE
data = wage[which, ]
sum(is.na(data))
colSums(is.na(data))
rownames(data) = 1:nrow(data)
#data = na.omit(data)
data$log_wage = log(data$salary / (data$weeks_worked * data$hours_worked))
## Prepare data and fit drf
## Define X and Y
X = data[,c(
'age',
'race',
'hispanic_origin',
'citizenship',
'nativity',
'marital',
'family_size',
'children',
'education_level',
'english_level',
'economic_region'
)]
X$occupation = unlist(lapply(as.character(data$occupation), function(s){return(substr(s, 1, 2))}))
X$occupation = as.factor(X$occupation)
X$industry = unlist(lapply(as.character(data$industry), function(s){return(substr(s, 1, 2))}))
X$industry[X$industry %in% c('32', '33', '3M')] = '31'
X$industry[X$industry %in% c('42')] = '41'
X$industry[X$industry %in% c('45', '4M')] = '44'
X$industry[X$industry %in% c('49')] = '48'
X$industry[X$industry %in% c('92')] = '91'
X$industry = as.factor(X$industry)
X=dummy_cols(X, remove_selected_columns = TRUE)
X = as.matrix(X)
Y = data[,c('sex', 'log_wage')]
Y$sex = (Y$sex == 'male')
Y = as.matrix(Y)
На самом деле это гораздо больше наблюдений, чем нам нужно, и вместо этого мы случайным образом отбираем 4000 точек обучающих данных для анализа.
train_idx = sample(1:nrow(data), 4000, replace = FALSE) ## Focus on training data Ytrain=Y[train_idx,] Xtrain=X[train_idx,]
Опять же, это просто иллюстрация — на самом деле вам нужно взять как можно больше точек данных. Предполагаемая плотность заработной платы обоих полов для этих 4000 точек данных представлена на рисунке 1 с этим кодом:
## Plot the test data without adjustment
plotdfunadj = data[train_idx, ]
plotdfunadj$weight=1
plotdfunadj$plotweight[plotdfunadj$sex=='female'] = plotdfunadj$weight[plotdfunadj$sex=='female']/sum(plotdfunadj$weight[plotdfunadj$sex=='female'])
plotdfunadj$plotweight[plotdfunadj$sex=='male'] = plotdfunadj$weight[plotdfunadj$sex=='male']/sum(plotdfunadj$weight[plotdfunadj$sex=='male'])
#pooled data
ggplot(plotdfunadj, aes(log_wage)) +
geom_density(adjust=2.5, alpha = 0.3, show.legend = TRUE, aes(fill=sex, weight=plotweight)) +
theme_light()+
scale_fill_discrete(name = "gender", labels = c('female', "male"))+
theme(legend.position = c(0.83, 0.66),
legend.text=element_text(size=18),
legend.title=element_text(size=20),
legend.background = element_rect(fill=alpha('white', 0.5)),
axis.text.x = element_text(size=14),
axis.text.y = element_text(size=14),
axis.title.x = element_text(size=19),
axis.title.y = element_text(size=19))+
labs(x='log(hourly_wage)')
Вычисление процентной медианной разницы между двумя зарплатами, т.
(медианная заработная плата мужчин – средняя заработная плата женщин)/(медианная заработная плата женщин) *100,
мы получаем около 18 процентов. То есть медианная заработная плата мужчин на 18 процентов выше, чем у женщин по нескорректированным данным (!)
## Median Difference before adjustment! quantile_maleunadj = wtd.quantile(x=plotdfunadj$log_wage, weights=plotdfunadj$plotweight*(plotdfunadj$sex=='male'), normwt=TRUE, probs=0.5) quantile_femaleunadj = wtd.quantile(x=plotdfunadj$log_wage, weights=plotdfunadj$plotweight*(plotdfunadj$sex=='female'), normwt=TRUE, probs=0.5) (1-exp(quantile_femaleunadj)/exp(quantile_maleunadj))
Анализ
Теперь возникает вопрос, действительно ли это «несправедливо». То есть мы предполагаем причинно-следственный график выше, где пол (G) влияет как на заработную плату (W), так и на ковариаты X, которые, в свою очередь, влияют на W. Мы хотели бы знать, влияет ли гендер напрямую на заработную плату (жирная стрелка). То есть, если женщина и мужчина с одинаковыми характеристиками X=x получают одинаковую заработную плату, или она получает меньше, просто из-за ее пола.
Мы изучим это в двух условиях. Первый заключается в том, чтобы действительно зафиксировать X=x и использовать механизмы, описанные в предыдущей статье. Интуитивно, если мы зафиксируем все другие ковариаты, которые могут влиять на заработную плату помимо пола, и теперь сравним два распределения заработной платы, то любое наблюдаемое различие должно быть связано только с заработной платой.
Второй пытается количественно определить эту разницу по всем возможным значениям X. Это делается путем вычисления гипотетического распределения
W(мужчина, X(женщина)).
Эта величина представляет собой предполагаемую заработную плату, которую получает мужчина, если он обладает в точности женскими характеристиками. То есть мы просим заработную плату, которую получает женщина, если к ней относятся как к мужчине.
Обратите внимание, что это предполагает, что приведенный выше причинный график верен. В частности, предполагается, что Xучитывает все релевантные факторы, кроме пола, которые определяют заработную плату. Вполне может быть, что это не так, поэтому отказ от ответственности в начале этой статьи.
Изучение условных распределительных различий
Далее мы фиксируем x в произвольной точке:
i<-47 # Important: Test point needs to be a matrix test_point<-X[i,, drop=F]
На следующем рисунке показаны некоторые значения, содержащиеся в этой контрольной точке x — мы рассматриваем воспитателей с дипломом средней школы, которые состоят в браке и имеют 1 ребенка. С помощью DRF мы можем оценить и построить график плотности для двух групп при условии X=x:
# predict with the new framework
DRF = predictdrf(drf_fit, x=x)
weights <- DRF$weights
## Conditional Density Plotting
plotdfx = data[train_idx, ]
propensity = sum(weights[plotdfx$sex=='female'])
plotdfx$plotweight = 0
plotdfx$plotweight[plotdfx$sex=='female'] = weights[plotdfx$sex=='female']/propensity
plotdfx$plotweight[plotdfx$sex=='male'] = weights[plotdfx$sex=='male']/(1-propensity)
gg = ggplot(plotdfx, aes(log_wage)) +
geom_density(adjust=5, alpha = 0.3, show.legend=TRUE, aes(fill=sex, weight=plotweight)) +
labs(x='log(hourly wage)')+
theme_light()+
scale_fill_discrete(name = "gender", labels = c(sprintf("F: %g%%", round(100*propensity, 1)), sprintf("M: %g%%", round(100*(1-propensity), 1))))+
theme(legend.position = c(0.9, 0.65),
legend.text=element_text(size=18),
legend.title=element_text(size=20),
legend.background = element_rect(fill=alpha('white', 0)),
axis.text.x = element_text(size=14),
axis.text.y = element_text(size=14),
axis.title.x = element_text(size=19),
axis.title.y = element_text(size=19))+
annotate("text", x=-1, y=Inf, hjust=0, vjust=1, size=5, label = point_description(data[i,]))
plot(gg)
На этом графике видна явная разница в заработной плате, даже в случае этого фиксированного x (помните, что в этом случае все предполагаемые вмешивающиеся факторы фиксированы, поэтому мы просто напрямую сравниваем заработную плату). Теперь с помощью DRF мы оцениваем и тестируем среднюю разницу.
## Getting the respective weights
weightsmale<-weights*(Ytrain[, "sex"]==1)/sum(weights*(Ytrain[, "sex"]==1))
weightsfemale<-weights*(Ytrain[, "sex"]==0)/sum(weights*(Ytrain[, "sex"]==0))
## Choosing alpha:
alpha<-0.05
# Step 1: Doing Median comparison for fixed x
quantile_male = wtd.quantile(x=data$log_wage[train_idx], weights=matrix(weightsmale), normwt=TRUE, probs=0.5)
quantile_female = wtd.quantile(x=data$log_wage[train_idx], weights=matrix(weightsfemale), normwt=TRUE, probs=0.5)
(medianx<-unname(1-exp(quantile_female)/exp(quantile_male)))
mediandist <- sapply(DRF$weightsb, function(wb) {
wbmale<-wb*(Ytrain[, "sex"]==1)/sum(wb*(Ytrain[, "sex"]==1))
wbfemale<-wb*(Ytrain[, "sex"]==0)/sum(wb*(Ytrain[, "sex"]==0))
quantile_maleb = wtd.quantile(x=data$log_wage[train_idx], weights=matrix(wbmale), normwt=TRUE, probs=0.5)
quantile_femaleb = wtd.quantile(x=data$log_wage[train_idx], weights=matrix(wbfemale), normwt=TRUE, probs=0.5)
return( unname(1-exp(quantile_femaleb)/exp(quantile_maleb)) )
})
varx<-var(mediandist)
## Use Gaussian CI:
(upper<-medianx + qnorm(1-alpha/2)*sqrt(varx))
(lower<-medianx - qnorm(1-alpha/2)*sqrt(varx))
Это дает доверительный интервал средней разницы
(0.06, 0.40) or (6%, 40%)
Совершенно очевидно, что этот интервал не содержит нуля, и, таким образом, медиана разницы действительно значительна.
Используя функцию Witobj, мы можем сделать эту разницу более заметной.
Witobj<-Witdrf(drf_fit, x=test_point, groupingvar="sex", alpha=0.05)
hatmun<-function(y,Witobj){
c<-Witobj$c
k_Y<-Witobj$k_Y
Y<-Witobj$Y
weightsall1<-Witobj$weightsall1
weightsall0<-Witobj$weightsall0
Ky=t(kernelMatrix(k_Y, Y , y = y))
out<-list()
out$val <- tcrossprod(Ky, weightsall1 ) - tcrossprod(Ky, weightsall0 )
out$upper<- out$val+sqrt(c)
out$lower<- out$val-sqrt(c)
return( out )
}
all<-hatmun(sort(Witobj$Y),Witobj)
plot(sort(Witobj$Y),all$val , type="l", col="darkblue", lwd=2, ylim=c(min(all$lower), max(all$upper)),
xlab="log(wage)", ylab="witness function")
lines(sort(Witobj$Y),all$upper , type="l", col="darkgreen", lwd=2 )
lines(sort(Witobj$Y),all$lower , type="l", col="darkgreen", lwd=2 )
abline(h=0)
что приводит к рисунку:
Мы ссылаемся на сопутствующую статью для более подробного объяснения этой концепции. По существу его можно рассматривать как
условная плотность логарифма(зарплаты) мужчин с учетом x — условная плотность логарифма(зарплаты) женщин с учетом x
То есть функция условного свидетеля показывает, где плотность одной группы больше, чем другой, без необходимости фактической оценки плотности. В этом примере отрицательные значения означают, что плотность заработной платы женщин выше, чем у мужчин при условии x, а положительные значения означают, что плотность заработной платы женщин ниже. Поскольку мы уже оценили условные плотности выше, одна только условная функция-свидетель ничего не добавляет. Но это хорошо для иллюстрации. Действительно, мы видим, что вначале она отрицательна для значений, где условная плотность заработной платы женщин выше, чем условная плотность заработной платы мужчин. Наоборот, она становится положительной для больших величин, при которых условная плотность заработной платы мужчин выше, чем женщин. Таким образом, релевантная информация о двух плотностях суммируется на графике свидетельской функции: мы видим, что плотность заработной платы женщин выше при более низких значениях заработной платы и ниже при более высоких значениях, что указывает на то, что плотность смещена влево и женщины зарабатывают меньше. ! Кроме того, мы также можем предоставить 95% доверительные интервалы зеленого цвета, которые включают истинную функцию с 95%, равномерно по всем y. (Хотя действительно нужно много данных, чтобы сделать это достоверным). Поскольку эта однородная доверительная полоса не содержит нулевой линии между примерно 2 и 2,5, мы снова видим, что разница между двумя распределениями статистически значима.
Обусловленность определенным x позволяет нам детально изучить индивидуальные эффекты и с учетом неопределенности. Тем не менее, также может быть интересно изучить общий эффект. Мы делаем это, оценивая гипотетическое распределение в следующем разделе.
Оценка контрфактического распределения
Используя расчетные законы причинности на нашем предполагаемом причинно-следственном графике, можно вывести, что:
т.е. искомое распределение контрфактического значения получается путем усреднения условного распределения W | G=мужчина, X=x, больше x пола жен.
Поскольку распределения задаются как простые веса, это легко сделать с помощью DRF следующим образом:
## Add code ## Male is 1, Female is 0 # obtain all X from the female test population Xtestf<-Xtest[Ytest[,"sex"]==0,] # Obtain the conditional distribution of W | G=male, X=x, for x in the female # population. # These weights correspond to P(W, G=male | X=x ) weightsf<-predictdrf(drf_fit, x=Xtestf)$weights*(Ytrain[, "sex"]==1) weightsf<-weightsf/rowSums(weightsf) # The counterfactual distribution is the average over those weights/distributions counterfactualw<-colMeans(weightsf)
что приводит к следующей гипотетической оценке плотности:
plotdfc<-rbind(plotdfc, plotdfunadj[plotdfunadj$sex=='female',])
plotdfc$sex2<-c(rep(1, length(train_idx)), rep(0,nrow(plotdfunadj[plotdfunadj$sex=='female',])))
plotdfc$sex2<-factor(plotdfc$sex2)
#interventional distribution
ggplot(plotdfc, aes(log_wage)) +
geom_density(adjust=2.5, alpha = 0.3, show.legend=TRUE, aes(fill=sex2, weight=plotweight)) +
theme_light()+
scale_fill_discrete(name = "", labels = c("observed women's wages", "wages if treated as men"))+
theme(legend.position = c(0.2, 0.98),
legend.text=element_text(size=16),
legend.title=element_text(size=20),
legend.background = element_rect(fill=alpha('white', 0)),
axis.text.x = element_text(size=14),
axis.text.y = element_text(size=14),
axis.title.x = element_text(size=19),
axis.title.y = element_text(size=19))+
labs(x='log(hourly wage)')
Эти две плотности теперь представляют собой плотность заработной платы женщин, выделенную красным цветом, и плотность женщин, если бы с ними обращались как с мужчинами за установление заработной платы, выделенную зелено-голубоватым цветом. Очевидно, что плотности теперь ближе друг к другу, чем раньше — с поправкой на искажающие факторы разница в оплате труда мужчин и женщин стала меньше. Однако средняя разница по-прежнему
quantile_male = wtd.quantile(x=plotdfc$log_wage[plotdfc$sex2==1], weights=counterfactualw, normwt=TRUE, probs=0.5) quantile_female = wtd.quantile(x=plotdfunadj$log_wage, weights=plotdfunadj$plotweight*(plotdfunadj$sex=='female'), normwt=TRUE, probs=0.5) (1-exp(quantile_female)/exp(quantile_male))
0,11 или 11 процентов!
Таким образом, если наш анализ верен, 11% разницы в оплате труда можно отнести только к полу. Другими словами, несмотря на то, что нам удалось сократить разницу с 18 % медианного дохода в нескорректированных данных до 11 %, по-прежнему остается существенная разница, указывающая на «несправедливый» разрыв в оплате труда между полами (по крайней мере, если X действительно улавливает соответствующие помехи).
Заключение
В этой статье мы рассмотрели пример того, как можно использовать DRF для реального анализа данных. Мы исследовали как случай фиксированного x, для которого рассмотренные в статье методы позволяют строить меры неопределенности, так и распределение контрфактической величины. В обоих случаях мы увидели существенную, а в случае фиксированного x существенную разницу при корректировке доступных искажающих факторов.
Хотя я не проверял, может быть интересно посмотреть, как результаты этого маленького эксперимента соотносятся с более серьезным анализом. В любом случае, я надеюсь, что эта статья продемонстрировала, как можно использовать DRF в реальном анализе данных.
Дополнительный код для создания графиков
## Step 0: Choosing x
point_description = function(test_point){
out = ''
out = paste(out, 'job: ', test_point$occupation_description[1], sep='')
out = paste(out, '\nindustry: ', test_point$industry_description[1], sep='')
out = paste(out, '\neducation: ', test_point$education[1], sep='')
out = paste(out, '\nemployer: ', test_point$employer[1], sep='')
out = paste(out, '\nregion: ', test_point$economic_region[1], sep='')
out = paste(out, '\nmarital: ', test_point$marital[1], sep='')
out = paste(out, '\nfamily_size: ', test_point$family_size[1], sep='')
out = paste(out, '\nchildren: ', test_point$children[1], sep='')
out = paste(out, '\nnativity: ', test_point$nativity[1], sep='')
out = paste(out, '\nhispanic: ', test_point$hispanic_origin[1], sep='')
out = paste(out, '\nrace: ', test_point$race[1], sep='')
out = paste(out, '\nage: ', test_point$age[1], sep='')
return(out)
}
