How can I create a meaningful weighting for RMSE?

Question

Background

I should start off by saying I am not a mathematician and please excuse simple/stupid mistakes!

The goal of my exercise is to find the “best-fitting” model for the purpose of prediction. However, there are a couple of caveats:

1. I am only interested in predicting over the upper quartile of the response points;

2. under-predicting the value of the response is MUCH worse than over-predicting, and;

3. I have three specific values I wish to calculate – lets say 90, 99 and 99.99.

In this case I am working with support vector machines (SVMs) which my previous investigations have shown to be the best performing models out of a wide range of possible model types. For this particular exercise the array of possible models are SVMs with the same formula but different tuning parameters.

This is an example of some analogous data:

df <- structure(list(Covariate = c(100, 91.1242603550296, 86.9822485207101, 
100, 0, 0, 90.5325443786982, 95.8579881656805, 88.7573964497041, 
96.4497041420118, 82.2485207100592, 99.4082840236686, 99.4082840236686, 
98.8165680473373, 91.7159763313609, 59.1715976331361, 44.9704142011834, 
0, 100, 95.2662721893491, 100, 82.8402366863905, 7.69230769230769, 
81.6568047337278, 62.7218934911243, 97.6331360946746, 73.9644970414201, 
8.87573964497041, 0, 98.8165680473373, 78.1065088757396, 98.2248520710059, 
52.6627218934911, 96.4497041420118, 52.0710059171598, 0, 62.043795620438, 
84.6715328467153, 97.8102189781022, 4.37956204379562, 89.051094890511, 
99.2700729927007, 99.2700729927007, 97.0802919708029, 81.7518248175183, 
80.2919708029197, 90.5109489051095, 99.2700729927007, 96.3503649635037, 
0, 0, 94.8905109489051, 79.5620437956204, 67.8832116788321, 73.7226277372263, 
100, 97.0802919708029, 93.4306569343066, 86.8613138686131, 33.5766423357664, 
32.1167883211679, 46.7153284671533, 98.5401459854015, 95.6204379562044, 
86.1313868613139, 14.5985401459854, 92.7007299270073, 86.1313868613139, 
0, 77.3722627737226, 89.051094890511, 80.2919708029197, 98.1818181818182, 
96.3636363636364, 30.9090909090909, 0, 60.9090909090909, 100, 
0, 83.6363636363636, 88.1818181818182, 97.2727272727273, 0, 0, 
99.0909090909091, 100, 100, 91.8181818181818, 88.1818181818182, 
46.3636363636364, 50.9090909090909, 99.0909090909091, 97.2727272727273, 
100, 0, 92.7272727272727, 60.9090909090909, 90.9090909090909, 
57.2727272727273, 76.3636363636364, 94.5454545454545, 50, 98.1818181818182, 
16.3636363636364, 87.2727272727273, 92.7272727272727, 87.2727272727273, 
88.1818181818182, 10.7438016528926, 91.7355371900827, 98.3471074380165, 
60.3305785123967, 95.8677685950413, 0, 63.6363636363636, 71.900826446281, 
0, 74.3801652892562, 76.8595041322314, 0, 61.9834710743802, 0, 
0, 0, 84.297520661157, 47.1074380165289, 69.4214876033058, 97.5206611570248, 
100, 61.1570247933884, 90.0826446280992, 78.5123966942149, 10.7438016528926, 
100, 98.3471074380165, 100, 98.3471074380165, 93.3884297520661, 
90.9090909090909, 57.8512396694215, 57.8512396694215, 92.5619834710744, 
77.6859504132231, 69.4214876033058), Response = c(20, 14, 14, 
20, 0, 0, 14, 14, 14, 16, 10, 20, 20, 20, 16, 10, 10, 0, 16, 
16, 16, 10, 0, 12, 10, 12, 12, 0, 0, 20, 12, 16, 10, 12, 12, 
0, 14, 14, 16, 0, 14, 20, 16, 20, 14, 12, 12, 20, 20, 0, 0, 14, 
12, 10, 10, 20, 16, 16, 14, 10, 10, 10, 20, 16, 10, 0, 12, 12, 
0, 12, 16, 14, 16, 14, 0, 0, 12, 20, 0, 12, 14, 14, 0, 0, 20, 
20, 20, 14, 14, 10, 10, 20, 16, 16, 0, 12, 10, 10, 10, 16, 16, 
12, 20, 10, 12, 12, 16, 14, 0, 16, 20, 12, 14, 10, 10, 0, 0, 
12, 12, 10, 10, 0, 0, 0, 14, 12, 12, 20, 20, 14, 14, 14, 12, 
20, 20, 20, 16, 16, 14, 10, 10, 16, 16, 16)), row.names = 433:576, class = "data.frame")

To these ends @missuse has written some code to use leave-one-out cross-validation over the upper quartile of time points. At each iteration of the cross-validation I am currently calculating a weighted root mean square error (wRMSE), where I have made a simple (read crude) alteration which penalises model under-prediction 3× as harshly as over-predictions.

Here is the wRMSE function which I am feeding into the caret::train function:

wRMSESummary <- function (data,lev = NULL,model = NULL) {Res <- (data$obs - data$pred)
Res[Res < 0] <- Res[Res < 0] * 3 #Finishing wRMSE calculation
out <- mean(sqrt(Res^2), na.rm  = T)
names(out) <- "wRMSE"
out}

