用回归树进行亚组识别
subgroup identification with regression tree
最近一篇关于"Tutorial in biostatistics: data-driven subgroup identification and analysis in clinical trials "的论文（DOI: 10.1002/sim.7064）很感兴趣，想在"Performance of the tree-based regression approach"部分重现结果。但是，分区树好像没有得到我想要的结果。
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 上分裂。然而，它似乎并没有这样做。x1 主导了节点分区，而 x2 似乎不是分区过程的重要决定因素。代码有什么问题？
您指定的基于模型的 GLM 树尝试拟合以下模型：
logit(prob) = alpha_tree(x1-x4) * t，其中 alpha_tree(x1-x4) 是用于处理的树（基于分区变量 x1-x4）的子组特定系数指标 t。因此，该模型不包括 x1 和 x2 的截距或全局主效应的可能性 - 正如您的数据中模拟的那样。由于这些主要影响非常明显，因此该模型除了通过进一步拆分来捕获这些影响之外别无他法。因此，您示例中的大树。
在经典的 MOB 框架中（Zeileis et al. 2008, Journal of Computational and Graphical Statistics, doi:10.1198/106186008X319331), the only option would be to include every relevant regressor in the model to be partitioned, i.e., logit(mu) = beta0_tree(x1-x4) + beta1_tree(x1_x4) * x1 + beta2_tree(x1-x4) * x2 + alpha_tree(x1-x4) * t. This works and will detect subgroups with respect to only the treatment effect alpha * t as well but, of course, loses some efficiency because it re-estimates the beta coefficients in every subgroup as well. A discussion of this approach specifically for subgroup analyses is available in Seibold et al. (2016a), The International Journal of Biostatistics, doi:10.1515/ijb-2015-0032.
最近，我们建议对 MOB 进行改编，我们将其称为 PALM 树，用于部分加性线性模型树（Seibold 等人 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 的 post 剪枝（结合较大的显着性水平），有时也可以发现与 x2 的交互。据推测，更大的样本也会更频繁地产生 "true" 树。因此，该模型原则上似乎能够适合您模拟的模型，只是不适合这种特定设置。
就个人而言，我总是让截距 和 治疗效果因亚组而异。因为如果有任何我忽略的主要影响，这可能会产生更好的模型。如果每个子组中都包含一个截距，那么最好将 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))
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 = as.factor(rbinom(n,1,pr))
dt <- data.frame(x1,x2,x3=as.factor(x3),
x4,t=as.factor(t),y)
在两个协变量 space:
下，治疗效果应该具有统计学意义
dt1 <- dt[x1>=-0.3&x2>=0.4,]
dt2 <- dt[x1< -0.3&x2<0.4,]
chisq.test(table(dt1$t,dt1$y))
Pearson's Chi-squared test with Yates' continuity correction
data: table(dt1$t, dt1$y)
X-squared = 11.684, df = 1, p-value = 0.0006305
prop.table(table(dt1$t,dt1$y),1)
chisq.test(table(dt2$t,dt2$y))
prop.table(table(dt2$t,dt2$y),1)
chisq.test(table(dt$t,dt$y))
我使用 glmtree() 函数，结果与您的输出不同，左侧正确识别了与 x2 的交互，但右侧未能捕获 x2 上的分区。
library(partykit)
tr1 <- glmtree(y~t+x1+x2 | x1+x2+x3+x4,
data=dt,
family = binomial())
plot(tr1,tp_args = list(which = 1))
最近一篇关于"Tutorial in biostatistics: data-driven subgroup identification and analysis in clinical trials "的论文（DOI: 10.1002/sim.7064）很感兴趣，想在"Performance of the tree-based regression approach"部分重现结果。但是，分区树好像没有得到我想要的结果。
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 上分裂。然而，它似乎并没有这样做。x1 主导了节点分区，而 x2 似乎不是分区过程的重要决定因素。代码有什么问题？
您指定的基于模型的 GLM 树尝试拟合以下模型： logit(prob) = alpha_tree(x1-x4) * t，其中 alpha_tree(x1-x4) 是用于处理的树（基于分区变量 x1-x4）的子组特定系数指标 t。因此，该模型不包括 x1 和 x2 的截距或全局主效应的可能性 - 正如您的数据中模拟的那样。由于这些主要影响非常明显，因此该模型除了通过进一步拆分来捕获这些影响之外别无他法。因此，您示例中的大树。
在经典的 MOB 框架中（Zeileis et al. 2008, Journal of Computational and Graphical Statistics, doi:10.1198/106186008X319331), the only option would be to include every relevant regressor in the model to be partitioned, i.e., logit(mu) = beta0_tree(x1-x4) + beta1_tree(x1_x4) * x1 + beta2_tree(x1-x4) * x2 + alpha_tree(x1-x4) * t. This works and will detect subgroups with respect to only the treatment effect alpha * t as well but, of course, loses some efficiency because it re-estimates the beta coefficients in every subgroup as well. A discussion of this approach specifically for subgroup analyses is available in Seibold et al. (2016a), The International Journal of Biostatistics, doi:10.1515/ijb-2015-0032.
最近，我们建议对 MOB 进行改编，我们将其称为 PALM 树，用于部分加性线性模型树（Seibold 等人 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 的 post 剪枝（结合较大的显着性水平），有时也可以发现与 x2 的交互。据推测，更大的样本也会更频繁地产生 "true" 树。因此，该模型原则上似乎能够适合您模拟的模型，只是不适合这种特定设置。
就个人而言，我总是让截距 和 治疗效果因亚组而异。因为如果有任何我忽略的主要影响，这可能会产生更好的模型。如果每个子组中都包含一个截距，那么最好将 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))
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 = as.factor(rbinom(n,1,pr))
dt <- data.frame(x1,x2,x3=as.factor(x3),
x4,t=as.factor(t),y)
在两个协变量 space:
下，治疗效果应该具有统计学意义dt1 <- dt[x1>=-0.3&x2>=0.4,]
dt2 <- dt[x1< -0.3&x2<0.4,]
chisq.test(table(dt1$t,dt1$y))
Pearson's Chi-squared test with Yates' continuity correction
data: table(dt1$t, dt1$y)
X-squared = 11.684, df = 1, p-value = 0.0006305
prop.table(table(dt1$t,dt1$y),1)
chisq.test(table(dt2$t,dt2$y))
prop.table(table(dt2$t,dt2$y),1)
chisq.test(table(dt$t,dt$y))
我使用 glmtree() 函数，结果与您的输出不同，左侧正确识别了与 x2 的交互，但右侧未能捕获 x2 上的分区。
library(partykit)
tr1 <- glmtree(y~t+x1+x2 | x1+x2+x3+x4,
data=dt,
family = binomial())
plot(tr1,tp_args = list(which = 1))