[R] party for prediction [REPOST]

[Achim Zeileis]

On Sun, 14 Oct 2012, Ed wrote:

> This was an exceptionally helpful answer, I can only thank you again.
> I have plenty of avenues ahead where I was worried before I was
> getting trapped in a dead end. If all else fails, the idea of using
> anova is brilliant. Thank you!

No prob. Just for the record, in case anybody else wonders: MOB is based 
on score-type tests (or Lagrange multiplier type tests) to decide which 
variable should be used for partitioning in the next step. In principle, 
likelihood-ratio-type tests could also be used but require fitting models 
under all conceivable splitted models. Especially for all two-way splits 
in numerical partitioning variables, this can be very costly and hence the 
computationally cheaper score tests are used (which only require model 
fitting once per node). In your case with only categorical partitioning 
variable, computing the likelihood ratio test is possible and only a 
little bit more coputationally expensive. (And you indicated that you do 
not worry about that.)

To build your own new tree function, it may be useful to employ the 
"partykit" infrastructure. This has very useful base infrastructure for 
representing the tree part but has not yet be customized for conveniently 
accessing model fits in the leaves. It's on our to-do list though...

Best,
Z

> Ed
>
> On 14 October 2012 18:36, Achim Zeileis <Achim.Zeileis at uibk.ac.at> wrote:
>> On Sun, 14 Oct 2012, Ed wrote:
>>
>>> First up, thanks hugely for your response. I've been beating my head
>>> against this!
>>>
>>> On 14 October 2012 16:51, Achim Zeileis <Achim.Zeileis at uibk.ac.at> wrote:
>>>>
>>>> I'm not sure what you mean by "integral vector". If you want to apply the
>>>> approach to hundreds of thousands of observations, I gues that these are
>>>> categorical (maybe even binary?) but maybe not...
>>>
>>>
>>> I'm sorry I can't go into the details of the data, I would if I could.
>>> z are categorical variables represented as integers, mostly ordered,
>>> but not all. I've tried fitting them as integers, as well as ordered,
>>> but O don't think it made a huge difference.
>>
>>
>> The tests performed for categorical partitioning variables are rather
>> different from the tests for numerical partitioning variables. If all of the
>> variables are categorical, this may not be immediately obvious, but the
>> factor coding should be more appropriate (especially if the number of levels
>> is small or moderate).
>>
>>
>>>> If I recall correctly, we kept linearModel as simple as we did to save as
>>>> much time as possible. This can be particularly important when one of the
>>>> partitioning variables has many possible splits and the linearModel has
>>>> to
>>>> be fitted thousands of times.
>>>
>>>
>>> I can appreciate that, but maybe having an alternative linearModel which
>>> will predict when the fit is degenerate would be worth including? I'm happy
>>> to contribute what I have, although it's pretty obvious stuff (and probably
>>> done suboptimally since I'm not much of an R coder at this point). For me at
>>> least, even with huge datasets, the speed of party is quite good; it's
>>> getting a better result that's the problem.
>>
>>
>> As I explained in my last e-mail. In your situation this does not solve the
>> problem completely because subsequent the tests are also not adapted to
>> this. Setting the empirical estimating functions to zero for non-identified
>> coefficients might alleviate the problem but is not really a clean solution.
>>
>>
>>>> Also, mob() assesses the stability of all coefficients of the model in
>>>> all
>>>> nodes during partitioning. If any of the coefficients is not identified,
>>>> this would have to be excluded from all subsequent parameter stability
>>>> tests
>>>> in that node (and its child nodes). This is currently not provided for in
>>>> mob().
>>>
>>>
>>> Would pretending the coefficients were fit at 0 fool mob into doing
>>> something moderately meaningful here?
>>
>>
>> The coefficients are not looked at during fitting, only the estfun(). This
>> would have to be set to 0.
>>
>>
>>> If not, I would try to hack the code, but I'm honestly at something of
>>> a loss as to how to modify it and feed the results back into my
>>> interpreter. I have bytecode installed; I downloaded the source, but I
>>> haven't squared the circle of modifying the source and installing the
>>> result. I will check out the docs on writing extensions you suggest.
>>
>>
>> Writing (or modifying) R packages and installing them under Windows is
>> pretty standard and well documented. The pointers I gave you should
>> hopefully get you started.
>>
>>
>>>>> The second problem I have is that I get "Cholesky not positive definite"
>>>>> errors at some nodes. I guess this is because of numerical error and
>>>>> degeneracy in the covariance matrix? Any thoughts on how to avoid having
>>>>> this happen would be welcome; it is ignorable though for now.
>>>>
>>>>
>>>> This comes from the parameter stability tests and might be a result of an
>>>> unidentified (or close to unidentified) model fit.
>>>
>>>
>>> This is a great help to know. I improved my results quite considerably
>>> with aggressive scaling of everything (scaling the response and all
>>> the predictors to lie between 0 an 1). That deepened my tree by a
>>> factor of two or so (say depth 3 to 7) and improved the quality of fit
>>> substantially. Is there any way I can engage a more numerically robust
>>> Cholesky in mob?
>>
>>
>> No, I don't think that this is conceivable with the way this is implemented
>> at the moment. Instead of the currently implemented tests, one could in
>> principle use the likelihood ratio test which is invariant against parameter
>> transformations and doesn't need the covariance matrix.
>>
>>
>>>> With hundreds of thousands of observations, you would need some
>>>> additional
>>>> pruning strategy anyway. Significance test-based splitting will probably
>>>> overfit because tiny differences in the coefficients will be picked up at
>>>> such large sample sizes.
>>>
>>>
>>> I'm okay with overfitting, honestly. At the moment it is underfitting
>>> by quite a large amount I think (the quality of the predictions on the
>>> training set is not very high). The problem really is there is so much
>>> going on the data, but the "noise" level is probably very low. I
>>> wouldn't be surprised if my data was accurate to 5 or 6 s.f.
>>
>>
>> OK. Maybe with so little noise some other splitting strategy (not based on
>> significance tests) would be better?
>>
>>
>>>> Furthermore, computationally the extensive search over all possible
>>>> splits
>>>> might be too burdensome with this many observations.
>>>
>>>
>>> I have plenty of compute time/power and RAM, though R seems to be
>>> running single threaded. But even on a few million observations, it's
>>> still pretty fast and doesn't use more than 30 or so gig of memory. If
>>> it takes a day and requires 150gig of RAM, that is absolutely fine,
>>> even over that would be viable though less optimal.
>>
>>
>> With this setup, you may consider writing your own partitioning algorithm
>> using the same type of ideas as MOB. Instead of using the parameter
>> stability tests, you could use plain likelihood ratio tests or ANOVAs to
>> decide about splitting further.
>>
>> If d is the data in the current node, y and x are response/regressors, and
>> f1 to fn are your categorical partitioning variables, you could do
>>
>> m0 <- lm(y ~ x, data = d)
>> m1 <- lm(y ~ f1 * x, data = d)
>> ...
>> mn <- lm(y ~ fn * x, data = d)
>>
>> and then
>>
>> anova(m0, m1)
>> ...
>> anova(m0, mn)
>>
>> and then choose the lowest p-value for splitting. You could then also
>> parallelize the computations in the daughter nodes.
>>
>> None of this is readily coded in party though...
>>
>>
>>>> Hence, using some subsampling strategy might not be the worst thing.
>>>
>>>
>>> I've tried this at various degrees, but the data is really very
>>> complicated with not a lot of error. I'm trying to encourage party to
>>> fit more closely, which I thought more data might encourage. At the
>>> moment I'm a long way from a clean fit. I have subsampled at various
>>> levels down to 1%, and although that increases the depth of the tree
>>> and quality of fit, it still doesn't give a very good quality fit and
>>> can encourage it to overlook obvious aspects of the training set.
>>>
>>>>> I guess what I really want to know is:
>>>>> (a) has anyone else had this problem, and if so how did they overcome
>>>>> it?
>>>>
>>>>
>>>>
>>>> We have had non-identified model fits in binary GLMs (with quasi-complete
>>>> separation) where we then set estfun() to all zero so that partitioning
>>>> stops. But I don't think that such a strategy helps here.
>>>
>>>
>>> I've considered using rpart() to partition into cells of constant
>>> gradient, then fitting linear models myself to the cells. This is my
>>> next thought. I'm pretty sure partitioning over linear regression is
>>> the way forward for the data we have. I tried mars and glm but there
>>> are good reasons to think they're less reasonable, even though the fit
>>> wasn't particularly poor. I'm not particularly wedded to party's
>>> approach except that it looked like it immediately returned what we
>>> needed, and with some degree of "optimality" into the bargain.
>>
>>
>> Possibly boosting linear models may be another route to go.
>>
>>
>>>>> (b) is there any way to get a line or stack trace out of a try()
>>>>> without source modification?
>>>>
>>>>
>>>> Not sure, I don't know any off the top off my head.
>>>
>>>
>>> I guess I really will have to bite the bullet and try to figure out
>>> how to install modified libraries. Thanks.
>>
>>
>> I think that's easier compared to solving the conceptual problems (thinking
>> about how to best partition degenerate linear models etc.).
>> ;-)
>>
>> Good luck!
>> Best wishes,
>>
>> Z
>>
>>>>> (c) failing all of that, does anyone know of an alternative to mob
>>>>> that does the same thing; for better or worse I'm now committed to
>>>>> recursive partitioning over linear models, as per mob?
>>>>
>>>>
>>>> If your partitioning variables are particularly simple (e.g., all binary)
>>>> you could exploit that and it may be easier to write a custom function
>>>> for
>>>> your particular data. Then likelihood-ratio tests (rather than LM-type
>>>> tests) would also be easier to apply in case of unidentified parameters.
>>>>
>>>> But if there are partitioning variables with different measurement
>>>> scales,
>>>> then this will not be that simple...
>>>
>>>
>>> Unfortunately each partitioning variable is essentially a state
>>> indicator, taking values say 0,...,R where R is different for each
>>> component. I'm not a stats expert either; I've spent some time with
>>> the party manuals and papers, but I wouldn't be confident of
>>> implementing something like it in the time available to me (though if
>>> I have to I will, but that wouldn't be a good situation to be in).
>>>
>>>>> (d) failing all of this, does anyone have a link to a way to rebuild, or
>>>>> locally modify, an R package (preferably windows, but anything would
>>>>> do)?
>>>>
>>>>
>>>>
>>>> Have a look at the "Writing R Extensions" manual and the R for Windows
>>>> FAQ.
>>>
>>>
>>> Will do.
>>>
>>> Thank you very much for your responses, I really appreciate it.
>>>
>>> Best wishes,
>>>
>>> Ed
>>>
>>
>