And here is the code I am using to carry out LOOCV over the upper range of the response which @missuse wrote:

`#Creating folds for upper quartile LOO-CV
    df %>%
    mutate(id = row_number()) %>% #create a column called id which will hold the row numbers
    filter(TrtTime >= unname(summary(df$Response)[5])) %>% #subset data frame according to your description 
  split(.$id)  %>% #split the data frame into lists by id (row number)
  map(~ .x %>% select(id) %>% #clean up so it works with indexOut argument in trainControl
        unlist %>%
        unname) -> dfFolds #Turning test sets into training sets
dfFolds <- lapply(dfFolds, function(x) setdiff(1:nrow(df), x))`

and finally the code for feeding both these into the caret::train function via caret::traincontrol:

#Fitting traincontrol
dfcontrol <- trainControl(index = dfFolds, method = "cv", summaryFunction = wRMSESummary)

#Fitting SVM
Poly.SVR2 <- train(log(Response + 0.0001) ~ Covariate+ I(Covariate^2), 
                   data = df,
                   tuneLength = 5, method = "svmPoly", 
                   trControl = dfcontrol)

Problem

I want to find a less arbitrary value than just multiplying the error of under-predictions by 3. Hopefully one that gives me some idea of the chance my model is under-predicting (if possible). By multiplying by a scalar as I outlined above, I believe (please correct me if I am wrong) that I will have diminishing returns as the value of the scalar increases- i.e. the PE will tend towards a particular model as the scalar increases. Or put another way the point estimates that I am interested in calculating (see point 3 in the background) will converge as my scalar tends towards infinity. If so, I was wondering if it would be possible to calculate the receiver-operator curve (ROC) that would characterise the loss in change and from here work out the area under the curve (AUC) which could in turn be used to determine what the optimal scalar would be, by say trying to account for 95% of the change in the fitted values across the upper quartile.

Does that even make sense? This is really out of my comfort zone but any guidance would be greatly appreciated!

I am completely open to alternative ideas and really looking for critiques as to whether or not this makes sense! If it does, would anybody be able to show me how I might go about determining this optimal value, preferably in R.

Answer 1

What you are doing is quantile prediction (or "forecasting", in a time series context, but the differences to the non-time-series case are irrelevant here): you are looking for a single number for each instance (a "point prediction") that will be larger than the actual outcome with a prespecified probability, e.g., 90%, 95% or 99%.

Thus, you are trying to predict a quantile of the (unknown) probability distribution of a future outcome. Since your target probabilities are typically larger than 50%, you are mainly interested in upper quantiles, or the right tail of the distribution.

Note that this unknown distribution may depend on covariates, which makes life a bit harder.

Now, the problem is that most common statistical or machine learning algorithms try to fit and predict a central tendency of this unknown distribution - often the expected value, but frequently the median. (The two may be different if the distribution is asymmetric.) They do this by minimizing an objective function that has its minimum at the expectation, e.g., the mean squared error (MSE). An additional problem is that this is rarely or never explicitly stated, so people are left wondering what to do in case they are looking for a quantile, as you are.

As far as I understand, your approach is two-pronged:

1. You penalize under-predictions more strongly. Unfortunately, as long as the predictions aim at the expected value, there is no connection between a "good" penalization scheme and getting at the quantile you actually want.
2. You take the point prediction your SVM (or whatever tool) outputs and multiply it by a factor. Again, this won't work very well, because in general, there is no simple multiplicative relationship between the expectation (which your tool outputs) and high quantiles.

There are two standard ways of getting at quantiles.

1. The first one is quantile regression and variants thereof. The idea is to use a target function other than the MSE, one that is minimized by the quantile you are looking for. Then a straightforward linear regression is run, but using this quantile target function. We have a quantile-regression tag. One consequence is that you will get a different model for each target probability (and you may run into inconsistencies - the model for a 90% target may output a larger quantile prediction than the model for a 95% target). In principle, you may be able to feed these target functions into other minimizers and get a quantile version of non-regression methods.

   Roger Koenker is pretty much the quantile regression guru. His 2005 textbook is called Quantile Regression, he coauthored an introductory article with Hallock (Koenker & Hallock, 2001, Journal of Economic Perspectives), and he wrote the quantreg package for R.

2. The alternative is to use a tool that can output a full predictive density, then take the appropriate quantile of this predicted density. For instance, straightforward linear regression via ordinary least squares will output an expectation point prediction, but also give you a predictive density based on the $t$ distribution and an estimate of parameter and residual uncertainty. You can therefore use the [prediction-interval] from such a model for your quantile. (Be careful not to use the confidence-interval - this is a different thing, and the two are very often confused.)

   The advantage is that your quantiles will certainly not be inconsistent, since you will derive the 90% and the 95% quantile prediction from the same predictive density. The disadvantage is that you need full predictive densities. SVM, in particular, does not output these, so even if you found it to perform well in other applications, it may not be the best tool for quantile predictions.

Finally, you may be interested in Gneiting (2011, IJF), who discusses objective functions that are optimized by quantiles (see above), as well as ways of actually evaluating quantile predictions.