@@ -1,9 +1,9 @@

 Package: limma

-Version: 3.61.7

-Date: 2024-07-30

+Version: 3.61.8

+Date: 2024-08-02

 Title: Linear Models for Microarray Data

 Description: Data analysis, linear models and differential expression for microarray and RNA-seq data.

-Author: Gordon Smyth [cre,aut], Yifang Hu [ctb], Matthew Ritchie [ctb], Jeremy Silver [ctb], James Wettenhall [ctb], Davis McCarthy [ctb], Di Wu [ctb], Wei Shi [ctb], Belinda Phipson [ctb], Aaron Lun [ctb], Natalie Thorne [ctb], Alicia Oshlack [ctb], Carolyn de Graaf [ctb], Yunshun Chen [ctb], Goknur Giner [ctb], Mette Langaas [ctb], Egil Ferkingstad [ctb], Marcus Davy [ctb], Francois Pepin [ctb], Dongseok Choi [ctb], Charity Law [ctb], Mengbo Li [ctb]

+Author: Gordon Smyth [cre,aut], Yifang Hu [ctb], Matthew Ritchie [ctb], Jeremy Silver [ctb], James Wettenhall [ctb], Davis McCarthy [ctb], Di Wu [ctb], Wei Shi [ctb], Belinda Phipson [ctb], Aaron Lun [ctb], Natalie Thorne [ctb], Alicia Oshlack [ctb], Carolyn de Graaf [ctb], Yunshun Chen [ctb], Goknur Giner [ctb], Mette Langaas [ctb], Egil Ferkingstad [ctb], Marcus Davy [ctb], Francois Pepin [ctb], Dongseok Choi [ctb], Charity Law [ctb], Mengbo Li [ctb], Lizhong Chen [ctb]

 Maintainer: Gordon Smyth <[email protected]>

 License: GPL (>=2)

 Depends: R (>= 3.6.0)

@@ -15,13 +15,13 @@ importFrom("stats", "approx", "approxfun", "as.dendrogram", "as.dist",

            "dist", "density", "df", "dhyper", "dnorm", "filter", "fitted",

            "hat", "hclust", "integrate", "isoreg", "lm", "lm.fit", "lm.wfit", "loess",

            "lowess", "median", "model.matrix", "na.exclude", "na.omit",

-           "nlminb", "nls", "nls.control", "optim", "p.adjust",

+           "nlminb", "nls", "nls.control", "optim", "optimize", "p.adjust",

            "pbeta", "pchisq", "pexp", "pf", "pgamma", "phyper",

            "pnorm", "ppoints", "predict", "pt", "qchisq", "qf",

            "qnorm", "qt", "quantile", "residuals", "rnorm", "sd",

            "uniroot", "var", "weighted.mean")

 importFrom("utils", "getFromNamespace", "head", "read.delim", "read.table", "tail", "write.table")

-importFrom("statmod", "logmdigamma", "gauss.quad.prob", "mixedModel2Fit", "remlscore")

+importFrom("statmod", "gauss.quad.prob", "logmdigamma", "mixedModel2Fit", "remlscore")

 if( tools:::.OStype() == "windows" ) importFrom("utils", "winMenuAddItem")

 S3method("[",MAList)

