r - 用回归树进行亚组识别

Question

我对最近一篇关于“生物统计学教程:临床试验中数据驱动的亚组识别和分析”(DOI:10.1002/sim.7064)的论文非常感兴趣，我想在“Performance of基于树的回归方法”。然而，分区树似乎没有得到我想要的结果。

set.seed(123)
n <- 1000
x1 <- rnorm(n)
x2 <- runif(n)
t <- rbinom(n,1,0.5)
x3 <- rbinom(n,1,0.3)
x4 <- rnorm(n)
z <-1+ 2*x1-1.8*x2+(x1>=-0.3)*(x2>=0.4)*t-(x1< -0.3)*(x2<0.4)*t
pr = 1/(1+exp(-z))
y = rbinom(n,1,pr)
dt <- data.frame(x1,x2,x3,x4,t,y)

library(party)
mb <- mob(y~t-1|x1+x2+x3+x4,
  data=dt,
  model = glinearModel,
  family = binomial(),
  control=mob_control(minsplit=100))
plot(mb)

如上图所示。它应该在模拟中定义的截止值 -0.3 和 0.4 处在 x1 和 x2 上 split 。然而，它似乎并没有这样做。x1 主导了节点分区，而 x2 似乎不是分区过程的重要决定因素。代码有什么问题？

Answer 1

您指定的基于模型的 GLM 树尝试拟合以下模型:logit(prob) = alpha_tree(x1-x4) * t，其中 alpha_tree(x1-x4) 是处理指标 t 的树(基于分区变量 x1-x4)的子组特定系数。因此，该模型不包括 x1 和 x2 的截距或全局主效应的可能性 - 正如您的数据中模拟的那样。由于这些主要影响非常明显，因此该模型除了通过进一步拆分来捕获这些影响之外别无他法。因此，您的示例中的大树。

在经典的 MOB 框架中(Zeileis et al. 2008，Journal of Computational and Graphical Statistics，doi:10.1198/106186008X319331)，唯一的选择是将每个相关回归变量包含在要划分的模型中，即 logit(mu) = beta0_tree(x1-x4) + beta1_tree(x1_x4) * x1 + beta2_tree(x1-x4) * x2 + alpha_tree(x1-x4) * t。这行得通并且也将仅针对处理效果 alpha * t 检测子组，但是当然会失去一些效率，因为它也会重新估计每个子组中的 beta 系数。 Seibold 等人 (2016a), The International Journal of Biostatistics, doi:10.1515/ijb-2015-0032 中讨论了这种专门用于亚组分析的方法。 .

最近，我们建议对 MOB 进行改编，我们将其称为 PALM 树，用于部分加性线性模型树(Seibold et al. 2016b，http://arxiv.org/abs/1612.07498)。这允许拟合类型为 logit(mu) = beta0 + beta1 * x1 + beta2 * x2 + alpha_tree(x1-x4) * t 的模型，正如您在问题中模拟的那样。

经典的基于 GLM 的 MOB 树和 PALM 树的实现分别作为 glmtree() 和 palmtree() 在 partykit 推荐使用旧的 party 实现。将这些应用于上面的模拟数据，会产生以下结果。首先，重要的是所有分类分区变量也被标记为 factor 变量(以便选择正确的参数不稳定性测试):

dt <- transform(dt, x3 = factor(x3))

然后，我们可以仅使用特定于子组的处理效果、x1 和 x2 的全局主效果以及基于 x1、x2、x3、x4 的分区来拟合您模拟的模型。

library("partykit")
tr1 <- palmtree(y ~ t - 1 | x1 + x2 | x1 + x2 + x3 + x4, data = dt,
  family = binomial, minsize = 100)
tr1
## Partially additive generalized linear model tree (family: binomial)
## 
## Model formula:
## y ~ t - 1 | x1 + x2 + x3 + x4
## 
## Fitted party:
## [1] root
## |   [2] x1 <= -0.21797: n = 399
## |                t 
## |       -0.9245345 
## |   [3] x1 > -0.21797: n = 601
## |               t 
## |       0.6033979 
## 
## Number of inner nodes:    1
## Number of terminal nodes: 2
## Number of parameters per node: 1
## Objective function (negative log-likelihood): 432.5717
## 
## Linear fixed effects (from palm model):
## (Intercept)          x1          x2 
##   0.7140397   1.7589675  -1.1335779 

因此，这涵盖了数据生成过程中最重要的部分，但未能检测到与 x2 的交互。

plot(tr1)

我试过设置其他种子或使用基于 BIC 的后修剪(结合较大的显着性水平)，有时也可以发现与 x2 的交互。据推测，更大的样本也会更频繁地产生“真实”树。因此，该模型原则上似乎能够适合您模拟的模型，只是不适合这种特定设置。

就个人而言，我总是让截距和治疗效果因亚组而异。因为如果有任何我忽略的主要影响，这可能会产生更好的模型。如果每个子组中都包含一个截距，那么最好将 y 和 t 都编码为因子，从而在可视化中产生更好的图:

dt <- transform(dt, y = factor(y), t = factor(t))
tr2 <- palmtree(y ~ t | x1 + x2 | x1 + x2 + x3 + x4, data = dt,
  family = binomial, minsize = 100)
tr2
## Partially additive generalized linear model tree (family: binomial)
## 
## Model formula:
## y ~ t | x1 + x2 + x3 + x4
## 
## Fitted party:
## [1] root
## |   [2] x1 <= -0.26515: n = 382
## |       (Intercept)          t1 
## |         0.9839353  -1.1376981 
## |   [3] x1 > -0.26515: n = 618
## |       (Intercept)          t1 
## |         0.5331386   0.6566111 
## 
## Number of inner nodes:    1
## Number of terminal nodes: 2
## Number of parameters per node: 2
## Objective function (negative log-likelihood): 432.3303
## 
## Linear fixed effects (from palm model):
##        x1        x2 
##  1.964397 -1.078958 

对于此数据，这与上面的模型几乎吻合。但是显示更容易阅读，显示了两组之间绝对治疗效果的巨大差异:

plot(tr2)

最后，如果使用不可能包含附加主效应的普通旧 MOB，我们应该在每个子组中包含回归变量 x1 和 x2:

tr3 <- glmtree(y ~ t + x1 + x2 | x1 + x2 + x3 + x4, data = dt,
  family = binomial, minsize = 100)
tr3
## Generalized linear model tree (family: binomial)
## 
## Model formula:
## y ~ t + x1 + x2 | x1 + x2 + x3 + x4
## 
## Fitted party:
## [1] root
## |   [2] x1 <= -0.26515: n = 382
## |       (Intercept)          t1          x1          x2 
## |         0.5570219  -1.0511317   1.2533975  -1.3899679 
## |   [3] x1 > -0.26515: n = 618
## |       (Intercept)          t1          x1          x2 
## |         0.3573041   0.6943206   2.2910053  -0.9570403 
## 
## Number of inner nodes:    1
## Number of terminal nodes: 2
## Number of parameters per node: 4
## Objective function (negative log-likelihood): 429.2406

同样，这基本上找到了与以前相同的子组。但是，它会损失一点效率，因为必须在每个子组中估计更多的回归系数，而只有 t 系数在子组之间真正发生变化。

plot(tr3, tp_args = list(which = 1))