Set alpha=Inf as the default for testEmptyDrops. (#118)

This is because a finite alpha is not universally safer; despite
representing an overdispersed multinomial, using it to compute p-values
for a multinomial-distributed vector will actually be anticonservative!

This means that an inaccurate alpha cannot be brushed under the carpet.
Previously, I was trusting that any finite alpha would be safer than
alpha=Inf, but it seems that this is not true. We can't just ignore poor
estimates of alpha that we get from the assumed-ambient counts.

So, we keep it simple and set the default to alpha=Inf, which is exactly
correct for the multinomial case, i.e., independent sampling of molecules into
droplet (which should represent most scenarios, hopefully). This also reduces
runtime and aligns with the CellRanger defaults.

Author: LTLA

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: DropletUtils
-Version: 1.27.1
-Date: 2024-12-04
+Version: 1.27.2
+Date: 2024-12-10
 Title: Utilities for Handling Single-Cell Droplet Data
 Authors@R: c(
     person("Aaron", "Lun", role = "aut"),

R/emptyDrops.R

@@ -75,16 +75,20 @@
 #' The latter is a \dQuote{no frills} version that is largely intended for use within other functions.
 #' 
 #' @section Handling overdispersion:
+#' By default, \code{alpha=Inf} which means that the count vector is modelled with a multinomial distribution.
+#' This is appropriate when molecules are independently sampled into each droplet. 
+#'
 #' If \code{alpha} is set to a positive number, sampling is assumed to follow a Dirichlet-multinomial (DM) distribution.
 #' The parameter vector of the DM distribution is defined as the estimated ambient profile scaled by \code{alpha}.
-#' Smaller values of \code{alpha} model overdispersion in the counts, due to dependencies in sampling between molecules.
+#' Smaller values of \code{alpha} model overdispersion in the counts, due to dependencies in sampling between molecules (e.g., aggregates, PCR duplication).
 #' If \code{alpha=NULL}, a maximum likelihood estimate is obtained from the count profiles for all barcodes with totals less than or equal to \code{lower}.
-#' If \code{alpha=Inf}, the sampling of molecules is modelled with a multinomial distribution.
 #' 
 #' Users can check whether the model is suitable by extracting the p-values for all barcodes with \code{test.ambient=TRUE}.
 #' Under the null hypothesis, the p-values for presumed ambient barcodes (i.e., with total counts less than or equal to \code{lower}) should be uniformly distributed.
-#' Skews in the p-value distribution are indicative of an inaccuracy in the model and/or its estimates (of \code{alpha} or the ambient profile).
-#' 
+#' Skews in the p-value distribution are indicative of an inaccuracy in the model.
+#' For example, an inaccurate \code{alpha} or ambient profile will manifest in the overenrichment of low p-values.
+#' Conversely, very sparse data will often exhibit in a enrichment of p-values at 1 as the Good-Turing probabilities in the ambient profile cannot be zero.
+#'
 #' @section \code{NA} values in the results:
 #' We assume that barcodes with total UMI counts less than or equal to \code{lower} correspond to empty droplets.
 #' These are used to estimate the ambient expression profile against which the remaining barcodes are tested.
@@ -181,7 +185,7 @@ NULL
 
 #' @export
 #' @rdname emptyDrops
-testEmptyDrops <- function(m, lower=100, niters=10000, test.ambient=FALSE, ignore=NULL, alpha=NULL, round=TRUE, by.rank=NULL, known.empty=NULL, BPPARAM=SerialParam()) {
+testEmptyDrops <- function(m, lower=100, niters=10000, test.ambient=FALSE, ignore=NULL, alpha=Inf, round=TRUE, by.rank=NULL, known.empty=NULL, BPPARAM=SerialParam()) {
     ambfun <- function(mat, totals) {
         assumed.empty <- .get_putative_empty(totals, lower, by.rank, known.empty)
         astats <- .compute_ambient_stats(mat, totals, assumed.empty)

inst/NEWS.Rd

@@ -6,6 +6,10 @@
 \item Use a more stable algorithm for identifying the knee point in \code{barcodeRanks()}.
 The new algorithm is based on maximizing the distance from a line between the plateau and the inflection point.
 Previously, we tried to minimize the signed curvature but this was susceptible to many local minima due to the instability of the empirical second derivative, even after smoothing.
+
+\item Set \code{alpha=Inf} as the default for \code{testEmptyDrops()}.
+This is motivated by the realization that an underestimated \code{alpha} can still yield anticonservative p-values and is not universally safer than \code{alpha=Inf}.
+Defaulting \code{alpha=Inf} is preferable as it is at least correct in the expected case of multinomial sampling.
 }}
 
 \section{Version 1.18.0}{\itemize{