置信区间:使用带有 logistic() 模型的 mob() 树(partykit 包)
Confidence intervals: Using mob() trees (partykit package) with logistic() model
如何为我的 MOB 对象终端节点中的模型计算 CI?
或者,我可以提取终端节点的每个模型,以便在我的例子中,它是一个逻辑回归模型并计算每个模型的 CI?
谢谢
基础知识
使用 refit.modelparty()
,您可以重新拟合与 modelparty
对象的节点关联的模型,如 mob()
或 glmtree()
返回的那样。然后您可以将常用的函数应用于这些以提取您感兴趣的信息。
注意事项
学习树后的推理不再“诚实”,即在您的情况下,置信区间的覆盖范围可能不会达到标称水平。为了获得真实的结果,您必须拆分数据并在一个子样本上学习树,然后在另一个子样本的终端节点中重新拟合模型。
如果使用glmtree()
那么每个节点内的模型都不是完整的glm
对象,以便更有效地拟合树。因此,来自 confint()
方法的概况似然置信区间不是开箱即用的。您可以求助于 Wald 置信区间(例如,在 lmtest
的 coefci()
中可用)或重新拟合模型(类似于评论 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
如何为我的 MOB 对象终端节点中的模型计算 CI? 或者,我可以提取终端节点的每个模型,以便在我的例子中,它是一个逻辑回归模型并计算每个模型的 CI?
谢谢
基础知识
使用 refit.modelparty()
,您可以重新拟合与 modelparty
对象的节点关联的模型,如 mob()
或 glmtree()
返回的那样。然后您可以将常用的函数应用于这些以提取您感兴趣的信息。
注意事项
学习树后的推理不再“诚实”,即在您的情况下,置信区间的覆盖范围可能不会达到标称水平。为了获得真实的结果,您必须拆分数据并在一个子样本上学习树,然后在另一个子样本的终端节点中重新拟合模型。
如果使用
glmtree()
那么每个节点内的模型都不是完整的glm
对象,以便更有效地拟合树。因此,来自confint()
方法的概况似然置信区间不是开箱即用的。您可以求助于 Wald 置信区间(例如,在lmtest
的coefci()
中可用)或重新拟合模型(类似于评论 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