# How can regression trees be fit in WinBUGS/OpenBUGS/JAGS?

There is an R package called BayesTree which can fit regression trees in Bayesian environment. However, this way only simple regression is possible. I would like to use regression trees as a part of a bigger hierarchical model instead of the simple GLM formulas (CARTs are proven to give better results in certain applications, see e.g. Hu et al 2011).

Is it possible (and how) to fit regression trees as a part of the model in WinBUGS/OpenBUGS/JAGS? Are there any such packages for these pieces of software?

• In Hu et al (2011) you can read: "A fully Bayesian simulation from the posterior distribution could have been implemented via a greedy search algorithm. However, currently this is computationally infeasible because the parameter space is large and has an inflexible hierarchical structure", so the answer is no. If you ask about other software like BUGS, there is Stan. However in your case you probably have to program the algorithm by hand from the scratch. – Tim Nov 19 '14 at 11:01
• @Tim, thanks. You mentioned Stan, do you think it can do it? – Curious Nov 19 '14 at 19:11
• Stan enables you to use control statements, for loops and do much more programming (BUGS does not), so if I would to try to implement what Tom Minka suggests below, I wouldn't even bother with BUGS and use Stan from the beginning. However, I didn't ever try this kind of models in Stan so sorry I am unable to tell you more. Try the Stan users mailing list: groups.google.com/forum/#!forum/stan-users – Tim Nov 20 '14 at 7:57

Here is a suggestion of how the Bayesian CART model from Hu et al (2011) might be implemented in WinBUGS. This method requires you to fix the tree depth in advance. Given a complete binary tree, any node at depth 1 (i.e. children of the root) can be identified by an indicator $b_1$. Any node at depth 2 can be identified by a pair of indicators $(b_1,b_2)$. And so on. Let $x[i,v]$ be the $v$th predictor variable in the $i$th instance. We can encode the path through the tree by thresholding $x[i,v]$ according to the parameters of the chosen node at each depth. This code works for trees of depth 3:

for(i in 1:n) {
b1[i] <- 1 + step(x[i,variable1] - split1)
b2[i] <- 1 + step(x[i,variable2[b1[i]]] - split2[b1[i]])
b3[i] <- 1 + step(x[i,variable3[b1[i],b2[i]]] - split3[b1[i],b2[i]])
y[i] ~ dnorm(mean[b1[i],b2[i],b3[i]], prec[b1[i],b2[i],b3[i]])
}


You just need to add appropriate priors for $variable$ (the splitting variable at each node), $split$ (the splitting threshold at each node), and $(mean, prec)$ of the leaves. There are 8 possible leaves, but they need not all be used. For example, if one of the splits is outside the range of the data then it effectively prunes away part of the tree. You could encourage such pruning with an appropriate prior on splits.

• I am afraid this wouldn't work. Are you sure you understand how tree regression works? The method will first split using a variable with highest impact - i.e. the most important variable. The method will automatically select the order of variables by importance. There is no fixed order, like "this is variable1, this is variable2...". Random forests are even more complicated, they build plenty of such trees. – Curious Nov 20 '14 at 10:37
• The model I wrote will automatically select the order of variables by importance, since it is inferring the values of variable1, split1, etc. – Tom Minka Nov 20 '14 at 18:14