18
$\begingroup$

Some of you might have read this nice paper:

O’Hara RB, Kotze DJ (2010) Do not log-transform count data. Methods in Ecology and Evolution 1:118–122. klick.

In my field of research (ecotoxicology) we're dealing with poorly replicated experiments and GLMs are not widely used. So I did a similar simulation as O’Hara & Kotze (2010), but mimicked ecotoxicological data.

Power simulations:

I simulated data from a factorial design with one control group ($\mu_c$) and 5 treatment groups ($\mu_{1-5}$). Abundance in treatment 1 was identical to the control ($\mu_1 = \mu_c$), abundances in treatments 2-5 was half of the abundance in the control ($\mu_{2-5} = 0.5 \mu_c$). For the simulations I varied the sample size (3,6,9,12) and the abundance in the control group (2, 4 ,8, ... , 1024). Abundances were drawn from a negative binomial distributions with fixed dispersion parameter ($\theta = 3.91$). 100 datasets were generated and analysed using a negative binomial GLM and a gaussian GLM + log-transformed data.

The results are as expected: the GLM has greater power, especially when not many animals were sampled. enter image description here Code is here.

Type I Error:

Next I looked at type one error. Simulations were done as above, however all groups had same abundance ($\mu_c = \mu_{1-5}$).

However, the results are not as expected: enter image description here Negative binomial GLM showed a greater Type-I error compared to LM + transformation. As expected the difference vanished with increasing sample size. Code is here.

Question:

Why is there an increased Type-I Error compared to lm+transformation?

If we have poor data (small sample size, low abundance (many zeros)), should we then use lm+transformation? Small sample sizes (2-4 per treatment) are typical for such experiments and cannot be increased easily.

Although, the neg. bin. GLM can be justified as being appropriate for this data, lm + transformation may prevent us from type 1 errors.

Improve this question
$\endgroup$
6
  • 1
    $\begingroup$ Not an answer to your main question, but something for readers to note: unless you make the actual type I error equivalent for the two procedures, comparing power doesn't make sense; I can always make the power higher for the lower one (in this case the take logs and fit normal one) by lifting its type I error. On the other hand, if you specify the particular situation (sample size, abundance), you can get the type I error rate (e.g. by simulation), and so work out what nominal rate to test at to achieve the desired type I error rate, so their power becomes comparable. $\endgroup$
    – Glen_b
    Commented Sep 9, 2014 at 21:16
  • $\begingroup$ Are the y-axis values in your plots averaged across the 100 data sets? $\endgroup$
    – shadowtalker
    Commented Sep 10, 2014 at 2:17
  • $\begingroup$ I should clarify my comment: in the case where statistics are inherently discrete you don't have perfect control of type I error rate, but you can generally make the type I error rates quite close. In situations where you can't get them close enough together to be comparable, about the only way to make them comparable is with randomized tests. $\endgroup$
    – Glen_b
    Commented Sep 10, 2014 at 2:37
  • $\begingroup$ @ssdecontrol: No, it's just the proportion of datasets (out of the 100) where p < $\alpha$ $\endgroup$
    – EDi
    Commented Sep 10, 2014 at 8:05
  • 1
    $\begingroup$ There are two issues: (i) is that the approximations are asymptotic but $n$ is not infinity, so the approximation is just that, an approximation -- this would be an issue whether there was discreteness or not, and would lead to a significance level other than the nominal one (but if it's continuous it's something you can adjust for); (ii) there's the issue of discreteness, which prevents you from getting an exact significance level if if you do adjust for it. $\endgroup$
    – Glen_b
    Commented Sep 10, 2014 at 8:59

3 Answers 3

Reset to default
17
$\begingroup$

This is an extremely interesting problem. I reviewed your code and can find no immediately obvious typo.

I would like to see you redo this simulation but use the maximum likelihood test to make inference about the heterogeneity between groups. This would involve refitting a null model so that you can get estimates of the $\theta$s under the null hypothesis of homogeneity in rates between groups. I think this is necessary because the negative binomial model is not a linear model (the rate is parameterized linearly, but the $\theta$s are not). Therefore I am not convinced drop1 argument provides correct inference.

Most tests for linear models do not require you to recompute the model under the null hypothesis. This is because you can calculate the geometric slope (score test) and approximate the width (Wald test) using parameter estimates and estimated covariance under the alternative hypothesis alone.

