Dispersion Modeling

Let Y_ij denote an RNA-Seq read count for the i^th gene (i = 1, ⋯, m) of the j^th experimental or observational unit (j = 1, ⋯, n), and X_j the associated p-dimensional explanatory variable. Suppose Y_ij ∼ NB(μ_ij, ϕ_ij) where μ_ij is the mean and ϕ_ij is the dispersion parameter in the NB2 parameterization. Suppose also that with ${\pi }_{ij}=\exp (X_j^{\prime }{\beta }_i)$, where s_j is the library size (the number of RNA-Seq reads in the biological sample from unit j), and R_j is an optional normalization factor estimated beforehand [6, 10, 11] and treated as known. In this formulation, π_ij is the mean relative frequency of occurrence of RNA-Seq reads associated with gene i, which is taken to be the expression level of gene i associated with observational or experimental unit j.

We label some of the ways to model the nuisance parameters ϕ_ij as follows:

1. Genewise: ϕ_ij = ϕ_i (constant within each gene i across all conditions j), with m parameters for NB dispersion.
2. Common: ϕ_ij = ϕ (constant for all gene/condition combinations), with one parameter for NB dispersion.
3. NBP: log(ϕ_ij) = α₀ + α₁ log(π_ij), equivalent to assuming NBP response distribution discussed in Di et al. [7], with two parameters for NB dispersion.

We also introduce here a new approach, in which the dispersion parameter trend is quadratic on the log scale:

1. 4. NBQ: log(ϕ_ij) = α₀ + α₁ log(π_ij) + α₂ [log(π_ij)]², with three parameters for NB dispersion.

An important related method estimates the ϕ_ij’s via non-parametric regression:

1. 5. Non-parameteric: ϕ_ij is estimated in a first step as a smooth function of $\log \left({\widehat{\varphi }}_{ij}\right)$ on $\log \left({\widehat{\mu }}_{ij}\right)$, and then treated as known in the second step of regression coefficient inference.

In addition, there are variants that use an average of trend and individually-estimated dispersion parameters, based on empirical Bayes considerations [12]:

1. 6. Tagwise-common: ϕ_ij is estimated as a weighted average of the common and genewise estimates, based on empirical Bayes calculations.
2. 7. Tagwise-trend: ϕ_ij is estimated as a weighted average of the non-parametric and genewise estimates, based on empirical Bayes calculations.

Methods for inference from the genewise, common, non-parametric, tagwise-common, and tagwise-trend approaches are available in the edgeR Bioconductor package [13, 14]. The non-parametric method is also available in the DESeq and NBPSeq packages [6, 15]. The NBP and NBQ approaches are implemented in NBPSeq [15, 16].

The details of estimation for these methods are important but are not relevant to the proposed diagnostic tools and so are not discussed here. The adequacy of the models for RNA-Seq data is not yet well understood. We wish to use the model diagnostic tools proposed in this article to judge the degree of fit of the various models on real RNA-Seq data—particularly the fit of simple parametric models for the trend of log(ϕ) as a function of π and the degree of noise, if any, about this trend, so that realistic robustness and power studies can follow.

To further clarify this point, Fig. 1 shows a log-log scatter plot of method-of-moments-like estimated NB2 dispersion parameters, $\widehat{\varphi }$, versus estimated mean relative frequencies, $\widehat{\pi }$, for each of 19,623 genes from a single sample of size three of a pilot Arabidopsis RNA-Seq study examined in [7]. The curves on the plot are estimated dispersion trends based on the models described above. Polynomial gamma log-linear regression models of $\widehat{\varphi }$ on $\log (\widehat{\pi })$ were used for quick-and-dirty testing and quantification of the trend, as follows. The linear model explains 24.1% of the variability in logged dispersion parameter estimates. A quadratic term (with p-value < 0.0001) explains an additional 7.2% of variability. A cubic term (p-value < 0.0001) in a full cubic model explains less than 0.1% additional variability. This plot and informal analysis suggest that the common ϕ model is inadequate; the trend is primarily, but not entirely, linear; and that a quadratic model captures essentially all of the trend in this particular dataset.

