Существует очень интригующее практическое решение
За исключением некоторых чрезмерно очищенных наборов данных, которые можно найти в Интернете, пропущенные значения есть везде. На самом деле, чем сложнее и больше набор данных, тем больше вероятность наличия пропущенных значений. Пропущенные значения — увлекательная область статистических исследований, но на практике они часто доставляют неудобства.
Если вы имеете дело с проблемой прогнозирования, когда вы хотите предсказать переменную Y по p-мерным ковариатам X=(X_1,… ,X_p) и вы сталкиваетесь с пропущенными значениями в X, есть интересное решение для древовидных методов. Этот метод на самом деле довольно старый, но, по-видимому, он замечательно работает в широком диапазоне наборов данных. Речь идет о отсутствующем критерии заложенных в признаках (МИА; [1]). Хотя есть много хороших статей об отсутствующих значениях (таких как эта), этот мощный подход кажется несколько недооцененным. В частности, вам не нужно каким-либо образом вводить, удалять или прогнозировать ваши пропущенные значения, а вместо этого можно просто запустить свой прогноз, как если бы у вас были полностью наблюдаемые данные.
Я быстро объясню, как работает сам метод, а затем приведу пример с распределенным случайным лесом (DRF), который объясняется здесь. Я выбрал DRF, потому что это очень общая версия Random Forest (в частности, его также можно использовать для предсказания случайного вектора Y) и потому, что я несколько предвзят. . MIA фактически реализован для обобщенного случайного леса (GRF), который охватывает широкий спектр реализаций леса. В частности, поскольку реализация DRF на CRAN основана на GRF, после небольшой модификации он может использовать и метод MIA.
Конечно, имейте в виду, что это быстрое решение, которое (насколько мне известно) не имеет теоретических гарантий. В зависимости от механизма пропуска, это может сильно исказить анализ. С другой стороны, наиболее часто используемые методы работы с отсутствующими значениями не имеют никаких теоретических гарантий или прямо известно, что они искажают анализ, и, по крайней мере, эмпирически, MIA работает хорошо и
Как это работает
Напомним, что в RF разбиения строятся в форме X_j ‹ Sили X_j ≥ S для измерения j=1,…,p. Чтобы найти это значение разделения S, он оптимизирует какой-либо критерий по Y, например критерий CART. Таким образом, наблюдения последовательно делятся с помощью правил принятия решений, зависящих от X.
В исходной статье это объясняется немного запутанно, но, насколько я понимаю, MIA работает следующим образом: рассмотрим образец (Y_1, X_1),… , (Y_n, X_n), с
X_i=(X_i1,…,X_ip)’.
Разделение без пропущенных значений — это просто поиск значения S, как указано выше, а затем выбрасывание всех Y_i с X_ij ‹ S в узле 1 и всех Y_i с X_ij ≥ S в узле 2. Вычисляя целевой критерий, такой как CART, для каждого значения S, мы можем выбрать лучший из них. С отсутствующими значениями вместо этого есть 3 варианта для каждого возможного значения разделения S для рассмотрения:
- Используйте обычное правило для всех наблюдений i, таких как X_ij, и отправьте i на узел 1, если X_ij равно отсутствующий.
- Используйте обычное правило для всех наблюдений i, таких как X_ij, и отправьте i на узел 2, если X_ij отсутствует.
- Игнорируйте обычное правило и просто отправьте i на узел 1, если X_ij отсутствует, и на узел 2, если он наблюдается.
Какому из этих правил следовать, снова решается в соответствии с используемым нами критерием Y_i.
Маленький пример
Здесь необходимо отметить, что пакет drf на CRAN еще не обновлен новейшей методологией. В будущем будет момент, когда все это будет реализовано в одном пакете на CRAN(!) Однако на данный момент существует две версии:
Если вы хотите использовать быструю реализацию drf с отсутствующими значениями (без доверительных интервалов), вы можете использовать функцию «drfown», прикрепленную в конце этой статьи. Этот код адаптирован из
lorismichel/drf: распределенные случайные леса (Cevid et al., 2020) (github.com)
С другой стороны, если вам нужны доверительные интервалы для ваших параметров, используйте этот (более медленный) код
drfinference/drf-foo.R на главной · JeffNaef/drfinference (github.com)
В частности, drf-foo.R содержит все необходимое в последнем случае.
Мы сосредоточимся на более медленном коде с доверительными интервалами, как описано в этой статье, а также рассмотрим тот же пример, что и в указанной статье:
set.seed(2) n<-2000 beta1<-1 beta2<--1.8 # Model Simulation X<-mvrnorm(n = n, mu=c(0,0), Sigma=matrix(c(1,0.7,0.7,1), nrow=2,ncol=2)) u<-rnorm(n=n, sd = sqrt(exp(X[,1]))) Y<- matrix(beta1*X[,1] + beta2*X[,2] + u, ncol=1)
Обратите внимание, что это гетероскедастическая линейная модель с p=2 и с дисперсией ошибки, зависящей от значений X_1. Теперь мы также добавляем отсутствующие значения в X_1 способом случайного пропущенного состояния (MAR):
prob_na <- 0.3 X[, 1] <- ifelse(X[, 2] <= -0.2 & runif(n) < prob_na, NA, X[, 1])
Это означает, что X_1 отсутствует с вероятностью 0,3, если X_2 имеет значение меньше -0,2. Таким образом, вероятность отсутствия X_1 зависит от X_2, что называется «случайным отсутствием». Это уже сложная ситуация, и есть информация, которую можно получить, изучив структуру пропущенных значений. То есть отсутствие не является «полностью случайным отсутствием (MCAR)», потому что отсутствие X_1 зависит от значения X_2. Это, в свою очередь, означает, что распределение X_2, из которого мы получаем данные, отличается в зависимости от того, отсутствует ли X_1 или нет. Это, в частности, означает, что удаление строк с отсутствующими значениями может серьезно повлиять на анализ.
Теперь мы фиксируем x и оцениваем условное математическое ожидание и дисперсию с учетом X=x,точно как в предыдущей статье.
# Choose an x that is not too far out x<-matrix(c(1,1),ncol=2) # Choose alpha for CIs alpha<-0.05
Затем мы также подгоняем DRF и прогнозируем веса для контрольной точки x (что соответствует прогнозированию условного распределения Y|X =x):
## Fit the new DRF framework drf_fit <- drfCI(X=X, Y=Y, min.node.size = 5, splitting.rule='FourierMMD', num.features=10, B=100) ## predict weights DRF = predictdrf(drf_fit, x=x) weights <- DRF$weights[1,]
Пример 1: условное ожидание
Сначала мы оцениваем условное математическое ожидание Y|X=x.
# Estimate the conditional expectation at x: condexpest<- sum(weights*Y) # Use the distribution of weights, see below distofcondexpest<-unlist(lapply(DRF$weightsb, function(wb) sum(wb[1,]*Y) )) # Can either use the above directly to build confidence interval, or can use the normal approximation. # We will use the latter varest<-var(distofcondexpest-condexpest) # build 95%-CI lower<-condexpest - qnorm(1-alpha/2)*sqrt(varest) upper<-condexpest + qnorm(1-alpha/2)*sqrt(varest) round(c(lower, condexpest, upper),2) # without NAs: (-1.00, -0.69 -0.37) # with NAs: (-1.15, -0.67, -0.19)
Примечательно, что значения, полученные с НС, очень близки к значениям из первого анализа без НС в предыдущей статье! Это действительно поразительно для меня, так как с этим отсутствующим механизмом нелегко иметь дело. Интересно, что расчетная дисперсия оценщика также удваивается: примерно с 0,025 без пропущенных значений до примерно 0,06 с пропущенными значениями.
Истина представлена как:
так что у нас есть небольшая ошибка, но доверительные интервалы содержат правду, как и должны.
Результат выглядит аналогично для более сложной цели, такой как условная дисперсия:
# Estimate the conditional expectation at x:
condvarest<- sum(weights*Y^2) - condexpest^2
distofcondvarest<-unlist(lapply(DRF$weightsb, function(wb) {
sum(wb[1,]*Y^2) - sum(wb[1,]*Y)^2
} ))
# Can either use the above directly to build confidence interval, or can use the normal approximation.
# We will use the latter
varest<-var(distofcondvarest-condvarest)
# build 95%-CI
lower<-condvarest - qnorm(1-alpha/2)*sqrt(varest)
upper<-condvarest + qnorm(1-alpha/2)*sqrt(varest)
c(lower, condvarest, upper)
# without NAs: (1.89, 2.65, 3.42)
# with NAs: (1.79, 2.74, 3.69)
Здесь разница в оценочных значениях несколько больше. Поскольку истина дается как
оценка с NA даже немного точнее (хотя, конечно, это, скорее всего, просто случайность). Снова оценка дисперсии оценщика (дисперсии) увеличивается с пропущенными значениями с 0,15 (отсутствие пропущенных значений) до 0,23.
Заключение
В этой статье мы обсудили MIA, который представляет собой адаптацию метода разбиения в Random Forest для работы с пропущенными значениями. Поскольку он реализован в GRF и DRF, его можно использовать широко, и рассмотренный нами небольшой пример показывает, что он работает замечательно хорошо.
Тем не менее, я хотел бы еще раз отметить, что нет никакой теоретической гарантии согласованности или того, что доверительные интервалы имеют смысл даже для очень большого количества точек данных. Причины отсутствия значений многочисленны, и нужно быть очень осторожным, чтобы не исказить анализ из-за небрежного обращения с этой проблемой. Метод MIA ни в коем случае не является хорошо понятным решением этой проблемы. Тем не менее, на данный момент это кажется разумным быстрым решением, которое, по-видимому, может в некоторой степени использовать шаблон отсутствующих данных. Если кто-то делает/имеет более обширный анализ моделирования, мне было бы любопытно узнать о результатах.
Код
require(drf)
drfown <- function(X, Y,
num.trees = 500,
splitting.rule = "FourierMMD",
num.features = 10,
bandwidth = NULL,
response.scaling = TRUE,
node.scaling = FALSE,
sample.weights = NULL,
sample.fraction = 0.5,
mtry = min(ceiling(sqrt(ncol(X)) + 20), ncol(X)),
min.node.size = 15,
honesty = TRUE,
honesty.fraction = 0.5,
honesty.prune.leaves = TRUE,
alpha = 0.05,
imbalance.penalty = 0,
compute.oob.predictions = TRUE,
num.threads = NULL,
seed = stats::runif(1, 0, .Machine$integer.max),
compute.variable.importance = FALSE) {
# initial checks for X and Y
if (is.data.frame(X)) {
if (is.null(names(X))) {
stop("the regressor should be named if provided under data.frame format.")
}
if (any(apply(X, 2, class) %in% c("factor", "character"))) {
any.factor.or.character <- TRUE
X.mat <- as.matrix(fastDummies::dummy_cols(X, remove_selected_columns = TRUE))
} else {
any.factor.or.character <- FALSE
X.mat <- as.matrix(X)
}
mat.col.names.df <- names(X)
mat.col.names <- colnames(X.mat)
} else {
X.mat <- X
mat.col.names <- NULL
mat.col.names.df <- NULL
any.factor.or.character <- FALSE
}
if (is.data.frame(Y)) {
if (any(apply(Y, 2, class) %in% c("factor", "character"))) {
stop("Y should only contain numeric variables.")
}
Y <- as.matrix(Y)
}
if (is.vector(Y)) {
Y <- matrix(Y,ncol=1)
}
#validate_X(X.mat)
if (inherits(X, "Matrix") && !(inherits(X, "dgCMatrix"))) {
stop("Currently only sparse data of class 'dgCMatrix' is supported.")
}
drf:::validate_sample_weights(sample.weights, X.mat)
#Y <- validate_observations(Y, X)
# set legacy GRF parameters
clusters <- vector(mode = "numeric", length = 0)
samples.per.cluster <- 0
equalize.cluster.weights <- FALSE
ci.group.size <- 1
num.threads <- drf:::validate_num_threads(num.threads)
all.tunable.params <- c("sample.fraction", "mtry", "min.node.size", "honesty.fraction",
"honesty.prune.leaves", "alpha", "imbalance.penalty")
# should we scale or not the data
if (response.scaling) {
Y.transformed <- scale(Y)
} else {
Y.transformed <- Y
}
data <- drf:::create_data_matrices(X.mat, outcome = Y.transformed, sample.weights = sample.weights)
# bandwidth using median heuristic by default
if (is.null(bandwidth)) {
bandwidth <- drf:::medianHeuristic(Y.transformed)
}
args <- list(num.trees = num.trees,
clusters = clusters,
samples.per.cluster = samples.per.cluster,
sample.fraction = sample.fraction,
mtry = mtry,
min.node.size = min.node.size,
honesty = honesty,
honesty.fraction = honesty.fraction,
honesty.prune.leaves = honesty.prune.leaves,
alpha = alpha,
imbalance.penalty = imbalance.penalty,
ci.group.size = ci.group.size,
compute.oob.predictions = compute.oob.predictions,
num.threads = num.threads,
seed = seed,
num_features = num.features,
bandwidth = bandwidth,
node_scaling = ifelse(node.scaling, 1, 0))
if (splitting.rule == "CART") {
##forest <- do.call(gini_train, c(data, args))
forest <- drf:::do.call.rcpp(drf:::gini_train, c(data, args))
##forest <- do.call(gini_train, c(data, args))
} else if (splitting.rule == "FourierMMD") {
forest <- drf:::do.call.rcpp(drf:::fourier_train, c(data, args))
} else {
stop("splitting rule not available.")
}
class(forest) <- c("drf")
forest[["ci.group.size"]] <- ci.group.size
forest[["X.orig"]] <- X.mat
forest[["is.df.X"]] <- is.data.frame(X)
forest[["Y.orig"]] <- Y
forest[["sample.weights"]] <- sample.weights
forest[["clusters"]] <- clusters
forest[["equalize.cluster.weights"]] <- equalize.cluster.weights
forest[["tunable.params"]] <- args[all.tunable.params]
forest[["mat.col.names"]] <- mat.col.names
forest[["mat.col.names.df"]] <- mat.col.names.df
forest[["any.factor.or.character"]] <- any.factor.or.character
if (compute.variable.importance) {
forest[['variable.importance']] <- variableImportance(forest, h = bandwidth)
}
forest
}
Цитаты
[1] Twala, B.E.T.H., M.C. Jones и David J. Hand. Хорошие методы работы с отсутствующими данными в деревьях решений. Письма о распознавании образов 29, 2008 г.
