Montecarlo模拟使用mob()算法(partykit包)恢复正确识别模型的计数
我正在使用来自mob()包的partykit函数，在解析所获得的模型时会遇到一些问题。
我的主要目的是检查一个样本大小需要有多大，以便在出现中断时检测数据生成过程(DGP)的实际结构。
下面的代码对带有中断的数据执行Montecarlo模拟，并试图确定该中断是否被M-波动测试捕获。
更具体地说，我希望计算模型实际捕获DGP的模拟总数(nreps)的次数，这取决于固定的样本大小(N)，以获得捕获实际DGP需要多少数据的感觉。
在代码的末尾，您可以看到我从模拟中得到的列表。问题是我无法恢复显示在控制台上的信息。
此外，我对如何计算“正确识别的模型”有一些疑问。现在，我要考虑的是，如果模型在指定的区域(z1>0)有一个正确的变量(z1>0)，并且对中断区域有一定的容忍度--例如，如果中断位于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)} - 函数，该函数生成数据并使用暴徒算法拟合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的结构，其中我找不到显示在控制台上的信息(即节点数、参数、阈值等)。
回答 1
Stack Overflow用户
发布于 2020-11-23 21:31:33
首先，我建议使用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根partysplit中的$split (如果有的话)对应的$node对象可能最容易手动提取：
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系数的均方误差。
更新:根节点之外的，您可以使用nodeapply()从各个节点提取信息。然而，我不建议评估所有的分裂，因为它变得越来越难匹配哪一个估计分裂匹配真正的分裂。相反，我们经常评估拟合的分区与真正的分区相比有多相似，例如，使用调整后的兰德指数。同样，您可以在上面提到的arXiv论文中找到一个示例。
https://stackoverflow.com/questions/64970510
复制相似问题
腾讯云开发者