A simple model for trend in NB dispersion parameter ϕ as a function of mean relative frequency π is a good starting point for reducing the number of nuisance parameters, but the evidence of a trend does not imply that the ϕ’s fall exactly on the trend; there may be additional variability in ϕ for genes with the same value of π. The main questions we wish to address with diagnostic tools are the following: (1) Does the NB assumption hold for a very rich model (for both regression and dispersion)? (2) What relatively simple models are adequate for describing ϕ as a function of π? (3) Is there evidence of additional biological variability in ϕ between genes having the same value of π?

Other Related Work on Model Diagnostics

Best et al. [17] extended Anscombe’s tests of fit for the NB distribution by using fourth order smooth tests, but these tests don’t extend in an obvious way to regression models for non-exponential family response distributions. The test we propose in this paper gives similar results to theirs for independent and identically distributed samples and can also be used for the procedures that involve non-parametric trend fitting and empirical Bayes averaging. Esnaola et al. [18] proposed a larger family of response distributions for RNA-Seq analysis, which permits the testing of NB as a special case; but we do not believe the approach (validated under extensively replicated experiments) is suitable for the small sample sizes we have in mind here.

Similar regression diagnostic approaches that use Monte Carlo or resampling to derive null sampling distributions of diagnostic quantities have been previously proposed for several situations. For ordinary linear regression, Atkinson [19] proposed half normal plots of jackknife residuals. For logistic regression, Landwehr et al. [20] proposed an “empirical probability plot” in which ordered residuals from the observed data are plotted against their expected values (or median values), as computed by Monte Carlo simulations. Their simulation procedure, which resembles parametric bootstrapping, is based on the estimated parameters from the fitted model. Similar graphical displays were adopted as informal checks of various count models. For example, Svetliza et al. [21] considered normal probability plots for log-linear Poisson, log-linear NB and non-linear NB models. Garay et al. [22] evaluated GOF between zero-inflated Poisson (ZIP) and zero-inflated NB (ZINB) models. Both of these used simulated envelopes in their plots, but with standard normal quantiles (instead of quantiles from simulations) on the x-axis. None of the aforementioned papers provided statistical tests for evaluating model lack-of-fit.

GOF Tests for Univariate NB Regression and the Empirical Probability Plot

We first consider univariate NB regression in this section and then return to the RNA-Seq problem of NB regression for each of many genes later. For regression data with counted response, we wish to determine whether any NBP model fits and, because of the convenience of NB2 estimation programs, whether the NB2 model fits in particular. We propose two GOF tests and an associated residual plot. The null hypothesis is that the count data follow the assumed NB regression model. In particular, the means follow the log-linear regression model, and the dispersion parameters follow the specified dispersion model (i.e., NB2, NBP, etc.). The test p-values provide an overall assessment of fit and the plot shows whether a small GOF p-value might be due to a small portion of the data. We use the same notation as in the “Background/Dispersion Modeling” subsection, but without the subscript i. In the RNA-Seq context, the methods of this section apply to a single gene. We start with Pearson residuals: $r_j=\left(y_j-{\widehat{\mu }}_j\right)/{\widehat{s}}_j$ , where μ̂_j is the estimated NB mean and ŝ_j is the estimated NB standard deviation of y_j from the particular model being tested, for j = 1, ⋯, n.

We first propose an empirical probability plot of the ordered Pearson residuals r_(j) versus the sampling distribution medians for each ordered Pearson residual, Med[r_(j)], assuming the proposed NB regression model is correct. To approximate the medians, we simulate a large number of NB regression datasets of the same size and form as the observed one, using the data estimates as parameters for simulation; fit the same NB regression model to each simulated dataset; extract the ordered Pearson residuals; and retain the sample medians for each ordered residual. This is exactly the Landwehr et al.’s approach [20] applied to NB regression. A 95% pointwise prediction envelope (in dashed blue lines) can be formed from the similarly estimated 2.5^th and 97.5^th percentiles of the ordered residuals. Fig. 2 from the Results section shows an empirical probability plot for an earthquake event dataset.

