From 51738d9807da32447038ec0253c2413addfbc448 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Mon, 26 May 2025 14:51:30 -0700 Subject: [PATCH 01/13] msqrobLm(): rename data arg to colsMetadata as that's what it is --- R/msqrob.R | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/R/msqrob.R b/R/msqrob.R index 2247aa7..28c561f 100644 --- a/R/msqrob.R +++ b/R/msqrob.R @@ -57,9 +57,10 @@ msqrobLm <- function(y, robust = TRUE, maxitRob = 5) { myDesign <- model.matrix(formula, data) + # apply the model to each protein models <- apply(y, 1, function(y, design) { - ## computatability check + ## computability check obs <- is.finite(y) type <- "fitError" model <- list( From f6c0cf2cb28e378def303e407e6d6de00c5b6f97 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Mon, 26 May 2025 14:53:22 -0700 Subject: [PATCH 02/13] .matchQuantColsOrder() match colsMetadata order to the data columns --- R/msqrob-utils.R | 24 ++++++++++++++++++++++++ R/msqrob.R | 1 + 2 files changed, 25 insertions(+) diff --git a/R/msqrob-utils.R b/R/msqrob-utils.R index 9449e32..cdd119d 100644 --- a/R/msqrob-utils.R +++ b/R/msqrob-utils.R @@ -668,3 +668,27 @@ makeContrast <- function(contrasts, parameterNames) { ex, " to character" ) } + +# Match Quantitative Columns Order for msqrobLmxxx() methods +.matchQuantColsOrder <- function(colsMetadata, data) { + if (!"quantCols" %in% names(colsMetadata)) { + warning("colsMetadata does not contain a 'quantCols' column. Skipping reordering.") + return(colsMetadata) + } + + data_cols <- colnames(data) + quant_cols_values <- colsMetadata$quantCols + + missing_cols <- setdiff(data_cols, quant_cols_values) + if (length(missing_cols) > 0) { + stop( + "Column(s) '", paste(missing_cols, collapse = "', '"), + "' from data not found in colsMetadata$quantCols." + ) + } + + match_indices <- match(data_cols, quant_cols_values) + reordered_metadata <- colsMetadata[match_indices, , drop = FALSE] + + return(reordered_metadata) +} diff --git a/R/msqrob.R b/R/msqrob.R index 28c561f..9834e86 100644 --- a/R/msqrob.R +++ b/R/msqrob.R @@ -56,6 +56,7 @@ msqrobLm <- function(y, data, robust = TRUE, maxitRob = 5) { + data <- .matchQuantColsOrder(data, y) myDesign <- model.matrix(formula, data) # apply the model to each protein models <- apply(y, 1, From 3df378505c1dd811e017f9916fd56a3e605b5ff8 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Tue, 27 May 2025 14:36:40 -0700 Subject: [PATCH 03/13] use [[ ]] for getting single element --- R/msqrob.R | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/R/msqrob.R b/R/msqrob.R index 9834e86..92adc16 100644 --- a/R/msqrob.R +++ b/R/msqrob.R @@ -147,7 +147,7 @@ msqrobLm <- function(y, ## Put variance and degrees of freedom in appropriate slots for (i in seq_len(length(models))) { mydf <- hlp$df.prior + getDF(models[[i]]) - models[[i]]@varPosterior <- as.numeric(hlp$var.post[i]) + models[[i]]@varPosterior <- as.numeric(hlp$var.post[[i]]) models[[i]]@dfPosterior <- as.numeric(mydf) } @@ -335,7 +335,7 @@ msqrobLmer <- function(y, ) for (i in seq_len(length(models))) { - models[[i]]@varPosterior <- as.numeric(hlp$var.post[i]) + models[[i]]@varPosterior <- as.numeric(hlp$var.post[[i]]) models[[i]]@dfPosterior <- as.numeric(hlp$df.prior + getDF(models[[i]])) } return(models) From 1c55018e33205c57ed19baa8b4c5d3faad702cda Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Tue, 27 May 2025 14:37:30 -0700 Subject: [PATCH 04/13] fix whitespace --- R/msqrob.R | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/R/msqrob.R b/R/msqrob.R index 92adc16..78dd7ae 100644 --- a/R/msqrob.R +++ b/R/msqrob.R @@ -342,7 +342,8 @@ msqrobLmer <- function(y, } ## Fit the mixed models with ridge regression -.ridge_msqrobLmer <- function(y,rowdata=NULL,formula,coldata, doQR, robust,maxitRob=1,tol = 1e-06){ +.ridge_msqrobLmer <- function(y, rowdata=NULL, formula, coldata, doQR, + robust, maxitRob=1, tol = 1e-06){ #Create the matrix containing the variable information data <- .create_data(y,rowdata,coldata) @@ -470,7 +471,8 @@ msqrobLmer <- function(y, } ## Fit the mixed models without ridge regression -.noridge_msqrobLmer <- function(y,rowdata=NULL,formula,coldata, robust,maxitRob=0, tol = 1e-06 ){ +.noridge_msqrobLmer <- function(y, rowdata=NULL, formula, coldata, + robust, maxitRob=0, tol = 1e-06){ #Create the matrix containing the variable information data <- .create_data(y,rowdata,coldata) From b68d4bff34edf6e8feb6cb9e6eca9166714b60d5 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Tue, 27 May 2025 14:40:52 -0700 Subject: [PATCH 05/13] msqrobLmer(): reduce code duplication --- R/msqrob.R | 59 ++++++++++++++++-------------------------------------- 1 file changed, 17 insertions(+), 42 deletions(-) diff --git a/R/msqrob.R b/R/msqrob.R index 78dd7ae..241d5a4 100644 --- a/R/msqrob.R +++ b/R/msqrob.R @@ -285,49 +285,24 @@ msqrobLmer <- function(y, y <- split.data.frame(y, featureGroups) - if (ridge == TRUE){ - if(is.null(rowdata)){ - models <- bplapply(y, - FUN = .ridge_msqrobLmer, - "formula" = formula, - "coldata" = data, - "doQR" = doQR, - "robust"=robust, - "maxitRob" = maxitRob, - "tol" =tol) - } else{ - models <- bpmapply(FUN = .ridge_msqrobLmer, - y, rowdata, - MoreArgs = list("formula" = formula, - "coldata" = data, - "doQR" = doQR, - "robust"=robust, - "maxitRob" = maxitRob, - "tol" =tol)) - } - - }else{ - - if(is.null(rowdata)){ - models <- bplapply(y, - FUN = .noridge_msqrobLmer, - "formula" = formula, - "coldata" = data, - "robust"=robust, - "maxitRob" = maxitRob, - "tol" =tol) - } else{ - models <- bpmapply(FUN = .noridge_msqrobLmer, - y, rowdata, - MoreArgs = list("formula" = formula, - "coldata" = data, - "robust"=robust, - "maxitRob" = maxitRob, - "tol" =tol)) - } + fit_args <- list( + formula = formula, + coldata = data, + robust = robust, + maxitRob = maxitRob, + tol = tol + ) + if (ridge) { + fit_func <- .ridge_msqrobLmer + fit_args$doQR <- doQR + } else { + fit_func <- .noridge_msqrobLmer + } + models <- if(is.null(rowdata)){ + do.call(bplapply, append(list(y, FUN=fit_func), fit_args)) + } else { + bpmapply(FUN = fit_func, y, rowdata, MoreArgs = fit_args) } - - hlp <- limma::squeezeVar( var = vapply(models, getVar, numeric(1)), From 164622852bc15f866bd1465874a712d2225146ac Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Tue, 27 May 2025 14:43:45 -0700 Subject: [PATCH 06/13] msqrobLmer(): skip NA when calling squeezeVar --- R/msqrob.R | 11 ++++++++--- 1 file changed, 8 insertions(+), 3 deletions(-) diff --git a/R/msqrob.R b/R/msqrob.R index 241d5a4..47573e2 100644 --- a/R/msqrob.R +++ b/R/msqrob.R @@ -304,13 +304,18 @@ msqrobLmer <- function(y, bpmapply(FUN = fit_func, y, rowdata, MoreArgs = fit_args) } + modelVars <- vapply(models, getVar, numeric(1)) + modelDfs <- vapply(models, getDF, numeric(1)) + modelMask <- !is.na(modelVars) & !is.na(modelDfs) + varPosterior <- rep_along(models, NA_real_) hlp <- limma::squeezeVar( - var = vapply(models, getVar, numeric(1)), - df = vapply(models, getDF, numeric(1)) + var = modelVars[modelMask], + df = modelDfs[modelMask], ) + varPosterior[modelMask] <- hlp$var.post for (i in seq_len(length(models))) { - models[[i]]@varPosterior <- as.numeric(hlp$var.post[[i]]) + models[[i]]@varPosterior <- as.numeric(varPosterior[[i]]) models[[i]]@dfPosterior <- as.numeric(hlp$df.prior + getDF(models[[i]])) } return(models) From 9aba3de8669f48f94ae2da6c8f4bae02aa2636c5 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Tue, 27 May 2025 14:44:53 -0700 Subject: [PATCH 07/13] msqrobLmer(): keep.model argument --- R/msqrob.R | 45 +++++++++++++++++++++++++++------------------ 1 file changed, 27 insertions(+), 18 deletions(-) diff --git a/R/msqrob.R b/R/msqrob.R index 47573e2..29c6e56 100644 --- a/R/msqrob.R +++ b/R/msqrob.R @@ -267,7 +267,8 @@ msqrobLmer <- function(y, maxitRob = 1, doQR = TRUE, featureGroups=NULL, - lmerArgs = list(control = lmerControl(calc.derivs = FALSE))){ + lmerArgs = list(control = lmerControl(calc.derivs = FALSE)), + keep.model = FALSE) { #Get the featureGroups variable if (is.null(featureGroups)){ @@ -290,7 +291,8 @@ msqrobLmer <- function(y, coldata = data, robust = robust, maxitRob = maxitRob, - tol = tol + tol = tol, + keep.model = keep.model ) if (ridge) { fit_func <- .ridge_msqrobLmer @@ -323,7 +325,8 @@ msqrobLmer <- function(y, ## Fit the mixed models with ridge regression .ridge_msqrobLmer <- function(y, rowdata=NULL, formula, coldata, doQR, - robust, maxitRob=1, tol = 1e-06){ + robust, maxitRob=1, tol = 1e-06, + keep.model = FALSE){ #Create the matrix containing the variable information data <- .create_data(y,rowdata,coldata) @@ -441,7 +444,7 @@ msqrobLmer <- function(y, } }, silent = TRUE) - model <- .create_model(betas, vcovUnscaled, sigma, df.residual, w, model) + model <- .create_model(betas, vcovUnscaled, sigma, df.residual, w, model, keep.model = keep.model) } return(StatModel(type = type, @@ -452,7 +455,8 @@ msqrobLmer <- function(y, ## Fit the mixed models without ridge regression .noridge_msqrobLmer <- function(y, rowdata=NULL, formula, coldata, - robust, maxitRob=0, tol = 1e-06){ + robust, maxitRob=0, tol = 1e-06, + keep.model = FALSE){ #Create the matrix containing the variable information data <- .create_data(y,rowdata,coldata) @@ -498,7 +502,8 @@ msqrobLmer <- function(y, } }, silent = TRUE) - model <- .create_model(betas, vcovUnscaled, sigma, df.residual, w, model) + model <- .create_model(betas, vcovUnscaled, sigma, df.residual, w, + model, keep.model = keep.model) } return(StatModel(type = type, @@ -574,6 +579,7 @@ msqrobLmer <- function(y, betaB } +#' Calculate the weighted REML residual degrees of freedom #' @importFrom stats resid .getDfLmer <- function(object) { w <- object@frame$"(weights)" @@ -595,21 +601,24 @@ msqrobLmer <- function(y, return(model) } -.create_model <- function(betas, vcovUnscaled, sigma, df.residual, w, model){ +.create_model <- function(betas, vcovUnscaled, sigma, df.residual, w, model, keep.model = FALSE){ if (df.residual<2L){ - model <- list(coefficients = NA, - vcovUnscaled = NA, - sigma = NA, - df.residual = NA, - w = NA) + res <- list(coefficients = NA, + vcovUnscaled = NA, + sigma = NA, + df.residual = NA, + w = NA) } else { - model <- list(coefficients = betas, - vcovUnscaled = vcovUnscaled, - sigma = sigma, - df.residual = df.residual, - w = model@frame$`(weights)`) + res <- list(coefficients = betas, + vcovUnscaled = vcovUnscaled, + sigma = sigma, + df.residual = df.residual, + w = model@frame$`(weights)`) } - return(model) + if (keep.model) { + res$model <- model + } + return(res) } From 0b9e870bcf2b9f27a9e9ddd3b3d6247eb2c5449c Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Tue, 24 Jun 2025 14:20:47 -0700 Subject: [PATCH 08/13] fix typo --- R/msqrob.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/msqrob.R b/R/msqrob.R index 29c6e56..3075c4a 100644 --- a/R/msqrob.R +++ b/R/msqrob.R @@ -352,7 +352,7 @@ msqrobLmer <- function(y, formula <- formula(y ~ (1|ridge)) } else { if (nobars(formula)[[2]] != ~1){ - #udpate formula to remove any fixed effect variables and replace with ridge + # update formula to remove any fixed effect variables and replace with ridge formula <- formula( paste0("y ~ (1|ridge) + ", paste0("(",paste(findbars(formula), collapse=")+("),")"))) } else { From 5faa82cb94641fad1c3b8717c8a73217a782c1d2 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Tue, 24 Jun 2025 14:21:03 -0700 Subject: [PATCH 09/13] use NA_real_ --- R/msqrob.R | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/R/msqrob.R b/R/msqrob.R index 3075c4a..9d8fa9e 100644 --- a/R/msqrob.R +++ b/R/msqrob.R @@ -449,8 +449,8 @@ msqrobLmer <- function(y, return(StatModel(type = type, params = model, - varPosterior = as.numeric(NA), - dfPosterior = as.numeric(NA))) + varPosterior = NA_real_, + dfPosterior = NA_real_)) } ## Fit the mixed models without ridge regression @@ -508,8 +508,8 @@ msqrobLmer <- function(y, return(StatModel(type = type, params = model, - varPosterior = as.numeric(NA), - dfPosterior = as.numeric(NA))) + varPosterior = NA_real_, + dfPosterior = NA_real_)) } From 2343a953e1d5da61393b0d5b9ef9baf351fa75a2 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Tue, 27 May 2025 14:45:35 -0700 Subject: [PATCH 10/13] topFeatures(): support alt. methods for DDoF calc add ddf.method parameter to topFeatures() to support using dof_xxx() methods from the parameters package for DDoF calculation instead of using the residual DoF. --- DESCRIPTION | 1 + R/topFeatures.R | 59 +++++++++++++++++++++++++++++++++++++++++++++---- 2 files changed, 56 insertions(+), 4 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index ad098aa..8babaff 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -60,6 +60,7 @@ Imports: codetools Suggests: multcomp, + parameters, gridExtra, knitr, BiocStyle, diff --git a/R/topFeatures.R b/R/topFeatures.R index e4a6c51..934849d 100644 --- a/R/topFeatures.R +++ b/R/topFeatures.R @@ -13,7 +13,23 @@ #' Options, in increasing conservatism, include ‘"none"’, #' ‘"BH"’, ‘"BY"’ and ‘"holm"’. See ‘p.adjust’ for the complete #' list of options. Default is "BH" the Benjamini-Hochberg method -#' to controle the False Discovery Rate (FDR). +#' to control the *False Discovery Rate* (FDR). +#' @param ddf.method `character` specifying the method to calculate the +#' *denominator degrees of freedom* (DoF) for the t-test. +#' Options are: +#' - "residual": (default) use the posterior residual DoF (\code{getDfPosterior}), +#' which was the only DDoF method for *msqrob2* before v2.XX. +#' It can significantly overestimate the denominator DoF, +#' resulting in overly significant p-values. +#' - "ML1": use the `dof_ml1` function from the *parameters* package for +#' heuristic approximation of the DDoF, which is much faster than +#' the Kenward-Roger method, while still providing an adequate approximation. +#' - "KenwardRoger": use the `dof_kenward` function from the *parameters* package for +#' *Kenward-Roger approximation* of the DDoF. This is the most accurate, but +#' also most computationally intensive method. +#' - "Satterthwaite": use the `dof_satterthwaite` function from the *parameters* package for +#' *Satterthwaite approximation* of the DDoF. This is a faster, but less accurate +#' alternative to *Kenward-Roger* method, which may not work well for small sample sizes. #' @param sort `boolean(1)` to indicate if the features have to be sorted according #' to statistical significance. #' @param alpha `numeric` specifying the cutoff value for adjusted p-values. @@ -42,8 +58,14 @@ #' @rdname topTable #' @author Lieven Clement #' @export - -topFeatures <- function(models, contrast, adjust.method = "BH", sort = TRUE, alpha = 1) { +topFeatures <- function( + models, + contrast, + adjust.method = "BH", + ddf.method = c("residual", "ML1", "KenwardRoger", "Satterthwaite"), + sort = TRUE, + alpha = 1 +) { if (!is(contrast, "matrix")) { contrast <- as.matrix(contrast) } @@ -63,11 +85,40 @@ topFeatures <- function(models, contrast, adjust.method = "BH", sort = TRUE, alp numeric(1), L = contrast )) - df <- vapply(models, getDfPosterior, numeric(1)) + ddf.method <- match.arg(ddf.method) + if (ddf.method == "residual") { + df <- vapply(models, getDfPosterior, numeric(1)) + } else if (ddf.method %in% c("KenwardRoger", "ML1", "Satterthwaite")) { + if (!requireNamespace("parameters", quietly = TRUE)) { + stop("parameters package is required to calculate the denominator DoF using ", ddf.method, " method.") + } + ddf_func <- if (ddf.method == "KenwardRoger") { + parameters::dof_kenward + } else if (ddf.method == "ML1") { + parameters::dof_ml1 + } else if (ddf.method == "Satterthwaite") { + parameters::dof_satterthwaite + } else { + stop("Unsupported ddf.method=", ddf.method) + } + df <- bplapply(models, function(model) { + if ("model" %in% names(model@params) && is(model@params$model, "lmerMod")) { + lmm <- model@params$model + tryCatch(min(ddf_func(lmm)[rownames(contrast)]), + error = function(e) NA_real_) + } else { + NA_real_ + } + }) + df <- unlist(df, recursive = FALSE, use.names = FALSE) + } else { + stop("Unsupported ddf.method=", ddf.method) + } t <- logFC / se pval <- pt(-abs(t), df) * 2 adjPval <- p.adjust(pval, method = adjust.method) out <- data.frame(logFC, se, df, t, pval, adjPval) + rownames(out) <- names(models) if (alpha < 1) { signif <- adjPval < alpha out <- na.exclude(out[signif, ]) From 8411f296d6c7c64c9ca2933bbc0d3effb604f61a Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Tue, 24 Jun 2025 16:13:48 -0700 Subject: [PATCH 11/13] .create.model(): remove unused w arg Calling .create.model() with undefined w causes CI failures. For weights, .create.model() uses model@frame$`(weights)`. --- R/msqrob.R | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/R/msqrob.R b/R/msqrob.R index 9d8fa9e..bd4b6d6 100644 --- a/R/msqrob.R +++ b/R/msqrob.R @@ -444,7 +444,8 @@ msqrobLmer <- function(y, } }, silent = TRUE) - model <- .create_model(betas, vcovUnscaled, sigma, df.residual, w, model, keep.model = keep.model) + model <- .create_model(betas, vcovUnscaled, sigma, df.residual, + model, keep.model = keep.model) } return(StatModel(type = type, @@ -502,7 +503,7 @@ msqrobLmer <- function(y, } }, silent = TRUE) - model <- .create_model(betas, vcovUnscaled, sigma, df.residual, w, + model <- .create_model(betas, vcovUnscaled, sigma, df.residual, model, keep.model = keep.model) } @@ -601,7 +602,7 @@ msqrobLmer <- function(y, return(model) } -.create_model <- function(betas, vcovUnscaled, sigma, df.residual, w, model, keep.model = FALSE){ +.create_model <- function(betas, vcovUnscaled, sigma, df.residual, model, keep.model = FALSE){ if (df.residual<2L){ res <- list(coefficients = NA, vcovUnscaled = NA, From 4b75f986b6a67e98f94248e4e6fc426ea58197a9 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Tue, 24 Jun 2025 16:25:31 -0700 Subject: [PATCH 12/13] msqrobLmer(): document keep.model --- R/msqrob.R | 3 +++ man/msqrobLmer.Rd | 6 +++++- 2 files changed, 8 insertions(+), 1 deletion(-) diff --git a/R/msqrob.R b/R/msqrob.R index bd4b6d6..a15f294 100644 --- a/R/msqrob.R +++ b/R/msqrob.R @@ -211,6 +211,9 @@ msqrobLm <- function(y, #' ‘lmerControl’ documentation of the lme4 package for more details. #' Default is `list(control = lmerControl(calc.derivs = FALSE))` #' +#' @param keep.model `boolean(1)` to specify if the `lme4` fitted model should be kept in the +#' `$model` field of the output. Default is `FALSE`. +#' #' @examples #' #' # Load example data diff --git a/man/msqrobLmer.Rd b/man/msqrobLmer.Rd index 7d42f81..cb6d6bd 100644 --- a/man/msqrobLmer.Rd +++ b/man/msqrobLmer.Rd @@ -15,7 +15,8 @@ msqrobLmer( maxitRob = 1, doQR = TRUE, featureGroups = NULL, - lmerArgs = list(control = lmerControl(calc.derivs = FALSE)) + lmerArgs = list(control = lmerControl(calc.derivs = FALSE)), + keep.model = FALSE ) } \arguments{ @@ -61,6 +62,9 @@ containing control parameters, including the nonlinear optimizer to be used and parameters to be passed through to the nonlinear optimizer, see the ‘lmerControl’ documentation of the lme4 package for more details. Default is \code{list(control = lmerControl(calc.derivs = FALSE))}} + +\item{keep.model}{\code{boolean(1)} to specify if the \code{lme4} fitted model should be kept in the +\verb{$model} field of the output. Default is \code{FALSE}.} } \value{ A list of objects of the \code{StatModel} class. From a1506a0152121252b65fe351c82e81081979b467 Mon Sep 17 00:00:00 2001 From: Alexey Stukalov Date: Tue, 24 Jun 2025 16:26:29 -0700 Subject: [PATCH 13/13] topFeatures(): regenerate docs --- man/topTable.Rd | 28 ++++++++++++++++++++++++++-- 1 file changed, 26 insertions(+), 2 deletions(-) diff --git a/man/topTable.Rd b/man/topTable.Rd index 8d3cc76..77c44b8 100644 --- a/man/topTable.Rd +++ b/man/topTable.Rd @@ -4,7 +4,14 @@ \alias{topFeatures} \title{Toplist of DE proteins, peptides or features} \usage{ -topFeatures(models, contrast, adjust.method = "BH", sort = TRUE, alpha = 1) +topFeatures( + models, + contrast, + adjust.method = "BH", + ddf.method = c("residual", "ML1", "KenwardRoger", "Satterthwaite"), + sort = TRUE, + alpha = 1 +) } \arguments{ \item{models}{A list with elements of the class \code{StatModel} that are @@ -20,7 +27,24 @@ the p-values for multiple testing. Options, in increasing conservatism, include ‘"none"’, ‘"BH"’, ‘"BY"’ and ‘"holm"’. See ‘p.adjust’ for the complete list of options. Default is "BH" the Benjamini-Hochberg method -to controle the False Discovery Rate (FDR).} +to control the \emph{False Discovery Rate} (FDR).} + +\item{ddf.method}{\code{character} specifying the method to calculate the +\emph{denominator degrees of freedom} (DoF) for the t-test. +Options are: +- "residual": (default) use the posterior residual DoF (\code{getDfPosterior}), +which was the only DDoF method for \emph{msqrob2} before v2.XX. +It can significantly overestimate the denominator DoF, +resulting in overly significant p-values. +- "ML1": use the \code{dof_ml1} function from the \emph{parameters} package for +heuristic approximation of the DDoF, which is much faster than +the Kenward-Roger method, while still providing an adequate approximation. +- "KenwardRoger": use the \code{dof_kenward} function from the \emph{parameters} package for +\emph{Kenward-Roger approximation} of the DDoF. This is the most accurate, but +also most computationally intensive method. +- "Satterthwaite": use the \code{dof_satterthwaite} function from the \emph{parameters} package for +\emph{Satterthwaite approximation} of the DDoF. This is a faster, but less accurate +alternative to \emph{Kenward-Roger} method, which may not work well for small sample sizes.} \item{sort}{\code{boolean(1)} to indicate if the features have to be sorted according to statistical significance.}