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(2+1) ## Check that this is about right
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.
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
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:
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,
xlab = "Base Damage (millions of 2014 US$); log transformed scale",
ylab="Deaths; log transformed; offset=1",
gph2 <- gph + layer(panel.text(x[c(13,84)], y[c(13,84)],
labels=hurricNamed[c(13,84), "Name"], pos=3,
labels=hurricNamed[c(13,84), "Year"], pos=1,
ab <- coef(MASS::rlm(funy(deaths) ~ log(BaseDam2014), data=hurricNamed))
gph3 <- gph2+layer(panel.abline(ab, b=ab, 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,