r - 使用partykit自定义拆分规则

Question

这篇文章遵循这个问题:https://stackoverflow.com/questions/31234329/rpart-user-defined-implementation

我对可以使用自定义标准处理树木生长的工具非常感兴趣，这样我就可以测试不同的模型。

我尝试使用 partykit R 包来生长一棵树，其拆分规则由 Cox 模型的负对数似然给出(在 Cox 模型的情况下为对数准似然)并拟合了 Cox 模型在每一片叶子里。

正如我在阅读有关 MOB 函数的小插图时所理解的，有两种方法可以实现我自己的拆分标准，即让 fit 函数返回列表或模型对象。

出于我的目的，我尝试了这两种解决方案，但未能使其发挥作用。

解决方案 1:返回一个列表对象:

我以“mob”小插图中的“乳腺癌数据集”为例。

我试过这个:

cox1 = function(y,x, start = NULL, weights = NULL, offset = NULL, ...,
           estfun = FALSE, object = TRUE){
  res_cox = coxph(formula = y ~ x )
list(
  coefficients = res_cox$coefficients,
  objfun = - res_cox$loglik[2],
  object = res_cox)
}


mob(formula = Surv(time, cens) ~ horTh + pnodes - 1 | age + tsize + tgrade + progrec +
  estrec + menostat , 
    data = GBSG2 ,
    fit = cox1,
    control = mob_control(alpha = 0.0001) )

有关于 X 矩阵奇异性的警告，并且 mob 函数是具有单个节点的树(即使 alpha 值较小)。

请注意，运行 coxph 函数时 X 矩阵没有奇异性问题:

res_cox = coxph( formula = Surv(time, cens) ~ horTh + pnodes  ,
             data = GBSG2 )

解决方案 2:返回 coxph.object :

我试过这个:

cox2 = function(y,x, start = NULL, weights = NULL, offset = NULL, ... ){
  res_cox = coxph(formula = y ~ x )
}

logLik.cox2 <- function(object, ...)
  structure( - object$loglik[2], class = "logLik")

mob(formula = Surv(time, cens) ~ horTh + pnodes - 1 | age + tsize + tgrade + progrec +
  estrec + menostat , 
    data = GBSG2 ,
    fit = cox2,
    control = mob_control(alpha = 0.0001 ) )

所以这次我沿着“progrec”变量进行了拆分:

Model-based recursive partitioning (cox2)

Model formula:
Surv(time, cens) ~ horTh + pnodes - 1 | age + tsize + tgrade + 
progrec + estrec + menostat

Fitted party:
[1] root
|   [2] progrec <= 21: n = 281
|         xhorThno  xhorThyes    xpnodes 
|       0.19306661         NA 0.07832756 
|   [3] progrec > 21: n = 405
|         xhorThno  xhorThyes    xpnodes 
|       0.64810352         NA 0.04482348 

Number of inner nodes:    1
Number of terminal nodes: 2
Number of parameters per node: 3
Objective function: 1531.132
Warning message:
In coxph(formula = y ~ x) : X matrix deemed to be singular; variable 2

我想知道我的解决方案 1 有什么问题。

我也为回归问题尝试了类似的方法并得到相同的结果，以单个叶子结尾:

data("BostonHousing", package = "mlbench")

BostonHousing <- transform(BostonHousing,
                       chas = factor(chas, levels = 0:1, labels = c("no", "yes")),
                       rad = factor(rad, ordered = TRUE))


linear_reg = function(y,x, start = NULL, weights = NULL, offset = NULL, ...,
                  estfun = FALSE, object = TRUE){
  res_lm = glm(formula = y ~ x , family = "gaussian")
  list(
    coefficients = res_lm$coefficients,
    objfun = res_lm$deviance,
    object = res_lm )
}

mob( formula = medv ~ log(lstat) + I(rm^2) | zn + indus + chas + nox +
   + age + dis + rad + tax + crim + b + ptratio, 
     data = BostonHousing ,
     fit = linear_reg)

另外，我想知道将变量用于“在节点中拟合模型”和“进行拆分”是否没有问题。

先感谢您。

我可能还有其他关于partykit 功能的问题。

Answer 1

cox1() 的问题和 linear_reg()您设置的功能是您不提供估计功能又名分数贡献。由于这些是选择 split 变量的推理的基础，如果没有提供这些，算法根本不会 split 。最近看到这个 answer对于这个问题的一些讨论。

但是对于 coxph()对象(与上面链接的讨论中的 fitdistr() 示例不同)很容易获得这些估计函数或分数，因为有一个 estfun()可用的方法。所以你的 cox2()方法是去这里更容易的路线。

后者不能正常工作的原因是由于 coxph() 中对拦截的特殊处理所致。 .在内部，这总是强制截距进入模型，但随后会从设计矩阵中省略第一列。通过 mob() 连接时你需要小心不要搞砸，因为 mob()建立自己的模型矩阵。因为你排除了拦截，mob()认为它可以估计horTh的两个水平.但事实并非如此，因为在 Cox-PH 模型中未识别截距。

在这种情况下 (IMO) 的最佳解决方案如下:您让 mob()设置一个截距，然后在将模型矩阵传递给 coxph() 时再次排除它.因为有coef() , logLik() , 和 estfun()结果对象的方法，您可以使用 cox2() 的简单设置功能。

包和数据:

library("partykit")
library("survival")
data("GBSG2", package = "TH.data")

拟合功能:

cox <- function(y, x, start = NULL, weights = NULL, offset = NULL, ... ) {
  x <- x[, -1]
  coxph(formula = y ~ 0 + x)
}

将 MOB 树拟合到 GBSG2数据:

mb <- mob(formula = Surv(time, cens) ~ horTh + pnodes | age + tsize + tgrade + progrec + estrec + menostat, 
  data = GBSG2, fit = cox)
mb
## Model-based recursive partitioning (cox)
## 
## Model formula:
## Surv(time, cens) ~ horTh + pnodes | age + tsize + tgrade + progrec + 
##     estrec + menostat
## 
## Fitted party:
## [1] root: n = 686
##       xhorThyes     xpnodes 
##     -0.35701115  0.05768026  
## 
## Number of inner nodes:    0
## Number of terminal nodes: 1
## Number of parameters per node: 2
## Objective function: 1758.86