# Choosing between LM and GLM for a log-transformed response variable

I'm trying to understand the philosophy behind using a Generalized Linear Model (GLM) vs a Linear Model (LM). I've created an example data set below where:

$$\log(y) = x + \varepsilon$$

The example does not have the error $\varepsilon$ as a function of the magnitude of $y$, so I would assume that a linear model of the log-transformed y would be the best. In the example below, this is indeed the case (I think) - since the AIC of the LM on the log-transformed data is lowest. The AIC of the Gamma distribution GLM with a log-link function has a lower sum of squares (SS), but the additional degrees of freedom result in a slightly higher AIC. I was surprised that the Gaussian distribution AIC is so much higher (even though the SS is the lowest of the models).

I am hoping to get some advice on when one should approach GLM models - i.e. is there something I should look for in my LM model fit residuals to tell me that another distribution is more appropriate? Also, how should one proceed in selecting an appropriate distribution family.

[EDIT]: I have now adjusted the summary statistics so that the SS of the log-transformed linear model is comparable to the GLM models with the log-link function. A graph of the statistics is now shown.

### Example

set.seed(1111)
n <- 1000
y <- rnorm(n, mean=0, sd=1)
y <- exp(y)
hist(y, n=20)
hist(log(y), n=20)

x <- log(y) - rnorm(n, mean=0, sd=1)
hist(x, n=20)

df  <- data.frame(y=y, x=x)
df2 <- data.frame(x=seq(from=min(df$x), to=max(df$x),,100))

#models
mod.name <- "LM"
assign(mod.name, lm(y ~ x, df))
summary(get(mod.name))
plot(y ~ x, df)
lines(predict(get(mod.name), newdata=df2) ~ df2$x, col=2) mod.name <- "LOG.LM" assign(mod.name, lm(log(y) ~ x, df)) summary(get(mod.name)) plot(y ~ x, df) lines(exp(predict(get(mod.name), newdata=df2)) ~ df2$x, col=2)

mod.name <- "LOG.GAUSS.GLM"
assign(mod.name, glm(y ~ x, df, family=gaussian(link="log")))
summary(get(mod.name))
plot(y ~ x, df)
lines(predict(get(mod.name), newdata=df2, type="response") ~ df2$x, col=2) mod.name <- "LOG.GAMMA.GLM" assign(mod.name, glm(y ~ x, df, family=Gamma(link="log"))) summary(get(mod.name)) plot(y ~ x, df) lines(predict(get(mod.name), newdata=df2, type="response") ~ df2$x, col=2)

#Results
model.names <- list("LM", "LOG.LM", "LOG.GAUSS.GLM", "LOG.GAMMA.GLM")

plot(y ~ x, df, log="y", pch=".", cex=3, col=8)
lines(predict(LM, newdata=df2) ~ df2$x, col=1, lwd=2) lines(exp(predict(LOG.LM, newdata=df2)) ~ df2$x, col=2, lwd=2)
lines(predict(LOG.GAUSS.GLM, newdata=df2, type="response") ~ df2$x, col=3, lwd=2) lines(predict(LOG.GAMMA.GLM, newdata=df2, type="response") ~ df2$x, col=4, lwd=2)
legend("topleft", legend=model.names, col=1:4, lwd=2, bty="n")

res.AIC <- as.matrix(
data.frame(
LM=AIC(LM),
LOG.LM=AIC(LOG.LM),
LOG.GAUSS.GLM=AIC(LOG.GAUSS.GLM),
LOG.GAMMA.GLM=AIC(LOG.GAMMA.GLM)
)
)

res.SS <- as.matrix(
data.frame(
LM=sum((predict(LM)-y)^2),
LOG.LM=sum((exp(predict(LOG.LM))-y)^2),
LOG.GAUSS.GLM=sum((predict(LOG.GAUSS.GLM, type="response")-y)^2),
LOG.GAMMA.GLM=sum((predict(LOG.GAMMA.GLM, type="response")-y)^2)
)
)

res.RMS <- as.matrix(
data.frame(
LM=sqrt(mean((predict(LM)-y)^2)),
LOG.LM=sqrt(mean((exp(predict(LOG.LM))-y)^2)),
LOG.GAUSS.GLM=sqrt(mean((predict(LOG.GAUSS.GLM, type="response")-y)^2)),
LOG.GAMMA.GLM=sqrt(mean((predict(LOG.GAMMA.GLM, type="response")-y)^2))
)
)

