# Which gamma regression model to use for extrapolation?

I'm looking for a regression model which would satify these requirements:

1. My target variable follows the exponential distribution, so to my understanding I should use gamma loss function. I have tried this feature from XGBoost, LightGBM, pyGAM, and StatsModels GLM and they all seem to be working well.
2. The model should be able to at least some extrapolation in respect of features $$x_e$$. In my training data one feature has only values $$x_e$$ = 1, 2, 3..., 10, and I'm really interested about the prediction at $$x_e$$=0. Linear models seem to perform well, but tree based models (XGBoost, LightGBM) are very bad at this.
3. Some non-linearity or other types of flexibility to fit complex relationships between the input variables. For example with BayesianRidge I get almost the same performance in training and test sets, which tells me that I should be able to build a better model by adding complexity.

My data set contains around one million samples and I'm using 10 features (both discrete and continuous). Potentially I might get 10-20 new features in the future.

Here are the results I have so far. Since I want the model to perform well on the extrapolation in respect of input x, my test set consists of all the samples with $$x_e$$=1.

• Linear regression softplus link Training loss: 0.949 Test loss: 0.673
• Bayesin Ridge regression Training loss: 0.931 Test loss: 0.653
• Bayesin Ridge with softplus link Training loss: 0.949 Test loss: 0.673
• Elastic Net Training loss: 0.930 Test loss: 0.660
• Elastic Net with softplus link Training loss: 0.946 Test loss: 0.655
• Elastic Net CV Training loss: 0.931 Test loss: 0.654
• Elastic Net CV with log link Training loss: 0.984 Test loss: 0.688
• Elastic Net CV with softplus link Training loss: 0.947 Test loss: 0.664
• GLM with log link and gamma loss Training loss: 0.927 Test loss: 0.657
• GAM with log link gamma loss Training loss: 0.927 Test loss: 0.657
• LGBMRegressor (max_depth 3) Training loss: 0.899 Test loss: 0.678

All models seem to have much better loss at the test set than the training set. Normally this would indicate over fitting, but I believe that the reason is that the values at $$x_e$$=1 are lower than those at other values of $$x_e$$.

Initially I had a completely randomized train-test split. Based on those results LGBM was model for interpolation, but now it is clear that the other models outperform it in extrapolation. Based on those results, none of the current models were overfit. I think that the model should have only linear or other very simple dependency to $$x_e$$, but there is probably potential to increase the model complexity in respect of the other features.

• Just to state the obvious: Why not a GAM with Gamma response? (In R what would be mgcv::gam with family=Gamma(link=log)) It will naturally work with a Gamma distributed response function as well as account for non-linear relations between the response and the explanatory variables. – usεr11852 Mar 7 '19 at 23:52
• Thanks for the advice! I tried this Python library pygam.readthedocs.io/en/latest/api/gammagam.html which seems to have similar features as mgcv::gam. When I fit the model with 20 samples I get pretty nice results. However, when I fit it with all 900000 samples in my training set I get error "PIRLS optimization has diverged". Is it even possible to fit a GAM to a large data set like this? – Heikki Pulkkinen Mar 8 '19 at 8:50
• Yes, it is possible. I commonly use dataset having 1M+ instances (in R). The obvious thing to suggest in your case is to limit the number of basis functions used; this should make the optimisation task solved by PIRLS easier (as it will have smaller dimensionality). This can be done by setting s(..., n_splines=X) where X is something smaller than the default $20$ in pyGAM. – usεr11852 Mar 8 '19 at 9:24
• Sorry, is there a typo at the "LGBMRegressor with gamma objective and max_depth 3"? cause it seems to clearly provide the lowest test loss (0.928) but the text mentions that BRR as the best model where the test loss is higher (0.960). To that extent, GAM with gamma distribution has a lower test loss (0.955) than all other regression routines... (Also, did the n_splines worked or you used something else?) – usεr11852 Mar 11 '19 at 14:15
• Yes, with your help on n_splines I got the GAM work properly. I used n_splines=10. And yes, while both LGBM and GAM give better loss than BRR, the extrapolation does not look correct (see requirement 2). – Heikki Pulkkinen Mar 11 '19 at 14:31

