GLM: effect of link function on choice of transformation of covariate

Question

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))

, 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?

Answer 1

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}