1
$\begingroup$

It struck me that if I have data of the form below,

library(data.table)
dt = data.table(xxx = rep(1:50, 10))
dt[, y := rgamma(500, 1.1 * xxx + runif(500)/2,)]
with(dt, plot(xxx, y))

data

, that before I would have though "ah this is a linear effect, so I should include a linear effect in my Gamma-regression", but upon realizing that in Gamma-regression one uses the log-link, and that any prediction from my model is obviously transformed via the exponential function $\hat y_i= \exp(\beta_0 + \beta_1xxx_i)$,

such that the best model would have a logarithmic transformation of the covariate $xxx$. When I check the Gamma deviance, this seems to be correct:

library(mgcv)
model = mgcv::gam(y ~ xxx, family = Gamma(link = "log"), data = dt)
model2 = mgcv::gam(y ~ log(xxx), family = Gamma(link = "log"), data = dt)

dt[, pred := predict(model, type ="response", newdata=dt)]
dt[, pred2 := predict(model2, type ="response", newdata=dt)]

dt[, dev := (y - pred)/pred - log(y/pred)]
dt[, dev2 := (y - pred2)/pred2 - log(y/pred2)]

dt[, list(sum(dev), sum(dev2))]

Deviance for model1: 43.9819, and for model2: 19.53396.

Am I loosing my mind after a long summer break, or am I thinking correctly?

$\endgroup$
3
$\begingroup$

If you believe the conditional expectation of $Y$ given $xxx$ is truly linear, your approach is completely reasonable. Your approach is equivalent to using the identity link function. Let $g()$ be an inverse link function:

\begin{align} E[Y_i|X_i] = g(X_i\beta) \end{align} Then if $g_1(X_i\beta) = exp(X_i\beta)$, i.e under the log-link and $g_2(X_i\beta) = X_i\beta$ under the indentity. Then

\begin{align} g_1(log(X_i)\beta) &= exp(log(X_i)\beta)\\ &= X_i\beta\\ &= g_2(X_i\beta) \end{align}

Holding as long as $log(X_i)$ is defined, i.e has positive support.

In general, the typical objection to such a link function is that you cannot guarantee that the predicted outcome lies in the support implied by a gamma random variable, i.e $\hat{Y} \in [0, \infty)$. Since the random variable $xxx$ is a positive random variable, this is a very reasonable approach, because there is no need to extrapolate to cases which might force your predicted values outside of the target support. In practice we can run into real data sets where our sample data has these properties but the underlying random variables do not in which case finding ways to reasonably extrapolate from our model might prove difficult or problematic.

In your case it is justified and $E[Y|X]$ can be shown to be in fact linear in $X$. $Y \sim gamma(\alpha(X,Z),1)$, where $X$ is $xxx$ and $Z$ is a uniform random variable $Z \sim U(0, 0.5)$.

$\alpha(X,Z) = 1.1*X + Z$ Using standard result that for the expectation gamma random variables $E[Y|X,Z] = \frac{\alpha(X,Z)}{\beta}$ and the fact that $X$ and $Z$ are independent \begin{align} E[Y|X] &= \int_0^{0.5} E[Y|X,Z]f(z)dz \text{ [independence of } X \text{ and } Z]\\ &= \int_0^{0.5}\frac{\alpha(X,z)}{1}f(z)dz\\ &= 1.1X + \int_0^{0.5}z*2dz\\ &= 1.1X + 0.25 \end{align} Which is linear in $X$.

In practice, in the case where we did not simulate the data, but rather just saw this relationship from the data, it is likely that we would just use linear regression since it implies linear conditional expectations and has nice robustness properties with respect to error distributions, but gamma with the identity link will give us a near identical answer in most cases.

Proof of first line: \begin{align} \int E[Y|X,Z] f(z)dz &= \int\int y f(y|x,z)f(z)dydz\\ &= \int\int y \frac{f(y,x,z)}{f(x,z)}f(z)dydz\\ &= \int\int y \frac{f(y,x,z)}{f(x)f(z)}f(z)dydz \text{ [Independence] }\\ &= \int\int y \frac{f(y,x,z)}{f(x)}dydz\\ &= \int y \frac{f(y,x)}{f(x)} dy\\ &= \int yf(y|x)dy\\ &= E[Y|X] \end{align}

| cite | improve this answer | |
$\endgroup$
  • $\begingroup$ Thanks for answering. So, I know that my data are gamma-distributed, and I would never use a linear regression. I could have stated this explicitly. You could imagine I had more covariates, factors and continuous variables, such that the case does not simplify down to using the identity link, and my question would then be: if one knows that the link function is the logarithm, and you observe some linear relationship for one of your covariates, then should you now always model this as y ~ log(covariate with linear effect) + covariate_2 + covariate_3 .... $\endgroup$ – Erosennin Aug 6 at 18:56
  • $\begingroup$ As a statistician, I'm terrified to say the word always, but the way I understand your situation I would say the following. If you are thinking about your model in a regression sense (not a generative probabilistic model for example) then implicitly we are usually decomposing the outcome into a conditional expectation and then the remaining error (you are confident is gamma, in general we may turn to more semi-parametric methods which relax such assumptions). Then yeah if you know the conditional expectation function is partially linear and partially multiplicative, then model it that way) $\endgroup$ – Tyrel Stokes Aug 6 at 19:26
  • $\begingroup$ Typically we care more about the conditional expectation part and it's the part we want to get most right. Notice that you could use the identity link and then model the non-linear components inside an exponential function equivalently. The point then isn't the link but to get the conditional expectation correct, whatever that means. You may find it useful to look at partially linear regression models as they are a more general case of your problem. We assume one part of the conditional expectation is linear, another part is not, and the error is typically semi-parametric. $\endgroup$ – Tyrel Stokes Aug 6 at 19:28

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.