A function called within `QL.fit`

to fit a negative binomial GLM of each gene for given design matrix

1 2 |

`counts` |
RNA-seq data matrix of integer expression counts. Each row contains observations from a single gene. Each column contains observations from a single experimental unit. |

`design` |
A single element from the |

`log.offset` |
A vector of log-scale, additive factors used to adjust estimated log-scale means for differences in library sizes across samples. Commonly used offsets include, |

`nb.disp` |
estimated negative binomial dispersion parameters obtained from either |

`print.progress` |
logical. If TRUE, the function will provide an update on what gene (row number) is being analyzed. Updates occur frequently to start then eventually occur every 5000 genes. |

`bias.fold.tolerance` |
A numerical value no smaller than 1. If the bias reduction of maximum likelihood estimates of (log) fold change is likely to result in a ratio of fold changes greater than this value, then bias reduction will be performed on such genes. Setting |

This functions fits, for each row of `counts`

, a negative binomial log-linear model through GLM framework with the over-dispersion parameter fixed. Asymptotic biases of regression coefficients (i.e., log fold changes) are then estimated by a plug-in estimate [eqn. (15.4) of McCullagh and Nelder, 1989] from the last iteration of iteratively re-weighted least squares (IWLS) procedure. The fitted response values are then compared with or without such a bias term. If the ratio of fitted response values are larger than `bias.fold.tolerance`

for any observation, then the bias-reduction (not bias-correction) procedure according to Firth (1993) and Kosmidis & Firth (2009) is applied to such rows of `counts`

, by adjusting the score equation with a term based on the observed information. Such bias-reduced estimates are more advantageous than directly subtracting the estimated bias from the maximum likelihood estimates as the latter may not exist (e.g., when all counts in one treatment group are zeros).

list containing:

`dev` |
vector containing the deviance for each gene under a negative binomial model fit to design matrix specified by |

`means` |
matrix of fitted mean values for each gene |

`parms` |
matrix of estimated coefficients for each gene. Note that these are given on the log scale. (i.e. intercept coefficient reports log(average count) and non-intercept coefficients report estimated (and bias reduced, when appropriate) log fold-changes.) |

Steve Lund (lundsp@gmail.com), Long Qu (long.qu@wright.edu)

Firth (1993) "Bias reduction of maximum likelihood estimates" *Biometrika*, **80**, 27–38.

Kosmidis and Firth (2009) "Bias reduction in exponential family nonlinear models" *Biometrika*, **96**, 793–804.

Lund, Nettleton, McCarthy and Smyth (2012) "Detecting differential expression in RNA-sequence data using quasi-likelihood with shrunken dispersion estimates" emphSAGMB, **11**(5).

McCullagh and Nelder (1989) "Generalized Linear Models", 2nd edition. London: Chapman and Hall.

Questions? Problems? Suggestions? Tweet to @rdrrHQ or email at ian@mutexlabs.com.

All documentation is copyright its authors; we didn't write any of that.