The OP has done a great job exploring a variety of different techniques. As commented, given that the response variable is Gamma-distributed it makes sense to consider a GLM and/or a GAM for Gamma distributed variables. Particularly for the use of GAMs, if the computation burden appears too much we might want to consider restricted the basis functions used by the GAM (in the case of pyGAM used here that being achieved by setting s(..., n_splines=X) where X is something smaller than the default $$20$$.

The main point to rectify is the use of evaluating the error of each method. Simple random resampling by cross-validation is providing us an indication on "interpolation" rather than "extrapolation" performance. Here, given $$x_e = \{1, \dots, 10\}$$, we focus on predicting $$x_e =0$$; therefore it is more reasonable to use instances where $$x_e = 1$$ in our validation set and instances where $$x_e =\{2, \dots, 10\}$$ in our training set. Note that in-sample errors are rather misleading for a extrapolation task; there is no "overfitting" perse because the validation and training set do not refer to the same sample/population. On that matter, the fact we get simple models (Elastic Net regression and Bayesian Ridge regression) as our top-performing routines is not too surprising. When extrapolating most bets are off (e.g. see the CV thread: What is wrong with extrapolation?) and commonly simple methods outperform complex ones (e.g. see the CV thread Best method for short time-series).

As a final note, it is always prudent to get estimates of the variability of our performance metric. If possible we should set aside a number of observations, fit our candidate models to the remaining data, and evaluate the models in the data we set aside. This should be repeated multiple times. In effect, what is described is nested cross-validation for model selection; only particular will be that the for each loop the hold-out set is such that $$x_e = 1$$. (Once again) CV has a great thread on the matter: Nested cross validation for model selection and Model selection and cross-validation: The right way. In short, the outer loop will be used to assess the performance of the particular model (e.g. Ridge regression), and the inner loop will be used to select the best model (the regularisation parameter $$\lambda$$ for the case of Ridge). A simplified and very succinct Python example of nested CV would be as follows:

myS = cross_val_score(GridSearchCV(linear_model.ElasticNet(), param_grid, cv=5),
myX, myY, cv=5)
print("CV scores: ", myS)
print("Mean CV scores & Std. Dev.: {:.3f} {:.3f}".format(myS.mean(), myS.std()))


Linear models are great for extrapolation, since they don't make too complex assumptions which might not generalize outside of the training set. Unfortunately that also means that they are not able to use the full potential of other features which are not changed when extrapolating.

Since extrapolation is done only in respect of input $$x_e$$ we can split all inputs $$X$$ into vector $$x_e$$ and matrix $$X_{other}$$.

One approach could be to add more features based on the values of $$X_{other}$$. These could be interactions between these features, squared values, logarithms etc. This could work well especially if you have a good idea on how each feature should contribute to the output.

One could also build a custom model, which fits only a linear relationship in respect of $$x_e$$ but uses more complex relations for $$X_{other}$$. Here is an example model built with keras which fits a neural network with 4 hidden layers for $$X_{other}$$, then concatenates the last hidden layer with $$x_e$$ and uses one final layer with softplus activation (to ensure that the predictions are positive). With keras

from keras.layers import Input, Dense, Lambda, concatenate
from keras.models import Model
import keras.backend as K
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import StandardScaler

def loss(y_true, y_pred):
return -K.log(1/y_pred) + y_true/y_pred

all_inputs = Input(shape=(n_inputs,))

# Only the first feature will be used for extrapolation
x_e = Lambda(lambda x: x[:,:1], output_shape=(1,))(all_inputs)

x_other = Lambda(lambda x: x[:,1:], output_shape=(n_inputs-1,))(all_inputs)
hidden_1 = Dense(20, activation='relu')(x_other)
hidden_2 = Dense(20, activation='relu')(hidden_1)
hidden_3 = Dense(10, activation='relu')(hidden_2)
hidden_4 = Dense(20, activation='relu')(hidden_3)

merged = concatenate([hidden_3, x_e])
preds = Dense(1,activation='softplus')(merged)

model = Model(inputs=all_inputs, outputs=preds)