@@ -6,18 +6,21 @@ backgroundCorrect <- function(RG, method="auto", offset=0, printer=RG$printer, n

 {

 #	Apply background correction to microarray data

 #	Gordon Smyth

-#	12 April 2003.  Last modified 30 August 2010.

+#	12 April 2003.  Last modified 1 August 2024.

 	if(is.null(method)) method <- "auto"

+#	Single-channel background correction

 	if(is(RG,"EListRaw")) {

 		RG$E <- backgroundCorrect.matrix(RG$E,RG$Eb,method=method,offset=offset,normexp.method=normexp.method,verbose=verbose)

 		RG$Eb <- NULL

 		return(RG)

 	}

-	if(class(RG)=="list" && !is.null(RG$R)) RG <- new("RGList",RG)

+#	If RG is an unclassed list with a red channel, then convert to RGList

+	if(!isS4(RG) && is.list(RG) && !is.null(RG$R)) RG <- new("RGList",RG)

+#	Two-color background correction

 	if(is(RG,"RGList")) {

 		RG$R <- backgroundCorrect.matrix(RG$R,RG$Rb,method=method,offset=offset,normexp.method=normexp.method,verbose=verbose)

 		RG$Rb <- NULL

@@ -26,6 +29,7 @@ backgroundCorrect <- function(RG, method="auto", offset=0, printer=RG$printer, n

 		return(RG)

 	}

+#	Failing all other classes, coerce RG to a matrix

 	RG <- as.matrix(RG)

 	if(!(method %in% c("none","normexp"))) stop(method,"correction requires background intensities")

 	RG <- backgroundCorrect.matrix(RG,method=method,offset=offset,normexp.method=normexp.method,verbose=verbose)

@@ -1,14 +1,14 @@

 #  EMPIRICAL BAYES FUNCTIONS

-eBayes <- function(fit,proportion=0.01,stdev.coef.lim=c(0.1,4),trend=FALSE,robust=FALSE,winsor.tail.p=c(0.05,0.1))

+eBayes <- function(fit,proportion=0.01,stdev.coef.lim=c(0.1,4),trend=FALSE,robust=FALSE,winsor.tail.p=c(0.05,0.1),legacy=NULL)

 #	Empirical Bayes statistics to select differentially expressed genes.

 #	Accepts and returns an MArrayLM object.

 #	Gordon Smyth

-#	4 August 2003.  Last modified 29 Nov 2022.

+#	Created 4 Aug 2003.  Last modified 2 Aug 2024.

 {

 	if(!is.list(fit)) stop("fit is not a valid MArrayLM object")

 	if(is.logical(trend) && trend && is.null(fit$Amean)) stop("Need Amean component in fit to estimate trend")

-	eb <- .ebayes(fit=fit,proportion=proportion,stdev.coef.lim=stdev.coef.lim,trend=trend,robust=robust,winsor.tail.p=winsor.tail.p)

+	eb <- .ebayes(fit=fit,proportion=proportion,stdev.coef.lim=stdev.coef.lim,trend=trend,robust=robust,winsor.tail.p=winsor.tail.p,legacy=legacy)

 	fit$df.prior <- eb$df.prior

 	fit$s2.prior <- eb$s2.prior

 	fit$var.prior <- eb$var.prior

@@ -28,11 +28,10 @@ eBayes <- function(fit,proportion=0.01,stdev.coef.lim=c(0.1,4),trend=FALSE,robus

 	fit

 }

-.ebayes <- function(fit,proportion=0.01,stdev.coef.lim=c(0.1,4),trend=FALSE,robust=FALSE,winsor.tail.p=c(0.05,0.1))

+.ebayes <- function(fit,proportion=0.01,stdev.coef.lim=c(0.1,4),trend=FALSE,robust=FALSE,winsor.tail.p=c(0.05,0.1),legacy=NULL)

 #	Empirical Bayes statistics to select differentially expressed genes

 #	Gordon Smyth

-#	Created 8 Sept 2002.  Last revised 10 Feb 2022.

-#	Made a non-exported function 18 Feb 2018.

+#	Created 8 Sep 2002. Made non-exported function 18 Feb 2018. Last revised 2 Aug 2024.

 {

 	coefficients <- fit$coefficients

 	stdev.unscaled <- fit$stdev.unscaled

@@ -55,7 +54,7 @@ eBayes <- function(fit,proportion=0.01,stdev.coef.lim=c(0.1,4),trend=FALSE,robus

 	}

 #	Moderated t-statistic

-	out <- squeezeVar(sigma^2, df.residual, covariate=covariate, robust=robust, winsor.tail.p=winsor.tail.p)

+	out <- squeezeVar(sigma^2, df.residual, covariate=covariate, robust=robust, winsor.tail.p=winsor.tail.p, legacy=legacy)

 	out$s2.prior <- out$var.prior

 	out$s2.post <- out$var.post

 	out$var.prior <- out$var.post <- NULL

@@ -1,9 +1,9 @@

 #	EMPIRICAL BAYES SQUEEZING OF VARIANCES

-squeezeVar <- function(var, df, covariate=NULL, robust=FALSE, winsor.tail.p=c(0.05,0.1))

+squeezeVar <- function(var, df, covariate=NULL, robust=FALSE, winsor.tail.p=c(0.05,0.1), legacy=NULL)

 #	Empirical Bayes posterior variances

 #	Gordon Smyth

-#	Created 2 March 2004.  Last modified 29 July 2024.

+#	Created 2 March 2004.  Last modified 1 August 2024.

 {

 	n <- length(var)

@@ -14,13 +14,25 @@ squeezeVar <- function(var, df, covariate=NULL, robust=FALSE, winsor.tail.p=c(0.

 #	When df==0, guard against missing or infinite values in var

 	if(length(df)>1L) var[df==0] <- 0

+#	Choose legacy or new method depending whether df are unequal

+	if(is.null(legacy)) {

+		dfp <- df[df>0]

+		legacy <- identical(min(dfp),max(dfp))

+	}

+

 #	Estimate hyperparameters

-	if(robust) {

-		fit <- fitFDistRobustly(var, df1=df, covariate=covariate, winsor.tail.p=winsor.tail.p)

-		df.prior <- fit$df2.shrunk

+	if(legacy) {

+		if(robust) {

+			fit <- fitFDistRobustly(var, df1=df, covariate=covariate, winsor.tail.p=winsor.tail.p)

+			df.prior <- fit$df2.shrunk

+		} else {

+			fit <- fitFDist(var, df1=df, covariate=covariate)

+			df.prior <- fit$df2

+		}

 	} else {

-		fit <- fitFDist(var, df1=df, covariate=covariate)

-		df.prior <- fit$df2

+		fit <- fitFDistUnequalDF1(var, df1=df, covariate=covariate, robust=robust)

+		df.prior <- fit$df.shrunk

+		if(is.null(df.prior)) df.prior <- fit$df2

 	}

 	if(anyNA(df.prior)) stop("Could not estimate prior df")

@@ -34,28 +46,29 @@ squeezeVar <- function(var, df, covariate=NULL, robust=FALSE, winsor.tail.p=c(0.

 #	Squeeze posterior variances given hyperparameters

 #	NAs not allowed in df.prior

 #	Gordon Smyth

-#	Created 5 May 2016. Last modified 29 July 2024.

+#	Created 5 May 2016. Last modified 1 August 2024.

 {

-	n <- length(var)

-	isfin <- is.finite(df.prior)

-	if(all(isfin)) return( (df*var + df.prior*var.prior) / (df+df.prior) )

-

-#	From here, at least some df.prior are infinite

+#	If df.prior are all finite, use canonical formula

+	m <- max(df.prior)

+	if(is.finite(m)) return( (df*var + df.prior*var.prior) / (df+df.prior) )

-#	For infinite df.prior, return var.prior

+#	Set var.post to var.prior of length n

+	n <- length(var)

 	if(identical(length(var.prior),n)) {

 		var.post <- var.prior

 	} else {

 		var.post <- rep_len(var.prior, length.out=n)

 	}

-#	Maybe some df.prior are finite

-	if(any(isfin)) {

-		i <- which(isfin)

-		if(length(df)>1) df <- df[i]

-		df.prior <- df.prior[i]

-		var.post[i] <- (df*var[i] + df.prior*var.post[i]) / (df+df.prior)

-	}

+#	Check if df.prior all Inf

+	m <- min(df.prior)

+	if(m > 1e100) return(var.post)

+

+#	Only some df.prior are finite

+	i <- which(is.finite(df.prior))

+	if(length(df)>1L) df <- df[i]

+	df.prior <- df.prior[i]

+	var.post[i] <- (df*var[i] + df.prior*var.post[i]) / (df+df.prior)

 	var.post

 }

@@ -1,9 +1,9 @@

 ###  treat.R

-treat <- function(fit, fc=1.2, lfc=NULL, trend=FALSE, robust=FALSE, winsor.tail.p=c(0.05,0.1))

+treat <- function(fit, fc=1.2, lfc=NULL, trend=FALSE, robust=FALSE, winsor.tail.p=c(0.05,0.1), legacy=NULL)

 #	Moderated t-statistics relative to a logFC threshold.

 #	Davis McCarthy, Gordon Smyth

-#	25 July 2008.  Last revised 6 May 2023.

+#	25 July 2008.  Last revised 2 Aug 2024.

 {

 #	Check fit

 	if(!is(fit,"MArrayLM")) stop("fit must be an MArrayLM object")

@@ -28,7 +28,7 @@ treat <- function(fit, fc=1.2, lfc=NULL, trend=FALSE, robust=FALSE, winsor.tail.

 	} else {

 		covariate <- NULL

 	}

-	sv <- squeezeVar(sigma^2, df.residual, covariate=covariate, robust=robust, winsor.tail.p=winsor.tail.p)

+	sv <- squeezeVar(sigma^2, df.residual, covariate=covariate, robust=robust, winsor.tail.p=winsor.tail.p, legacy=legacy)

 	fit$df.prior <- sv$df.prior

 	fit$s2.prior <- sv$var.prior

 	fit$s2.post <- sv$var.post

@@ -3,6 +3,24 @@

 \encoding{UTF-8}

+\section{Version 3.61.8}{\itemize{

+

+  \item

+  New function fitFDistUnequalDF1() for improved empirical Bayes

+  hyperparameter estimation when the residual degrees of freedom are

+  not all equal between genes. The new hyperpameter estimation will be

+  used by eBayes(), treat() and squeezeVar() by default whenever the

+  residual degrees of freedom are unequal, but the legacy behavior

+  can be preserved for backward compatibility by setting

+  `legacy=TRUE`.

+

+  \item

+  New argument `keep.EList` for voomaLmFit(). The observation weights

+  generated by voomaLmFit() are now stored only if `keep.EList=TRUE`.

+

+}}

+

+

 \section{Version 3.60.0}{\itemize{

   \item

@@ -1,3 +1,24 @@

+ 1 Aug 2024: limma 3.61.8

+

+- New function fitFDistUnequalDF1(), which improves on fitFDist() and

+  fitFDistRobustly(), with special attention to contexts where the

+  df1 argument is unequal between cases and can include small

+  fractional values, such as those resulting from the edgeR 4.0

+  quasi-likelihood pipeline.

+

+- New argument `legacy` for squeezeVar(). If `legacy=FALSE`, then it

+  uses fitFDistUnequalDF1() instead of fitFDist() or

+  fitFDistRobustly(). If `legacy=NULL` then fitFDistUnequalDF1() will

+  be used when the degrees of freedom `df` are not all equal.

+

+- New argument `legacy` for eBayes() and treat(). If `legacy=FALSE`,

+  then the new hyperparameter estimation implemented via

+  fitFDistUnequalDF1() will be used, otherwise the legacy method will

+  be used. If `legacy=NULL` then the new method will be used when the

+  residual degrees of freedom are not all equal.

+

+- Improved checking for the class of the input in backgroundCorrect.R.

+

 30 Jul 2024: limma 3.61.7

 - Clean up some code in fitFDistRobustly().

@@ -58,7 +58,7 @@ After fitting a linear model, the standard errors are moderated using a simple e

 A moderated t-statistic and a log-odds of differential expression is computed for each contrast for each gene.

 \code{treat} tests whether log-fold-changes are greater than a threshold rather than merely different to zero.

-\code{\link{eBayes}} and \code{\link{treat}} use internal functions \code{\link{squeezeVar}}, \code{\link{fitFDist}}, \code{\link{tmixture.matrix}} and \code{\link{tmixture.vector}}.

+\code{\link{eBayes}} and \code{\link{treat}} use internal functions \code{\link{squeezeVar}}, \code{\link{fitFDist}}, \code{\link{fitFDistRobustly}}, \code{\link{fitFDistUnequalDF1}}, \code{\link{tmixture.matrix}} and \code{\link{tmixture.vector}}.

 }

 \section{Summarizing Model Fits}{

@@ -104,7 +104,7 @@ Smyth, G. K. (2004). Linear models and empirical Bayes methods for assessing dif

 \emph{Statistical Applications in Genetics and Molecular Biology}, \bold{3}, No. 1, Article 3.

 \url{https://gksmyth.github.io/pubs/ebayes.pdf}

-Smyth, G. K., Michaud, J., and Scott, H. (2005). The use of within-array replicate spots for assessing differential expression in microarray experiments. Bioinformatics 21(9), 2067-2075.

+Smyth, G. K., Michaud, J., and Scott, H. (2005). The use of within-array replicate spots for assessing differential expression in microarray experiments. \emph{Bioinformatics} 21(9), 2067-2075.

 }

 \seealso{

@@ -7,9 +7,9 @@

 \usage{

 eBayes(fit, proportion = 0.01, stdev.coef.lim = c(0.1,4),

-       trend = FALSE, robust = FALSE, winsor.tail.p = c(0.05,0.1))

+       trend = FALSE, robust = FALSE, winsor.tail.p = c(0.05,0.1), legacy = NULL)

 treat(fit, fc = 1.2, lfc = NULL, trend = FALSE,

-      robust = FALSE, winsor.tail.p = c(0.05,0.1))

+      robust = FALSE, winsor.tail.p = c(0.05,0.1), legacy = NULL)

 }

 \arguments{

@@ -19,7 +19,8 @@ treat(fit, fc = 1.2, lfc = NULL, trend = FALSE,

   \item{stdev.coef.lim}{numeric vector of length 2, assumed lower and upper limits for the standard deviation of log2-fold-changes for differentially expressed genes}

   \item{trend}{logical, should an intensity-dependent trend be allowed for the prior variance? If \code{FALSE} then the prior variance is constant. Alternatively, \code{trend} can be a row-wise numeric vector, which will be used as the covariate for the prior variance.}

   \item{robust}{logical, should the estimation of \code{df.prior} and \code{var.prior} be robustified against outlier sample variances?}

-  \item{winsor.tail.p}{numeric vector of length 1 or 2, giving left and right tail proportions of \code{x} to Winsorize. Used only when \code{robust=TRUE}.}

+  \item{winsor.tail.p}{numeric vector of length 1 or 2, giving left and right tail proportions of \code{x} to Winsorize (if \code{robust=TRUE} and \code{legacy=TRUE}).}

+  \item{legacy}{logical. If \code{FALSE} then the new hyperparameter estimation (introduced in limma 3.61.8) will be used, otherwise the legacy methods will be used. If \code{NULL}, then the new method will be used whenever the residual degrees of freedom are not all equal.}

   \item{fc}{a minimum fold-change below which changes are not considered scientifically meaningful.}

   \item{lfc}{a minimum log2-fold-change below which changes not considered scientifically meaningful. Defaults to \code{log2(fc)}.}

 }

@@ -60,7 +61,7 @@ For each gene (row), the moderated F-statistic tests whether all the contrasts a

 The F-statistic is an overall test computed from the set of t-statistics for that probe.

 This is exactly analogous the relationship between t-tests and F-statistics in conventional anova, except that the residual mean squares have been moderated between genes.

-The estimates \code{s2.prior} and \code{df.prior} are computed by \code{fitFDist}.

+The estimates \code{s2.prior} and \code{df.prior} are computed by one of \code{fitFDist}, \code{fitFDistRobustly} or \code{fitFDistUnequalDF1} (depending on settings for \code{robust} and \code{legacy}).

 \code{s2.post} is the weighted average of \code{s2.prior} and \code{sigma^2} with weights proportional to \code{df.prior} and \code{df.residual} respectively.

 The log-odds of differential expression \code{lods} was called the \emph{B-statistic} by Loennstedt and Speed (2002).

 The F-statistics \code{F} are computed by \code{classifyTestsF} with \code{fstat.only=TRUE}.

@@ -102,6 +103,12 @@ limma-trend is recommended for RNA-seq analysis when the library sizes are reaso

 If \code{robust=TRUE} then the robust empirical Bayes procedure of Phipson et al (2016) is used.

 This is frequently useful to protect the empirical Bayes procedure against hyper-variable or hypo-variable genes, especially when analysing RNA-seq data.

 See \code{\link{squeezeVar}} for more details.

+

+In limma 3.61.8 (August 2024), the new function \code{\link{fitFDistUnequalDF1}} was introduced to improve estimation of the hyperparameters \code{s2.prior} and \code{df.prior}, especially when not all genes have the same residual degrees of freedom.

+\code{fitFDistUnequalDF1} is a potential replacement for the original functions \code{fitFDist} and \code{fitFDistRobustly} and the argument \code{legacy} is provided to control backward compatibility.

+The new hyperparameter estimation will be used if \code{legacy=FALSE} and the original methods will be used if \code{legacy=TRUE}.

+If \code{legacy=NULL}, then the new method will be used if the residual degrees of freedom are unequal and the original methods otherwise.

+Unequal residual degrees of freedom arise in limma pipelines when the expression matrix includes missing values or from the quasi-likelihood pipeline in edgeR v4.

 }

 \note{

@@ -1,14 +1,17 @@

 \name{fitFDist}

 \alias{fitFDist}

 \alias{fitFDistRobustly}

+\alias{fitFDistUnequalDF1}

+

 \title{Moment Estimation of Scaled F-Distribution}

 \description{

 Moment estimation of the parameters of a scaled F-distribution given one of the degrees of freedom.

-This function is called internally by \code{eBayes} and \code{squeezeVar} and is not usually called directly by a user.

+These functions are called internally by \code{eBayes} and \code{squeezeVar} and is not usually called directly by a user.

 }

 \usage{

-fitFDist(x, df1, covariate=NULL)

-fitFDistRobustly(x, df1, covariate=NULL, winsor.tail.p=c(0.05,0.1), trace=FALSE)

+fitFDist(x, df1, covariate = NULL)

+fitFDistRobustly(x, df1, covariate = NULL, winsor.tail.p = c(0.05,0.1), trace = FALSE)

+fitFDistUnequalDF1(x, df1, covariate = NULL, robust = FALSE, prior.weights = NULL)

 }

 \arguments{

   \item{x}{numeric vector or array of positive values representing a sample from a scaled F-distribution.}

@@ -16,9 +19,11 @@ fitFDistRobustly(x, df1, covariate=NULL, winsor.tail.p=c(0.05,0.1), trace=FALSE)

   \item{covariate}{if non-\code{NULL}, the estimated scale value will depend on this numeric covariate.}

   \item{winsor.tail.p}{numeric vector of length 1 or 2, giving left and right tail proportions of \code{x} to Winsorize.}

   \item{trace}{logical value indicating whether a trace of the iteration progress should be printed.}

+  \item{robust}{logical. Should outlier values of \code{x} be down-weighted with results similar to \code{fitFDistRobustly}?}

+  \item{prior.weights}{numeric vector of (non-negative) prior weights.}

 }

 \details{

-\code{fitFDist} implements an algorithm proposed by Smyth (2004) and Phipson et al (2016).

+\code{fitFDist()} implements an algorithm proposed by Smyth (2004) and Phipson et al (2016).

 It estimates \code{scale} and \code{df2} under the assumption that \code{x} is distributed as \code{scale} times an F-distributed random variable on \code{df1} and \code{df2} degrees of freedom.

 The parameters are estimated using the method of moments, specifically from the mean and variance of the \code{x} values on the log-scale.

@@ -30,13 +35,27 @@ When \code{covariate} is supplied, a loess trend is estimated for the \code{x} v

 The robust method is described by Phipson et al (2016).

 As well as estimating the F-distribution for the bulk of the cases, i.e., with outliers discounted, \code{fitFDistRobustly} also returns an estimated F-distribution with reduced df2 that might be appropriate for each outlier case.

+

+\code{fitFDistUnequalDF1} was introduced in limma 3.61.8 and gives special attention to the possibility that the degrees of freedom \code{df1} might be unequal and might include non-integer values.

+The most important innovation of \code{fitFDistUnequalDF1} is downweighting of observations with lower degrees of freedom, to give more precise estimation overall.

+It also allows the possibility of prior weights, which can be used to downweight unreliable \code{x} values for reasons other than small \code{df1}.

+\code{fitFDistUnequalDF1} implements a different robust estimation strategy to \code{fitFDistRobustly}.

+Instead of Winsorizing the \code{x} values, potential outliers are instead downweighted using the prior weights.

+Whereas \code{fitFDist} and \code{fitFDistRobustly} use unweighted moment estimation for both \code{scale} and \code{df2},

+\code{fitFDistUnequalDF1} uses weighted moment estimation for \code{scale} and profile maximum likelihood for \code{df2}.

+

+\code{fitFDistUnequalDF1} gives improved performance over \code{fitFDist} and \code{fitFDistRobustly}, especially when the degrees of freedom are unequal but also to a lesser extent when the degrees of freedom are equal.

+Unequal residual degrees of freedom arise in limma pipelines when the expression matrix includes missing values or from the quasi-likelihood pipeline in edgeR v4 (Chen et al 2024).

+The edgeR v4 pipeline produces fractional degrees of freedom including, potentially, degrees of freedom less than 1.

 }

+

 \note{

 The algorithm used by \code{fitFDistRobustly} was revised slightly in limma 3.27.6.

 The \code{prob.outlier} value, which is the lower bound for \code{df2.shrunk}, may be slightly smaller than previously.

 }

+

 \value{

-\code{fitFDist} produces a list with the following components:

+\code{fitFDist} or \code{fitFDistUnequalDF1} with \code{robust=FALSE} produces a list with the following components:

   \item{scale}{scale factor for F-distribution. A vector if \code{covariate} is non-\code{NULL}, otherwise a scalar.}

   \item{df2}{the second degrees of freedom of the fitted F-distribution.}

@@ -47,19 +66,24 @@ The \code{prob.outlier} value, which is the lower bound for \code{df2.shrunk}, m

   \item{df2.shrunk}{numeric vector of values for the second degrees of freedom, with shrunk values for outliers. Most values are equal to \code{df2}, but outliers have reduced values depending on how extreme each case is. All values lie between \code{df2.outlier} and \code{df2}.}

 }

-\author{Gordon Smyth and Belinda Phipson}

+\author{Gordon Smyth, Belinda Phipson (\code{fitFDistRobustly}) and Lizhong Chen (\code{fitFDistUnequalDF1}).}

 \references{

-Smyth, G. K. (2004).

+Smyth GK (2004).

 Linear models and empirical Bayes methods for assessing differential expression in microarray experiments.

 \emph{Statistical Applications in Genetics and Molecular Biology} Volume 3, Issue 1, Article 3.

 \doi{10.2202/1544-6115.1027}

 \url{https://gksmyth.github.io/pubs/ebayes.pdf}

-Phipson, B, Lee, S, Majewski, IJ, Alexander, WS, and Smyth, GK (2016).

+Phipson B, Lee S, Majewski IJ, Alexander WS, Smyth GK (2016).

 Robust hyperparameter estimation protects against hypervariable genes and improves power to detect differential expression.

 \emph{Annals of Applied Statistics} 10, 946-963.

 \doi{10.1214/16-AOAS920}

+

+Chen Y, Chen L, Lun ATL, Baldoni PL, Smyth GK (2024).

+edgeR 4.0: powerful differential analysis of sequencing data with expanded functionality and improved support for small counts and larger datasets.

+\emph{bioRxiv} 2024.01.21.576131.

+\doi{10.1101/2024.01.21.576131}

 }

 \seealso{

@@ -5,7 +5,8 @@

 Normalize the columns of a matrix, cyclicly applying loess normalization to normalize each pair of columns to each other.

 }

 \usage{

-normalizeCyclicLoess(x, weights = NULL, span = 0.7, adaptive.span = FALSE, iterations = 3, method = "fast")

+normalizeCyclicLoess(x, weights = NULL, span = 0.7, adaptive.span = FALSE,

+                     iterations = 3, method = "fast")

 }

 \arguments{

   \item{x}{numeric matrix, or object which can be coerced to a numeric matrix, containing log-expression values.}

@@ -47,7 +48,7 @@ A matrix of the same dimensions as \code{x} containing the normalized values.

 \references{

 Bolstad BM, Irizarry RA, Astrand M, Speed TP (2003).

 A comparison of normalization methods for high density oligonucleotide array data based on bias and variance.

-\emph{Bioinformatics} \bold{19}, 185-193.

+\emph{Bioinformatics} 19, 185-193.

 Ballman KV, Grill DE, Oberg AL, Therneau TM (2004).

 Faster cyclic loess: normalizing RNA arrays via linear models.

@@ -7,15 +7,17 @@ Squeeze a set of sample variances together by computing empirical Bayes posterio

 }

 \usage{

-squeezeVar(var, df, covariate=NULL, robust=FALSE, winsor.tail.p=c(0.05,0.1))

+squeezeVar(var, df, covariate = NULL, robust = FALSE, winsor.tail.p = c(0.05,0.1),

+           legacy = NULL)

 }

 \arguments{

 \item{var}{numeric vector of independent sample variances.}

-\item{df}{numeric vector of degrees of freedom for the sample variances.}

-\item{covariate}{if non-\code{NULL}, \code{var.prior} will depend on this numeric covariate. Otherwise, \code{var.prior} is constant.}

+\item{df}{numeric vector of degrees of freedom for the sample variances. Can be a unit vector or of same length as \code{var}.}

+\item{covariate}{numeric covariate of same length as \code{var} for estimating a trended prior variance. If \code{NULL}, then the prior variance \code{var.prior} is constant.}

 \item{robust}{logical, should the estimation of \code{df.prior} and \code{var.prior} be robustified against outlier sample variances?}

-\item{winsor.tail.p}{numeric vector of length 1 or 2, giving left and right tail proportions of \code{x} to Winsorize. Used only when \code{robust=TRUE}.}

+\item{winsor.tail.p}{numeric vector of length 1 or 2, giving left and right tail proportions of \code{x} to Winsorize when \code{robust=TRUE}.}

+\item{legacy}{logical. If \code{FALSE} then the new function \code{fitFDistUnequalDF1} will be called internally, otherwise the legacy functions \code{fitFDist} or \code{fitFDistRobustly} will be used. If \code{NULL}, then \code{fitFDistUnequalDF1} will be used whenever the degrees of freedom \code{df} are not all equal.}

 }

 \details{

@@ -32,10 +34,17 @@ The scale and degrees of freedom of this prior distribution are estimated from t

 The effect of this function is to squeeze the variances towards a common value, or to a global trend if a \code{covariate} is provided.

 The squeezed variances have a smaller expected mean square error to the true variances than do the sample variances themselves.

+The amount of squeezing is controlled by the \code{prior.df}.

+Both the global trend and the prior df are estimated internally but fitting an F-distribution to the sample variances, using either \code{fitFDist()} or \code{fitFDistRobustly()} or \code{fitFDistUnequalDF1()}.

+

 If \code{covariate} is non-null, then the scale parameter of the prior distribution is assumed to depend on the covariate.

 If the covariate is average log-expression, then the effect is an intensity-dependent trend similar to that in Sartor et al (2006).

-\code{robust=TRUE} implements the robust empirical Bayes procedure of Phipson et al (2016) which allows some of the \code{var} values to be outliers.

+\code{robust=TRUE} implements the robust empirical Bayes procedure of Phipson et al (2016), which allows some of the \code{var} values to be outliers.

+

+The \code{legacy} argument was added in limma version 3.61.8 (August 2024).

+If \code{legacy=FALSE}, then the new function \code{fitFDistUnequalDF1()} provides improved estimation of the global trend and prior df hyperparameters, especially when the \code{df} values are unequal.

+\code{legacy=TRUE} provides legacy behavior for backward compatibility.

 }

 \note{

@@ -74,7 +83,7 @@ incorporating corrections to 30 June 2009.

 \seealso{

 This function is called by \code{\link{eBayes}}.

-This function calls \code{\link{fitFDist}}.

+This function calls \code{\link{fitFDist}}, \code{\link{fitFDistRobustly}} or \code{\link{fitFDistUnequalDF1}} .

 An overview of linear model functions in limma is given by \link{06.LinearModels}.

 }

@@ -6,7 +6,8 @@ Combine voom observational-level precision weights with sample-specific quality

 }

 \usage{

 voomWithQualityWeights(counts, design = NULL, lib.size = NULL, normalize.method = "none",

-             plot = FALSE, span = 0.5, adaptive.span = FALSE, var.design = NULL, var.group = NULL,

+             plot = FALSE, span = 0.5, adaptive.span = FALSE,

+             var.design = NULL, var.group = NULL,

              method = "genebygene", maxiter = 50, tol = 1e-5, trace = FALSE,

              col = NULL, \dots)

 }

@@ -12,8 +12,8 @@ Equivalent to calling vooma(), arrayWeights(), duplicateCorrelation() and lmFit(

 \usage{

 voomaLmFit(y, design = NULL, prior.weights = NULL, block = NULL,

            sample.weights = FALSE, var.design = NULL, var.group = NULL, prior.n = 10,

-           predictor = NULL, span = NULL, legacy.span = FALSE, plot = FALSE, save.plot = FALSE,

-           keep.EList = TRUE)

+           predictor = NULL, span = NULL, legacy.span = FALSE,

+           plot = FALSE, save.plot = FALSE, keep.EList = TRUE)

 }

 \arguments{

@@ -1,5 +1,5 @@

-R version 4.4.0 RC (2024-04-16 r86458 ucrt) -- "Puppy Cup"

+R version 4.4.1 (2024-06-14 ucrt) -- "Race for Your Life"

 Copyright (C) 2024 The R Foundation for Statistical Computing

 Platform: x86_64-w64-mingw32/x64

@@ -917,52 +917,52 @@ alpha          1        1       0

 beta          -1        0       1

 $df.prior

-[1] 9.306153

+[1] 8134.845

 $s2.prior

-[1] 0.923179

+[1] 1.021387

 $var.prior

-[1] 17.33142 17.33142 12.26855

+[1] 15.664973 11.397823  9.122785

 $proportion

 [1] 0.01

 $s2.post

-[1] 1.2677996 0.6459499 1.1251558 0.9097727 0.9124980

+[1] 1.021963 1.020793 1.021694 1.021289 1.021321

 $t

      alpha-beta mu+alpha   mu+beta

-[1,]   3.847656 2.580411 -2.860996

-[2,]   6.637308 5.113018 -4.273553

-[3,]   3.692066 1.720376 -3.500994

-[4,]   3.464003 2.926234 -1.972606

-[5,]   3.175181 2.566881 -1.814221

+[1,]   4.285525 2.874066 -3.186582

+[2,]   5.279861 4.067315 -3.399536

+[3,]   3.874497 1.805382 -3.673984

+[4,]   3.269417 2.761856 -1.861797

+[5,]   3.001258 2.426278 -1.714845

 $df.total

-[1] 15.30615 15.30615 15.30615 15.30615 13.30615

+[1] 28 28 28 28 28

 $p.value

        alpha-beta     mu+alpha      mu+beta

-[1,] 1.529450e-03 0.0206493481 0.0117123495

-[2,] 7.144893e-06 0.0001195844 0.0006385076

-[3,] 2.109270e-03 0.1055117477 0.0031325769

-[4,] 3.381970e-03 0.0102514264 0.0668844448

-[5,] 7.124839e-03 0.0230888584 0.0922478630

+[1,] 1.945874e-04 0.0076518793 0.0035226472

+[2,] 1.290875e-05 0.0003507232 0.0020449743

+[3,] 5.877177e-04 0.0817788496 0.0009997966

+[4,] 2.854794e-03 0.0100337538 0.0731588447

+[5,] 5.599946e-03 0.0219470179 0.0974229293

 $lods

-     alpha-beta  mu+alpha    mu+beta

-[1,]  -1.013417 -3.702133 -3.0332393

-[2,]   3.981496  1.283349 -0.2615911

-[3,]  -1.315036 -5.168621 -1.7864101

-[4,]  -1.757103 -3.043209 -4.6191869

-[5,]  -2.257358 -3.478267 -4.5683738

+     alpha-beta   mu+alpha    mu+beta

+[1,]  0.7356274 -2.7480831 -1.9651516

+[2,]  3.2466921  0.1146964 -1.4669585

+[3,] -0.2839280 -4.8267525 -0.8071472

+[4,] -1.7300021 -2.9939736 -4.6385674

+[5,] -2.1848828 -3.4282344 -4.5756693

 $F

-[1]  7.421911 22.203107  7.608327  6.227010  5.060579

+[1]  9.207280 14.049948  8.378781  5.547069  4.521367

 $F.p.value

-[1] 5.581800e-03 2.988923e-05 5.080726e-03 1.050148e-02 2.320274e-02

+[1] 1.013549e-04 8.103854e-07 2.316764e-04 3.913618e-03 1.090148e-02

 > 

 > ### uniquegenelist

@@ -1273,4 +1273,4 @@ GO:0043009 chordate embryonic developm...  BP  3  0    2 1.0000000 0.025788497

 > 

 > proc.time()

    user  system elapsed 

-   3.06    0.31    3.36 

+   1.32    0.15    2.62