# Loss function in for gamma objective function in regression in XGBoost?

Suppose I want to predict $$y$$ from a set of predictors $$x$$. $$y$$ is gamma distributed, so I want to use gamma regression with XGBoost. The help page of XGBoost specifies, for the objective parameter (loss function):

reg:gamma: gamma regression with log-link. Output is a mean of gamma distribution. It might be useful, e.g., for modeling insurance claims severity, or for any outcome that might be gamma-distributed.

What is the explicit formula for this loss function? I haven't been able to find it. Also, Why is the log-link required? As far as I understand, log is usually used as the link function when using Generalized Linear Models, but in gradient boosting it would make any sense...

I'm not actually sure whether "gamma regression" is officially defined (it doesn't appear to have a wikipedia page for example), but if I were to define it (and some googling around suggests I'm not alone here), I would define it as setting up my regression problem so that for a given input vector $$\underline{x}$$, I predict a value y. When I do this, I meant that I think that the true value of the target will be gamma distributed with mean y.

How does this differ from least squares? One setup which leads to OLS is that you say that for given $$\underline{x}$$, the target variable will be normally distributed and your prediction $$y(\underline{x})$$ is the mean of that distribution. Of course a normal distribution is parametrised through its mean and variance, but it turns out that you don't need to know the variance in order to calculate the cost function you need to optimise, and thus this parameter doesn't need to be passed to xgboost.

For the gamma distribution however, this is different. Let's go through the maths. For the gamma distribution parametrised as $$\frac{1}{\Gamma (k)\theta ^{k}}x^{k-1}e^{-\frac{x}{\theta}}$$, the mean is given by $$k\theta$$ and the variance by $$k \theta ^{2}$$

Thus let's reparametrise in terms of $$\mu$$ and $$\theta$$ so that the distribution is given by $$\frac{1}{\Gamma (\frac{\mu}{\theta})\theta ^{\frac{\mu}{\theta}}}x^{\frac{\mu}{\theta}-1}e^{-\frac{x}{\theta}}$$

So for a given dataset, if you predict a bunch of $$\hat{y}_{i}$$ for target values $$y_{i}$$, the likelihood is given by

$$\prod _{i=1}^{N} \frac{1}{\Gamma (\frac{\hat{y}_{i}}{\theta})\theta ^{\frac{\hat{y}_{i}}{\theta}}}y_{i}^{\frac{\hat{y}_{i}}{\theta} -1}e^{-\frac{y_{i}}{\theta}}$$

and thus the negative (xgboost assumes a cost function you're trying to minimise) log-likelihood is

$$\sum _{i=1}^{N} \ln \Gamma (\frac{\hat{y}_{i}}{\theta}) + \frac{\hat{y}_{i}}{\theta}\ln \theta - (\frac{\hat{y}_{i}}{\theta}-1) \ln y_{i} + \frac{y_{i}}{\theta}$$

Compare this to the Gaussian Regression case where the negative log likelihood is given by $$\frac{1}{\sigma ^{2}}\sum _{i=1}^{N} \left(y_{i} - \hat{y}_{i}\right)^{2}$$

The the latter case, the $$\frac{1}{\sigma ^{2}}$$ is a constant term out the front. If you were doing linear regression or even xgboost without regularisation, this would mean that no matter what value you changed $$\sigma$$ to, the linear regressor/xgboost you trained would turn out to be exactly the same, so "Gaussian regression with $$\sigma = 10$$ and Gaussian regression with $$\sigma = 1$$ lead to the same predictions". This is no longer true when you have a regulariser, but you can always suck the value of $$\sigma$$ into the definition of the regulariser to get around this, and this is why the OLS formula never includes a $$\sigma$$ in it.

In the gamma case however, because of the $$\theta$$ factor contained in the $$\Gamma$$ function and the $$\ln \theta$$, you can't just pull the factor of $$\theta$$ outside of the summation.

For xgboost, you now need to pass it the elementwise first and second derivatives of the cost function wrt $$\hat{y}_{i}$$. This is where basic calculus doesn't get you all the way, you'll likely need to look up that the derivative of the logarithm of the gamma function is given by the digramma function $$\psi (z)$$.

The (elementwise) first derivative of the loss will be given by (by the xgboost definition which is $$G_{i}=\frac{\partial L}{\partial \hat{y}_{i}}$$):

$$G_{i} = \frac{1}{\theta}\psi (\frac{\hat{y}_{i}}{\theta}) + \frac{1}{\theta}\ln \theta - \frac{1}{\theta}\ln y_{i}$$

The second derivative will require derivatives of the digamma function, I don't know much about this but some googling tells me you need the trigamma function $$\psi _{1}(z)$$ which is the derivative of the digamma function, thus

$$H_{i}=\frac{1}{\theta ^{2}}\psi_{1}(\frac{\hat{y}_{i}}{\theta})$$

Again, note that you still have to supply $$\theta$$ up front as a hyperparameter, pass this to xgboost and then train a new xgboost model every time you wish to investigate another $$\theta$$

Finally, it's worth noting that I just did this derivation myself today, I haven't lifted anything other than the definition of the gamma distribution from elsewhere, so there could easily be minor algebra error, I'd feel more comfortable if somebody else independently verified my workings.

Edit: Alternately, you could parametrise the other way around:

You could use k as your free parameter and $$\theta = \frac{\mu}{k}$$, thus your gamma distribution is $$\frac{1}{\Gamma(k)(\frac{\mu}{k})^{k}}x^{k-1}e^{-\frac{xk}{\mu}}$$

and thus your negative log-likelihood is given by

$$L = \sum_{i=1}^{N}\left[\ln \Gamma (k) + k\ln \hat{y}_{i} - k\ln k - (k-1)\ln y_{i} + \frac{y_{i}k}{\hat{y}_{i}}\right]$$

this parametrisation is easier to differentiate, you get

$$G_{i}=\frac{k}{\hat{y}_{i}} - \frac{y_{i}k}{\hat{y}_{i}^{2}}$$

and

$$H_{i}=-\frac{k}{\hat{y}_{i}^{2}}+ 2\frac{y_{i}k}{\hat{y}_{i}^{3}}$$

but you still need to pass k as a hyperparameter

• Thanks for your answer. The derivation of the negative likelihood seems right but, in XGBoost you can definitely use objective reg::gamma and no $\theta$ parameter has to be inputed. There is also a gamma-nloglik eval_metric available, but I have not been able to obtain a formula for the first and second derivatives of this loss function.
– D1X
Aug 27 '20 at 9:41
• My c++ isn't good enough to fully understand the source code (github.com/dmlc/xgboost/blob/…) but it looks to me like it's hard-coded to be 1. I would suggest using the learning API and coding this up yourself, that will allow you to pass $\theta$ manually. Last I used it, the syntax (in python) was awkward and the cost function didn't allow for other kwargs so you need to make theta a global variable. Aug 27 '20 at 9:56
• Okay, I may be wrong (please somebody correct me if I am), but this is how I understand it. In regression we are interested in predicting $E(y \ | \ x) = \mu$. For a gamma distributed value, this would be $\mu = k \theta$, for $k, \theta > 0$ shape and scale respectively. Thus, for predicting $\mu$ we can just fix $\theta = 1$ and then adjust just the $k$, as this can adjust for any parameter of $\mu$. This means the distribution of $y \ | \ x$ will probably not be very reliable, but still $\mu$ will be, which suffices for deterministic prediction.
– D1X
Aug 27 '20 at 10:49
• (continuation) Setting $\theta = 1$ also simplifies the log likelihood, which would become $\sum _{i=1}^{N} \ln \Gamma (\hat{y}_{i}) - (\hat{y}_{i}-1) \ln y_{i} + y_{i}$.
– D1X
Aug 27 '20 at 11:02
• I think you could parametrise the other way around, let k be your free parameter, then your cost function is: $\sum_{i=1}^{N}\left[\ln \Gamma(k) + k \ln \hat{y}_{i} -k\ln k - (k-1)\ln y_{i} +\frac{y_{i}k}{\hat{y}_{i}}\right]$, but you still have to arbitrarily choose k, or run for different k values and see which gives you the best training loss. I think this parametrisation looks more like the source code and also doesn't require di- and tr-gamma functions. Will update my original post in due course to include this parametrisation. Aug 27 '20 at 11:56