Since negative binomial is not linear, I think you will need to fit a null model.

EDIT:

I edited the code and got the following: enter image description here

Edited code here: https://github.com/aomidpanah/simulations/blob/master/negativeBinomialML.r

Improve this answer
$\endgroup$
7
  • $\begingroup$ But I think that drop1() does internally re-fit the null model ... $\endgroup$
    – Ben Bolker
    Commented Sep 9, 2014 at 18:21
  • 4
    $\begingroup$ Nope! Look at how complicated the negbinom fitter is in the glm.nb code, the $\theta$s are estimated using the EM algorithm. drop1 has no default method for negative binomial. However, the logLik function does (see getS3method('logLik', 'negbin'). I have fixed the OP's code and shown that doing this gives inference that's almost identical to the Poisson model under the null (which still has a significant calibration issue due to other problems). $\endgroup$
    – AdamO
    Commented Sep 9, 2014 at 18:25
  • $\begingroup$ would like to +1 again but I can't. Nice. $\endgroup$
    – Ben Bolker
    Commented Sep 9, 2014 at 18:31
  • $\begingroup$ Thanks! I just looked the code of both drop1() and lrtest(). You're right, drop1.glm uses glm.fit which gives the wrong deviance. Wasn't aware that we can't use drop1() with glm.nb()! $\endgroup$
    – EDi
    Commented Sep 10, 2014 at 9:02
  • $\begingroup$ So the typical score and Wald tests are invalid in the negative binomial model? $\endgroup$
    – shadowtalker
    Commented Sep 10, 2014 at 11:29
10
+50
$\begingroup$

The O'Hara and Kotze paper (Methods in Ecology and Evolution 1:118–122) is not a good starting point for discussion. My most serious concern is the claim in point 4 of the summary:

We found that the transformations performed poorly, except . . .. The quasi-Poisson and negative binomial models ... [showed] little bias.

The mean $\lambda$ for a Poisson or negative binomial distribution is for a distribution that, for values of $\theta$ <= 2 and for the range of values of the mean $\lambda$ that was investigated, is highly positively skew. The means of the fitted normal distributions are on a scale of log(y+c) (c is the offset), and estimate E(log(y+c)]. This distribution is much closer to symmetric than is the distribution of y.

O'Hara and Kotze's simulations compare E(log(y+c)], as estimated by mean(log(y+c)), with log(E[y+c]). They can be, and in the cases noted are, very different. Their graphs do not compare a negative binomial with a log(y+c) fit, but rather compare mean(log(y+c)] with log(E[y+c]). On the log($\lambda$) scale shown in their graphs, it is actually the negative binomial fits that are more biased!

The following R code illustrates the point:

x <- rnbinom(10000, 0.5, mu=2)  
## NB: Above, this 'mu' was our lambda. Confusing, is'nt it?
log(mean(x+1))
[1] 1.09631
log(2+1)  ## Check that this is about right
[1] 1.098612

mean(log(x+1))
[1] 0.7317908

Or try

log(mean(x+.5))
[1] 0.9135269
mean(log(x+.5))
[1] 0.3270837

The scale on which the parameters are estimated matters a great deal!

If one samples from a Poisson, of course one expects the Poisson to do better, if judged by the criteria used to fit the Poisson. Ditto for a negative binomial. The difference may not be all that great, if the comparison is fair. Real data (e.g., maybe, in some genetic contexts) may sometimes be very close to Poisson. When they depart from Poisson, the negative binomial may or may not work well. Likewise, especially if $\lambda$ is of the order of maybe 10 or more, for modeling log(y+1) using standard normal theory.

Note that standard diagnostics work better on a scale of log(x+c). The choice of c may not matter too much; often 0.5 or 1.0 make sense. Also it is a better starting point for investigating Box-Cox transformations, or the Yeo-Johnson variant of Box-Cox. [Yeo, I. and Johnson, R. (2000)]. See further the help page for powerTransform() in R’s car package. R's gamlss package makes it possible to fit negative binomial types I (the common variety) or II, or other distributions that model the dispersion as well as the mean, with power transform links of 0 (=log, i.e., log link) or more. Fits may not always converge.

Example: Deaths vs Base Damage Data are for named Atlantic hurricanes that reached the US mainland. Data are available (name hurricNamed) from a recent release of the DAAG package for R. The help page for the data has details.

Robust loglinear vs negative binomial fit

The graph compares a fitted line obtained using a robust linear model fit, with the curve obtained by transforming a negative binomial fit with log link onto the log(count+1) scale used for the y-axis on the graph. (Note that one has to use something akin to a log(count+c) scale, with positive c, to show the points and the fitted "line" from the negative binomial fit on the same graph.) Note the large bias that is evident for the negative binomial fit on the log scale. The robust linear model fit is much less biased on this scale, if one assumes a negative binomial distribution for the counts. A linear model fit would be unbiased under the classical normal theory assumptions. I found the bias astonishing when I first created what was essentially the above graph! A curve would fit the data better, but the difference is within the bounds of the usual standards of statistical variability. The robust linear model fit does a poor job for counts at the low end of the scale.

Note --- Studies with RNA-Seq Data: Comparison of the two styles of model has been of interest for analysis of count data from gene expression experiments. The following paper compares the use of a robust linear model, working with log(count+1), with the use of negative binomial fits (as in the Bioconductor package edgeR). Most counts, in the RNA-Seq application that is primarily in mind, are large enough that suitably weighed log-linear model fits work extremely well.

Law, CW, Chen, Y, Shi, W, Smyth, GK (2014). Voom: precision weights unlock linear model analysis tools for RNA-seq read counts. Genome Biology 15, R29. http://genomebiology.com/2014/15/2/R29

NB also the recent paper:

Schurch NJ, Schofield P, Gierliński M, Cole C, Sherstnev A, Singh V, Wrobel N, Gharbi K, Simpson GG, Owen-Hughes T, Blaxter M, Barton GJ (2016). How many biological replicates are needed in an RNA-seq experiment and which differential expression tool should you use? RNA http://www.rnajournal.org/cgi/doi/10.1261/rna.053959.115

It is interesting that the linear model fits using the limma package (like edgeR, from the WEHI group) stand up extremely well (in the sense of showing little evidence of bias), relative to results with many replicates, as the number of replicates is reduced.

R code for the graph above:

library(latticeExtra, quietly=TRUE)
hurricNamed <- DAAG::hurricNamed
ytxt <- c(0, 1, 3, 10, 30, 100, 300, 1000)
xtxt <- c(1,10, 100, 1000, 10000, 100000, 1000000 )
funy <- function(y)log(y+1)
gph <- xyplot(funy(deaths) ~ log(BaseDam2014), groups= mf, data=hurricNamed,
   scales=list(y=list(at=funy(ytxt), labels=paste(ytxt)),
           x=list(at=log(xtxt), labels=paste(xtxt))),
   xlab = "Base Damage (millions of 2014 US$); log transformed scale",
   ylab="Deaths; log transformed; offset=1",
   auto.key=list(columns=2),
   par.settings=simpleTheme(col=c("red","blue"), pch=16))
gph2 <- gph + layer(panel.text(x[c(13,84)], y[c(13,84)],
           labels=hurricNamed[c(13,84), "Name"], pos=3,
           col="gray30", cex=0.8),
        panel.text(x[c(13,84)], y[c(13,84)],
           labels=hurricNamed[c(13,84), "Year"], pos=1, 
           col="gray30", cex=0.8))
ab <- coef(MASS::rlm(funy(deaths) ~ log(BaseDam2014), data=hurricNamed))

gph3 <- gph2+layer(panel.abline(ab[1], b=ab[2], col="gray30", alpha=0.4))
## 100 points that are evenly spread on a log(BaseDam2014) scale
x <- with(hurricNamed, pretty(log(BaseDam2014),100))
df <- data.frame(BaseDam2014=exp(x[x>0])) 
hurr.nb <- MASS::glm.nb(deaths~log(BaseDam2014), data=hurricNamed[-c(13,84),])
df[,'hatnb'] <- funy(predict(hurr.nb, newdata=df, type='response'))
gph3 + latticeExtra::layer(data=df,
       panel.lines(log(BaseDam2014), hatnb, lwd=2, lty=2, 
           alpha=0.5, col="gray30"))
Improve this answer
$\endgroup$
5
  • 2
    $\begingroup$ Thanks for your comment Mr. Maindonald. In the last two years there were also some more papers (focusing more on hypothesis testing, then bias): Ives 2015, Warton et al 2016, Szöcs 2015. $\endgroup$
    – EDi
    Commented May 28, 2016 at 9:47
  • $\begingroup$ maybe it's a good starting point for discussion, even if this particular point is problematic? (I would argue more generally that this is a reason not to focus on bias too much, but rather to consider something like RMSE ... [disclaimer, I haven't reread these papers lately, and I've only read the abstract of the Warton paper ...] $\endgroup$
    – Ben Bolker
    Commented May 29, 2016 at 23:30
  • 1
    $\begingroup$ Warton et al's (2016) point, that data properties should be the grounds for choice, is crucial. Quantile-quantile plots are a good way to compare the details of the fits. In particular the fit at one or other or both extremes can be important for some applications. Zero-inflated or hurdle models can be an effective refinement for getting zero counts right. At the upper extreme, any of the models under discussion may be severely compromised. Warton et al do, commendably, have one example. I'd like to see comparisons across a wide range of ecological data sets. $\endgroup$
    – John Maindonald
    Commented May 30, 2016 at 23:10
  • $\begingroup$ But aren't in ecological datasets the species in the lower part (=rare species) interesting? Shouldn't be too hard to compile some ecological datasets and compare... $\endgroup$
    – EDi
    Commented Jun 5, 2016 at 22:07
  • $\begingroup$ Actually, it is for the low end of the damage category that the negative binomial model seems, for the hurricane deaths data, to be least satisfactory. R's gamlss package has a function that makes it easy to compare centiles of the fitted distribution with centiles of the data: $\endgroup$
    – John Maindonald
    Commented Aug 4, 2016 at 4:04
7
$\begingroup$

The original post reflects Tony Ives' paper: Ives (2015). It's clear that significance testing gives different results to parameter estimation.

John Maindonald explains why the estimates are biased, but his ignorance of the background is annoying - he criticises us for showing that a method we all agree is flawed is flawed. A lot of ecologists do blindly log transform, and we were trying to point out the problems with doing that.

There's a more nuanced discussion here: Warton (2016)

Ives, A. R. (2015), For testing the significance of regression coefficients, go ahead and log-transform count data. Methods Ecol Evol, 6: 828–835. doi:10.1111/2041-210X.12386

Warton, D. I., Lyons, M., Stoklosa, J. and Ives, A. R. (2016), Three points to consider when choosing a LM or GLM test for count data. Methods Ecol Evol. doi:10.1111/2041-210X.12552

Improve this answer
$\endgroup$
5
  • $\begingroup$ Welcome to CV. While helpful, this response is mostly of the "link-only" type answer. Links do change and de-link. It would be more helpful to the CV if you were to elaborate on the key points in each one. $\endgroup$
    – user78229
    Commented May 28, 2016 at 11:57
  • $\begingroup$ Thanks for the response. I think the paper of Warton et al. coins the current state of discussion. $\endgroup$
    – EDi
    Commented May 28, 2016 at 12:10
  • $\begingroup$ Thanks & welcome! I've taken the liberty of adding the references in full. $\endgroup$
    – Scortchi ♦
    Commented May 28, 2016 at 22:21
  • 1
    $\begingroup$ Please outline the main points being made in the new references, and where it makes sense to, also relate them back to the original question. This is a valuable contribution but at present is closer to a comment on another answer than an answer to the question (which should provide context for links, for example). A few additional sentences of context would help the post substantially. $\endgroup$
    – Glen_b
    Commented May 28, 2016 at 23:08
  • 3
    $\begingroup$ Specifically, my comments address point 4 in the O'Hara and Kotze paper: "We found that the transformations performed poorly, except . . .. The quasi-Poisson and negative binomial models ... [showed] little bias." The simulations are a comment on the comparison between the expected mean on a scale of y (the counts), versus the expected mean on a scale of log(y+c), for a highly positively skew distribution, nothing more. The negative binomial parameter lambda is unbiased on the scale of y, while the log-normal mean is unbiased (under normality on that scale) on a scale of log(y+c) . $\endgroup$
    – John Maindonald
    Commented May 29, 2016 at 22:35

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.