How to model bias in dependent variable?

Question

I am trying to model the odds of soccermatches in play, based on the odds at start of the match and possesion during the game. My dataset contains:

Start_odd (x1)  Possesion (x2) Market_odd_observed (y)
0.67            80             0.90
0.45            75             0.63 etc

Start_odd is on a scale of 0-1. Possesion is on a scale of 0-100. Market_odd is on a scale of 0-1.

The GAM-model is fitted using mgcv:

Family: gaussian 
    Link function: identity 

    Formula:
    Market_odd_observed ~ s(Start_odd , k = 20) + s(Possesion , k = 20) + ti(Start_odd , 
        Possesion , k = c(10, 10))

    Parametric coefficients:
                 Estimate Std. Error t value Pr(>|t|)    
    (Intercept) 7.394e-01  4.609e-05   16043   <2e-16 ***
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

    Approximate significance of smooth terms:
                              edf Ref.df      F p-value    
    s(Start_odd )             18.87  19.00 288685  <2e-16 ***
    s(Possesion )             18.95  19.00 190429  <2e-16 ***
    ti(Start_odd ,Possesion ) 69.69  75.33  12433  <2e-16 ***
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

    R-sq.(adj) =  0.947   Deviance explained = 94.7%
    -REML = -1.134e+06  Scale est. = 0.0012332  n = 587663

If I plot the residuals by the dependent variable I still see a pattern in the data with a upward slope:

I conclude that there is some bias in the model. The problem is that I can not include the dependent variable as an interaction term since this is the outcome I try to predict. Is it unusual to look at the residuals grouped by the dependent variable?

I have tried to fit a second gam-model with the predictions from the model above as the input. Unfortunately the RMSE is exacly the same and the pattern is still there.

I have also plotted the residuals by the predictions. In that case the bias is not there as can be seen in this plot:

Is there an alternative method to improve the model?

Next I have fitted a catagorical GAM on the winflag of the match (0 or 1). The results are the same as above.

Next I have plotted the observed odds vs predictions:

And transparant:

Next I have grouped_by the errors by observed market odd:

And grouped by prediction:

I expect it is not related to the use of the GAM since there is similar pattern using a neural network. What could be the explanation that the models do not fit this pattern?

Thanks a lot!

I have added an example to illustrate the answer from Aksakal:

library(tidyverse)
library(ggplot2)
library(mgcv)
library(mlbench)

data("BostonHousing")


gam_y <-
  gam(
    medv ~ s(nox) + s(rm) + s(dis) ++s(tax) + s(ptratio) + s(lstat) ,
    method = "REML",
    data = BostonHousing
  )

y_pred <- predict(gam_y)
predictions <-
  cbind(BostonHousing$medv, y_pred, resi = BostonHousing$medv - y_pred)
predictions <- as.data.frame(predictions)
colnames(predictions)[1] <- "medv"

ggplot(predictions, mapping = aes(x = medv, y = resi)) +
  geom_point(alpha = 100 / 100) +
  geom_smooth(method = lm) +
  labs(y = "residual", x = "house price observed (y)") +
  ggtitle("residuals by y")

ggplot(predictions, mapping = aes(x = y_pred, y = resi)) +
  geom_point(alpha = 100 / 100) +
  geom_smooth(method = lm) +
  labs(y = "residual", x = "house price predicted (y)") +
  ggtitle("residuals by y")

Answer 1

when you work with probabilities consider cross-entropy loss instead of fitting the least squares of deviations (residuals). the most straightforward application of this is logit regression.

For instance, consider a logit link function in your code to combine GLM with GAM.

On residuals $e=y-\hat y$ (observed minus prdicted not the other way!), when you plot them vs the observed $y$, they'll exhibit negative bias. The reason is that when $y\to 1$, then predictions will tend to be below observed. In an extreme case where observed is 1 unless you have an absolutely perfect predictive power you must have $\hat y<y=1$, therefore on the right end you must have positive bias, then with similar consideration for $y\to 0$ on the left end you must have negative bias in the plot.

Generally, in any model $y=f(x)+\varepsilon$ you have the errors baked in the dependent variables, hence the plots $\varepsilon\sim y=f(x)+\varepsilon$ have the erros both in x- and y-axes, thus making the plots correlated (sloped). It is, therefore, preferable to plot $\varepsilon\sim f(x)$, i.e. error vs predicted not observed.