We then wish to provide a global GOF test to accompany the empirical probability plot of Pearson residuals. A natural starting point is a test based on the Pearson statistic, i.e., the sum of squared Pearson residuals. The classical use of the χ² reference distribution is not appropriate here, but the null sampling distribution may be approximated by the Monte Carlo estimate. A p-value can be obtained as the proportion of simulated samples that produce a Pearson statistic as extreme or more extreme than the observed one. This is similar to the approach of Best et al. [17] in its application of parametric bootstrap to obtain a GOF p-value. We have found that our procedure and theirs give very similar results for samples with a common mean, but the approach discussed in Best et al. [17] requires estimation of higher order moments, which is difficult for the regression models we have in mind for RNA-Seq analysis.

Since the simulations provide estimated sampling distributions for each ordered residual, a finer test statistic is available as the sum of squared differences of the ordered residuals from their sampling distribution medians. We believed this was a worthwhile test statistic to consider given that we had already obtained approximate sampling distributions for each ordered residual to obtain the type of diagnostic plot discussed in Landwehr et al. [20]. The test is also related to the SAM graphical procedure of Tusher et al. [23] for identifying differential gene expression from microarray.

The following algorithm defines the diagnostic empirical probability plot of residuals and the Monte Carlo GOF test p-value based on the second test statistic.

1. #1. Fit an NB regression model from the data Y⁽⁰⁾ = (Y₁, ⋯, Y_n)^T; estimate all unknown dispersion parameters, e.g., ${\widehat{\alpha }}_0$ , ${\widehat{\alpha }}_1$ , ⋯, and regression coefficients ${\widehat{\beta }}^{(0)}$; calculate Pearson residuals ${\mathbf{\text{r}}}^{(0)}=(r_1^{(0)},\cdots ,r_n^{(0)})$ and the mean vector ${\widehat{\mu }}^{(0)}$.
2. #2. For h = 1, ⋯, R:
   1. Simulate a random vector Y^(h) from $\mathrm{\text{NB}}({\widehat{\mu }}^{(0)},{\widehat{\alpha }}_0,{\widehat{\alpha }}_1,\cdots )$.
   2. Compute and retain Pearson residuals r^(h).
3. #3. Find the median, 2.5^th and 97.5^th percentiles of the Monte Carlo sampling distribution for each ordered residual, denoted by ${\tilde{r}}_{\left(j\right)}^{50}$ , ${\tilde{r}}_{\left(j\right)}^{2.5}$ and ${\tilde{r}}_{\left(j\right)}^{97.5}$ , respectively; plot the ordered residuals from the observed data against the Monte Carlo medians; draw a 95% pointwise prediction envelope on the plot by connecting the ${\tilde{r}}_{\left(j\right)}^{2.5}$ ’s (for lower bound) and the ${\tilde{r}}_{\left(j\right)}^{97.5}$ ’s (for upper bound).
4. #4. Compute the sum of squared deviations of ordered residuals from the medians of their sampling distributions $d^{(h)}={\displaystyle \sum }_{j=1}^n{\left(r_{(j)}^{(h)}-{\tilde{r}}_{(j)}^{50}\right)}^2$ for the observed data (h = 0) and for the simulated samples (h = 1, ⋯, R); compute a Monte Carlo GOF test p-value by: (1) where 𝟙(A) is the indicator function equal to 1 if the event A is true and 0 otherwise [24].

The Pearson GOF p-value is computed in the same way, but using the sum of squared residuals as the test statistic. The two statistics are visualized on the empirical probability plots (see, e.g., Fig. 2) as the sum of squared deviations about the y = 0 line and the sum of squared deviations about the y = x dotted line. For this reason, we call the latter statistic the sum of squared vertical distances.

In using Monte Carlo simulation in lieu of theory to obtain the sampling distributions of the ordered residuals and test statistics, it would be ideal, but impossible, to use the true parameter values rather than their estimates from the data. Using estimated parameters in the simulation may lead to conservative test results, especially in small sample situations. As the sample sizes increase, we expect the sampling distribution of the residuals to be about the same.

We use R = 999 Monte Carlo samples so that the (binomial) standard error in p-value estimation is 0.016 for p-values near 0.5, 0.007 for p-values near 0.05, and 0.003 for p-values near 0.01. As pointed out in North et al. [25], adding 1 to both the numerator and the denominator in Equation (1) produces a slightly biased estimate of the true p-value but with the correct Type-I error rate, in contrast to the unbiased but anti-conservative p-value obtained without adding the 1.

