【问题标题】:How to get terminal nodes for a new observation from an rpart object?如何从 rpart 对象获取新观察的终端节点?
【发布时间】:2015-05-31 23:56:28
【问题描述】:

说我有

head(kyphosis)
inTrain <- sample(1:nrow(kyphosis), 45, replace = F)
TRAIN_KYPHOSIS <- kyphosis[inTrain,]
TEST_KYPHOSIS <- kyphosis[-inTrain,]

(kyph_tree <- rpart(Number ~ ., data = TRAIN_KYPHOSIS))

对于TEST_KYPHOSIS中的每个观察,如何从拟合对象中获取终端节点?

如何从每个测试观察映射到的终端节点获取摘要,例如偏差和预测值?

【问题讨论】:

    标签: r rpart


    【解决方案1】:

    rpart 实际上有这个功能,但它没有暴露(奇怪的是,这是一个相当明显的要求)。

    predict_nodes <-
        function (object, newdata, na.action = na.pass) {
            where <-
                if (missing(newdata)) 
                    object$where
                else {
                    if (is.null(attr(newdata, "terms"))) {
                        Terms <- delete.response(object$terms)
                        newdata <- model.frame(Terms, newdata, na.action = na.action, 
                                               xlev = attr(object, "xlevels"))
                        if (!is.null(cl <- attr(Terms, "dataClasses"))) 
                            .checkMFClasses(cl, newdata, TRUE)
                    }
                    rpart:::pred.rpart(object, rpart:::rpart.matrix(newdata))
                }
            as.integer(row.names(object$frame))[where]
        }
    

    然后:

    > predict_nodes(kyph_tree, TEST_KYPHOSIS)
     [1] 5 3 4 3 3 5 5 3 3 3 3 5 5 4 3 5 4 3 3 3 3 4 3 4 4 5 5 3 4 4 3 5 3 5 5 5
    

    【讨论】:

    • 为什么rpart:::pred.rpart(object, rpart:::rpart.matrix(newdata))会导致预测节点?
    • @goldisfine,因为这是 rpart 在内部计算预测节点的方式。此功能在内部使用,但未公开。
    • @VitoshKa 感谢发布解决方案。这是树的基本部分!没有这部分几乎无法使用。
    【解决方案2】:

    一种选择是将rpart 对象转换为party 包中party 类的对象。这提供了一个用于处理递归聚会的通用工具包。转换很简单:

    library("partykit")
    (kyph_party <- as.party(kyph_tree))
    
    Model formula:
    Number ~ Kyphosis + Age + Start
    
    Fitted party:
    [1] root
    |   [2] Start >= 15.5: 2.933 (n = 15, err = 10.9)
    |   [3] Start < 15.5
    |   |   [4] Age >= 112.5: 3.714 (n = 14, err = 18.9)
    |   |   [5] Age < 112.5: 5.125 (n = 16, err = 29.8)
    
    Number of inner nodes:    2
    Number of terminal nodes: 3
    

    (为了获得精确的可重复性,请在运行我的代码之前使用 set.seed(1) 运行您的问题中的代码。)

    对于此类的对象,plot()predict()fitted() 等有一些更灵活的方法。例如,plot(kyph_party) 比默认的plot(kyph_tree) 产生更多信息显示。 fitted() 方法提取两列 data.frame,其中包含拟合的节点编号和在训练数据上观察到的响应。

    kyph_fit <- fitted(kyph_party)
    head(kyph_fit, 3)
    
      (fitted) (response)
    1        5          6
    2        2          2
    3        4          3
    

    通过它,您可以轻松计算您感兴趣的任何数量,例如,每个节点内的均值、中位数或残差平方和。

    tapply(kyph_fit[,2], kyph_fit[,1], mean)
    
           2        4        5 
    2.933333 3.714286 5.125000 
    
    tapply(kyph_fit[,2], kyph_fit[,1], median)
    
    2 4 5 
    3 4 5 
    
    tapply(kyph_fit[,2], kyph_fit[,1], function(x) sum((x - mean(x))^2))
    
           2        4        5 
    10.93333 18.85714 29.75000 
    

    除了简单的tapply(),您还可以使用您选择的任何其他函数来计算分组统计表。

    现在要从测试数据TEST_KYPHOSIS 中了解哪个观察到树中的哪个节点,您可以简单地使用predict(..., type = "node") 方法:

    kyph_pred <- predict(kyph_party, newdata = TEST_KYPHOSIS, type = "node")
    head(kyph_pred)
    
     2  3  4  6  7 10 
     4  4  5  2  2  5 
    

    【讨论】:

    • 您的解决方案产生的结果与kyph_tree$where 相同,但与下面的 VitoshKa 解决方案获得的结果不同。
    • VitoshKa 的predict_nodes() 解决方案和partykit 中的predict(..., type = "node") 解决方案不会产生完全相同的标签,因为节点ID 的分配方式略有不同。但信息实际上是等价的。退房:table(predict_nodes(kyph_tree, TEST_KYPHOSIS), predict(kyph_party, newdata = TEST_KYPHOSIS, type = "node"))。由于标签不同,它可能无法诊断,但存在 1:1 匹配。但这是因为partykit 为递归分区提供了通用解决方案,而不是特定于rpart
    • 感谢您的回复。我很难理解kyph_tree$where 返回的内容。它看起来不像返回终端节点标签
    • 其实就是训练数据上的终端节点标签/IDTRAIN_KYPHOSIS,比如看看table(kyph_tree$where)。如果您将其与 table(fitted(kyph_party)[,1])table(predict(kyph_party, type = "node")) 进行比较,您会看到相同的绝对频率,但可能会看到不同的标签(取决于树的结构)。
    • 好点。 $where 中的 ID 与 $frame 中的行相关,即,与打印/绘图中使用的标签不同,但对应于汇总表中的 ID。如果您使用该汇总表中的行名,您将获得与打印/绘图相同的标签:table(as.numeric(rownames(kyph_tree$frame))[kyph_tree$where]).
    猜你喜欢
    • 2016-08-13
    • 2013-01-18
    • 2017-08-09
    • 2022-11-12
    • 2016-06-10
    • 1970-01-01
    • 1970-01-01
    • 2017-06-03
    相关资源
    最近更新 更多