这篇文章遵循这个问题: https://stackoverflow.com/questions/31234329/rpart-user-defined-implementation
我对可以使用自定义标准处理树木生长的工具非常感兴趣,这样我就可以测试不同的模型。
我尝试使用 partykit R 包来生长一棵树,其拆分规则由 Cox 模型的负对数似然给出(在 Cox 模型的情况下是对数准似然),并且每个叶子中都拟合了一个 Cox 模型。
正如我理解阅读有关 MOB 函数的小插图一样,有两种方法可以实现我自己的拆分条件,即让 fit 函数返回列表或模型对象。
出于我的目的,我尝试了两种解决方案,但未能使其正常工作。
解决方案 1 : 返回一个列表对象:
我以"暴民"小插曲中的"乳腺癌数据集"为例。
我试过这个:
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矩阵奇点的警告,并且暴民函数具有单个节点的树(即使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)
我还想知道使用变量来"在节点中拟合模型"和"进行拆分"是否没有问题。
提前谢谢你。
我可能会对派对套件功能有其他问题。
您设置的cox1()
和linear_reg()
函数的问题在于您没有提供估计函数(即分数贡献)。由于这些是选择拆分变量的推理的基础,因此如果未提供这些,则算法根本不会拆分。有关此问题的一些讨论,请参阅最近的答案。
但是对于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