Diagnostic Tools for RNA-Seq Modeling

A major step for improving RNA-Seq analysis is the comparative evaluation of the models and methods for incorporating commonalities of NB dispersion parameters within and across genes, as described in the Introduction section and displayed in Fig. 1. For studying these models on a given RNA-Seq dataset, we propose fitting them, calculating the squared vertical distance NB GOF p-value for each of a randomly selected sample of genes (i.e., testing the univariate NB regression model fit for each gene individually, using the NB2 dispersion parameter estimated according to the particular dispersion model), drawing a uniform QQ plot of the p-values, and calculating a single p-value using the Fisher’s method.

Let p_i be the GOF p-value for gene i, based on the parameter estimates from the global model being tested. Fisher’s method produces a single GOF p-value by testing the conformity of the p_i’s from m genes to a standard uniform distribution. If the m single-gene p-values are independent and all follow standard uniform distribution, the test statistic $X^2=-2{\sum }_{i=1}^m\log (p_i)$ follows a ${\chi }_{(2m)}^2$ distribution. In this multiple genes case, the null hypothesis is that the NB counts for each gene follow the assumed NB regression model (the same as in the univariate case) and the NB dispersion parameters across all genes follow the specified dispersion model.

Although it is possible to base this on all genes, we elect to reduce the computational burden by selecting a random sample m* genes, and use m* = 1,000 as a computationally tolerable value. We do not have a direct way to study the suitability of this sample size for testing whether the p-values follow a uniform(0,1) distribution with Fisher’s method; but we do have an indirect approach that helps. Let P be the proportion of genes with p-values less than 0.05. A binomial 95% confidence interval for P from a sample of 1,000 genes has half-width 0.0135, so we would be likely to detect lack-of-fit to the uniform(0,1) if the actual proportion is 0.0635 or greater. Although we use Fisher’s method rather than this (arbitrary) binomial test, the binomial calculation provides some clarification of the type of departure from the uniform(0,1) that we are likely to detect with a sample of 1,000 genes. The m p-values are not exactly independent, as required for the theory of Fisher’s combined test, but are approximately so because the global parameter estimates are based on such a large number of genes.

As we noted earlier, the Pearson statistic can give misleading results if there are combinations of under- and over-dispersion relative to the response distributional model being tested. We have found this problem to be exacerbated in the RNA-Seq setting and so focus only on the squared vertical distance estimator, which does not suffer from the same problem.

An extremely small (or large) X² value can be due to either a small number of extreme single-gene p-values or a large number of moderately small (or large) p-values. The Fisher’s method itself cannot distinguish the different possibilities. Along with the Fisher’s combination of p-values, we suggest a uniform QQ plot of individual p-values to help reveal the nature of any lack-of-fit, indicated by a higher than expected proportion of small p-values. The proportion of genes with small p-values may have some effect on the thinking about appropriate models.

The Earthquake Event Dataset

