Which gamma regression model to use for extrapolation? I'm looking for a regression model which would satify these requirements:


*

*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.

*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.

*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.
 A: 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()))

A: 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)
model.compile(optimizer='adam', loss=loss, metrics=['mse'])

model = Pipeline([
        ('scale', StandardScaler()),
        ('keras', model),
])

With a model similar to this (I added some customization to my specific problem) I was able to get a test loss of 0.653. This model also seems to be best so far for interpolation within the training set.
