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



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