We consider an earthquake event dataset when illustrating the empirical probability plot and the GOF tests for univariate NB regression. The dataset (provided in the NBGOF package) contains the frequencies of all earthquakes of a given magnitude (reported to one decimal place) for magnitudes from 4.5 to 9.1, that occurred between January 1, 1964 to December 31, 2012 (Source: Composite Earthquake Catalog, Advanced National Seismic System, Northern California Earthquake Data Center (NCEDC), http://quake.geo.berkeley.edu/cnss/). The empirical probability plots with GOF test results, based on a log-linear regression of mean number of earthquakes on magnitude, are shown in Fig. 2. Neither the NB2 nor NBP model shows lack-of-fit.

The Arabidopsis Study

We use an Arabidopsis dataset when demonstrating the diagnostic tools for RNA-Seq modeling. Arabidopsis thaliana has been intensively studied as a model organism in plant biology. The Arabidopsis data discussed in Di et al. [7] contain RNA-Seq reads that aligned to more than 25,000 genes from two groups of Arabidopsis samples of size three each. The two groups of size three each were derived from plants inoculated with ΔhrcC of Pseudomonas syringae pv tomato DC3000 or 10 mM MgCl₂ (mock). The dataset used in this article comes from Di et al. [7], which is a subset of the data described in Cumbie et al. [26].

Simulation Parameter Specifications

We specify the parameters in the simulation studies in the subsection “Results/Error Rates of GOF Tests in Simulations” (Tables 1 and 2) as follows: the mean is determined by μ = exp(X′ β) with the coefficient β = (15,−1.5). The design matrix X takes an intercept and a covariate equally spaced from 4 to 8 of length n = 5,10,50 and 100. The resulting mean levels approximately range from 20 to 8,100.

For the NB2 model fit: the NB2 responses are simulated under α₀ = log(0.1), α₁ = 0 and ϕ = 0.1. The NB1 responses are obtained by simulating NB2 with α₀ = log(0.5), α₁ = −1, and ϕ = 0.5/μ. The “NB2 plus outliers” responses are simulated with the same dispersions as in NB2, except we randomly double 20% of responses (as outliers). The “NB2 plus noise” responses are simulated with the same dispersions as in NB2, except ϕ is specified as 0.1 ⋅ exp(G), where G ∼ 𝒩(0,1). The α in the variance μ + ϕμ^α is determined by α = α₁ + 2.

For the NBP model fit: the specifications are almost the same as in the NB2 model fit above, except for the simulated NB2 data, we use α₀ = log(0.05) so that ϕ = 0.05.

GOF Tests for Univariate NB Regression and the Empirical Probability Plot

Diagnostic Tools for RNA-Seq Modeling

In this section, we apply our proposed GOF test and diagnostic graphics to the analysis of an Arabidopsis RNA-Seq dataset.

The following are the GOF p-values (in parentheses) from fitting the seven dispersion models (described in the “Background/Dispersion Modeling” subsection) to a random sample of 1,000 genes: common (< 0.0001); NBP (0.04); NBQ (0.94); trended (0.21); genewise (> 0.9999); tagwise-common (> 0.9999) and tagwise-trend (> 0.9999). The corresponding uniform QQ plots of p-values are shown in Fig. 4.

The p-values are greater than 0.9999 for the genewise and tagwise models. These unusually large p-values indicate that our GOF test may be conservative for these models. We note that for these models, the number of dispersion parameters to be estimated increases with the number of genes. In this sense, it is inherently more challenging to judge goodness-of-fit for these models. As can be seen in the uniform QQ plots of single-gene p-values, there are fewer small p-values than expected from a uniform(0,1) distribution. The extra large Fisher’s p-value is most likely due to conservativeness in the tests from small sample sizes, as evident in Table 1, when individual (genewise and tagwise) dispersion parameters are estimated from the small sample. Even a slight degree of conservativeness in the individual NB GOF tests can produce a very small Fisher combination statistic when there are so many p-values being combined. The evidence from Table 1, and our experience with simulations and other RNA-Seq datasets, suggest that this conservativeness diminishes with increasing sample size.

The consistency of the data with the NBQ trend model and the non-parametric “trended” approach suggests that models without noise about the trend may be adequate. We need to exercise caution with this conclusion, though—since the lack of evidence for lack-of-fit does not prove “fit” and that the test may be conservative at this sample size. Nevertheless, the apparent fit of NBQ in conjunction with the evidence of lack-of-fit for NBP is intriguing.

Fig. 1 shows the mean-dispersion plot (log-log scale) with six fitted dispersion models (common, NBP, NBQ, trended, tagwise-common and tagwise-trend) based on the mock treatment group alone. The genewise estimates are not included since there is no implied trend associated with that method.

In addition to the real RNA-Seq data analyses, we also performed simulation studies where the datasets are generated according to the “NB2 + noise” model. We then tested for GOF of the simple models and the tagwise approaches. Details on the simulation specifications and the results are provided in Supporting Information S1 File. The simulation results are as expected from the simulation setting.

Software Information

The proposed approach is implemented as an R package named NBGOF (version 0.1.6, available in the Unix-like platforms) released at the first author’s github page: https://github.com/gu-mi/NBGOF, under GPL-2 License. The package also includes all datasets analyzed in this article. The R codes for reproducing all results in this article are available at the first author’s github page.