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

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.

• Could you describe your models more clearly? If your models do not already incorporate the uncertainty in your dependent variables, it is difficult to reconcile the following goals: (1) capturing uncertainty in the underlying independent variables in your final predictions, and (2) avoiding refitting the models. Your bootstrap approach is generally sound, though. – mkt - Reinstate Monica May 18 '18 at 10:02
• @mkt Hi, I added last paragraph. Thanks, indeed I understand your point. I guess I'll have to make some assumptions about the prediction error at some point.? – agenis May 18 '18 at 15:01
• @mkt let me add that bootstrap has no meaning for, like, physical models. In the sense that their predictions do not depend on the calibration data... – agenis May 22 '18 at 12:55

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.

• Thanks @david Dale, for putting this into some formal equation. To answer remark, some models used are statistical, but some are purely physical relying on heavy/long computation so refitting a great number of times is excluded. And for those models, there is no guarantee about the bias (they are probably biased that's why i'm more interested in comparing predictions than the precise prediction value). I guess the hypothesis p(e|A)=Norm is OK but P(e|A)=P(e) is wrong for me (heteroscedastic). Can you explain more the convolution part, i've a hard time getting this. thanks – agenis May 18 '18 at 8:55
• About the convolution part, see Wikipedia: en.wikipedia.org/wiki/Convolution_of_probability_distributions. I made a mistake in the formula, let me fix it. – David Dale May 18 '18 at 9:16

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.

• thanks, the reference you give seems quite interesting, i will read it. Unfortunately, as i specified in the last paragraph, i want to have a result that does not require to fit again the model.. – agenis May 17 '18 at 19:02
• I edited my answer. If you have few data points and can only fit the model once then you may be in an unworkable situation. – Nat May 18 '18 at 16:40

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.

• Hi thanks for your answer! Indeed bootstrap could be an idea, but i'm afraid i'm forced to exclude it for 2 reasons: first, some of the databases contains only a few points, and bootstrapping might sometimees give 0% or 100% and that's obviously wrong. Secondly, operationnal platform where it will be implemented does not allow to recompute the model (which is the basis of bootstrap) – agenis May 7 '18 at 15:20
• Tricky. In that case, I think you need to estimate the conditional distributions of growth given A and growth given B. Then, find P(Y|A - Y|B > 0)  – timwiz May 7 '18 at 22:47
• Then again, you run into the same problem. Having only a few data points probably won't lead to very reliable distribution function estimates. – timwiz May 7 '18 at 22:49