png("stats.png", height=7, width=10, units="in", res=300)
#x11(height=7, width=10)
par(mar=c(10,5,2,1), mfcol=c(1,3), cex=1, ps=12)
barplot(res.AIC, main="AIC", las=2)
barplot(res.SS, main="SS", las=2)
barplot(res.RMS, main="RMS", las=2)
dev.off()


-

Good effort for thinking through this issue. Here's an incomplete answer, but some starters for the next steps.

First, the AIC scores - based on likelihoods - are on different scales because of the different distributions and link functions, so aren't comparable. Your sum of squares and mean sum of squares have been calculated on the original scale and hence are on the same scale, so can be compared, although whether this is a good criterion for model selection is another question (it might be, or might not - search the cross validated archives on model selection for some good discussion of this).

For your more general question, a good way of focusing on the problem is to consider the difference between LOG.LM (your linear model with the response as log(y)); and LOG.GAUSS.GLM, the glm with the response as y and a log link function. In the first case the model you are fitting is:

$\log(y)=X\beta+\epsilon$;

and in the glm() case it is:

$\log(y+\epsilon)=X\beta$

and in both cases $\epsilon$ is distributed $\mathcal{N}(0,\sigma^2)$.

-
The characterization of the glm does not look correct: on the left hand side is a random variable $\epsilon$ while the right hand side contains only data and parameters but no random variables. –  whuber Oct 30 at 22:08
It's an odd way of putting it I know @whuber but comes from $E(Y) = g^{-1}(X\beta)$ becomes $g(E(Y)) = X\beta$. The point being that the link function goes around the $E(Y)$ –  Peter Ellis Nov 6 at 22:59

In a more general way, E[ln(Y|x)] and ln([E(Y|X]) are not the same. Also the variance assumptions made by GLM are more flexible than in OLS, and for certain modeling situation as counts variance can be different taking distinct distribution families.

About the distribution family in my opinion is a question about the variance and its relation with the mean. For example in a gaussian family we have constant variance. In a gamma family we have variance as a quadratic function of the mean. Plot your standarized residuals vs the fitted values and see how they are.

-
+1 for actually relating to the question of how to choose the right family (and I'd say there's room for some more elaboration here) –  etov Jun 18 at 10:44

Unfortunately, your R code does not lead to an example where $\log(y) = x + \varepsilon$. Instead, your example is $x = \log(y) + \varepsilon$. The errors here are horizontal, not vertical; they are errors in $x$, not errors in $y$. Intuitively, it seems like this shouldn't make a difference, but it does. You may want to read my answer here: What is the difference between linear regression on y with x and x with y? Your setup complicates the issue of what the "right" model is. Strictly, the right model is reverse regression:

ly = log(y)
REVERSE.REGRESSION = lm(x~ly)
summary(REVERSE.REGRESSION)
# Call:
# lm(formula = x ~ ly)
#
# Residuals:
#      Min       1Q   Median       3Q      Max
# -2.93996 -0.64547 -0.01351  0.63133  2.92991
#
# Coefficients:
#             Estimate Std. Error t value Pr(>|t|)
# (Intercept)  0.01563    0.03113   0.502    0.616
# ly           1.01519    0.03138  32.350   <2e-16 ***
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
# Residual standard error: 0.984 on 998 degrees of freedom
# Multiple R-squared:  0.5119,    Adjusted R-squared:  0.5114
# F-statistic:  1047 on 1 and 998 DF,  p-value: < 2.2e-16


Metrics for this model (like the AIC) won't be comparable to your models. However, we know that this is the right model based on the data generating process, and notice that the estimated coefficients are right on target.

-
Thanks for your comment. I admit, the example data could have been better, but I believe it to be correct in how it generated errors. In the example, there is no intercept and the slope is 1. If you turn around the line x = log(y) - rnorm(n, mean=0, sd=1), you get log(y) = x + rnorm(n, mean=0, sd=1). If @whuber's comment spawned your answer (I believe it did), then I believe he is not referring to the data generation, but rather the GLM model formulation by @peterellis . –  Marc in the box Oct 31 at 7:51