8

我目前正在使用该rpart软件包将回归树拟合到具有相对较少的观察和数千个具有两个可能值的分类预测变量的数据。

通过在较小的数据上测试包，我知道在这种情况下，我是否将回归变量声明为分类（即因子）或保持原样（它们被编码为 +/-1）并不重要。

但是，我仍然想了解为什么将我的解释变量作为因素传递会显着减慢算法的速度（尤其是因为我很快就会得到新数据，其中响应需要 3 个不同的值并将它们视为连续的将不再是一种选择）。当然应该反过来吗？

这是一个模拟我的数据的示例代码：

library(rpart)

x <- as.data.frame(matrix(sample(c(-1, +1), 50 * 3000, replace = T), nrow = 50))
y <- rnorm(50)

x.fac <- as.data.frame(lapply(x, factor))

现在比较：

system.time(rpart( y ~ ., data = x, method = 'anova'))

   user  system elapsed 
   1.62    0.21    1.85 

system.time(rpart( y ~ ., data = x.fac, method = 'anova'))

   user  system elapsed 
   246.87  165.91  412.92

处理每个变量（因子）只有一个潜在的分裂可能性比处理整个范围的潜在分裂（对于连续变量）更简单、更快捷，所以我对这种rpart行为最感到困惑。任何澄清/建议将不胜感激。

4

2 回答 2

6

只是以@gavin simpson 的上述发现为基础......我决定修改一下rpart.matrix，看看我是否可以对过多的执行时间做点什么。

问题归结为for循环的使用。通常，与;for相比，我是不可知论者。[sl]apply后者通常被认为更优雅，但我不会for在它工作正常时替换它，只是为了这个。特别是，我认为*apply有时会夸大其性能优势；for与旧的 S-Plus 时代相比，在速度和内存使用方面有了显着提高。

但在这种情况下不是。只需将 替换for为即可lapply将此示例的运行时间缩短 2 个数量级。看看其他人是否可以证实这一点会很好。

m <- model.frame(x.fac)


# call rpart.matrix
system.time(mm <- rpart:::rpart.matrix(m))
   user  system elapsed 
 208.25   88.03  296.99 


# exactly the same as rpart.matrix, but with for replaced by lapply
f <- function(frame)
{
    if (!inherits(frame, "data.frame") || is.null(attr(frame, 
        "terms"))) 
        return(as.matrix(frame))
    frame[] <- lapply(frame, function(x) {
        if (is.character(x))
            as.numeric(factor(x))
        else if(!is.numeric(x))
            as.numeric(x)
        else x
    })
    X <- model.matrix(attr(frame, "terms"), frame)[, -1L, drop = FALSE]
    colnames(X) <- sub("^`(.*)`", "\\1", colnames(X))
    class(X) <- c("rpart.matrix", class(X))
    X
}

system.time(mm2 <- f(m))
   user  system elapsed 
   0.65    0.04    0.70 


identical(mm, mm2)
[1] TRUE
于 2013-06-20T09:59:37.933 回答
6

您需要对代码进行分析以确定，但如果时间差异不是来自 R 在准备模型矩阵时必须将每个因子变量转换为两个二进制变量，我会感到惊讶。

尝试

Rprof("rpartProfile.Rprof")
rpart( y ~ ., data = x.fac, method = 'anova')
Rprof()

summaryRprof("rpartProfile.Rprof")

看看时间都花在了哪里。我现在已经完成了：

> summaryRprof("rpartProfile.Rprof")
$by.self
                          self.time self.pct total.time total.pct
"[[<-.data.frame"            786.46    72.45     786.56     72.46
"rpart.matrix"               294.26    27.11    1081.78     99.66
"model.frame.default"          1.04     0.10       3.00      0.28
"terms.formula"                0.96     0.09       0.96      0.09
"as.list.data.frame"           0.46     0.04       0.46      0.04
"makepredictcall.default"      0.46     0.04       0.46      0.04
"rpart"                        0.44     0.04    1085.38     99.99
"[[.data.frame"                0.16     0.01       0.42      0.04
"<Anonymous>"                  0.16     0.01       0.18      0.02
"match"                        0.14     0.01       0.22      0.02
"print"                        0.12     0.01       0.12      0.01
"model.matrix.default"         0.10     0.01       0.44      0.04
....

$by.total
                          total.time total.pct self.time self.pct
"rpart"                      1085.38     99.99      0.44     0.04
"rpart.matrix"               1081.78     99.66    294.26    27.11
"[[<-"                        786.62     72.47      0.06     0.01
"[[<-.data.frame"             786.56     72.46    786.46    72.45
"model.frame.default"           3.00      0.28      1.04     0.10
"eval"                          3.00      0.28      0.04     0.00
"eval.parent"                   3.00      0.28      0.00     0.00
"model.frame"                   3.00      0.28      0.00     0.00
"terms.formula"                 0.96      0.09      0.96     0.09
"terms"                         0.96      0.09      0.00     0.00
"makepredictcall"               0.50      0.05      0.04     0.00
"as.list.data.frame"            0.46      0.04      0.46     0.04
"makepredictcall.default"       0.46      0.04      0.46     0.04
"as.list"                       0.46      0.04      0.00     0.00
"vapply"                        0.46      0.04      0.00     0.00
"model.matrix.default"          0.44      0.04      0.10     0.01
"[["                            0.44      0.04      0.02     0.00
"model.matrix"                  0.44      0.04      0.00     0.00
....

$sample.interval
[1] 0.02

$sampling.time
[1] 1085.5

从上面注意到，大部分时间都花在了 function 上rpart.matrix：

> rpart:::rpart.matrix
function (frame) 
{
    if (!inherits(frame, "data.frame") || is.null(attr(frame, 
        "terms"))) 
        return(as.matrix(frame))
    for (i in 1:ncol(frame)) {
        if (is.character(frame[[i]])) 
            frame[[i]] <- as.numeric(factor(frame[[i]]))
        else if (!is.numeric(frame[[i]])) 
            frame[[i]] <- as.numeric(frame[[i]])
    }
    X <- model.matrix(attr(frame, "terms"), frame)[, -1L, drop = FALSE]
    colnames(X) <- sub("^`(.*)`", "\\1", colnames(X))
    class(X) <- c("rpart.matrix", class(X))
    X
}

但for大部分时间都花在了该函数中的循环上，本质上是转换每一列并将它们添加回数据框中。

于 2013-06-19T15:37:25.540 回答