63
$\begingroup$

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.

Many thanks in advance for your help.

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

enter image description here

enter image description here

Improve this question
$\endgroup$
2
  • $\begingroup$ The formula for predicted value in log.lm is incorrect. $exp(X \hat{beta}) $ gives the median of $y$. To get expected value, add $1/2 \times sigma^2$ in exponent $\endgroup$
    – pauljohn32
    Commented Dec 27, 2018 at 4:35
  • 1
    $\begingroup$ Another model, for which R does not offer a family, is a lognormal distribution. SAS will fit that, I don't know why R glm doesn't. Some suggest R package gamlss for tgat, but it never works understandably for me. Maybe you will have better luck. $\endgroup$
    – pauljohn32
    Commented Dec 27, 2018 at 4:43

4 Answers 4

Reset to default
28
$\begingroup$

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

Improve this answer
$\endgroup$
3
  • 6
    $\begingroup$ 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. $\endgroup$
    – whuber ♦
    Commented Oct 30, 2014 at 22:08
  • 4
    $\begingroup$ 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)$ $\endgroup$
    – Peter Ellis
    Commented Nov 6, 2014 at 22:59
  • $\begingroup$ I found this very helpful: christoph-scherber.de/content/PDF%20Files/… $\endgroup$
    – Aditya
    Commented Jul 9, 2017 at 7:10
19
$\begingroup$

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.

Improve this answer
$\endgroup$
2
  • 1
    $\begingroup$ +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) $\endgroup$
    – etov
    Commented Jun 18, 2014 at 10:44
  • $\begingroup$ In LM we have a g(Y) = XB. In GLM we have a g(E[Y]) = XB It means that in LM we model a function of Y, and in GLM we model a function of the mean of Y. $\endgroup$
    – igorkf
    Commented Apr 12, 2020 at 19:12
8
$\begingroup$

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.

Improve this answer
$\endgroup$
1
  • $\begingroup$ 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 . $\endgroup$
    – Marc in the box
    Commented Oct 31, 2014 at 7:51
-1
$\begingroup$

The choice is based on your hypothesis on your variable.

the log transformation is based on $$\frac{\sqrt{\mathrm{Var}(X_t} }{\mathrm{E}(X_t)} = \mathrm{constant}$$

the gamma distribution is based on

$$\frac{\mathrm{Var}(X_t) }{\mathrm{E}(X_t)} = \mathrm{constant}$$

The log transformation rests on the hypothesis that,

$$\sqrt{\mathrm{Var}(X_t} = \mathrm{E}(X_t) \sigma$$

In this way,

$$\begin{alignat}{2} X_t &= X_t\\ & = \mathrm{E}(X_t) \cdot \frac{X_t}{\mathrm{E}(X_t)} \\ & = \mathrm{E}(X_t) \cdot \frac{X_t - \mathrm{E}(X_t) + \mathrm{E}(X_t)}{\mathrm{E}(X_t)} \\ & = \mathrm{E}(X_t) \cdot (1 + \frac{X_t - \mathrm{E}(X_t)}{\mathrm{E}(X_t)}) \\ \end{alignat}$$

Based on Taylor rule,

$$\log(1+x) \approx x$$

We get

$$\log(1 + \frac{X_t - \mathrm{E}(X_t)}{\mathrm{E}(X_t)}) = \frac{X_t - \mathrm{E}(X_t)}{\mathrm{E}(X_t)}$$

Thus,

$$\begin{alignat}{2} X_t &= \mathrm{E}(X_t) \cdot (1 + \frac{X_t - \mathrm{E}(X_t)}{\mathrm{E}(X_t)}) \\ \log X_t &= \log \mathrm{E}(X_t) + \log (1 + \frac{X_t - \mathrm{E}(X_t)}{\mathrm{E}(X_t)}) \\ &= \log \mathrm{E}(X_t) + \frac{X_t - \mathrm{E}(X_t)}{\mathrm{E}(X_t)} \\ \mathrm{E}(\log X_t) & \approx \log \mathrm{E}(X_t) \end{alignat}$$

However, the gamma distribution rests on the hypothesis that,

$$Y \sim \Gamma(\alpha, \beta)$$

$$\begin{cases} E(y_i) = \alpha_i \cdot \beta_i \\ Var(y_i) = \alpha_i \cdot \beta_i^2 \\ \end{cases} \to \frac{Var(y_i)}{E(y_i)} = \beta_i$$

Improve this answer
$\endgroup$
1
  • $\begingroup$ The gamma distribution used in GLMs fixes the shape parameter (it's the same for all observations, though unspecified); the model is parameterized in proportion to the scale parameter. In that case it has variance proportional to mean squared (s.d proportional to mean, constant coefficient of variation). $\endgroup$
    – Glen_b
    Commented Jan 23, 2020 at 23:59

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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