Dear Drew and Liz,

Thank you both very much for your help in explaining this to me. I was
pleasantly surprised to have found your answers in my mailbox after having
thought that I made a mistake in posting a question right before many go on
vacation! I'm happy that you found this to be an interesting example. You
have indeed helped me better understand the realm of GLMs.

I believe that I have made one error in my description of the problem -
that my concentrations (Y) may not actually have increasing measurement
error with concentration. My confusion comes from the fact that the error
of estimated chlorophyll via remote sensing algorithms is often given as a
percent (e.g. 50%). This relationship is usually derived via a log
concentration plot of in situ versus predicted values, where the error is
pretty much normally distributed when observed on this log-scale (see the
Globcolour validation report Appendices for an example:
http://www.globcolour.info/validation/report/GlobCOLOUR_FVR_v1.1.pdf). So
this is different in that the error of the model's prediction is plotted
versus an independent data set.

So, I believe for the reasons you mention, I will stick with the log-link
GLM, but rather assume a Gaussian distribution for my in situ estimated
(Y). This solves the problem of predicting negative values and actually
gives a better fit over that of the LM.

I unfortunately do not have a large data set for my model (nor an
independent data set for validation), so I am planning to present the
confidence intervals around my predicted values. Interestingly, this is not
out of line with the expectation that predictions of higher concentrations
have larger confidence intervals, and thus, if an independent data set were
available, the predicted vs in situ plot might closely resemble the more
typical remote sensing validation plots.

pred <- predict(GAUSS.LOG.GLM, type="response", se.fit=TRUE)

LIM <- range(c(pred$fit-(pred$se.fit*1.96), pred$fit+(pred$se.fit*1.96),
DF$Y))
plot(pred$fit ~ DF$Y, t="n", xlim=LIM, ylim=LIM)
abline(0,1,lty=1, col=8, lwd=2)
segments(x0=DF$Y, x1=DF$Y, y0=pred$fit-(pred$se.fit*1.96),
y1=pred$fit+(pred$se.fit*1.96), lwd=2)
points(pred$fit ~ DF$Y, pch=1, lwd=2)


Again thanks for your help,
Marc




On Wed, Dec 26, 2012 at 1:55 PM, Liz Pryde <elizabethpryde@gmail.com> wrote:

> Hi Marc,
>
> The main difference between the 2 methods is that with linear models (lm)
> you are trying to shape the data to fit the model's assumptions. In GLMs
> you're trying to fit a model to the data. The family you choose in order to
> estimate the error distribution of your data (account for the mean-variance
> relationship, in your case the increasing variance with mean). The link
> aims to linearise the relationship between your predictor(s) and your
> response variable. You can usually decide what this is by looking at the
> distribution of the data and the way it was collected. What is 'better'
> about the GLM is that it gives reasonable results (ie non-negative where
> negative values are impossible) because it captures the error structure
> better leading to more accurate predictions.
>
> Of course, large amounts of zeroes are hard to model, but there are a
> number of ways around this - using negative binomial distribution,
> zero-inflated models, or positive distributions with a mass at zero (e.g.
> tweedie). With GLMs and their extensions you have more flexibility for such
> options.
>
> Really, the best way to check is to do your model diagnostics. Which
> approach gives less or no pattern to the residuals? That approach should be
> your answer.
>
> A good paper on the topic is Bob O'Hara's: "do not log transform count
> data" in Methods in Ecology & Evolution.
> You can get it here:
> http://www.mendeley.com/profiles/bob-ohara/
>
> Hope that helps,
> Liz
>
> On 26/12/2012, at 4:35 PM, Drew Tyre <atyre2@unl.edu> wrote:
>
> I can see why this would be frustrating - you come to different
> conclusions about X2 depending on which model you fit.
>
> Using a GLM with a log link and Gaussian observation distribution is
> not the same as using lm on log transformed data. In both cases the
> observations are assumed to be normally distributed with constant
> variance; this implies that the observations can come from -infinity
> to +infinity. In the log link glm however, the expected value of y is
> constrained to be larger than zero. That is, even though the
> observations y can be less than zero, the expected value of the y is
> strictly greater than zero because of the log link.
>
> In contrast, log transforming the y maps the domain [0,infinity) into
> (-infinity, infinity), which is exactly what is required for the
> validity of the assumed normal distribution of the errors. Thus the
> distribution of the observations on the original scale is log-normal,
> and all is well.
>
> The test of this idea is to see if a glm with log-link is allowed to
> have values of y <0:
> # Try a parametric bootstrap to see what the assumed models really look
> like
> mu <- fitted(GAUSS.LOG.GLM)
>
> # the dispersion parameter is the residual variance, roughly.
> disp = summary(GAUSS.LOG.GLM)$dispersion
> # and this generates negative numbers
> y.test <- rnorm(length(mu),mu,varfunc(mu)*sqrt(disp))
> # but you can still fit a log-link glm, as long as you provide starting
> values
> glm(y.test~X1+X2,data=DF,start=c(-5,12,-8),family=gaussian(link="log"))
>
> However, even if that explains why you get radically different
> answers, it doesn't really help you decide which one is better! The
> log(Y) lm and the Gamma glm are making the correct assumptions about
> the distribution of Y, that is, that Y > 0, and that var(Y) ~ Y. From
> a pure "prediction" perspective the Gaussian glm is doing better by
> making incorrect assumptions, but that extra performance is focused
> entirely on those few very large values. One way to see this is to
> look at plots of the residuals
>
> par(mfrow=c(2,2))
> plot(LOG.LM)
> plot(GAUSS.LOG.GLM)
>
> In the first plot, you can see that the residuals are pretty well
> distributed, and there is no heteroscedasticity. In the 2nd plot, the
> heteroscedasticity (dependence of var(Y) on Y) is clear in the
> scale-location plot, and you can see the strong effect of the large
> observations in the plot of residuals vs. leverage (points 15 and 25
> have cooks distance > 1) and are the 2 largest observations in the
> dataset. You mentioned that the measurement error increases with
> concentration - that suggests that these large points *should* be
> downweighted with respect to smaller observations - they are less
> accurate. This is happening naturally in models that allow for
> increasing variance with increasing Y (log-normal and gamma)
>
> I would go with the Gamma GLM, I think. The error distribution meets
> the assumptions, the residuals look ok, and the prediction performance
> is better than the LOG.LM.
>
> Thanks for posting an interesting example.
>
> hth
>
> On Dec 22, 2012, at 5:53 AM, Marc Taylor <marchtaylor@gmail.com> wrote:
>
> 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]]
>
>
> _______________________________________________
>
> R-sig-ecology mailing list
>
> R-sig-ecology@r-project.org
>
> https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
>
>
>
>
>
> --
> Drew Tyre
>
> School of Natural Resources
> University of Nebraska-Lincoln
> 416 Hardin Hall, East Campus
> 3310 Holdrege Street
> Lincoln, NE 68583-0974
>
> phone: +1 402 472 4054
> fax: +1 402 472 2946
> email: atyre2@unl.edu
> http://snr.unl.edu/tyre
> http://aminpractice.blogspot.com
> http://www.flickr.com/photos/atiretoo
>
>    [[alternative HTML version deleted]]
>
> _______________________________________________
> R-sig-ecology mailing list
> R-sig-ecology@r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-ecology
>
>

	[[alternative HTML version deleted]]

