Random forest prediction intervals - sum of new observations

Question

I am relatively new to random forests but have been experimenting with the various R packages available. So far so good. But for my application I am modelling a response variable with number of predictors and that is fine, but what I really need to do is use the resulting model to predict values for new data as it comes in. This is fine, I receive a new set of predictor variables and predict the values for the new unobserved response variable.

However, the problem is that I am not really interested in each of these individual new variables, but the sum of them. That is also fine, I just sum up each of the predicted values and that is the parameter estimate I am interested in. However, I would really also like to estimate a prediction interval associated with this. I can see that there are various methods and R packages to do this on the individual variables e.g. the quantile regression forest approach of Meinhausen (and package quantregForest), alternative approaches based on bootstrapping or jackknife. But have yet to find an approach that allows the prediction interval to be made on a derived parameter such as the sum over the individual variables?

I have found some very interesting comments on previous questions about random forest prediction intervals but no obvious information on how to do it for a derived parameter. Perhaps I am being naive and have missed it?

The closest I have found is this one: Sum of Random Forest prediction intervals?

which is basically my exact question, but it was not actually answered adequately.

If anyone has any ideas or could point me towards some discussion of this matter, I would be very appreciative.

Answer 1

Your problem is that you would need to know the correlation between the predictions. That doesn't seem like an easy thing to get. If you had it, you could use that when have a mean vector and its covariance matrix, then you can trivially calculate the Variance or SD of the sum of the components of the vector. Using the delta method you can calculate an approximate Variance or SD for more general transformations.

A rough hack that might work would be to empirically look at this on new data on which the RF was not trained and there might be some hope a single correlation matrix could approximately describe it. The problem is that the correlation matrix might depend on the data constellation (so it might be some complex function and vary substantially based on the input data) or that while it looks straightforward for the example data, you might later get new data where this is not a good approximation.

Most of the alternatives I can think of would involve moving to a different type of model. E.g. if you fit some Bayesian regression of some form (interactions and splines can often handle many of the things that make people think they need a machine learning model), it would be vary easy to form all sorts of transformations of predictive distributions. If it must be some complex machine learning model multi-output neural networks would be able to output multiple things and one could then do quantile regression approach to getting uncertainty on those things.

Answer 2

If you have enough data, I would use a conformal prediction-inspired approach for this task.

Use time series cross-validation to backtest how large your error typically is when predicting out N days, aggregating the predictions, and comparing against the aggregated actuals.

Let's say you're making daily predictions and interested in 30-day aggregated predictions and prediction intervals. In detail, the process would be:

1. Set aside a 30-day validation period completely after a training dataset.
2. Train your RF on the training data.
3. Predict on the 30-day validation period, aggregate to a single value.
4. Compare to the aggregated actual values and calculate the error
5. Step forward k days and repeat steps 1-4 as often as your data allows. k should ideally be 30, but if you don't have a lot of data, perhaps allowing some overlap and using 10 or 15 would be fine.
6. Take the absolute values of all the errors, and use quantiles to calculate prediction intervals. For example, if you wanted a 95% prediction interval, you would take the 0.95 quantile of the absolute errors, add and subtract that from the predicted value to get upper and lower prediction interval values.

If the target variable is skewed, it may work better to not take the absolute value of the errors, and instead take the 0.025 and 0.975 quantiles of the signed errors to produce a 95% prediction interval.

If you want to read more about conformal prediction generally, I found this resource helpful: https://www.stat.berkeley.edu/~ryantibs/statlearn-s23/lectures/conformal.pdf