Distribution of the difference between two predictions made by the same model

Question

I have a database of 200 trees grown in different conditions, since it's very expensive to get more points I want to predict with a machine learning model the growth potential of tree seedlings:

• X variables: type of soil, exposure, species, etc.
• Y variable: growth potential

to help decision-making on which of two species should I plant on a new location: I need to compare the predictions returned by my model, given that X and Y are both subject to uncertainty (I have a known error in the measurement of my X and Y variables (potentially a different error attached to each data point)). For instance, if Y(A)=60 and Y(B)=59, I would like to compute the probability that species A will indeed give a better growth than B, the probability that it will give the inverse outcome, or a similar outcome (maybe the distribution of the difference of prediction?)

(What I have done so far: I thought of doing some bootstrap on the points, compute again the model, and see in how many cases Y(A)>Y(B). But still, this doesn't integrate the uncertainty of the X variables so the result is far from correct (too confident imo). I also tried to plot the differences between all couples of points, but I don't know if this is valid. Eventually, I would prefer if the method does not involve any refitting of the model because for operational use it won't be possible.)

Any help on a method would be appreciated. Thanks!

EDIT precisions on the model type and why I don't want to refit it. The model can be either of 2 models: statistical model (e.g.: ML) or physical model based on simulation and very heavy to compute. The second one can be obviously subject to bias because it's a simulation. Let's assume I can have a measure of uncertainty for each X point, at least.

Answer 1

You say "I want to estimate distribution". But you don't have enough data to estimate this distribution. In this case, you have to make simplifying prior assumptions.

You say that prediction $\hat{Y}(A)$ of a model is uncertain: that is, the real growh potential $Y(A)$ might be different by $\varepsilon=Y(A)-\hat{Y}(A)$. Conditional distribution $P(\varepsilon|A)$ is what you are trying to estimate. Having this distribution, you can predict distribution of difference $Y(A)-Y(B)=x$ as convolution $p(x) \sim \int p(z|A)p(z-x|B)dz$. This already requires assumption that different trees are independent of each other, but this assumption is weak enough.

Now the stronger assumtptions which can make it as simple as possible:

• assume $P(\varepsilon|A)$ is normal - then it will be easy to integrate it. For many natural proccesses, this assumption is realistic
• assume $P(\varepsilon|A)$ does not depend on $A$ at all - then you have only to estimate two parameters, its mean and variance. If your data is limimted, this is enough.
• assume your model is on average unbiased (most of them are on IID data), then you have to estimate only variance of $\varepsilon$ . If you think that the model is biased, you can estimate the bias as well, but in the distribution of difference it will cancel out.

This variance should ideally be estimated out-of-the-train-sample. Cross validation is ideal for this, but if you can fit the model only once (why?), then you can split your 200 trees into train and test sets e.g. of 140 and 60 trees, fit the model on the first set, and estimate variance of the error on the second one. Let it be $\sigma^2$

Now when you have your $P(\varepsilon)\sim\mathcal{N}(0, \sigma^2)$, you can show that distribution of difference is also normal, $P(Y(A)-Y(B))\sim\mathcal{N}(\hat{Y}(A)-\hat{Y}(B), 2\sigma^2)$. Voila!

If you feel that the homoschedasticity assumption ($Var(\varepsilon|A)=\sigma^2$ is too strong, you have to somehow estimate this variance. It can be done non-parametrically (e.g. with bootstrap), but you say it is not an option. In this case, you can come up with a simple parametric model for $\sigma(A)$ - for example, $\sigma(A)=\alpha+\beta\hat{Y}(A)$, or $\sigma(A)=\alpha\hat{Y}(A)^\beta$. Then you can estimate parameters of this model instead of estimating a single $\sigma$. If you still believe that $(\hat{Y}(A)-Y(A))\sim\mathcal{N}(0, \sigma(A)$), then $Y(A)-Y(B)$ has the same mean as before and variance $\sigma^2(A)+\sigma^2(B)$ and is still normal.

Answer 2

EDIT:

For linear regression, if the variability in the errors of observation of X are small relative to the range of X then it is not much of a concern when fitting model. Therefore, you may not have to adjust your fitting procedure to account for error in inputs.

You are interested in prediction error, but you have few data points and can only fit the model once. You can't divide your dataset into training and test sets because you do not have enough data. You can't perform cross-validation because you can't refit the model several times. Therefore, I believe you can not estimate the prediction error.

If you can estimate the prediction error through cross validation, then you can use a MCS to get the distribution of error at two prediction locations and therefore the distribution of differences in predictions.

ORIGINAL:

I think you are interested in the prediction error. The simplest and most widely used method for estimating prediction error is cross-validation. You should be able to use cross-validation to get an idea of how accurate your model is for species A and B.

To incorporate input uncertainty you may consider using a Monte-Carlo Simulation where you generate model predictions with uncertainty estimates for a distribution of possible inputs. Then you can summarize the distributions of outputs for species A and B and compare.

However, the above approach will not account for model error that may arise due to not accounting for error in model inputs when fitting your model. You might consider "training with jitter" where you add some noise to your training data to improve generalization of model predictions. I am not very familiar with this approach. See "Machine Learning with Known Input Data Uncertainty Measure" by Czarnecki and Podolak for more information.

Answer 3

Cool problem!

You had a nice idea with the Bootstrap: $\frac{Y(A) > Y(B)}{N_{bootstrap}}$

I would argue that, yes, it does capture the uncertainty of the variables. If growth potential is all over the place for species A, and very consistent for B, then this should be reflected in the bootstrapped predictions.