动动发财的小手,点个赞吧!

除了在网上找到的一些适度清理的数据集之外,缺失值无处不在。事实上,数据集越简单、越大,呈现缺失值的可能性就越大。缺失值是统计钻研的一个令人着迷的畛域,但在实践中它们往往很麻烦。

如果您解决一个预测问题,想要从 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<-2000beta1<-1beta2<--1.8# Model SimulationX<-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.3X[, 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 outx<-matrix(c(1,1),ncol=2)# Choose alpha for CIsalpha<-0.05

而后咱们还拟合 DRF 并预测测试点 x 的权重(对应于预测 Y|X=x 的条件散布):

## Fit the new DRF frameworkdrf_fit <- drfCI(X=X, Y=Y, min.node.size = 5, splitting.rule='FourierMMD', num.features=10, B=100)## predict weightsDRF = 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 belowdistofcondexpest<-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 lattervarest<-var(distofcondexpest-condexpest)# build 95%-CIlower<-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^2distofcondvarest<-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 lattervarest<-var(distofcondvarest-condvarest)# build 95%-CIlower<-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,它是随机森林中决裂办法的一种改良,用于解决缺失值。因为它是在 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}

本文由mdnice多平台公布