Description Usage Arguments Details Value Author(s) References Examples

A function called within `QL.fit`

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

1 |

`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 |

`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 ([email protected]), Long Qu ([email protected])

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.

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 | ```
## Not run:
## no bias reduction
noReduction =
with(mockRNASeqData,
NBDev(counts, design.matrix, log(estimated.normalization),
estimated.nbdisp, bias.fold.tolerance=Inf)
)
## bias reduction for all genes
allReduction =
with(mockRNASeqData,
NBDev(counts, design.matrix, log(estimated.normalization),
estimated.nbdisp, bias.fold.tolerance=1)
)
## default: bias reduction for genes showing large bias
someReduction =
with(mockRNASeqData,
NBDev(counts, design.matrix, log(estimated.normalization),
estimated.nbdisp, bias.fold.tolerance=1.1)
)
## End(Not run)
``` |

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.