置信区间:使用带有 logistic() 模型的 mob() 树(partykit 包)

Confidence intervals: Using mob() trees (partykit package) with logistic() model

如何为我的 MOB 对象终端节点中的模型计算 CI? 或者,我可以提取终端节点的每个模型,以便在我的例子中,它是一个逻辑回归模型并计算每个模型的 CI?

谢谢

基础知识

使用 refit.modelparty(),您可以重新拟合与 modelparty 对象的节点关联的模型,如 mob()glmtree() 返回的那样。然后您可以将常用的函数应用于这些以提取您感兴趣的信息。

注意事项

  1. 学习树后的推理不再“诚实”,即在您的情况下,置信区间的覆盖范围可能不会达到标称水平。为了获得真实的结果,您必须拆分数据并在一个子样本上学习树,然后在另一个子样本的终端节点中重新拟合模型。

  2. 如果使用glmtree()那么每个节点内的模型都不是完整的glm对象,以便更有效地拟合树。因此,来自 confint() 方法的概况似然置信区间不是开箱即用的。您可以求助于 Wald 置信区间(例如,在 lmtestcoefci() 中可用)或重新拟合模型(类似于评论 1)。

插图

作为起点,让我们首先考虑一个基于 glmtree()coefci() 的简单示例:

## Pima Indians diabetes data
data("PimaIndiansDiabetes", package = "mlbench")
 
## recursive partitioning of a logistic regression model
library("partykit")
pid_tree <- glmtree(diabetes ~ glucose | pregnant +
  pressure + triceps + insulin + mass + pedigree + age,
  data = PimaIndiansDiabetes, family = binomial)

## extract fitted models in terminal nodes
pid_glm <- refit.modelparty(pid_tree,
  node = nodeids(pid_tree, terminal = TRUE))

## compute Wald confidence intervals
library("lmtest")
lapply(pid_glm, coefci)
## $`2`
##                    2.5 %      97.5 %
## (Intercept) -13.36226424 -6.54075507
## glucose       0.03496439  0.08245134
## 
## $`4`
##                   2.5 %      97.5 %
## (Intercept) -8.27393574 -5.13723535
## glucose      0.03466962  0.05900533
## 
## $`5`
##                   2.5 %      97.5 %
## (Intercept) -3.84548740 -1.69642032
## glucose      0.01529946  0.03177217

对于诚实的做法,我们可以这样做:

## id for sample splitting
n <- nrow(PimaIndiansDiabetes)
id <- sample(1:n, round(n/2))

## estimate tree on learning sample
pid_tree <- glmtree(diabetes ~ glucose | pregnant +
  pressure + triceps + insulin + mass + pedigree + age,
  data = PimaIndiansDiabetes[id,], family = binomial)

## out-of-sample prediction
pid_new <- PimaIndiansDiabetes[-id,]
pid_new$node <- predict(pid_tree, newdata = pid_new, type = "node")

## fit separate models on each subset, splitted by predicted node
pid_glm <- lapply(
  split(pid_new, pid_new$node),
  function(d) glm(diabetes ~ glucose, data = d, family = binomial)
)

## obtain profile likelihood confidence intervals
lapply(pid_glm, confint)
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## Waiting for profiling to be done...
## $`2`
##                   2.5 %     97.5 %
## (Intercept) -8.25379802 -4.5031005
## glucose      0.02743191  0.0569824
## 
## $`4`
##                    2.5 %     97.5 %
## (Intercept) -25.32777607 -4.5078931
## glucose       0.02090339  0.1617026
## 
## $`5`
##                   2.5 %      97.5 %
## (Intercept) -6.38422374 -2.66969443
## glucose      0.02222873  0.05060984

在基于逻辑 GLM 等线性预测器的模型中,也可以拟合表示整棵树的单个模型。您只需要包括所有回归变量与节点指示器的交互。要获得相同的参数和置信区间,您需要使用嵌套编码 (/) 而不是默认交互编码 (*):

pid_glm2 <- glm(diabetes ~ 0 + factor(node)/glucose,
  data = pid_new, family = binomial)
confint(pid_glm2)
## Waiting for profiling to be done...
##                              2.5 %      97.5 %
## factor(node)2          -8.25379802 -4.50310050
## factor(node)4         -25.32777607 -4.50789313
## factor(node)5          -6.38422332 -2.66969490
## factor(node)2:glucose   0.02743191  0.05698240
## factor(node)4:glucose   0.02090339  0.16170262
## factor(node)5:glucose   0.02222874  0.05060983