<p>First, I would recommend to use the <code>lmtree()</code> function and not vanilla <code>mob()</code>. The former is faster, comes with better defaults for printing and plotting, and has more options for predictions.</p>
<p>Second, I recommend that you consult the <code>vignette(&quot;partykit&quot;, package = &quot;partykit&quot;)</code> which explains how <code>party</code> objects are constructed and which classes and methods are involved.</p>
<p>Third, to determine which variable (if any) was used for splitting in the root node it is probably of interest to extract the results from all parameter instability tests. There is a dedicated <code>sctest()</code> (structural change test) method for obtaining this:</p>
<pre><code>library(&quot;strucchange&quot;)
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
</code></pre>
<p>The corresponding <code>partysplit</code> object for the <code>$split</code> (if any) in the root <code>$node</code> is probably easiest to extract manually:</p>
<pre><code>results[[1]]$node$split
## $varid
## [1] 4
## 
## $breaks
## [1] 0.0002853492
## 
## $index
## NULL
## 
## $right
## [1] TRUE
## 
## $prob
## NULL
## 
## $info
## NULL
## 
## attr(,&quot;class&quot;)
## [1] &quot;partysplit&quot;
results[[2]]$node$split
## NULL
</code></pre>
<p>The variable id pertains to the order of the variables in:</p>
<pre><code>names(results[[1]]$data)
## [1] &quot;y&quot;  &quot;x1&quot; &quot;x2&quot; &quot;z1&quot; &quot;z2&quot;
</code></pre>
<p>Finally, as for the question <em>what</em> to evaluate: Everything depends on identifying the correct variable for splitting. If this is done correctly, then the split point estimates converge fast to the true values, and then the parameter estimates also converge. See for example our recent arXiv paper (<a href="https://arxiv.org/abs/1906.10179" rel="nofollow noreferrer">https://arxiv.org/abs/1906.10179</a>) which contains a larger simulation study and also provides replication code.</p>
<p>Therefore, typically, I either evaluate the probability of selecting the correct variable in the first split. Or alternatively I look at the RMSE of the estimated vs.true coefficients <em>for each observation</em>.</p>
<p><em><strong>Update:</strong></em> Beyond the root node you can use <code>nodeapply()</code> to extract information from various nodes. However, I do not recommend to evaluate all splits because it becomes increasingly difficult to match which estimated split matches which of the true splits. Instead, we often assess how similar the fitted partition is compared to the true partition, e.g., using the adjusted Rand Index. Again, you can find an example for the in the arXiv paper mentioned above.</p>
Montecarlo Simulation using mob() algorithm (partykit package) to recover the count of the correctly identified models
<p>I am using the <code>mob()</code> function from <code>partykit</code> package and I am getting some problems when parsing the obtained model.</p>
<p>My main aim is to check <em>approximately</em> how large a sample size needs to be in order to detect the real
structure of a data generating process (DGP) when breaks are present.</p>
<p>The code below performs a Montecarlo simulation of data with breaks and tries to identify if the break was captured by the M-fluctuation test or not.</p>
<p>More specifically, I want to make a count on the number of times over the total number of simulations (<code>nreps</code>) that the model actually captures the DGP, conditional on a fixed sample size (<code>N</code>) to get a feeling of how many data should I need to capture the real DGP.</p>
<p>At the end of the code, you can see the list that I get out of my simulations. The problem is that I cannot recover the information displayed on the console.</p>
<p>Additionally, I have some doubts about how to make the count of the &quot;correctly identified models&quot;. For now, what I have in mind is to count as positive if the model has a break into the correct variable (<code>z1</code>) at the specified region (<code>z1&gt;0</code>) with some tolerance into the break region for example if the break is at <code>z1&gt;0.1</code> or <code>z1&gt;-0.4</code> it is also valid as a positive for me. Therefore, is there any simple way of counting the models that have the stated characteristics?</p>
<p>I hope my example is clear enough for you to help me out. Thank you a lot in advance.</p>
<h3>Simulate model</h3>
<ol>
<li>Ingredients to generate the DGP.</li>
</ol>
<pre><code>library(&quot;partykit&quot;)
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 &lt;- 300
#coeficients
betas &lt;- list()

betas$b0 &lt;- 1

betas$b1_up   &lt;- 2.4
betas$b1_down &lt;- 2

betas$b2_up   &lt;- 2.4
betas$b2_down &lt;- 2

#mob() ingredients
ols_formula &lt;- y ~ x1 + x2 | z1 + z2
# the &quot;&quot;0 +&quot;&quot; ---&gt;  supress the 'double' interecept
ols &lt;- function(y, x, start = NULL, weights = NULL, offset = NULL, ...) {lm(y ~ 0 + x)} 
</code></pre>
<ol start="2">
<li>Function that generates the data and fit OLS using mob algorithm.</li>
</ol>
<pre><code>reg_simulation_mob &lt;- function(...){
  #data
  data &lt;- data.frame(
    x1 = rnorm(N),
    x2 = rnorm(N),
    z1 = rnorm(N),
    z2 = rnorm(N),
    e  = rnorm(N))
  
  #dependent variable
  data$y &lt;-   betas$b0 + with(data,  ifelse(z1&gt;0,
                         betas$b1_up   * x1 + betas$b2_up   * x2 ,
                         betas$b1_down * x1 + betas$b2_down * x2 ) 
                        + e )
  
  
    #Estimate mob()-OLS
  ols_mob &lt;- mob(ols_formula,
                 data = data, 
                 fit = ols)
  
  # return(ols$coefficients)
  return(ols_mob)
}
</code></pre>
<ol start="3">
<li>Montecarlo simulation (only 2 trials) of the above-described setup.</li>
</ol>
<pre><code># N repetitions
nreps &lt;- 2

## Parallel version
results  &lt;- future_lapply(1:nreps, reg_simulation_mob, future.seed = 1234L)

</code></pre>
<h3>Obtained result</h3>
<p>As you can see below in the first trial (<code>results[[1]]</code>), the model finds the correct break but in the second (<code>results[[2]]</code>) it fails to find it.</p>
<pre><code>&gt; results
[[1]]
Model-based recursive partitioning (ols)

Model formula:
y ~ x1 + x2 | z1 + z2

Fitted party:
[1] root
|   [2] z1 &lt;= 0.00029: n = 140
|       x(Intercept)          xx1          xx2 
|          0.9597894    1.7552122    2.1360788 
|   [3] z1 &gt; 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

</code></pre>
<p>In the picture below, you can observe the structure of the list <code>results</code>, where I cannot find the information displayed on the console (i.e. number of nodes, parameters, threshold values, etc..)
<a href="https://i.stack.imgur.com/lhz7Q.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/lhz7Q.png" alt="enter image description here" /></a></p>