動動發财的小手，點個贊吧！
除了在網上找到的一些過度清理的資料集之外，缺失值無處不在。事實上，資料集越複雜、越大，出現缺失值的可能性就越大。缺失值是統計研究的一個令人着迷的領域，但在實踐中它們往往很麻煩。
如果您處理一個預測問題，想要從 p 維協變量 X=(X_1,…,X_p) 預測變量 Y，并且面臨 X 中的缺失值，那麼基于樹的方法有一個有趣的解決方案。這種方法實際上相當古老，但在各種資料集中似乎都表現得非常好。我說的是“缺失的屬性标準”（MIA；[1]）。雖然有很多關于缺失值的好文章（例如這篇文章），但這種強大的方法似乎有些未得到充分利用。特别是，不需要以任何方式插補、删除或預測缺失值，而是可以像完全觀察到的資料一樣運作預測。
我将快速解釋該方法本身是如何工作的，然後提供一個示例以及此處解釋的分布式随機森林 (DRF)。我選擇 DRF 是因為它是随機森林的一個非常通用的版本（特别是，它也可以用來預測随機向量 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)，
不缺失值的分割就是像上面那樣尋找值S，然後将節點1中所有X_ij < S的Y_i和節點2中所有X_ij ≥ S的Y_i扔進去。計算每個值S的目标标準，例如CART，我們可以選擇最好的一個。對于缺失值，每個候選分割值 S 有 3 個選項需要考慮：
- 對所有觀測值 i 使用通常的規則，使得 X_ij 被觀測到，如果 X_ij 丢失，則将 i 發送到節點 1。
- 對所有觀測值 i 使用通常的規則，以便觀測到 X_ij，如果缺少 X_ij，則将 i 發送到節點 2。
- 忽略通常的規則，如果 X_ij 缺失，則将 i 發送到節點 1；如果觀察到 X_ij，則将 i 發送到節點 2。
遵循這些規則中的哪一個再次根據我們使用的 Y_i 的标準來決定。
例子
需要指出的是，CRAN 上的 drf 包尚未使用最新的方法進行更新。将來有一天，所有這些都将在 CRAN 上的一個包中實作。但是，目前有兩個版本：
如果您想使用缺失值（無置信區間）的快速 drf 實作，您可以使用本文末尾附帶的“drfown”函數。此代碼改編自
lorismichel/drf: Distributional Random Forests (Cevid et al., 2020) (github.com)
另一方面，如果您想要參數的置信區間，請使用此（較慢的）代碼
drfinference/drf-foo.R at main · 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 值。現在我們還以随機缺失 (MAR) 方式向 X_1 添加缺失值：
prob_na <- 0.3
X[, 1] <- ifelse(X[, 2] <= -0.2 & runif(n) < prob_na, NA, X[, 1])
這意味着每當 X_2 的值小于 -0.2 時，X_1 缺失的機率為 0.3。是以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,]
條件期望
我們首先估計 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)
值得注意的是，使用 NA 獲得的值與上一篇文章中未使用 NA 的第一次分析得到的值非常接近！這确實令我震驚，因為這個缺失的機制并不容易處理。有趣的是，估計器的估計方差也翻倍，從沒有缺失值的大約 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。
結論
在本文[1]中，我們讨論了 MIA，它是随機森林中分裂方法的一種改進，用于處理缺失值。由于它是在 GRF 和 DRF 中實作的，是以它可以被廣泛使用，我們看到的小例子表明它工作得非常好。
然而，我想再次指出，即使對于大量資料點，也沒有一緻性或置信區間有意義的理論保證。缺失值的原因有很多，必須非常小心，不要因粗心處理這一問題而使分析産生偏差。 MIA 方法對于這個問題來說決不是一個很好了解的解決方案。然而，目前這似乎是一個合理的快速解決方案，它似乎能夠利用資料缺失的模式。如果有人進行了更廣泛的模拟分析，我會對結果感到好奇。
Code
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
}
Reference
[1]Source: https://towardsdatascience.com/random-forests-and-missing-values-3daaea103db0