Dear r-sig-ecology group,

I am having some difficulties understanding the difference between results
obtained for a multiple regression linear model fit either by a General
Linear Model (LM) or a Generalized Linear Model (GLM).

I have attached an example data set in which I am trying to fit Y ~ X1 +
X2. Y is a concentration and the distribution of observed values is highly
skewed towards zero (with no negative numbers). I would expect that the
error in measuring Y increases with it's value. An analog would be in the
estimation of Chla concentrations via remote sensing - e.g. the typical
error of a SeaWifs estimate is usually given as a percent of the
concentration so that large values have wider confidence intervals
surrounding them.

I originally started out fitting a multiple LM to the data without
log-transforming Y. This is probably incorrect given that Y is so skewed
toward zero, and that model can predict negative Y values (which is not
possible). Now, I am hoping to fit a model to the log-transformed Y, but am
unsure if I should remain with a simple LM or look into using a GLM. I have
yet to find a satisfactory answer of when GLM is more appropriate.

My initial explorations seem to indicate that there may be advantages to
using a GLM - specifically, one that assumes a Gaussian distribution of the
expectation of Y and with a log-link function.

In the example below, the Sum of squared differences (SS) as calculated
between log-transformed observed vs predicted Y is lowest for the LM. Based
on the original Y scale, SS is much lower for the GLM
(family=gaussian(link=log)) as it does a better job at fitting the high
concentrations.

So, my question is two-fold:
1) What to look for when deciding to choose a log-link GLM over a LM with a
log-transformed response variable?
2) How to best assess whether a GLM fit is better than LM? Is the SS of the
original obs. vs. pred. a valid criteria?


I would really appreciate any advise on these topics as I have been
frustrated by these decision for some time now.

Cheers




######EXAMPLE##########
DF <- structure(list(Y = c(0.158507483, 0.008510161, 0.002684648,
0.009587276,
0.001803681, 0.010173461, 0.002273384, 0.00345826, 0.051424454,
0.029937484, 0.194813452, 0.042138033, 0.022944148, 0.729585218,
0.887009621, 0.008899131, 0.001588576, 0.0216036, 0.001409499,
0.161051383, 0.026504919, 0.001495132, 0.059066545, 0.008317594,
0.490868633, 0.057027831), X1 = c(0.0974369543591941, -0.11971810600977,
-0.168908964300336, -0.0011723143713434, -0.200018273737778,
0.0536459384966756, -0.188248143615029, -0.154736748196712,
0.0529959236206016,
-0.152396350558232, 0.103766445240172, -0.0693365907826557,
-0.114615555500542,
0.488829422819801, 0.561719898192691, -0.0469180067616361,
0.0631502939411764,
-0.135689617930714, 0.0343957489602316, -0.0749974069726867,
-0.107592097416425, 0.067741017650224, -0.167713403634508,
0.275062271178857,
0.276065626134302, -0.0926000525628916), X2 = c(-0.19192408577628,
0.116576354094024, 0.208731289320505, -0.138772290234524,
0.364065047213473,
-0.1574052089755, 0.285540178523006, 0.29343767019163, -0.203222931158516,
0.0835579872715545, -0.157325117354138, -0.0242157560597033,
-0.175123479037643, 0.174087353210292, 0.246559485637939,
-0.43074835446357,
-0.0181308378901971, 0.0525230701557242, -0.121813588478372,
-0.0549274842561502, -0.115591654073407, -0.0190993986035446,
0.124566313208749, 0.138414677580375, -0.0981459346380045,
-0.319191657096572
)), .Names = c("Y", "X1", "X2"), class = "data.frame", row.names = c(NA,
-26L))

hist(DF$Y, n=10)
hist(log(DF$Y), n=10)

fmla1 <- formula(Y ~ X1 + X2)
fmla2 <- formula(log(Y) ~ X1 + X2)

############
###MODELS###
############

#Linear Models
LM <- lm(fmla1, DF)
summary(LM)

LOG.LM <- lm(fmla2, DF)
summary(LOG.LM)

#Generalized Linear Models
GAUSS.LOG.GLM <- glm(fmla1, DF, family=gaussian(link="log"))
summary(GAUSS.LOG.GLM)

GAMMA.LOG.GLM <- glm(fmla1, DF, family=Gamma(link="log"))
summary(GAMMA.LOG.GLM)


#LOG - predicted vs observed
YLIM <- range(c(exp(predict(LOG.LM)), exp(predict(GAUSS.LOG.GLM)),
exp(predict(GAMMA.LOG.GLM))))
plot(1,1, t="n", xlim=YLIM, ylim=YLIM, log="xy", xlab="Observed Y",
ylab="Predicted Y")
abline(0,1, col=8, lwd=2)
points(exp(predict(LOG.LM)) ~ DF$Y, col=1, pch=1)
points(exp(predict(GAUSS.LOG.GLM)) ~ DF$Y, col=2, pch=2)
points(exp(predict(GAMMA.LOG.GLM)) ~ DF$Y, col=3, pch=3)
points(predict(LM) ~ DF$Y, col=4, pch=4)
legend("topleft", legend=c("LOG.LM", "GAUSS.LOG.GLM", "GAMMA.LOG.GLM",
"LM"), col=1:4, pch=1:4)
text(YLIM[1], YLIM[1], labels="1:1", pos=4, col=8)

#LOG sum of squared differences
LOG.SS <- data.frame(
    LOG.LM = sum((predict(LOG.LM) - log(DF$Y))^2),
    GAUSS.LOG.GLM = sum((predict(GAUSS.LOG.GLM) - log(DF$Y))^2),
    GAMMA.LOG.GLM = sum((predict(GAMMA.LOG.GLM) - log(DF$Y))^2)
)
LOG.SS


#ORIG - predicted vs observed
YLIM <- range(c(exp(predict(LOG.LM)), exp(predict(GAUSS.LOG.GLM)),
exp(predict(GAMMA.LOG.GLM)), predict(LM)))
plot(1,1, t="n", xlim=YLIM, ylim=YLIM, log="", xlab="Observed Y",
ylab="Predicted Y")
abline(0,1, col=8, lwd=2)
text(YLIM[1], YLIM[1], labels="1:1", pos=3, col=8)
points(exp(predict(LOG.LM)) ~ DF$Y, col=1, pch=1)
points(exp(predict(GAUSS.LOG.GLM)) ~ DF$Y, col=2, pch=2)
points(exp(predict(GAMMA.LOG.GLM)) ~ DF$Y, col=3, pch=3)
points(predict(LM) ~ DF$Y, col=4, pch=4)
legend("topleft", legend=c("LOG.LM", "GAUSS.LOG.GLM", "GAMMA.LOG.GLM",
"LM"), col=1:4, pch=1:4)

ORIG.SS <- data.frame(
    LM = sum((predict(LM) - DF$Y)^2),
    LOG.LM = sum((exp(predict(LOG.LM)) - DF$Y)^2),
    GAUSS.LOG.GLM = sum((predict(GAUSS.LOG.GLM, type="response") - DF$Y)^2),
    GAMMA.LOG.GLM = sum((predict(GAMMA.LOG.GLM, type="response") - DF$Y)^2)
)
ORIG.SS

	[[alternative HTML version deleted]]

