Montecarlo Simulation 使用 mob() 算法（partykit 包）恢复正确识别模型的计数
Montecarlo Simulation using mob() algorithm (partykit package) to recover the count of the correctly identified models
我正在使用 partykit 包中的 mob() 函数，在解析获得的模型时遇到一些问题。
我的主要目的是检查大约需要多大的样本量才能检测到真实的
存在中断时的数据生成过程 (DGP) 的结构。
下面的代码对带有中断的数据执行蒙特卡洛模拟，并尝试确定中断是否被 M 波动测试捕获。
更具体地说，我想计算模型实际捕获 DGP 的次数占模拟总数 (nreps) 的次数，条件是固定样本大小 (N) 来了解我需要多少数据才能捕获真正的 DGP。
在代码的最后，你可以看到我从模拟中得到的列表。问题是我无法恢复控制台显示的信息。
另外，我对如何计算“正确识别的模型”有一些疑问。现在，我想到的是，如果模型在指定区域 (z1>0) 进入正确的变量 (z1) 并且对中断区域有一定的容差，则算作正数如果休息时间在 z1>0.1 或 z1>-0.4，这对我来说也是有效的。那么，有没有什么简单的方法可以计算出具有上述特征的模型？
我希望我的示例足够清楚，可以帮助您解决问题。非常感谢您。
模拟模型
- 生成 DGP 的成分。
library("partykit")
library(data.table) ## optional, but what I'll use to coerce the list into a DT
library(future.apply) ## for parallel stuff
plan(multisession) ## use all available cores
#sample size
N <- 300
#coeficients
betas <- list()
betas$b0 <- 1
betas$b1_up <- 2.4
betas$b1_down <- 2
betas$b2_up <- 2.4
betas$b2_down <- 2
#mob() ingredients
ols_formula <- y ~ x1 + x2 | z1 + z2
# the ""0 +"" ---> supress the 'double' interecept
ols <- function(y, x, start = NULL, weights = NULL, offset = NULL, ...) {lm(y ~ 0 + x)}
- 生成数据并使用 mob 算法拟合 OLS 的函数。
reg_simulation_mob <- function(...){
#data
data <- data.frame(
x1 = rnorm(N),
x2 = rnorm(N),
z1 = rnorm(N),
z2 = rnorm(N),
e = rnorm(N))
#dependent variable
data$y <- betas$b0 + with(data, ifelse(z1>0,
betas$b1_up * x1 + betas$b2_up * x2 ,
betas$b1_down * x1 + betas$b2_down * x2 )
+ e )
#Estimate mob()-OLS
ols_mob <- mob(ols_formula,
data = data,
fit = ols)
# return(ols$coefficients)
return(ols_mob)
}
- 上述设置的 Montecarlo 模拟（仅 2 次试验）。
# N repetitions
nreps <- 2
## Parallel version
results <- future_lapply(1:nreps, reg_simulation_mob, future.seed = 1234L)
得到的结果
正如您在下面的第一次试验 (results[[1]]) 中看到的那样，模型找到了正确的中断点，但在第二次试验 (results[[2]]) 中却找不到它。
> results
[[1]]
Model-based recursive partitioning (ols)
Model formula:
y ~ x1 + x2 | z1 + z2
Fitted party:
[1] root
| [2] z1 <= 0.00029: n = 140
| x(Intercept) xx1 xx2
| 0.9597894 1.7552122 2.1360788
| [3] z1 > 0.00029: n = 160
| x(Intercept) xx1 xx2
| 0.9371795 2.4745728 2.5087608
Number of inner nodes: 1
Number of terminal nodes: 2
Number of parameters per node: 3
Objective function: 422.2329
[[2]]
Model-based recursive partitioning (ols)
Model formula:
y ~ x1 + x2 | z1 + z2
Fitted party:
[1] root: n = 300
x(Intercept) xx1 xx2
1.015224 2.175625 2.200746
Number of inner nodes: 0
Number of terminal nodes: 1
Number of parameters per node: 3
Objective function: 422.3085
在下图中，您可以观察列表的结构results，其中我找不到控制台上显示的信息（即节点数、参数、阈值等。）
首先，我建议使用 lmtree() 函数而不是普通的 mob()。前者更快，带有更好的打印和绘图默认值，并且有更多的预测选项。
其次，我建议您查阅 vignette("partykit", package = "partykit")，其中解释了 party 对象的构造方式以及涉及哪些 类 和方法。
第三，要确定哪个变量（如果有的话）用于在根节点中进行拆分，可能有兴趣从所有参数不稳定性测试中提取结果。有一个专门的sctest()（结构变化测试）方法来获得这个：
library("strucchange")
sctest(results[[1]], node = 1)
## z1 z2
## statistic 4.072483e+01 6.1762164
## p.value 5.953672e-07 0.9153013
sctest(results[[2]], node = 1)
## z1 z2
## statistic 12.2810548 10.1944484
## p.value 0.2165527 0.4179998
根 $node 中 $split（如果有）对应的 partysplit 对象可能最容易手动提取：
results[[1]]$node$split
## $varid
## [1] 4
##
## $breaks
## [1] 0.0002853492
##
## $index
## NULL
##
## $right
## [1] TRUE
##
## $prob
## NULL
##
## $info
## NULL
##
## attr(,"class")
## [1] "partysplit"
results[[2]]$node$split
## NULL
变量 id 与以下变量的顺序有关：
names(results[[1]]$data)
## [1] "y" "x1" "x2" "z1" "z2"
最后，关于评估什么的问题：一切都取决于确定正确的分裂变量。如果这样做正确，那么分裂点估计会快速收敛到真实值，然后参数估计也会收敛。例如，请参阅我们最近的 arXiv 论文 (https://arxiv.org/abs/1906.10179)，其中包含更大的模拟研究并提供了复制代码。
因此，通常，我要么评估在第一次拆分中选择正确变量的概率。或者，我查看每个观察.
的估计 vs.true 系数 的 RMSE
更新： 在根节点之外，您可以使用 nodeapply() 从各个节点提取信息。但是，我不建议评估所有拆分，因为匹配哪个估计拆分与哪个真实拆分匹配变得越来越困难。相反，我们经常评估拟合分区与真实分区的相似程度，例如，使用调整后的兰德指数。同样，您可以在上述 arXiv 论文中找到一个示例。
这个答案是基于@AchimZeileis 教授在他的文章（https://arxiv.org/abs/1906.10179）中提供的参考，它专门用于我原始问题的第二部分，即关于问题：如何计算正确指定的模型（树）？
简短的回答。
文章将问题分为两种不同类型的数据生成过程（DGP）。在第一个中，数据只有一个变量中断（“stump”案例），作者根据 M 出现的次数统计正确识别模型的数量- 波动测试正确地识别了产生中断的变量（一个变量中只有一个真正的中断，9 个嘈杂的候选变量将在没有中断的情况下分裂变量）。第二个 DGP 是一个有两个断点的模型（“树” 案例），他们使用 Adjusted Rand Index (ARI) 来评估模型的性能，作为真实树与预测树的相似度。
非常的回答。
让我们分解一下@AchimZeileis 推荐的ARI for 6 different illustrative possible trees that can be obtained at different sample sizes. The code used here is highly based on the supplemental material of the article。
数据生成过程：树结构
真正的dgp有2个断点，如下图所示。第一个由变量 z2 生成，第二个由 z1 生成。在下面的代码片段中，delta 等于 1。第一次突破的阈值（取决于 z2）等于 0.3，第二次突破的阈值（取决于 z1) 等于-0.3（值可以在对象中看到xi = c(-0.3, 0.3)）
#function from https://arxiv.org/src/1906.10179v1/anc
dgp_tree <- function(nobs = 1000, delta = 1, xi = c(-0.3, 0.3),
sigma = 1, seed = 7,
changetype = "abrupt",
variation = "all",
beta0 = NULL, beta1 = NULL)
{
# check input values
if(variation != "all") stop("variation can only be set to 'all' in dgp_tree")
if(changetype != "abrupt") stop("changetype can only be abrupt in dgp_tree")
if(!is.null(beta0) | !is.null(beta1)) warning("values for beta0 or beta1 are ignored since variation='all' for dgp_tree")
set.seed(seed)
if(length(xi)==1){
xi1 <- xi2 <- xi
} else {
xi1 <- xi[1]
xi2 <- xi[2]
}
z1 <- runif(nobs,-1,1)
z2 <- runif(nobs,-1,1)
z3 <- rnorm(nobs, 0, 1)
z4 <- rnorm(nobs, 0, 1)
z5 <- rnorm(nobs, 0, 1)
z6 <- rnorm(nobs, 0, 1)
z7 <- rnorm(nobs, 0, 1)
z8 <- runif(nobs, -1, 1)
z9 <- runif(nobs, -1, 1)
z10 <- runif(nobs, -1, 1)
id <- numeric(length(z1))
x <- runif(nobs, min = -1, max = 1)
beta0 <- delta * (-1)^(z1<xi1) * 0^(z2<xi2)
beta1 <- delta * (-1)^(z2>=xi2)
id <- 1 + (z2>=xi2) + (z2>=xi2)*(z1>=xi1)
mu <- beta0 + beta1 * x
y <- rnorm(nobs, mu, sigma)
d <- data.frame(y = y, x = x,
z1 = z1, z2 = z2, z3 = z3, z4 = z4, z5 = z5, z6 = z6, z7 = z7, z8 = z8, z9 = z9, z10 = z10,
beta0 = beta0, beta1 = beta1, mu = mu, sigma = rep.int(sigma, times = length(y)), id = id)
return(d)
}
Adjusted Rand Index
其中包含的函数article, there is one to compute the ARI, and it is listed below to be used in the following examples. It resembles almost exactly letter by letter the notation used here.
# function to compute adjusted Rand Index from https://arxiv.org/src/1906.10179v1/anc
adj_rand_index <- function(x, y) {
tab <- table(x, y)
a <- rowSums(tab)
b <- colSums(tab)
M <- sum(choose(tab, 2))
N <- choose(length(x), 2)
A <- sum(choose(a, 2))
B <- sum(choose(b, 2))
c(ARI = (M - (A * B) / N) / (0.5 * (A + B) - (A * B) / N))
}
模拟数据
library(partykit)
library(future.apply) ## for parallel stuff
plan(multisession) ## use all available cores
ols_formula <- y ~ x | z1 + z2 +z3 +z4 + z5 +z6 +z7+ z8 +z9 +z10
ols <- function(y, x, start = NULL, weights = NULL, offset = NULL, ...) {lm(y ~ 0 + x)}
sim_ari <- function(n){
tree_data <- dgp_tree(nobs = n)
ols_mob <- mob(ols_formula,
data = tree_data,
fit = ols)
prednode <- predict(ols_mob ,
type = "node")
cross_table <- table(prednode,tree_data$id)
ari <- adj_rand_index(prednode,
tree_data$id)
print(n)
print(ari)
return(
list(
ols_mob = ols_mob,
cross_table = cross_table,
ari=ari,
data = tree_data)
)
}
n_levels <- c(55, ## no break
87, ## only one break
123, ## Correct structure, but poor performance
199, ## Nested break in second leaf
667, ## Additional break in first leaf
5000 ## Perfect model
)
ari <- future_lapply(n_levels, sim_ari, future.seed = 1234L)
六棵树的故事。
下面就ARI如何准确地捕捉到正确树与估计树的相似度，分析了六种情况。比较树的关键是 id，它显示了每个观察值应该属于树中的哪个叶子。例如，如果一个观测值的 id 值为 1，则它满足分配给上图中节点编号 2 的要求。另一方面，如果 id 等于 2，则观察应分配给同一图片中的节点号 4。最后，如果 id 等于 3，则将其分配给节点号 5。您可以在以下行中检查此推理 id <- 1 + (z2>=xi2) + (z2>=xi2)*(z1>=xi1)
第一棵树 (n=55)：没有中断
分析的第一个案例对应于未识别到中断的情况。这里，在这种情况下，ARI 等于 0。
##### First Tree (n=55): No break ####
ari[[1]][[1]]
## Fitted party:
## [1] root: n = 55
## x(Intercept) xx
## -0.01309586 0.39291089
##
## Number of inner nodes: 0
## Number of terminal nodes: 1
## Number of parameters per node: 2
## Objective function: 95.58631
这里有趣的是，所有观察值都分配给了根节点。因此，当crosstable预测节点prednode_1时，我们看到所有可能的id值都属于预测树的根节点[1]（基本上，因为没有其他选择。）通过使用函数 adj_rand_index()，您可以检查这是否导致 ARI 等于 0。
#data first tree (n=55)
data_1 <- ari[[1]][[4]]
#predicted node first iteration
data_1$prednode_1 <- predict(ari[[1]][[1]], type = "node")
#Cross table
with(data_1, table(prednode_1 ,id))
## id
## prednode_1 1 2 3
## 1 37 7 11
#adj_rand_index
ari[[1]][[3]]
第二棵树 (n=87)：只发现一个断点
这个案例很有趣，因为它部分地识别了树的结构（即缺少 z1 上的中断）。
##### Second Tree (n=87): Extra partition in node[5] ####
ari[[2]][[1]]
# Fitted party:
# [1] root
# | [2] z2 <= 0.29288: n = 57
# | x(Intercept) xx
# | 0.133293 1.082701
# | [3] z2 > 0.29288: n = 30
# | x(Intercept) xx
# | 0.2598309 -1.8014133
#
# Number of inner nodes: 1
# Number of terminal nodes: 2
# Number of parameters per node: 2
# Objective function: 122.0116
此外，我们可以检查当交叉预测节点和真实节点时，我们看到即使在这棵不完美的树中，一些观察结果也符合标准。这意味着 57 个观测值正确分配给了第一个节点，9 个观测值被正确分配给了第二个分支。最后，30 因为最后一个节点根本没有被识别而被错误分配。这导致 ARI 等于 0.8577366，这是第一棵树的巨大改进。
#data second iteration (n=87)
data_2 <- ari[[2]][[4]]
#predicted node first iteration
data_2$prednode_2 <- predict(ari[[2]][[1]], type = "node")
#Cross table
with(data_2, table(prednode_2 ,id))
# id
# prednode_2 1 2 3
# 2 57 0 0
# 3 1 9 20
#adj_rand_index
ari[[2]][[3]]
# > ari[[2]][[3]]
# ARI
# 0.8577366
第三棵树（n=123）：结构正确但性能较差
这个案例很有趣，因为它确实恢复了真实的树结构，但它的性能比仅部分识别其结构的后三个更差。
##### Third Tree (n=123): Correct structure but poor performance ####
ari[[3]][[1]]
# Fitted party:
# [1] root
# | [2] z2 <= 0.07319: n = 60
# | x(Intercept) xx
# | -0.1723388 1.1071878
# | [3] z2 > 0.07319
# | | [4] z1 <= -0.35485: n = 22
# | | x(Intercept) xx
# | | -0.7166565 -0.6791717
# | | [5] z1 > -0.35485: n = 41
# | | x(Intercept) xx
# | | 0.7096033 -0.8605967
#
# Number of inner nodes: 2
# Number of terminal nodes: 3
# Number of parameters per node: 2
# Objective function: 156.4397
下面我们可以看到，当我们对预测节点和真实节点进行交叉表时，我们看到 16 (10 + 6) 个观测值被错误分类，这一事实导致 ARI 为 0.6117612.
#data third iteration (n=123)
data_3 <- ari[[3]][[4]]
#predicted node first iteration
data_3$prednode_3 <- predict(ari[[3]][[1]], type = "node")
#Cross table
with(data_3, table(prednode_3 ,id))
# id
# prednode_3 1 2 3
# 2 60 0 0
# 4 6 16 0
# 5 10 0 31
#adj_rand_index
ari[[3]][[3]]
# > ari[[3]][[3]]
# ARI
# 0.6117612
第四棵树：第四棵树 (n=199)：节点 [5]
处的额外叶子
这里识别的树与原来的有偏差，因为它多了一片node[5]的叶子，这在真实数据中是不存在的。
##### Forth Tree (n=199): Extra leaf at node[5] ####
ari[[4]][[1]]
# Fitted party:
# [1] root
# | [2] z2 <= -0.19806: n = 79
# | x(Intercept) xx
# | 0.06455217 1.51512672
# | [3] z2 > -0.19806
# | | [4] z1 <= -0.27127: n = 44
# | | x(Intercept) xx
# | | -0.4863122 -0.3860951
# | | [5] z1 > -0.27127
# | | | [6] z2 <= 0.17481: n = 23
# | | | x(Intercept) xx
# | | | -0.1335096 0.2046050
# | | | [7] z2 > 0.17481: n = 53
# | | | x(Intercept) xx
# | | | 1.0868488 -0.0290925
#
# Number of inner nodes: 3
# Number of terminal nodes: 4
# Number of parameters per node: 2
# Objective function: 282.6727
这里真实节点和预测节点的交叉表计数很有趣，因为节点 [6] 和 [7] 在真实数据中不存在，但它们得到的观察结果应该是，例如，分配给节点 [1]（分别为 23 和 7 观察值。）这种错误分配将 ARI 指数减少到 0.4649789。
#data forth iteration (n=199)
data_4 <- ari[[4]][[4]]
#predicted node first iteration
data_4$prednode_4 <- predict(ari[[4]][[1]], type = "node")
#Cross table
with(data_4, table(prednode_4 ,id))
# id
# prednode_4 1 2 3
# 2 79 0 0
# 4 16 27 1
# 6 23 0 0
# 7 7 0 46
#adj_rand_index
ari[[4]][[3]]
# ARI
# 0.4649789
第五棵树 (n=667)：节点 [2]
处的额外叶子
这是结构不正确的树的另一个示例，其中额外的叶子（基于 z5 上的分区是不正确的（！））附加到节点 [2]。
##### Fifth Tree (n=667): Extra leaf at node[2] ####
ari[[5]][[1]]
# Fitted party:
# [1] root
# | [2] z2 <= 0.28476
# | | [3] z5 <= 0.76285: n = 322
# | | x(Intercept) xx
# | | -0.1322881 0.9535337
# | | [4] z5 > 0.76285: n = 96
# | | x(Intercept) xx
# | | 0.1686863 1.3878776
# | [5] z2 > 0.28476
# | | [6] z1 <= -0.32001: n = 89
# | | x(Intercept) xx
# | | -0.9139858 -0.7957158
# | | [7] z1 > -0.32001: n = 160
# | | x(Intercept) xx
# | | 0.7661154 -0.8656553
#
# Number of inner nodes: 3
# Number of terminal nodes: 4
# Number of parameters per node: 2
# Objective function: 927.9088
来自预测节点和正确节点的交叉表向我们展示了实际上属于第一个节点 [1] 的大多数观察值 (322) 被分配给了预测节点 [3].最后，这种糟糕的结构导致 0.6932132.`
的 ARI
#data third iteration (n=667)
data_5 <- ari[[5]][[4]]
#predicted node first iteration
data_5$prednode_5 <- predict(ari[[5]][[1]], type = "node")
#Cross table
with(data_5, table(prednode_5 ,id))
# id
# prednode_5 1 2 3
# 3 322 0 0
# 4 96 0 0
# 6 0 89 0
# 7 3 3 154
#adj_rand_index
ari[[5]][[3]]
# ARI
# 0.6932132
第六棵树（n=5000）：黄金树！
这棵最终树完美地恢复了数据，无论是在树结构中还是在将观察值分配给每个叶子上。
##### Sixth Tree (n=5000): Extra leaf at node[2] ####
ari[[6]][[1]]
# Fitted party:
# [1] root
# | [2] z2 <= 0.29971: n = 3187
# | x(Intercept) xx
# | -0.008719923 1.022232280
# | [3] z2 > 0.29971
# | | [4] z1 <= -0.30286: n = 609
# | | x(Intercept) xx
# | | -0.9488846 -0.9813765
# | | [5] z1 > -0.30286: n = 1204
# | | x(Intercept) xx
# | | 1.0281410 -0.9565637
#
# Number of inner nodes: 2
# Number of terminal nodes: 3
# Number of parameters per node: 2
# Objective function: 6992.848
这里我们可以从预测节点和真实节点的交叉表中看到，它完美地分配了每个观测值，它确实属于它，导致 ARI 等于 1 .
#data sixt iteration (n=5000)
data_6 <- ari[[6]][[4]]
#predicted node first iteration
data_6$prednode_6 <- predict(ari[[6]][[1]], type = "node")
#Cross table
with(data_6, table(prednode_6 ,id))
# id
# prednode_6 1 2 3
# 2 3187 0 0
# 4 0 609 0
# 5 0 0 1204
#adj_rand_index
ari[[6]][[3]]
# ARI
# 1
结论。
可以从上图中恢复一些重要的要点。
1.- ARI 可用于评估预测树的相似程度，该预测树在数据生成过程中可能具有与真实树非常不同的结构。
2.- 恢复树的正确结构不会导致 ARI 等于一。
3.- 不正确的树不一定会使 ARI 等于零。
最终模拟。
总而言之，这是一个小型模拟，可以了解 ARI 指数在样本量增加时的表现。
### Final simulation
n_levels <-seq(from= 10, to= 2000, by= 5)
ari <- lapply(n_levels, sim_ari)
ari_models<- function(i){
ari <- ari_sim[[i]]$ari
n <- nobs(ari_sim[[i]]$ols_mob)
return(
list(ari = ari , n= n )
)
}
ari_n_list <- lapply(1:length(ari_sim), ari_models)
df <- data.frame(matrix(unlist(ari_n_list), nrow=length(ari_n_list), byrow=T))
colnames(df) <- c("ARI" , "N")
library(ggplot2)
ggplot(df, aes(N)) +
geom_line(aes(y = ARI, colour = "ARI"))
我正在使用 partykit 包中的 mob() 函数，在解析获得的模型时遇到一些问题。
我的主要目的是检查大约需要多大的样本量才能检测到真实的 存在中断时的数据生成过程 (DGP) 的结构。
下面的代码对带有中断的数据执行蒙特卡洛模拟，并尝试确定中断是否被 M 波动测试捕获。
更具体地说，我想计算模型实际捕获 DGP 的次数占模拟总数 (nreps) 的次数，条件是固定样本大小 (N) 来了解我需要多少数据才能捕获真正的 DGP。
在代码的最后，你可以看到我从模拟中得到的列表。问题是我无法恢复控制台显示的信息。
另外，我对如何计算“正确识别的模型”有一些疑问。现在，我想到的是，如果模型在指定区域 (z1>0) 进入正确的变量 (z1) 并且对中断区域有一定的容差，则算作正数如果休息时间在 z1>0.1 或 z1>-0.4，这对我来说也是有效的。那么，有没有什么简单的方法可以计算出具有上述特征的模型？
我希望我的示例足够清楚，可以帮助您解决问题。非常感谢您。
模拟模型
- 生成 DGP 的成分。
library("partykit")
library(data.table) ## optional, but what I'll use to coerce the list into a DT
library(future.apply) ## for parallel stuff
plan(multisession) ## use all available cores
#sample size
N <- 300
#coeficients
betas <- list()
betas$b0 <- 1
betas$b1_up <- 2.4
betas$b1_down <- 2
betas$b2_up <- 2.4
betas$b2_down <- 2
#mob() ingredients
ols_formula <- y ~ x1 + x2 | z1 + z2
# the ""0 +"" ---> supress the 'double' interecept
ols <- function(y, x, start = NULL, weights = NULL, offset = NULL, ...) {lm(y ~ 0 + x)}
- 生成数据并使用 mob 算法拟合 OLS 的函数。
reg_simulation_mob <- function(...){
#data
data <- data.frame(
x1 = rnorm(N),
x2 = rnorm(N),
z1 = rnorm(N),
z2 = rnorm(N),
e = rnorm(N))
#dependent variable
data$y <- betas$b0 + with(data, ifelse(z1>0,
betas$b1_up * x1 + betas$b2_up * x2 ,
betas$b1_down * x1 + betas$b2_down * x2 )
+ e )
#Estimate mob()-OLS
ols_mob <- mob(ols_formula,
data = data,
fit = ols)
# return(ols$coefficients)
return(ols_mob)
}
- 上述设置的 Montecarlo 模拟（仅 2 次试验）。
# N repetitions
nreps <- 2
## Parallel version
results <- future_lapply(1:nreps, reg_simulation_mob, future.seed = 1234L)
得到的结果
正如您在下面的第一次试验 (results[[1]]) 中看到的那样，模型找到了正确的中断点，但在第二次试验 (results[[2]]) 中却找不到它。
> results
[[1]]
Model-based recursive partitioning (ols)
Model formula:
y ~ x1 + x2 | z1 + z2
Fitted party:
[1] root
| [2] z1 <= 0.00029: n = 140
| x(Intercept) xx1 xx2
| 0.9597894 1.7552122 2.1360788
| [3] z1 > 0.00029: n = 160
| x(Intercept) xx1 xx2
| 0.9371795 2.4745728 2.5087608
Number of inner nodes: 1
Number of terminal nodes: 2
Number of parameters per node: 3
Objective function: 422.2329
[[2]]
Model-based recursive partitioning (ols)
Model formula:
y ~ x1 + x2 | z1 + z2
Fitted party:
[1] root: n = 300
x(Intercept) xx1 xx2
1.015224 2.175625 2.200746
Number of inner nodes: 0
Number of terminal nodes: 1
Number of parameters per node: 3
Objective function: 422.3085
在下图中，您可以观察列表的结构results，其中我找不到控制台上显示的信息（即节点数、参数、阈值等。）
首先，我建议使用 lmtree() 函数而不是普通的 mob()。前者更快，带有更好的打印和绘图默认值，并且有更多的预测选项。
其次，我建议您查阅 vignette("partykit", package = "partykit")，其中解释了 party 对象的构造方式以及涉及哪些 类 和方法。
第三，要确定哪个变量（如果有的话）用于在根节点中进行拆分，可能有兴趣从所有参数不稳定性测试中提取结果。有一个专门的sctest()（结构变化测试）方法来获得这个：
library("strucchange")
sctest(results[[1]], node = 1)
## z1 z2
## statistic 4.072483e+01 6.1762164
## p.value 5.953672e-07 0.9153013
sctest(results[[2]], node = 1)
## z1 z2
## statistic 12.2810548 10.1944484
## p.value 0.2165527 0.4179998
根 $node 中 $split（如果有）对应的 partysplit 对象可能最容易手动提取：
results[[1]]$node$split
## $varid
## [1] 4
##
## $breaks
## [1] 0.0002853492
##
## $index
## NULL
##
## $right
## [1] TRUE
##
## $prob
## NULL
##
## $info
## NULL
##
## attr(,"class")
## [1] "partysplit"
results[[2]]$node$split
## NULL
变量 id 与以下变量的顺序有关：
names(results[[1]]$data)
## [1] "y" "x1" "x2" "z1" "z2"
最后，关于评估什么的问题：一切都取决于确定正确的分裂变量。如果这样做正确，那么分裂点估计会快速收敛到真实值，然后参数估计也会收敛。例如，请参阅我们最近的 arXiv 论文 (https://arxiv.org/abs/1906.10179)，其中包含更大的模拟研究并提供了复制代码。
因此，通常，我要么评估在第一次拆分中选择正确变量的概率。或者，我查看每个观察.
的估计 vs.true 系数 的 RMSE更新： 在根节点之外，您可以使用 nodeapply() 从各个节点提取信息。但是，我不建议评估所有拆分，因为匹配哪个估计拆分与哪个真实拆分匹配变得越来越困难。相反，我们经常评估拟合分区与真实分区的相似程度，例如，使用调整后的兰德指数。同样，您可以在上述 arXiv 论文中找到一个示例。
这个答案是基于@AchimZeileis 教授在他的文章（https://arxiv.org/abs/1906.10179）中提供的参考，它专门用于我原始问题的第二部分，即关于问题：如何计算正确指定的模型（树）？
简短的回答。
文章将问题分为两种不同类型的数据生成过程（DGP）。在第一个中，数据只有一个变量中断（“stump”案例），作者根据 M 出现的次数统计正确识别模型的数量- 波动测试正确地识别了产生中断的变量（一个变量中只有一个真正的中断，9 个嘈杂的候选变量将在没有中断的情况下分裂变量）。第二个 DGP 是一个有两个断点的模型（“树” 案例），他们使用 Adjusted Rand Index (ARI) 来评估模型的性能，作为真实树与预测树的相似度。
非常的回答。
让我们分解一下@AchimZeileis 推荐的ARI for 6 different illustrative possible trees that can be obtained at different sample sizes. The code used here is highly based on the supplemental material of the article。
数据生成过程：树结构
真正的dgp有2个断点，如下图所示。第一个由变量 z2 生成，第二个由 z1 生成。在下面的代码片段中，delta 等于 1。第一次突破的阈值（取决于 z2）等于 0.3，第二次突破的阈值（取决于 z1) 等于-0.3（值可以在对象中看到xi = c(-0.3, 0.3)）
#function from https://arxiv.org/src/1906.10179v1/anc
dgp_tree <- function(nobs = 1000, delta = 1, xi = c(-0.3, 0.3),
sigma = 1, seed = 7,
changetype = "abrupt",
variation = "all",
beta0 = NULL, beta1 = NULL)
{
# check input values
if(variation != "all") stop("variation can only be set to 'all' in dgp_tree")
if(changetype != "abrupt") stop("changetype can only be abrupt in dgp_tree")
if(!is.null(beta0) | !is.null(beta1)) warning("values for beta0 or beta1 are ignored since variation='all' for dgp_tree")
set.seed(seed)
if(length(xi)==1){
xi1 <- xi2 <- xi
} else {
xi1 <- xi[1]
xi2 <- xi[2]
}
z1 <- runif(nobs,-1,1)
z2 <- runif(nobs,-1,1)
z3 <- rnorm(nobs, 0, 1)
z4 <- rnorm(nobs, 0, 1)
z5 <- rnorm(nobs, 0, 1)
z6 <- rnorm(nobs, 0, 1)
z7 <- rnorm(nobs, 0, 1)
z8 <- runif(nobs, -1, 1)
z9 <- runif(nobs, -1, 1)
z10 <- runif(nobs, -1, 1)
id <- numeric(length(z1))
x <- runif(nobs, min = -1, max = 1)
beta0 <- delta * (-1)^(z1<xi1) * 0^(z2<xi2)
beta1 <- delta * (-1)^(z2>=xi2)
id <- 1 + (z2>=xi2) + (z2>=xi2)*(z1>=xi1)
mu <- beta0 + beta1 * x
y <- rnorm(nobs, mu, sigma)
d <- data.frame(y = y, x = x,
z1 = z1, z2 = z2, z3 = z3, z4 = z4, z5 = z5, z6 = z6, z7 = z7, z8 = z8, z9 = z9, z10 = z10,
beta0 = beta0, beta1 = beta1, mu = mu, sigma = rep.int(sigma, times = length(y)), id = id)
return(d)
}
Adjusted Rand Index
其中包含的函数article, there is one to compute the ARI, and it is listed below to be used in the following examples. It resembles almost exactly letter by letter the notation used here.
# function to compute adjusted Rand Index from https://arxiv.org/src/1906.10179v1/anc
adj_rand_index <- function(x, y) {
tab <- table(x, y)
a <- rowSums(tab)
b <- colSums(tab)
M <- sum(choose(tab, 2))
N <- choose(length(x), 2)
A <- sum(choose(a, 2))
B <- sum(choose(b, 2))
c(ARI = (M - (A * B) / N) / (0.5 * (A + B) - (A * B) / N))
}
模拟数据
library(partykit)
library(future.apply) ## for parallel stuff
plan(multisession) ## use all available cores
ols_formula <- y ~ x | z1 + z2 +z3 +z4 + z5 +z6 +z7+ z8 +z9 +z10
ols <- function(y, x, start = NULL, weights = NULL, offset = NULL, ...) {lm(y ~ 0 + x)}
sim_ari <- function(n){
tree_data <- dgp_tree(nobs = n)
ols_mob <- mob(ols_formula,
data = tree_data,
fit = ols)
prednode <- predict(ols_mob ,
type = "node")
cross_table <- table(prednode,tree_data$id)
ari <- adj_rand_index(prednode,
tree_data$id)
print(n)
print(ari)
return(
list(
ols_mob = ols_mob,
cross_table = cross_table,
ari=ari,
data = tree_data)
)
}
n_levels <- c(55, ## no break
87, ## only one break
123, ## Correct structure, but poor performance
199, ## Nested break in second leaf
667, ## Additional break in first leaf
5000 ## Perfect model
)
ari <- future_lapply(n_levels, sim_ari, future.seed = 1234L)
六棵树的故事。
下面就ARI如何准确地捕捉到正确树与估计树的相似度，分析了六种情况。比较树的关键是 id，它显示了每个观察值应该属于树中的哪个叶子。例如，如果一个观测值的 id 值为 1，则它满足分配给上图中节点编号 2 的要求。另一方面，如果 id 等于 2，则观察应分配给同一图片中的节点号 4。最后，如果 id 等于 3，则将其分配给节点号 5。您可以在以下行中检查此推理 id <- 1 + (z2>=xi2) + (z2>=xi2)*(z1>=xi1)
第一棵树 (n=55)：没有中断
分析的第一个案例对应于未识别到中断的情况。这里，在这种情况下，ARI 等于 0。
##### First Tree (n=55): No break ####
ari[[1]][[1]]
## Fitted party:
## [1] root: n = 55
## x(Intercept) xx
## -0.01309586 0.39291089
##
## Number of inner nodes: 0
## Number of terminal nodes: 1
## Number of parameters per node: 2
## Objective function: 95.58631
这里有趣的是，所有观察值都分配给了根节点。因此，当crosstable预测节点prednode_1时，我们看到所有可能的id值都属于预测树的根节点[1]（基本上，因为没有其他选择。）通过使用函数 adj_rand_index()，您可以检查这是否导致 ARI 等于 0。
#data first tree (n=55)
data_1 <- ari[[1]][[4]]
#predicted node first iteration
data_1$prednode_1 <- predict(ari[[1]][[1]], type = "node")
#Cross table
with(data_1, table(prednode_1 ,id))
## id
## prednode_1 1 2 3
## 1 37 7 11
#adj_rand_index
ari[[1]][[3]]
第二棵树 (n=87)：只发现一个断点
这个案例很有趣，因为它部分地识别了树的结构（即缺少 z1 上的中断）。
##### Second Tree (n=87): Extra partition in node[5] ####
ari[[2]][[1]]
# Fitted party:
# [1] root
# | [2] z2 <= 0.29288: n = 57
# | x(Intercept) xx
# | 0.133293 1.082701
# | [3] z2 > 0.29288: n = 30
# | x(Intercept) xx
# | 0.2598309 -1.8014133
#
# Number of inner nodes: 1
# Number of terminal nodes: 2
# Number of parameters per node: 2
# Objective function: 122.0116
此外，我们可以检查当交叉预测节点和真实节点时，我们看到即使在这棵不完美的树中，一些观察结果也符合标准。这意味着 57 个观测值正确分配给了第一个节点，9 个观测值被正确分配给了第二个分支。最后，30 因为最后一个节点根本没有被识别而被错误分配。这导致 ARI 等于 0.8577366，这是第一棵树的巨大改进。
#data second iteration (n=87)
data_2 <- ari[[2]][[4]]
#predicted node first iteration
data_2$prednode_2 <- predict(ari[[2]][[1]], type = "node")
#Cross table
with(data_2, table(prednode_2 ,id))
# id
# prednode_2 1 2 3
# 2 57 0 0
# 3 1 9 20
#adj_rand_index
ari[[2]][[3]]
# > ari[[2]][[3]]
# ARI
# 0.8577366
第三棵树（n=123）：结构正确但性能较差
这个案例很有趣，因为它确实恢复了真实的树结构，但它的性能比仅部分识别其结构的后三个更差。
##### Third Tree (n=123): Correct structure but poor performance ####
ari[[3]][[1]]
# Fitted party:
# [1] root
# | [2] z2 <= 0.07319: n = 60
# | x(Intercept) xx
# | -0.1723388 1.1071878
# | [3] z2 > 0.07319
# | | [4] z1 <= -0.35485: n = 22
# | | x(Intercept) xx
# | | -0.7166565 -0.6791717
# | | [5] z1 > -0.35485: n = 41
# | | x(Intercept) xx
# | | 0.7096033 -0.8605967
#
# Number of inner nodes: 2
# Number of terminal nodes: 3
# Number of parameters per node: 2
# Objective function: 156.4397
下面我们可以看到，当我们对预测节点和真实节点进行交叉表时，我们看到 16 (10 + 6) 个观测值被错误分类，这一事实导致 ARI 为 0.6117612.
#data third iteration (n=123)
data_3 <- ari[[3]][[4]]
#predicted node first iteration
data_3$prednode_3 <- predict(ari[[3]][[1]], type = "node")
#Cross table
with(data_3, table(prednode_3 ,id))
# id
# prednode_3 1 2 3
# 2 60 0 0
# 4 6 16 0
# 5 10 0 31
#adj_rand_index
ari[[3]][[3]]
# > ari[[3]][[3]]
# ARI
# 0.6117612
第四棵树：第四棵树 (n=199)：节点 [5]
处的额外叶子
这里识别的树与原来的有偏差，因为它多了一片node[5]的叶子，这在真实数据中是不存在的。
##### Forth Tree (n=199): Extra leaf at node[5] ####
ari[[4]][[1]]
# Fitted party:
# [1] root
# | [2] z2 <= -0.19806: n = 79
# | x(Intercept) xx
# | 0.06455217 1.51512672
# | [3] z2 > -0.19806
# | | [4] z1 <= -0.27127: n = 44
# | | x(Intercept) xx
# | | -0.4863122 -0.3860951
# | | [5] z1 > -0.27127
# | | | [6] z2 <= 0.17481: n = 23
# | | | x(Intercept) xx
# | | | -0.1335096 0.2046050
# | | | [7] z2 > 0.17481: n = 53
# | | | x(Intercept) xx
# | | | 1.0868488 -0.0290925
#
# Number of inner nodes: 3
# Number of terminal nodes: 4
# Number of parameters per node: 2
# Objective function: 282.6727
这里真实节点和预测节点的交叉表计数很有趣，因为节点 [6] 和 [7] 在真实数据中不存在，但它们得到的观察结果应该是，例如，分配给节点 [1]（分别为 23 和 7 观察值。）这种错误分配将 ARI 指数减少到 0.4649789。
#data forth iteration (n=199)
data_4 <- ari[[4]][[4]]
#predicted node first iteration
data_4$prednode_4 <- predict(ari[[4]][[1]], type = "node")
#Cross table
with(data_4, table(prednode_4 ,id))
# id
# prednode_4 1 2 3
# 2 79 0 0
# 4 16 27 1
# 6 23 0 0
# 7 7 0 46
#adj_rand_index
ari[[4]][[3]]
# ARI
# 0.4649789
第五棵树 (n=667)：节点 [2]
处的额外叶子
这是结构不正确的树的另一个示例，其中额外的叶子（基于 z5 上的分区是不正确的（！））附加到节点 [2]。
##### Fifth Tree (n=667): Extra leaf at node[2] ####
ari[[5]][[1]]
# Fitted party:
# [1] root
# | [2] z2 <= 0.28476
# | | [3] z5 <= 0.76285: n = 322
# | | x(Intercept) xx
# | | -0.1322881 0.9535337
# | | [4] z5 > 0.76285: n = 96
# | | x(Intercept) xx
# | | 0.1686863 1.3878776
# | [5] z2 > 0.28476
# | | [6] z1 <= -0.32001: n = 89
# | | x(Intercept) xx
# | | -0.9139858 -0.7957158
# | | [7] z1 > -0.32001: n = 160
# | | x(Intercept) xx
# | | 0.7661154 -0.8656553
#
# Number of inner nodes: 3
# Number of terminal nodes: 4
# Number of parameters per node: 2
# Objective function: 927.9088
来自预测节点和正确节点的交叉表向我们展示了实际上属于第一个节点 [1] 的大多数观察值 (322) 被分配给了预测节点 [3].最后，这种糟糕的结构导致 0.6932132.`
#data third iteration (n=667)
data_5 <- ari[[5]][[4]]
#predicted node first iteration
data_5$prednode_5 <- predict(ari[[5]][[1]], type = "node")
#Cross table
with(data_5, table(prednode_5 ,id))
# id
# prednode_5 1 2 3
# 3 322 0 0
# 4 96 0 0
# 6 0 89 0
# 7 3 3 154
#adj_rand_index
ari[[5]][[3]]
# ARI
# 0.6932132
第六棵树（n=5000）：黄金树！
这棵最终树完美地恢复了数据，无论是在树结构中还是在将观察值分配给每个叶子上。
##### Sixth Tree (n=5000): Extra leaf at node[2] ####
ari[[6]][[1]]
# Fitted party:
# [1] root
# | [2] z2 <= 0.29971: n = 3187
# | x(Intercept) xx
# | -0.008719923 1.022232280
# | [3] z2 > 0.29971
# | | [4] z1 <= -0.30286: n = 609
# | | x(Intercept) xx
# | | -0.9488846 -0.9813765
# | | [5] z1 > -0.30286: n = 1204
# | | x(Intercept) xx
# | | 1.0281410 -0.9565637
#
# Number of inner nodes: 2
# Number of terminal nodes: 3
# Number of parameters per node: 2
# Objective function: 6992.848
这里我们可以从预测节点和真实节点的交叉表中看到，它完美地分配了每个观测值，它确实属于它，导致 ARI 等于 1 .
#data sixt iteration (n=5000)
data_6 <- ari[[6]][[4]]
#predicted node first iteration
data_6$prednode_6 <- predict(ari[[6]][[1]], type = "node")
#Cross table
with(data_6, table(prednode_6 ,id))
# id
# prednode_6 1 2 3
# 2 3187 0 0
# 4 0 609 0
# 5 0 0 1204
#adj_rand_index
ari[[6]][[3]]
# ARI
# 1
结论。
可以从上图中恢复一些重要的要点。
1.- ARI 可用于评估预测树的相似程度，该预测树在数据生成过程中可能具有与真实树非常不同的结构。
2.- 恢复树的正确结构不会导致 ARI 等于一。
3.- 不正确的树不一定会使 ARI 等于零。
最终模拟。
总而言之，这是一个小型模拟，可以了解 ARI 指数在样本量增加时的表现。
### Final simulation
n_levels <-seq(from= 10, to= 2000, by= 5)
ari <- lapply(n_levels, sim_ari)
ari_models<- function(i){
ari <- ari_sim[[i]]$ari
n <- nobs(ari_sim[[i]]$ols_mob)
return(
list(ari = ari , n= n )
)
}
ari_n_list <- lapply(1:length(ari_sim), ari_models)
df <- data.frame(matrix(unlist(ari_n_list), nrow=length(ari_n_list), byrow=T))
colnames(df) <- c("ARI" , "N")
library(ggplot2)
ggplot(df, aes(N)) +
geom_line(aes(y = ARI, colour = "ARI"))