diff --git a/DESCRIPTION b/DESCRIPTION index d890bee3..818a948b 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -40,7 +40,8 @@ Imports: parallel, rlang, matter, - BiocParallel + BiocParallel, + RhpcBLASctl Suggests: BiocStyle, knitr, @@ -62,4 +63,4 @@ Packaged: 2017-10-20 02:13:12 UTC; meenachoi LinkingTo: Rcpp, RcppArmadillo -Config/roxygen2/version: 8.0.0 +Config/roxygen2/version: 8.1.0 diff --git a/NAMESPACE b/NAMESPACE index ba6c1354..5f156811 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -55,95 +55,124 @@ import(data.table) import(ggplot2) import(limma) import(lme4) -importFrom(BiocParallel,bpisup) -importFrom(BiocParallel,bplapply) -importFrom(BiocParallel,bpnworkers) -importFrom(BiocParallel,bpprogressbar) -importFrom(BiocParallel,bpstart) -importFrom(BiocParallel,bpstop) +importFrom(BiocParallel, + bpisup, + bplapply, + bpnworkers, + bpprogressbar, + bpstart, + bpstop +) importFrom(MASS,rlm) -importFrom(MSstatsConvert,DIANNtoMSstatsFormat) -importFrom(MSstatsConvert,DIAUmpiretoMSstatsFormat) -importFrom(MSstatsConvert,FragPipetoMSstatsFormat) -importFrom(MSstatsConvert,MSstatsBalancedDesign) -importFrom(MSstatsConvert,MSstatsClean) -importFrom(MSstatsConvert,MSstatsImport) -importFrom(MSstatsConvert,MSstatsLogsSettings) -importFrom(MSstatsConvert,MSstatsMakeAnnotation) -importFrom(MSstatsConvert,MSstatsPreprocess) -importFrom(MSstatsConvert,MZMinetoMSstatsFormat) -importFrom(MSstatsConvert,MaxQtoMSstatsFormat) -importFrom(MSstatsConvert,OpenMStoMSstatsFormat) -importFrom(MSstatsConvert,OpenSWATHtoMSstatsFormat) -importFrom(MSstatsConvert,PDtoMSstatsFormat) -importFrom(MSstatsConvert,ProgenesistoMSstatsFormat) -importFrom(MSstatsConvert,SkylinetoMSstatsFormat) -importFrom(MSstatsConvert,SpectronauttoMSstatsFormat) +importFrom(MSstatsConvert, + DIANNtoMSstatsFormat, + DIAUmpiretoMSstatsFormat, + FragPipetoMSstatsFormat, + MSstatsBalancedDesign, + MSstatsClean, + MSstatsImport, + MSstatsLogsSettings, + MSstatsMakeAnnotation, + MSstatsPreprocess, + MZMinetoMSstatsFormat, + MaxQtoMSstatsFormat, + OpenMStoMSstatsFormat, + OpenSWATHtoMSstatsFormat, + PDtoMSstatsFormat, + ProgenesistoMSstatsFormat, + SkylinetoMSstatsFormat, + SpectronauttoMSstatsFormat +) importFrom(Rcpp,sourceCpp) -importFrom(data.table,as.data.table) -importFrom(data.table,data.table) -importFrom(data.table,fifelse) -importFrom(data.table,melt) -importFrom(data.table,rbindlist) -importFrom(data.table,setDT) -importFrom(data.table,setDTthreads) -importFrom(data.table,uniqueN) +importFrom(RhpcBLASctl,blas_set_num_threads) +importFrom(data.table, + as.data.table, + data.table, + fifelse, + melt, + rbindlist, + setDT, + setDTthreads, + uniqueN +) importFrom(ggrepel,geom_text_repel) importFrom(gplots,heatmap.2) -importFrom(grDevices,dev.off) -importFrom(grDevices,hcl) -importFrom(grDevices,pdf) -importFrom(graphics,axis) -importFrom(graphics,image) -importFrom(graphics,legend) -importFrom(graphics,mtext) -importFrom(graphics,par) -importFrom(graphics,plot) -importFrom(graphics,plot.new) -importFrom(graphics,title) -importFrom(htmltools,div) -importFrom(htmltools,save_html) -importFrom(htmltools,tagList) +importFrom(grDevices, + dev.off, + hcl, + pdf +) +importFrom(graphics, + axis, + image, + legend, + mtext, + par, + plot, + plot.new, + title +) +importFrom(htmltools, + div, + save_html, + tagList +) importFrom(limma,squeezeVar) importFrom(lme4,lmer) importFrom(marray,maPalette) importFrom(matter,SnowfastParam) importFrom(methods,is) -importFrom(parallel,clusterExport) -importFrom(parallel,makeCluster) -importFrom(parallel,parLapply) -importFrom(parallel,stopCluster) -importFrom(plotly,add_trace) -importFrom(plotly,ggplotly) -importFrom(plotly,layout) -importFrom(plotly,plot_ly) -importFrom(plotly,style) -importFrom(plotly,subplot) +importFrom(parallel, + clusterExport, + makeCluster, + parLapply, + stopCluster +) +importFrom(plotly, + add_trace, + ggplotly, + layout, + plot_ly, + style, + subplot +) importFrom(preprocessCore,normalize.quantiles) importFrom(rlang,.data) -importFrom(stats,dist) -importFrom(stats,fitted) -importFrom(stats,formula) -importFrom(stats,hclust) -importFrom(stats,lm) -importFrom(stats,loess) -importFrom(stats,median) -importFrom(stats,na.omit) -importFrom(stats,p.adjust) -importFrom(stats,predict) -importFrom(stats,qbinom) -importFrom(stats,qnorm) -importFrom(stats,qt) -importFrom(stats,quantile) -importFrom(stats,resid) -importFrom(stats,residuals) -importFrom(stats,sd) -importFrom(stats,vcov) -importFrom(stats,xtabs) -importFrom(survival,Surv) -importFrom(survival,survreg) -importFrom(utils,combn) -importFrom(utils,sessionInfo) -importFrom(utils,setTxtProgressBar) -importFrom(utils,txtProgressBar) +importFrom(stats, + dist, + dnorm, + fitted, + formula, + hclust, + lm, + lm.fit, + loess, + median, + model.frame, + model.matrix, + model.response, + na.omit, + p.adjust, + pnorm, + predict, + qbinom, + qnorm, + qt, + quantile, + resid, + residuals, + sd, + vcov, + xtabs +) +importFrom(survival, + Surv, + survreg +) +importFrom(utils, + combn, + sessionInfo, + setTxtProgressBar, + txtProgressBar +) useDynLib(MSstats, .registration=TRUE) diff --git a/R/MSstatsSummarizeWithMultipleCores.R b/R/MSstatsSummarizeWithMultipleCores.R index 1d4678fa..215771e1 100644 --- a/R/MSstatsSummarizeWithMultipleCores.R +++ b/R/MSstatsSummarizeWithMultipleCores.R @@ -278,7 +278,7 @@ #' @noRd .build_summarize_worker <- function( use_TMP, impute, censored_symbol, remove50missing, - aft_iterations, equal_variance + aft_iterations, equal_variance, aft_solver, aft_verbose ) { unpack_fn <- .unpack_protein_slot use_TMP_ <- use_TMP @@ -287,6 +287,8 @@ remove50missing_ <- remove50missing aft_iterations_ <- aft_iterations equal_variance_ <- equal_variance + aft_solver_ <- aft_solver + aft_verbose_ <- aft_verbose function(record) { meta <- record$meta @@ -294,12 +296,13 @@ result <- if (use_TMP_) { MSstatsSummarizeSingleTMP( protein_dt, impute_, censored_symbol_, - remove50missing_, aft_iterations_) + remove50missing_, aft_iterations_, aft_solver_, aft_verbose_) } else { MSstatsSummarizeSingleLinear( protein_dt, impute_, censored_symbol_, remove50missing_, aft_iterations_, - equal_variances = equal_variance_) + equal_variances = equal_variance_, + aft_solver = aft_solver_, aft_verbose = aft_verbose_) } result } @@ -313,9 +316,11 @@ list(worker = i, pid = Sys.getpid(), max_rss_mb = .max_rss_mb()) } +#' @importFrom RhpcBLASctl blas_set_num_threads .warmup_worker <- function(i) { library(MSstats, quietly = TRUE, warn.conflicts = FALSE) data.table::setDTthreads(1) + RhpcBLASctl::blas_set_num_threads(1) NULL } @@ -330,6 +335,14 @@ #' @param equal_variance only for method = "linear"; assume equal variance among feature intensities #' @param numberOfCores number of cores for parallel processing (Linux/Mac only) #' @param aft_iterations number of AFT model iterations +#' @param aft_solver only used when impute = TRUE; linear solve +#' used in the AFT imputation model's Newton-Raphson step: "cholesky" +#' (default, via \code{survival::survreg}), "cg" (conjugate gradient), or +#' "pcg" (conjugate gradient with a Jacobi/inverse-diagonal preconditioner). +#' "cg"/"pcg" are experimental. +#' @param aft_verbose if \code{TRUE}, \code{message()} AFT fitting diagnostics +#' (problem size, elapsed time, and for "cg"/"pcg" per-Newton-iteration +#' conjugate-gradient counts) for every protein fit. Default \code{FALSE}. #' @param verbose whether to print verbose output #' @param BPPARAM optional \code{BiocParallelParam} instance #' @param track_memory whether to report per-worker maximum RSS memory usage. @@ -361,15 +374,18 @@ MSstatsSummarizeWithMultipleCores <- function( equal_variance, numberOfCores = 1L, aft_iterations = 90L, + aft_solver = "cholesky", + aft_verbose = FALSE, verbose = FALSE, BPPARAM = NULL, track_memory = FALSE, max_proteins_per_worker = 50L ) { + .checkAFTSolver(aft_solver) if (numberOfCores <= 1L && is.null(BPPARAM)) { return(MSstatsSummarizeWithSingleCore( input, method, impute, censored_symbol, - remove50missing, equal_variance, aft_iterations)) + remove50missing, equal_variance, aft_iterations, aft_solver, aft_verbose)) } start_time <- proc.time()[["elapsed"]] @@ -419,7 +435,7 @@ MSstatsSummarizeWithMultipleCores <- function( worker_fn <- .build_summarize_worker( use_TMP, impute, censored_symbol, remove50missing, - aft_iterations, equal_variance) + aft_iterations, equal_variance, aft_solver, aft_verbose) if (is.null(BPPARAM)) { tasks <- if (max_proteins_per_worker > 0L) { diff --git a/R/dataProcess.R b/R/dataProcess.R index ef23f518..d9f9c1bf 100755 --- a/R/dataProcess.R +++ b/R/dataProcess.R @@ -62,6 +62,20 @@ #' a logfile named `MSstats_dataProcess_log_progress.log` is created to #' track progress. Only works for Linux & Mac OS. Default is 1. #' @param aft_iterations Number of iterations for AFT model fitting. Default is 90. +#' @param aft_solver Which linear solve to use for the AFT imputation +#' model's Newton-Raphson step: "cholesky" (default) delegates to +#' \code{survival::survreg}, which solves it via Cholesky factorization. +#' "cg" solves the same Newton step with a vendored conjugate-gradient +#' routine instead; "pcg" is the same conjugate-gradient routine with a +#' Jacobi (inverse-diagonal) preconditioner, which can reduce the number +#' of conjugate-gradient iterations needed. "cg"/"pcg" are experimental +#' alternatives, currently opt-in only. +#' @param aft_verbose If \code{TRUE}, \code{message()} diagnostics for +#' every protein fit: problem size and elapsed fitting time for all +#' solvers, plus per-Newton-iteration conjugate-gradient iteration counts +#' and timing when \code{aft_solver} is "cg" or "pcg" - useful for +#' evaluating solver time complexity, but produces one block of output +#' per protein, so leave at the default \code{FALSE} for routine runs. #' @inheritParams .documentFunction #' #' @importFrom utils sessionInfo @@ -130,7 +144,8 @@ dataProcess = function( equalFeatureVar = TRUE, censoredInt = "NA", MBimpute = TRUE, remove50missing = FALSE, fix_missing = NULL, maxQuantileforCensored = 0.999, use_log_file = TRUE, append = FALSE, verbose = TRUE, log_file_path = NULL, - numberOfCores = 1, aft_iterations=90 + numberOfCores = 1, aft_iterations=90, aft_solver = "cholesky", + aft_verbose = FALSE ) { MSstatsConvert::MSstatsLogsSettings(use_log_file, append, verbose, log_file_path, @@ -143,6 +158,7 @@ dataProcess = function( list(method = summaryMethod, equal_var = equalFeatureVar), list(symbol = censoredInt, MB = MBimpute), colnames(raw)) + .checkAFTSolver(aft_solver) peptides_dict = makePeptidesDictionary(as.data.table(unclass(raw)), normalization) input = MSstatsPrepareForDataProcess(raw, logTrans, fix_missing) @@ -164,9 +180,10 @@ dataProcess = function( input = MSstatsPrepareForSummarization(input, summaryMethod, MBimpute, censoredInt, remove_uninformative_feature_outlier) summarized = tryCatch(MSstatsSummarizeWithMultipleCores(input, summaryMethod, - MBimpute, censoredInt, - remove50missing, equalFeatureVar, - numberOfCores, aft_iterations), + MBimpute, censoredInt, + remove50missing, equalFeatureVar, + numberOfCores, aft_iterations, + aft_solver, aft_verbose), error = function(e) { print(e) NULL @@ -211,8 +228,9 @@ dataProcess = function( #' head(summarized[[1]][[1]]) # run-level summary #' MSstatsSummarizeWithSingleCore = function(input, method, impute, censored_symbol, - remove50missing, equal_variance, aft_iterations = 90) { - + remove50missing, equal_variance, aft_iterations = 90, + aft_solver = "cholesky", aft_verbose = FALSE) { + .checkAFTSolver(aft_solver) is_labeled_reference = "is_labeled_ref" %in% colnames(input) && any(input$is_labeled_ref, na.rm = TRUE) if (is_labeled_reference) { @@ -227,8 +245,9 @@ MSstatsSummarizeWithSingleCore = function(input, method, impute, censored_symbol for (protein_id in seq_len(num_proteins)) { single_protein = input[protein_indices[[protein_id]],] summarized_results[[protein_id]] = MSstatsSummarizeSingleTMP( - single_protein, impute, censored_symbol, remove50missing, - aft_iterations) + single_protein, impute, censored_symbol, remove50missing, + aft_iterations, aft_solver = aft_solver, + aft_verbose = aft_verbose) setTxtProgressBar(pb, protein_id) } close(pb) @@ -237,8 +256,9 @@ MSstatsSummarizeWithSingleCore = function(input, method, impute, censored_symbol for (protein_id in seq_len(num_proteins)) { single_protein = input[protein_indices[[protein_id]],] summarized_result = MSstatsSummarizeSingleLinear( - single_protein, impute, censored_symbol, - remove50missing, aft_iterations) + single_protein, impute, censored_symbol, + remove50missing, aft_iterations, aft_solver = aft_solver, + aft_verbose = aft_verbose) summarized_results[[protein_id]] = summarized_result setTxtProgressBar(pb, protein_id) @@ -256,9 +276,16 @@ MSstatsSummarizeWithSingleCore = function(input, method, impute, censored_symbol #' @param remove50missing if TRUE, proteins with more than 50\% missing values in each run are removed #' @param aft_iterations number of iterations for AFT model fitting #' @param equal_variances if TRUE, observation are assumed to be homoskedastic -#' +#' @param aft_solver Which linear solve to use for the AFT imputation +#' model's Newton-Raphson step: "cholesky" (default, via +#' \code{survival::survreg}), "cg" (conjugate gradient), or "pcg" +#' (conjugate gradient with a Jacobi/inverse-diagonal preconditioner). +#' @param aft_verbose If \code{TRUE}, log AFT fitting diagnostics for +#' every protein fit. See \code{.fitSurvival}'s and +#' \code{.fitSurvivalCG}'s \code{verbose}. +#' #' @return list with protein-level data -#' +#' #' @importFrom stats xtabs #' #' @export @@ -286,8 +313,11 @@ MSstatsSummarizeSingleLinear = function(single_protein, censored_symbol, remove50missing, aft_iterations = 90, - equal_variances = TRUE) { + equal_variances = TRUE, + aft_solver = "cholesky", + aft_verbose = FALSE) { ABUNDANCE = RUN = FEATURE = PROTEIN = LogIntensities = NULL + .checkAFTSolver(aft_solver) cols = intersect( colnames(single_protein), @@ -315,7 +345,8 @@ MSstatsSummarizeSingleLinear = function(single_protein, } else { single_protein[, cols, with = FALSE] } - survival_fit = .fitSurvival(fit_data, aft_iterations) + survival_fit = .fitAFTModel(fit_data, aft_iterations, aft_solver, + aft_verbose) sigma2 = survival_fit$scale^2 single_protein[, c("predicted", "imputation_var") := { @@ -437,9 +468,12 @@ MSstatsSummarizeSingleLinear = function(single_protein, #' head(single_protein_summary[[1]]) #' MSstatsSummarizeSingleTMP = function(single_protein, impute, censored_symbol, - remove50missing, aft_iterations = 90) { + remove50missing, aft_iterations = 90, + aft_solver = "cholesky", + aft_verbose = FALSE) { newABUNDANCE = n_obs = n_obs_run = RUN = FEATURE = LABEL = NULL predicted = censored = NULL + .checkAFTSolver(aft_solver) cols = intersect(colnames(single_protein), c("newABUNDANCE", "cen", "RUN", "FEATURE", "ref_covariate")) is_labeled_reference = "is_labeled_ref" %in% colnames(single_protein) && @@ -455,6 +489,7 @@ MSstatsSummarizeSingleTMP = function(single_protein, impute, censored_symbol, # Flag to track convergence warning converged = TRUE + convergence_messages = character(0) fit_data = if (is_labeled_reference) { single_protein[(!is_labeled_ref), cols, with = FALSE] @@ -464,14 +499,30 @@ MSstatsSummarizeSingleTMP = function(single_protein, impute, censored_symbol, # Try to fit survival model and catch convergence warnings survival_fit = withCallingHandlers({ - .fitSurvival(fit_data, aft_iterations) + .fitAFTModel(fit_data, aft_iterations, aft_solver, aft_verbose) }, warning = function(w) { - if (grepl("converge", conditionMessage(w), ignore.case = TRUE)) { - message("Convergence warning caught: ", conditionMessage(w)) + warning_message = conditionMessage(w) + if (grepl("converge", warning_message, ignore.case = TRUE)) { + convergence_messages <<- c(convergence_messages, + warning_message) converged <<- FALSE + invokeRestart("muffleWarning") } }) + protein_name = as.character(unique(single_protein$PROTEIN))[1] + log_fun = getOption("MSstatsLog") + if (!converged) { + msg = paste0("CONVERGENCE WARNING for protein: ", protein_name, + " (", length(convergence_messages), + " warning(s)) - ", + paste(unique(convergence_messages), collapse = " | ")) + message(msg) + if (is.function(log_fun)) { + log_fun("INFO", msg) + } + } + if (converged) { single_protein[, predicted := predict(survival_fit, newdata = .SD)] } else { diff --git a/R/utils_cgsolve.R b/R/utils_cgsolve.R new file mode 100644 index 00000000..3d8eefd8 --- /dev/null +++ b/R/utils_cgsolve.R @@ -0,0 +1,187 @@ +#' Solve a system of linear equations Ax = b with the conjugate gradient method. +#' +#' @section Conjugate gradient is similar to gradient descent, except with +#' how it picks its search direction. +#' +#' Conjugate gradient is an iterative method that solves for "x" in Ax=b. +#' The method reformulates the linear system as an optimization problem where +#' one attempts to minimize Ax^2 - bx. Similar to gradient descent, the solution +#' is initialized at some arbitrary point, then the method iteratively descends +#' toward the optimal solution. But as opposed to gradient descent, which +#' moves in the direction of steepest descent, conjugate gradient picks a +#' different search direction. +#' +#' @section Conjugate gradient search direction is determined based on a linear +#' transformation of the previous iteration's search direction. +#' +#' \strong{Pseudocode for one iteration, written as matrix operations.} +#' Starting from \code{x_0 = 0}, \code{r_0 = b}, \code{z_0 = M^-1 r_0}, +#' \code{p_0 = z_0}, each step computes +#' \preformatted{ +#' q_k = A p_k (the only matrix-vector +#' product per iteration) +#' alpha_k = (r_k' z_k) / (p_k' q_k) step length; p_k' A p_k is +#' the "curvature" +#' x_k+1 = x_k + alpha_k p_k +#' r_k+1 = r_k - alpha_k q_k (= b - A x_k+1, updated +#' without a new product) +#' z_k+1 = M^-1 r_k+1 (elementwise scaling for +#' Jacobi) +#' beta_k = (r_k+1' z_k+1) / (r_k' z_k) (guarantees p_k A p_k+1 = 0, +#' i.e. A-conjugacy) +#' p_k+1 = z_k+1 + beta_k p_k +#' } +#' +#' @section If A is a (p x p) matrix, conjugate gradient is guaranteed to +#' converge in p iterations +#' +#' \strong{Residuals are orthogonal to the step taken.} Given a direction +#' \code{p_k}, CG picks the step size that minimizes the objective along the +#' search direction by setting the derivative to zero with respect to alpha: +#' \preformatted{ +#' d/d alpha phi(x_k + alpha p_k) +#' = p_k' (A x_k + alpha A p_k - b) +#' = -p_k' r_k + alpha p_k' A p_k = 0 +#' => alpha_k = (p_k' r_k) / (p_k' A p_k) +#' +#' x_k+1 = x_k + alpha_k p_k +#' r_k+1 = b - A x_k+1 = r_k - alpha_k A p_k +#' p_k' r_k+1 = p_k' r_k - alpha_k p_k' A p_k +#' = p_k' r_k - (p_k' r_k / p_k' A p_k) p_k' A p_k +#' = p_k' r_k - p_k' r_k +#' = 0 +#' } +#' +#' \strong{Residuals are orthogonal to all other residuals.} +#' +#' Dotting \code{r_k+1 = r_k - alpha_k A p_k} with \code{p_k-1} instead: +#' +#' \preformatted{ +#' p_k-1' r_k = 0 +#' p_k-1' A p_k = 0 (by design of beta to ensure A-conjugacy) +#' p_k-1' r_k+1 = p_k-1' r_k - alpha_k p_k-1' A p_k = 0 +#' +#' r_k = p_k - beta_k-1 p_k-1 +#' r_k' r_k+1 = (p_k - beta_k-1 p_k-1)' r_k+1 +#' = p_k' r_k+1 - beta_k-1 (p_k-1' r_k+1) +#' = 0 - beta_k-1 (0) +#' = 0 +#' } +#' +#' So each new residual is orthogonal to the previous one as well. Since in +#' a p-dimensional system, there can only be p orthogonal residuals, CG +#' is guaranteed to converge in p-iterations. +#' +#' @section Time complexity is empirically linear with respect to the number +#' of entries in the coefficient matrix +#' +#' The time complexity of CG is O(kn), where n is the number of entries in +#' matrix A, and k is the number of iterations. As shown earlier, k is +#' guaranteed to be at most the number of rows in matrix A, i.e. n^0.5. +#' +#' But we can ensure k is small and near a constant given certain data +#' constraints. In the case of MSstats, the following intuition leads +#' k to be empirically small: +#' +#' 1. MSstats Hessian matrix contains zeros on +#' (run x run) or (feature x feature) entries +#' 2. Sparsity makes it more likely diagonals of a matrix dominate +#' 3. Because diagonals more likely dominate, scaling the matrix by the inverse +#' of the diagonals (Jacobi preconditioning) makes the eigenvalues of the +#' resultant matrix to cluster well +#' 4. Well-clustered eigenvalues reduces the number of iterations needed +#' +#' @param coefficient_matrix symmetric positive (semi-)definite matrix, +#' e.g. the negative Hessian from a Newton step. +#' @param right_hand_side vector the system is solved against, e.g. the +#' gradient/score vector from a Newton step. +#' @param use_jacobi_preconditioner if \code{TRUE}, precondition with the +#' inverse of \code{coefficient_matrix}'s own diagonal - cheap to apply, +#' and often enough to cut down the number of iterations needed when the +#' diagonal dominates (as it typically does for an AFT negative Hessian, +#' where each parameter's own curvature tends to be much larger than its +#' cross-terms with the other parameters). Defaults to \code{FALSE}, which +#' reduces exactly to plain (unpreconditioned) conjugate gradient. +#' +#' @return a list with: \code{solution}, the numeric vector solving +#' (approximately) \code{coefficient_matrix \%*\% solution = +#' right_hand_side}; \code{iterations}, how many conjugate-gradient steps +#' were actually taken; \code{converged}, whether the residual tolerance +#' was met; and \code{positive_definite}, whether \code{coefficient_matrix} +#' behaved as positive definite throughout (a caller can fall back to a +#' different matrix, e.g. a Gauss-Newton approximation, when this is +#' \code{FALSE}). +#' +#' @keywords internal +#' @noRd +.cgSolve = function(coefficient_matrix, right_hand_side, + use_jacobi_preconditioner = FALSE) { + number_of_unknowns = nrow(coefficient_matrix) + relative_tolerance = 1e-8 + max_iterations = 10 * number_of_unknowns + solution = rep(0, number_of_unknowns) + + apply_preconditioner = if (use_jacobi_preconditioner) { + diagonal = diag(coefficient_matrix) + inverse_diagonal = ifelse( + is.finite(diagonal) & diagonal > 0, 1 / diagonal, 1) + function(vector) inverse_diagonal * vector + } else { + identity + } + + residual = right_hand_side - drop(coefficient_matrix %*% solution) + preconditioned_residual = apply_preconditioner(residual) + search_direction = preconditioned_residual + residual_size = sum(residual * residual) + residual_dot_preconditioned_residual = + sum(residual * preconditioned_residual) + + convergence_threshold = + (relative_tolerance * max(sqrt(sum(right_hand_side^2)), 1))^2 + + positive_definite = TRUE + iterations_used = 0 + + for (iteration in seq_len(max_iterations)) { + if (residual_size <= convergence_threshold) { + break + } + iterations_used = iteration + matrix_times_search_direction = + drop(coefficient_matrix %*% search_direction) + curvature = sum(search_direction * matrix_times_search_direction) + if (!is.finite(curvature) || curvature <= 0) { + positive_definite = FALSE + warning(".cgSolve: coefficient_matrix is not positive definite ", + "along the current search direction; returning the ", + "best iterate found so far") + break + } + + step_length = residual_dot_preconditioned_residual / curvature + solution = solution + step_length * search_direction + residual = residual - step_length * matrix_times_search_direction + new_residual_size = sum(residual * residual) + + new_preconditioned_residual = apply_preconditioner(residual) + new_residual_dot_preconditioned_residual = + sum(residual * new_preconditioned_residual) + search_direction = new_preconditioned_residual + + (new_residual_dot_preconditioned_residual / + residual_dot_preconditioned_residual) * search_direction + residual_size = new_residual_size + residual_dot_preconditioned_residual = + new_residual_dot_preconditioned_residual + } + + converged = residual_size <= convergence_threshold + if (!converged && positive_definite) { + warning(".cgSolve: did not converge within max_iterations = ", + max_iterations, " iterations; returning the best iterate ", + "found so far") + } + + list(solution = solution, iterations = iterations_used, + converged = converged, positive_definite = positive_definite) +} diff --git a/R/utils_checks.R b/R/utils_checks.R index cf994a59..fb34ad99 100644 --- a/R/utils_checks.R +++ b/R/utils_checks.R @@ -79,6 +79,15 @@ MSstatsPrepareForDataProcess = function(input, log_base, fix_missing) { sink() } +#' Check that aft_solver is one of the supported AFT solvers +#' @param aft_solver string: "cholesky", "cg", or "pcg" +#' @keywords internal +#' @noRd +.checkAFTSolver = function(aft_solver) { + checkmate::assertChoice(aft_solver, c("cholesky", "cg", "pcg"), + .var.name = "aft_solver") +} + #' Check validity of parameters to dataProcess function #' @param log_base of logarithmic transformation #' @param normalization_method string: "quantile", "equalizemedians", "FALSE", diff --git a/R/utils_imputation.R b/R/utils_imputation.R index 99bffcf1..cf210ae2 100644 --- a/R/utils_imputation.R +++ b/R/utils_imputation.R @@ -1,59 +1,680 @@ +#' Decide which predictors go into a single protein's AFT imputation model +#' +#' MSstats fits an accelerated-failure-time (AFT) model per protein to +#' impute left-censored values, and predictors are chosen based on how much +#' information is actually available: whether this is a labeled (SRM) +#' experiment (\code{ref_covariate}), whether there is more than one feature +#' to estimate a \code{FEATURE} effect for, and whether there are enough +#' uncensored observations to estimate that effect at all. +#' +#' @param input data.table with columns \code{newABUNDANCE}, \code{cen}, +#' \code{RUN}, \code{FEATURE}, \code{LABEL}, and (for labeled experiments) +#' \code{ref_covariate}. +#' +#' @return a formula whose left side is +#' \code{Surv(newABUNDANCE, cen, type = "left")}. +#' #' @importFrom data.table uniqueN -#' @importFrom survival survreg Surv +#' @importFrom survival Surv #' @keywords internal -.fitSurvival = function(input, aft_iterations) { +#' @noRd +.buildAFTFormula = function(input) { FEATURE = RUN = NULL - + missingness_filter = is.finite(input$newABUNDANCE) n_total = nrow(input[missingness_filter, ]) n_features = data.table::uniqueN(input[missingness_filter, FEATURE]) n_runs = data.table::uniqueN(input[missingness_filter, RUN]) is_labeled = data.table::uniqueN(input$LABEL) > 1 - countdf = n_total < n_features + n_runs - 1 - # TODO: set.seed here? - set.seed(100) + not_enough_data_for_feature_effect = n_total < n_features + n_runs - 1 + if (is_labeled) { - if (length(unique(input$FEATURE)) == 1) { - # with single feature, not converge, wrong intercept - # need to check - fit = survreg(Surv(newABUNDANCE, cen, type='left') ~ RUN + ref_covariate, - data = input, dist = "gaussian", - control = list(maxiter=aft_iterations)) + if (length(unique(input$FEATURE)) == 1 || + not_enough_data_for_feature_effect) { + # with a single feature (or too little data), a FEATURE term + # either adds nothing or keeps the model from converging / + # gives it the wrong intercept - need to check + Surv(newABUNDANCE, cen, type = "left") ~ RUN + ref_covariate } else { - if (countdf) { - fit = survreg(Surv(newABUNDANCE, cen, type='left') ~ RUN + ref_covariate, - data = input, dist = "gaussian", - control = list(maxiter=aft_iterations)) - } else { - fit = survreg(Surv(newABUNDANCE, cen, type='left') ~ FEATURE + RUN + ref_covariate, - data = input, dist = "gaussian", - control = list(maxiter=aft_iterations)) - } + Surv(newABUNDANCE, cen, type = "left") ~ + FEATURE + RUN + ref_covariate } } else { - if (n_features == 1L) { - fit = survreg(Surv(newABUNDANCE, cen, type = "left") ~ RUN, - data = input, dist = "gaussian", - control = list(maxiter=aft_iterations)) + if (n_features == 1L || not_enough_data_for_feature_effect) { + Surv(newABUNDANCE, cen, type = "left") ~ RUN } else { - if (countdf) { - fit = survreg(Surv(newABUNDANCE, cen, type = "left") ~ RUN, - data = input, dist = "gaussian", - control = list(maxiter=aft_iterations)) - } else { - fit = survreg(Surv(newABUNDANCE, cen, type = "left") ~ FEATURE + RUN, - data = input, dist = "gaussian", - control = list(maxiter=aft_iterations)) - } + Surv(newABUNDANCE, cen, type = "left") ~ FEATURE + RUN } } +} + +#' Fit an AFT survival model with SurvReg dependency +#' +#' @param input data.table with the columns \code{.buildAFTFormula} needs. +#' @param aft_iterations maximum number of iterations for AFT model fitting. +#' @param verbose if \code{TRUE}, \code{message()} the problem size +#' (observations and parameters) before fitting and the wall time the fit +#' took afterwards, mirroring what \code{.fitSurvivalCG}'s \code{verbose} +#' reports. Meant for comparing solvers, not for routine use (this fits +#' one protein at a time). +#' +#' @importFrom stats model.frame model.matrix +#' @importFrom survival survreg +#' @keywords internal +#' @noRd +.fitSurvival = function(input, aft_iterations, verbose = FALSE) { + set.seed(100) + aft_formula = .buildAFTFormula(input) + if (verbose) { + model_frame = model.frame(aft_formula, data = input) + design_matrix = model.matrix(attr(model_frame, "terms"), model_frame) + message(sprintf( + "[AFT-Cholesky] starting fit: %d observations, %d parameters", + nrow(design_matrix), ncol(design_matrix) + 1)) + } + fit_start_time = Sys.time() + fit = survreg(aft_formula, data = input, dist = "gaussian", + control = list(maxiter = aft_iterations)) + if (verbose) { + message(sprintf( + "[AFT-Cholesky] finished: %d iterations, %.4f sec", + fit$iter[length(fit$iter)], + as.numeric(Sys.time() - fit_start_time, units = "secs"))) + } fit$y = NULL fit$linear.predictors = NULL fit } +#' Per-observation log-likelihood and derivatives for a Gaussian AFT model +#' +#' @param linear_predictor current linear predictor +#' (\code{model_matrix \%*\% coefficients}). +#' @param log_scale current log of the scale parameter. +#' @param observed_value observed value (or, for censored rows, the +#' detection-limit ceiling substituted in by +#' \code{.setCensoredByThreshold}). +#' @param exact_indicator \code{1} for an exact/uncensored observation, +#' \code{0} for one left-censored below \code{observed_value}. +#' +#' @return a list with the total \code{log_likelihood}, and +#' per-observation vectors \code{gradient_wrt_linear_predictor}, +#' \code{second_derivative_wrt_linear_predictor}, +#' \code{gradient_wrt_log_scale}, \code{second_derivative_wrt_log_scale}, +#' and \code{cross_derivative} +#' (d2 log_likelihood / d linear_predictor d log_scale). +#' +#' @importFrom stats dnorm pnorm +#' @keywords internal +#' @noRd +.aftGaussianDerivatives = function(linear_predictor, log_scale, + observed_value, exact_indicator) { + scale = exp(log_scale) + inverse_scale_squared = 1 / scale^2 + distance_from_prediction = observed_value - linear_predictor + standardized_distance = distance_from_prediction / scale + + density_at_standardized_distance = dnorm(standardized_distance) + cumulative_probability_at_standardized_distance = + pnorm(standardized_distance) + is_exact_observation = (exact_indicator == 1) + + exact_log_likelihood = + log(density_at_standardized_distance) - log_scale + exact_gradient_wrt_linear_predictor = standardized_distance / scale + exact_log_density_curvature = + (standardized_distance^2 - 1) * inverse_scale_squared + exact_second_derivative_wrt_linear_predictor = + exact_log_density_curvature - + exact_gradient_wrt_linear_predictor^2 + exact_gradient_wrt_log_scale_before_adjustment = + exact_gradient_wrt_linear_predictor * distance_from_prediction + exact_cross_derivative = + distance_from_prediction * exact_log_density_curvature - + exact_gradient_wrt_linear_predictor * + (exact_gradient_wrt_log_scale_before_adjustment + 1) + exact_second_derivative_wrt_log_scale = + distance_from_prediction^2 * exact_log_density_curvature - + exact_gradient_wrt_log_scale_before_adjustment * + (1 + exact_gradient_wrt_log_scale_before_adjustment) + exact_gradient_wrt_log_scale = + exact_gradient_wrt_log_scale_before_adjustment - 1 + + exact_density_underflowed = density_at_standardized_distance <= 0 + exact_log_likelihood = + ifelse(exact_density_underflowed, -200, exact_log_likelihood) + exact_gradient_wrt_linear_predictor = ifelse( + exact_density_underflowed, -standardized_distance / scale, + exact_gradient_wrt_linear_predictor) + exact_second_derivative_wrt_linear_predictor = ifelse( + exact_density_underflowed, -1 / scale, + exact_second_derivative_wrt_linear_predictor) + exact_gradient_wrt_log_scale = + ifelse(exact_density_underflowed, 0, exact_gradient_wrt_log_scale) + exact_cross_derivative = + ifelse(exact_density_underflowed, 0, exact_cross_derivative) + exact_second_derivative_wrt_log_scale = ifelse( + exact_density_underflowed, 0, + exact_second_derivative_wrt_log_scale) + + censored_log_likelihood = + log(cumulative_probability_at_standardized_distance) + censoring_hazard = density_at_standardized_distance / + (cumulative_probability_at_standardized_distance * scale) + censored_gradient_wrt_linear_predictor = -censoring_hazard + censored_log_density_curvature = + -standardized_distance * density_at_standardized_distance * + inverse_scale_squared / + cumulative_probability_at_standardized_distance + censored_second_derivative_wrt_linear_predictor = + censored_log_density_curvature - + censored_gradient_wrt_linear_predictor^2 + censored_gradient_wrt_log_scale = + censored_gradient_wrt_linear_predictor * distance_from_prediction + censored_cross_derivative = + distance_from_prediction * censored_log_density_curvature - + censored_gradient_wrt_linear_predictor * + (censored_gradient_wrt_log_scale + 1) + censored_second_derivative_wrt_log_scale = + distance_from_prediction^2 * censored_log_density_curvature - + censored_gradient_wrt_log_scale * (1 + censored_gradient_wrt_log_scale) + + censored_probability_underflowed = + cumulative_probability_at_standardized_distance <= 0 + censored_log_likelihood = ifelse( + censored_probability_underflowed, -200, censored_log_likelihood) + censored_gradient_wrt_linear_predictor = ifelse( + censored_probability_underflowed, -standardized_distance / scale, + censored_gradient_wrt_linear_predictor) + censored_second_derivative_wrt_linear_predictor = ifelse( + censored_probability_underflowed, 0, + censored_second_derivative_wrt_linear_predictor) + censored_gradient_wrt_log_scale = ifelse( + censored_probability_underflowed, 0, censored_gradient_wrt_log_scale) + censored_cross_derivative = ifelse( + censored_probability_underflowed, 0, censored_cross_derivative) + censored_second_derivative_wrt_log_scale = ifelse( + censored_probability_underflowed, 0, + censored_second_derivative_wrt_log_scale) + + list( + log_likelihood = sum(ifelse( + is_exact_observation, exact_log_likelihood, + censored_log_likelihood)), + gradient_wrt_linear_predictor = ifelse( + is_exact_observation, exact_gradient_wrt_linear_predictor, + censored_gradient_wrt_linear_predictor), + second_derivative_wrt_linear_predictor = ifelse( + is_exact_observation, + exact_second_derivative_wrt_linear_predictor, + censored_second_derivative_wrt_linear_predictor), + gradient_wrt_log_scale = ifelse( + is_exact_observation, exact_gradient_wrt_log_scale, + censored_gradient_wrt_log_scale), + second_derivative_wrt_log_scale = ifelse( + is_exact_observation, exact_second_derivative_wrt_log_scale, + censored_second_derivative_wrt_log_scale), + cross_derivative = ifelse( + is_exact_observation, exact_cross_derivative, + censored_cross_derivative) + ) +} + +#' Evaluate the Gaussian AFT log-likelihood and derivatives at a +#' parameter guess +#' +#' @param design_matrix model matrix of the AFT fit. +#' @param coefficients current regression coefficients. +#' @param log_scale current log of the scale parameter. +#' @param observed_value observed (or censoring-threshold) values. +#' @param exact_indicator \code{1} for exact rows, \code{0} for +#' left-censored rows. +#' +#' @return the list returned by \code{.aftGaussianDerivatives}. +#' +#' @keywords internal +#' @noRd +.evaluateAFTLogLikelihood = function(design_matrix, coefficients, + log_scale, observed_value, + exact_indicator) { + .aftGaussianDerivatives( + drop(design_matrix %*% coefficients), log_scale, + observed_value, exact_indicator) +} + +#' Assemble the AFT gradient vector +#' +#' @param design_matrix model matrix of the AFT fit. +#' @param derivatives output of \code{.aftGaussianDerivatives}. +#' +#' @return gradient of the log-likelihood with respect to the regression +#' coefficients followed by the log scale. +#' +#' @keywords internal +#' @noRd +.buildAFTGradient = function(design_matrix, derivatives) { + c(as.vector(crossprod( + design_matrix, derivatives$gradient_wrt_linear_predictor)), + sum(derivatives$gradient_wrt_log_scale)) +} + +#' Assemble the negative Hessian of the AFT log-likelihood +#' +#' @param design_matrix model matrix of the AFT fit. +#' @param derivatives output of \code{.aftGaussianDerivatives}. +#' +#' @return negative Hessian of the log-likelihood over the regression +#' coefficients and the log scale. +#' +#' @keywords internal +#' @noRd +.buildAFTNegativeHessian = function(design_matrix, derivatives) { + regression_block = -crossprod( + design_matrix, + design_matrix * derivatives$second_derivative_wrt_linear_predictor) + cross_block = -as.vector( + crossprod(design_matrix, derivatives$cross_derivative)) + scale_block = -sum(derivatives$second_derivative_wrt_log_scale) + rbind(cbind(regression_block, cross_block), + c(cross_block, scale_block)) +} + +#' Check that an AFT log-likelihood evaluation is finite +#' +#' @param derivatives output of \code{.aftGaussianDerivatives}. +#' +#' @return \code{TRUE} if the log-likelihood and all first/second +#' derivatives used by the Newton step are finite. +#' +#' @keywords internal +#' @noRd +.isFiniteAFTFit = function(derivatives) { + is.finite(derivatives$log_likelihood) && + all(is.finite(derivatives$gradient_wrt_linear_predictor)) && + all(is.finite(derivatives$gradient_wrt_log_scale)) && + all(is.finite(derivatives$second_derivative_wrt_linear_predictor)) && + all(is.finite(derivatives$second_derivative_wrt_log_scale)) +} + +#' Gauss-Newton (outer-product-of-gradients) approximation to the AFT +#' negative Hessian +#' +#' A fallback when the negative Hessian is not positive definite. +#' +#' @param design_matrix model matrix of the AFT fit. +#' @param derivatives output of \code{.aftGaussianDerivatives}. +#' +#' @return crossproduct of the per-observation gradient contributions. +#' +#' @keywords internal +#' @noRd +.buildGaussNewtonApproximation = function(design_matrix, derivatives) { + per_observation_gradient_contributions = cbind( + design_matrix * derivatives$gradient_wrt_linear_predictor, + derivatives$gradient_wrt_log_scale) + crossprod(per_observation_gradient_contributions) +} + +#' Run \code{.cgSolve} with its "not positive definite" warning muffled +#' +#' @param ... passed to \code{.cgSolve}. +#' +#' @return the output of \code{.cgSolve}. +#' +#' @keywords internal +#' @noRd +.cgSolveMufflingPDWarning = function(...) { + withCallingHandlers( + .cgSolve(...), + warning = function(w) { + if (grepl("not positive definite", conditionMessage(w))) { + invokeRestart("muffleWarning") + } + }) +} + +#' Solve for one AFT Newton-Raphson step with conjugate gradient +#' +#' @param design_matrix model matrix of the AFT fit. +#' @param negative_hessian output of \code{.buildAFTNegativeHessian}. +#' @param derivatives output of \code{.aftGaussianDerivatives}. +#' @param gradient output of \code{.buildAFTGradient}. +#' @param use_jacobi_preconditioner passed to \code{.cgSolve}. +#' +#' @return a list with the \code{step}, \code{primary_iterations}, +#' \code{used_fallback}, and \code{fallback_iterations}. +#' +#' @keywords internal +#' @noRd +.solveAFTNewtonStep = function(design_matrix, negative_hessian, + derivatives, gradient, + use_jacobi_preconditioner) { + primary_solve = .cgSolveMufflingPDWarning( + negative_hessian, gradient, + use_jacobi_preconditioner = use_jacobi_preconditioner) + if (primary_solve$positive_definite) { + list(step = primary_solve$solution, + primary_iterations = primary_solve$iterations, + used_fallback = FALSE, fallback_iterations = 0L) + } else { + fallback_solve = .cgSolve( + .buildGaussNewtonApproximation(design_matrix, derivatives), + gradient, + use_jacobi_preconditioner = use_jacobi_preconditioner) + list(step = fallback_solve$solution, + primary_iterations = primary_solve$iterations, + used_fallback = TRUE, + fallback_iterations = fallback_solve$iterations) + } +} + +#' Fit a Gaussian, left-censored AFT model with a conjugate-gradient +#' Newton step (rather than a cholesky solve). Maximum likelihood estimation +#' loop was written to match the survival package (survreg6.c) to ensure +#' results match with survreg. +#' +#' @section Under the hood, the AFT model is fit with maximum likelihood +#' estimation, where the objective is a Gaussian density for observed rows and +#' CDF for censored rows. +#' +#' \code{phi} for the standard normal density +#' (\code{dnorm}), and \code{Phi} for its CDF (\code{pnorm}). +#' Maximum likelihood picks the set of parameter values \code{theta} +#' under which the data we saw were most probable. What we "saw" differs by +#' whether a row is observed or censored. +#' +#' (\code{exact_indicator}): +#' \itemize{ +#' \item An \emph{observed} (exact) row has a known value, so it +#' contributes the normal density evaluated at that value: +#' \code{L_i = (1 / sigma) phi(z_i)}. +#' \item A \emph{censored} row is one whose intensity fell below the +#' detection limit. Its true value is unknown; all we know is that it lies +#' somewhere below the threshold \code{c_i}. +#' \code{L_i = P(Y_i <= c_i) = Phi((c_i - mu_i) / sigma)}. +#' } +#' +#' Taking logs and summing over rows gives the objective that is maximized: +#' \preformatted{ +#' l(theta) = sum_{observed} [ log phi(z_i) - log sigma ] +#' + sum_{censored} log Phi(z_i) +#' } +#' +#' The first sum is, up to a constant, ordinary least squares; the second +#' is what pulls \code{mu_i} and \code{sigma} toward values that make the +#' censored rows plausibly low. If a censored row's \code{mu_i} is well +#' above its threshold, \code{Phi(z_i)} is tiny and \code{l} is heavily +#' penalized. +#' +#' @section Gradient ascent is performed to maximize the log likelihood. +#' +#' The maximum log likelihood is where the gradient is zero. We compute +#' derivatives with respect to each parameter: +#' +#' \preformatted{ +#' observed: d l_i / d mu_i = z_i / sigma +#' censored: d l_i / d mu_i = -phi(z_i) / (sigma Phi(z_i)) +#' } +#' +#' The observed term is the usual least-squares residual pull; the censored +#' term always pushes \code{mu_i} down, strongly +#' when the prediction sits above the threshold and negligibly when it is +#' already well below. +#' +#' @section Step size is determined with the negative Hessian. +#' +#' Newton's method updates \code{theta} using the step that exactly +#' maximizes a second-order Taylor approximation of the +#' log-likelihood around the current estimate. +#' +#' \preformatted{ +#' l(theta) ~ l(theta_0) + g'(theta - theta_0) +#' + 1/2 (theta - theta_0)' H (theta - theta_0) +#' } +#' +#' Setting the derivative of this quadratic to zero and solving for +#' \code{theta} gives the update: +#' +#' \preformatted{ +#' theta_new = theta + (-H)^-1 * gradient +#' } +#' +#' Another way to think about this is that the Newton method rescales each +#' component of the gradient by an amount determined by local curvature, +#' rather than applying a single global step size. For example, if \code{-H} +#' were diagonal, this would reduce to an entry-specific +#' step size for each parameter: \code{theta_new_i = theta_i + +#' gradient_i / (-H_ii)}. Directions with sharp curvature (large +#' \code{|H_ii|}) get small steps, since the gradient there changes +#' quickly and is only locally reliable; directions with flat curvature +#' get large steps. +#' +#' In general \code{-H} is not diagonal, so \code{(-H)^-1} does not +#' just rescale each gradient entry independently. It captures how +#' curvature in one parameter's direction depends on the value of +#' another. This coupling is what makes Newton's method converge faster +#' than methods that rescale each coordinate independently. +#' +#' @param input data.table, the same shape \code{.fitSurvival} expects. +#' @param aft_iterations maximum number of log-likelihood evaluations the +#' fit may spend. Newton-Raphson iterations and the step-halvings used to +#' recover from an overshooting step share this one budget; once it is +#' exhausted fitting stops and the last accepted coefficients and scale +#' are returned (with a non-convergence warning), rather than failing. +#' @param convergence_tolerance stop once the change in log-likelihood +#' between iterations falls below this (matches the default +#' \code{rel.tolerance} in \code{survival::survreg.control}). +#' @param use_jacobi_preconditioner if \code{TRUE}, precondition every +#' conjugate-gradient solve with the inverse of the current negative +#' Hessian's own diagonal (see \code{.cgSolve}'s +#' \code{use_jacobi_preconditioner}). +#' @param verbose if \code{TRUE}, \code{message()} a line per +#' Newton-Raphson iteration - conjugate-gradient iterations used, whether +#' the Gauss-Newton fallback (see below) was needed, elapsed time, and the +#' resulting log-likelihood - plus a one-line summary once fitting +#' finishes. +#' +#' @return a fitted model of class \code{"survreg"}, with one added field: +#' \code{cg_diagnostics}, a data.frame with one row per Newton-Raphson +#' iteration recording the conjugate-gradient iteration counts and timing +#' described above +#' +#' @importFrom stats model.frame model.matrix model.response lm.fit sd +#' @keywords internal +#' @noRd +.fitSurvivalCG = function(input, aft_iterations, + convergence_tolerance = 1e-9, + use_jacobi_preconditioner = FALSE, + verbose = FALSE) { + model_frame = model.frame(.buildAFTFormula(input), data = input) + model_terms = attr(model_frame, "terms") + design_matrix = model.matrix(model_terms, model_frame) + number_of_coefficients = ncol(design_matrix) + number_of_parameters = number_of_coefficients + 1 + number_of_observations = nrow(design_matrix) + + if (verbose) { + message(sprintf( + "[AFT-CG] starting fit: %d observations, %d parameters, preconditioner = %s", + number_of_observations, number_of_parameters, + if (use_jacobi_preconditioner) "jacobi" else "none")) + } + + response = model.response(model_frame) + observed_value = response[, 1] + exact_indicator = response[, 2] + + initial_fit = lm.fit(design_matrix, observed_value) + coefficients = initial_fit$coefficients + coefficients[!is.finite(coefficients)] = 0 + residual_standard_deviation = sd(initial_fit$residuals) + log_scale = log(max(residual_standard_deviation, 1e-4)) + + current_fit = + .evaluateAFTLogLikelihood(design_matrix, coefficients, log_scale, + observed_value, exact_indicator) + current_log_likelihood = current_fit$log_likelihood + number_of_iterations_used = 0 + converged = FALSE + iterations_remaining = aft_iterations + cg_diagnostics = vector("list", aft_iterations) + + iteration = 0 + while (iterations_remaining > 0) { + iteration = iteration + 1 + iterations_remaining = iterations_remaining - 1 + number_of_iterations_used = iteration + iteration_start_time = Sys.time() + + gradient = .buildAFTGradient(design_matrix, current_fit) + negative_hessian = + .buildAFTNegativeHessian(design_matrix, current_fit) + newton_step = .solveAFTNewtonStep( + design_matrix, negative_hessian, current_fit, gradient, + use_jacobi_preconditioner) + + elapsed_seconds = + as.numeric(Sys.time() - iteration_start_time, units = "secs") + cg_diagnostics[[iteration]] = data.frame( + newton_iteration = iteration, + cg_iterations = newton_step$primary_iterations + + newton_step$fallback_iterations, + used_gauss_newton_fallback = newton_step$used_fallback, + elapsed_seconds = elapsed_seconds) + if (verbose) { + message(sprintf( + "[AFT-CG] newton iter %d: cg iterations = %d%s, %.4f sec", + iteration, + newton_step$primary_iterations + newton_step$fallback_iterations, + if (newton_step$used_fallback) " (Gauss-Newton fallback used)" else "", + elapsed_seconds)) + } + + candidate_coefficients = + coefficients + newton_step$step[seq_len(number_of_coefficients)] + candidate_log_scale = + log_scale + newton_step$step[number_of_coefficients + 1] + + number_of_halvings = 0 + repeat { + candidate_fit = .evaluateAFTLogLikelihood( + design_matrix, candidate_coefficients, candidate_log_scale, + observed_value, exact_indicator) + candidate_improves = .isFiniteAFTFit(candidate_fit) && + candidate_fit$log_likelihood >= current_log_likelihood + if (candidate_improves || iterations_remaining <= 0) { + break + } + iterations_remaining = iterations_remaining - 1 + number_of_halvings = number_of_halvings + 1 + if (number_of_halvings == 1 && + (log_scale - candidate_log_scale) > 1.1) { + candidate_log_scale = log_scale - 1.1 + } + candidate_coefficients = + (candidate_coefficients + 2 * coefficients) / 3 + candidate_log_scale = (candidate_log_scale + 2 * log_scale) / 3 + } + + if (!candidate_improves) { + break + } + + relative_change = + abs(1 - current_log_likelihood / candidate_fit$log_likelihood) + absolute_change = + abs(candidate_fit$log_likelihood - current_log_likelihood) + + coefficients = candidate_coefficients + log_scale = candidate_log_scale + current_fit = candidate_fit + current_log_likelihood = candidate_fit$log_likelihood + + if (relative_change <= convergence_tolerance || + absolute_change <= convergence_tolerance) { + converged = TRUE + break + } + } + + if (!converged) { + warning("AFT model (CG solver) did not converge within its ", + "iteration budget; returning the last accepted ", + "coefficients") + } + + cg_diagnostics = do.call( + rbind, cg_diagnostics[seq_len(number_of_iterations_used)]) + + if (verbose) { + message(sprintf( + paste0("[AFT-CG] finished: %d newton iterations, ", + "%d total cg iterations, %.4f sec total, converged = %s"), + number_of_iterations_used, sum(cg_diagnostics$cg_iterations), + sum(cg_diagnostics$elapsed_seconds), converged)) + } + + final_negative_hessian = + .buildAFTNegativeHessian(design_matrix, current_fit) + variance_covariance_matrix = tryCatch( + solve(final_negative_hessian), + error = function(e) MASS::ginv(final_negative_hessian)) + + fitted_coefficients = coefficients + names(fitted_coefficients) = colnames(design_matrix) + + is_factor_column = vapply(model_frame, is.factor, logical(1)) + + fit = list( + coefficients = fitted_coefficients, + var = variance_covariance_matrix, + scale = exp(log_scale), + terms = model_terms, + xlevels = lapply(model_frame[is_factor_column], levels), + dist = "gaussian", + iter = number_of_iterations_used, + loglik = current_log_likelihood, + cg_diagnostics = cg_diagnostics + ) + class(fit) = "survreg" + fit +} + +#' Fit the AFT imputation model with the requested solver +#' +#' Shared dispatch used by both \code{MSstatsSummarizeSingleLinear} and +#' \code{MSstatsSummarizeSingleTMP} so the \code{aft_solver}/ +#' \code{aft_verbose} logic lives in one place instead of being duplicated +#' at both call sites. +#' +#' @param input data.table, the same shape \code{.fitSurvival} expects. +#' @param aft_iterations maximum number of iterations for AFT model fitting. +#' @param aft_solver "cholesky" (default, via \code{survival::survreg}), +#' "cg" (conjugate gradient), or "pcg" (conjugate gradient with a +#' Jacobi/inverse-diagonal preconditioner). +#' @param aft_verbose passed through to the chosen solver's +#' \code{verbose}: \code{.fitSurvivalCG}'s for "cg"/"pcg", +#' \code{.fitSurvival}'s for "cholesky". +#' +#' @return a fitted model of class \code{"survreg"}. +#' +#' @keywords internal +#' @noRd +.fitAFTModel = function(input, aft_iterations, aft_solver = "cholesky", + aft_verbose = FALSE) { + .checkAFTSolver(aft_solver) + if (aft_solver == "pcg") { + .fitSurvivalCG(input, aft_iterations, + use_jacobi_preconditioner = TRUE, verbose = aft_verbose) + } else if (aft_solver == "cg") { + .fitSurvivalCG(input, aft_iterations, verbose = aft_verbose) + } else { + .fitSurvival(input, aft_iterations, verbose = aft_verbose) + } +} + #' Get predicted values from a survival model #' @param input data.table #' @return numeric vector of predictions diff --git a/inst/tinytest/test_dataProcess.R b/inst/tinytest/test_dataProcess.R index e0113862..e5a6f5dd 100644 --- a/inst/tinytest/test_dataProcess.R +++ b/inst/tinytest/test_dataProcess.R @@ -347,19 +347,14 @@ expect_equal( ) ) -# MSstatsSummarizeSingleTMP: SRM imputation — H rows must NOT be imputed ------ -# For SRM experiments, H is the normalization reference and must never be -# imputed. Only censored L rows (is_labeled_ref=FALSE) should receive a -# predicted value from the survival model. - -make_srm_impute_input <- function() { - runs <- paste0("R", 1:4) - levels_rc <- c("0", runs) - f1 <- data.table::data.table( +make_srm_imputation_input <- function() { + run_names <- paste0("R", 1:4) + reference_covariate_levels <- c("0", run_names) + feature_one_rows <- data.table::data.table( PROTEIN = "P1", FEATURE = "F1", LABEL = c("H","H","H","H", "L","L","L","L"), - RUN = c(runs, runs), + RUN = c(run_names, run_names), # F1-H-R1 censored (H reference — must NOT be imputed) # F1-L-R2 censored (light peptide — MUST be imputed) newABUNDANCE = c(NA, 10.5, 11.0, 11.5, 14.0, NA, 15.0, 15.5), @@ -367,59 +362,161 @@ make_srm_impute_input <- function() { cen = c(0L, 1L, 1L, 1L, 1L, 0L, 1L, 1L), is_labeled_ref = c(TRUE,TRUE,TRUE,TRUE, FALSE,FALSE,FALSE,FALSE) ) - f2 <- data.table::data.table( + feature_two_rows <- data.table::data.table( PROTEIN = "P1", FEATURE = "F2", LABEL = c("H","H","H","H", "L","L","L","L"), - RUN = c(runs, runs), + RUN = c(run_names, run_names), newABUNDANCE = c(10.0,10.5,11.0,11.5, 14.0,14.5,15.0,15.5), censored = rep(FALSE, 8), cen = rep(1L, 8), is_labeled_ref = c(TRUE,TRUE,TRUE,TRUE, FALSE,FALSE,FALSE,FALSE) ) - dt <- data.table::rbindlist(list(f1, f2)) - dt[, ref_covariate := factor( + srm_input <- data.table::rbindlist(list(feature_one_rows, feature_two_rows)) + srm_input[, ref_covariate := factor( ifelse(is_labeled_ref == FALSE, as.character(RUN), "0"), - levels = levels_rc + levels = reference_covariate_levels )] - dt[, FEATURE := factor(FEATURE)] - dt[, RUN := factor(RUN)] - dt[, n_obs := 4L] - dt[, n_obs_run := 2L] - dt[, ANOMALYSCORES := NA_real_] - dt + srm_input[, FEATURE := factor(FEATURE)] + srm_input[, RUN := factor(RUN)] + srm_input[, n_obs := 4L] + srm_input[, n_obs_run := 2L] + srm_input[, ANOMALYSCORES := NA_real_] + srm_input } -result_srm_imp <- MSstatsSummarizeSingleTMP( - make_srm_impute_input(), - impute = TRUE, - censored_symbol = "NA", - remove50missing = FALSE, - aft_iterations = 90 -) +make_srm_imputation_input_with_noise <- function(seed) { + input <- make_srm_imputation_input() + set.seed(seed) + input[cen == 1L, + newABUNDANCE := newABUNDANCE + rnorm(.N, sd = 0.01)] + input +} -survival_srm <- result_srm_imp[[2]] +get_censored_row_predictions <- function(input, aft_solver) { + survival_predictions <- MSstatsSummarizeSingleTMP( + input, + impute = TRUE, + censored_symbol = "NA", + remove50missing = FALSE, + aft_iterations = 90, + aft_solver = aft_solver + )[[2]] + get_feature_one_prediction <- function(label, run) { + survival_predictions[ + as.character(FEATURE) == "F1" & + as.character(LABEL) == label & + as.character(RUN) == run, + predicted + ] + } + list( + censored_heavy = get_feature_one_prediction("H", "R1"), + censored_light = get_feature_one_prediction("L", "R2") + ) +} -# Censored H reference row: predicted must remain NA (not imputed) -h_cens_pred <- survival_srm[ - as.character(FEATURE) == "F1" & - as.character(LABEL) == "H" & - as.character(RUN) == "R1", - predicted -] -expect_true( - length(h_cens_pred) > 0 && all(is.na(h_cens_pred)), - info = "MSstatsSummarizeSingleTMP SRM: censored H rows must NOT receive an imputed predicted value" +expect_heavy_not_imputed_and_light_imputed <- function(predictions, + description) { + expect_true( + length(predictions$censored_heavy) > 0 && + all(is.na(predictions$censored_heavy)), + info = sprintf("MSstatsSummarizeSingleTMP SRM (%s): censored H rows must NOT receive an imputed predicted value", description) + ) + expect_true( + length(predictions$censored_light) > 0 && + all(is.finite(predictions$censored_light)), + info = sprintf("MSstatsSummarizeSingleTMP SRM (%s): censored L rows must receive a finite imputed predicted value", description) + ) +} + +aft_solvers <- c("cholesky", "cg", "pcg") +noisy_input_predictions_by_solver <- lapply( + setNames(nm = aft_solvers), + function(solver) { + get_censored_row_predictions( + make_srm_imputation_input_with_noise(seed = 1), solver + ) + } ) -# Censored L row: predicted must be a finite imputed value -l_cens_pred <- survival_srm[ - as.character(FEATURE) == "F1" & - as.character(LABEL) == "L" & - as.character(RUN) == "R2", - predicted -] -expect_true( - length(l_cens_pred) > 0 && all(is.finite(l_cens_pred)), - info = "MSstatsSummarizeSingleTMP SRM: censored L rows must receive a finite imputed predicted value" +for (solver in aft_solvers) { + expect_heavy_not_imputed_and_light_imputed( + noisy_input_predictions_by_solver[[solver]], + sprintf("aft_solver = %s", solver) + ) +} + +for (solver in setdiff(aft_solvers, "cholesky")) { + expect_equal( + noisy_input_predictions_by_solver[[solver]]$censored_light, + noisy_input_predictions_by_solver[["cholesky"]]$censored_light, + tolerance = 1e-6, check.attributes = FALSE, + info = sprintf("MSstatsSummarizeSingleTMP SRM: aft_solver = %s should closely match aft_solver = cholesky", solver) + ) +} + +srm_input <- make_srm_imputation_input_with_noise(seed = 1) + +expect_error( + get_censored_row_predictions(srm_input, "cgp"), + pattern = "aft_solver", + info = "MSstatsSummarizeSingleTMP should reject an unsupported aft_solver" +) +expect_error( + MSstatsSummarizeSingleLinear(srm_input, impute = TRUE, + censored_symbol = "NA", + remove50missing = FALSE, + aft_solver = "cgp"), + pattern = "aft_solver", + info = "MSstatsSummarizeSingleLinear should reject an unsupported aft_solver" +) +expect_error( + MSstatsSummarizeWithSingleCore(srm_input, "TMP", TRUE, "NA", FALSE, TRUE, + aft_solver = "cgp"), + pattern = "aft_solver", + info = "MSstatsSummarizeWithSingleCore should reject an unsupported aft_solver" ) +expect_error( + MSstatsSummarizeWithMultipleCores(srm_input, "TMP", TRUE, "NA", FALSE, + TRUE, aft_solver = "cgp"), + pattern = "aft_solver", + info = "MSstatsSummarizeWithMultipleCores should reject an unsupported aft_solver" +) +expect_error( + dataProcess(DDARawData, aft_solver = "cgp", use_log_file = FALSE, + verbose = FALSE), + pattern = "aft_solver", + info = "dataProcess should reject an unsupported aft_solver before summarization" +) + +for (solver in aft_solvers) { + raw_warnings <- character(0) + convergence_messages <- character(0) + unconverged_summary <- withCallingHandlers( + MSstatsSummarizeSingleTMP(srm_input, impute = TRUE, + censored_symbol = "NA", + remove50missing = FALSE, + aft_iterations = 2, aft_solver = solver), + warning = function(w) { + raw_warnings <<- c(raw_warnings, conditionMessage(w)) + invokeRestart("muffleWarning") + }, + message = function(m) { + convergence_messages <<- c(convergence_messages, + conditionMessage(m)) + invokeRestart("muffleMessage") + }) + expect_equal( + sum(grepl("CONVERGENCE WARNING", convergence_messages)), 1L, + info = sprintf("MSstatsSummarizeSingleTMP (aft_solver = %s): a non-converged fit should emit one combined convergence message", solver) + ) + expect_false( + any(grepl("converge", raw_warnings, ignore.case = TRUE)), + info = sprintf("MSstatsSummarizeSingleTMP (aft_solver = %s): the solver's raw convergence warning should be muffled", solver) + ) + expect_true( + all(is.na(unconverged_summary[[2]]$predicted)), + info = sprintf("MSstatsSummarizeSingleTMP (aft_solver = %s): a non-converged fit should not be used for imputation", solver) + ) +} diff --git a/inst/tinytest/test_utils_cgsolve.R b/inst/tinytest/test_utils_cgsolve.R new file mode 100644 index 00000000..1b9ceaaa --- /dev/null +++ b/inst/tinytest/test_utils_cgsolve.R @@ -0,0 +1,108 @@ +make_random_solvable_matrix <- function(size, seed, diagonal_boost = 0.01) { + # Builds a random symmetric positive-definite matrix + set.seed(seed) + random_matrix <- matrix(rnorm(size * size), size, size) + random_matrix %*% t(random_matrix) + diag(size) * diagonal_boost +} + +for (size in c(2, 5, 10, 30, 80)) { + coefficient_matrix <- make_random_solvable_matrix(size, seed = size) + set.seed(size + 1000) + right_hand_side <- rnorm(size) + + iterative_result <- MSstats:::.cgSolve(coefficient_matrix, right_hand_side) + exact_solution <- solve(coefficient_matrix, right_hand_side) + + expect_equal( + iterative_result$solution, exact_solution, tolerance = 1e-6, + info = paste0(".cgSolve should give the same answer as solve() on a ", + "random solvable symmetric system of size ", size) + ) + expect_true( + iterative_result$converged && iterative_result$positive_definite, + info = paste0("A well-behaved symmetric system of size ", size, + " should report converged = TRUE and ", + "positive_definite = TRUE") + ) + expect_true( + iterative_result$iterations >= 1 && + iterative_result$iterations <= size * 10, + info = paste("The number of steps taken should be at least one and", + "should not go over the maximum allowed") + ) +} + +make_random_unsolvable_matrix <- function(size, seed) { + unsolvable_matrix <- make_random_solvable_matrix(size, seed) + unsolvable_matrix[1, ] <- 0 + unsolvable_matrix[, 1] <- 0 + unsolvable_matrix +} + +unsolvable_matrix <- make_random_unsolvable_matrix(10, seed = 42) +set.seed(43) +right_hand_side <- rnorm(10) + +expect_warning( + unsolvable_result <- MSstats:::.cgSolve(unsolvable_matrix, right_hand_side), + info = paste("An unsolvable matrix should produce a warning rather than", + "an error or an endless loop") +) +expect_true( + all(is.finite(unsolvable_result$solution)), + info = paste("An unsolvable system should still return a partial answer", + "made of ordinary finite numbers") +) +expect_false( + unsolvable_result$converged && unsolvable_result$positive_definite, + info = paste("An unsolvable matrix should be flagged by reporting", + "converged = FALSE, positive_definite = FALSE, or both") +) + +make_large_diagonal_matrix <- function(size, seed) { + set.seed(seed) + small_off_diagonal_entries <- + matrix(runif(size * size, -0.1, 0.1), size, size) + small_off_diagonal_entries <- + (small_off_diagonal_entries + t(small_off_diagonal_entries)) / 2 + diag(small_off_diagonal_entries) <- 0 + diag(size) * runif(size, 5, 10) + small_off_diagonal_entries +} + +large_diagonal_matrix <- make_large_diagonal_matrix(40, seed = 11) +set.seed(12) +large_diagonal_right_hand_side <- rnorm(40) +large_diagonal_exact_solution <- + solve(large_diagonal_matrix, large_diagonal_right_hand_side) + +result_without_scaling <- MSstats:::.cgSolve( + large_diagonal_matrix, large_diagonal_right_hand_side) +result_with_scaling <- MSstats:::.cgSolve( + large_diagonal_matrix, large_diagonal_right_hand_side, + use_jacobi_preconditioner = TRUE) + +expect_equal( + result_with_scaling$solution, large_diagonal_exact_solution, + tolerance = 1e-6, + info = paste("With diagonal scaling turned on, .cgSolve should still", + "give the same answer as solve()") +) +expect_true( + result_with_scaling$iterations <= result_without_scaling$iterations, + info = paste("Diagonal scaling should not need more steps than", + "no scaling when the diagonal entries are large (without", + "scaling =", result_without_scaling$iterations, + ", with scaling =", result_with_scaling$iterations, ")") +) + +zero_diagonal_matrix <- make_random_solvable_matrix(8, seed = 55) +zero_diagonal_matrix[3, 3] <- 0 +set.seed(56) +zero_diagonal_right_hand_side <- rnorm(8) +expect_true( + all(is.finite(suppressWarnings(MSstats:::.cgSolve( + zero_diagonal_matrix, zero_diagonal_right_hand_side, + use_jacobi_preconditioner = TRUE))$solution)), + info = paste("A zero on the diagonal should not cause diagonal scaling", + "to return infinite or missing values") +) diff --git a/inst/tinytest/test_utils_imputation.R b/inst/tinytest/test_utils_imputation.R index 03095cf5..2a441d20 100644 --- a/inst/tinytest/test_utils_imputation.R +++ b/inst/tinytest/test_utils_imputation.R @@ -1,133 +1,194 @@ -make_surv_labeled_single <- function() { - runs <- paste0("R", 1:3) - dt <- data.table::rbindlist(list( - data.table::data.table( - FEATURE = factor(rep("F1", 9)), - RUN = factor(rep(runs, each = 3)), - LABEL = "H", - newABUNDANCE = seq(10.1, by = 0.1, length.out = 9), - cen = 1L - ), - data.table::data.table( - FEATURE = factor(rep("F1", 9)), - RUN = factor(rep(runs, each = 3)), - LABEL = "L", - newABUNDANCE = seq(14.1, by = 0.1, length.out = 9), - cen = 1L - ) - )) - ref_vals <- ifelse(dt$LABEL == "L", as.character(dt$RUN), "0") - dt[["ref_covariate"]] <- factor(ref_vals, levels = c("0", runs)) - dt -} +number_of_iterations <- 90 -make_surv_labeled_multi_welldetermined <- function() { - features <- paste0("F", 1:3) - runs <- paste0("R", 1:4) - dt <- data.table::CJ(FEATURE = features, RUN = runs, LABEL = c("H", "L")) - dt[, FEATURE := factor(FEATURE)] - dt[, RUN := factor(RUN)] - dt[, newABUNDANCE := 10 + as.integer(FEATURE) + as.integer(RUN) * 0.5 + - ifelse(LABEL == "L", 4, 0)] - dt[, cen := 1L] - ref_vals <- ifelse(dt$LABEL == "L", as.character(dt$RUN), "0") - dt[["ref_covariate"]] <- factor(ref_vals, levels = c("0", runs)) - dt +add_reference_covariate <- function(input) { + input[, ref_covariate := factor( + ifelse(LABEL == "L", as.character(RUN), "0"), + levels = c("0", levels(RUN)))] } -make_surv_labeled_underdetermined <- function() { - runs <- c("R1", "R2", "R3") - dt_h <- data.table::data.table( - FEATURE = factor(paste0("F", 1:8)), - RUN = factor(rep_len(runs, 8)), - LABEL = "H", - newABUNDANCE = seq(10, by = 0.5, length.out = 8), - cen = 1L +make_survival_input <- function(number_of_features, number_of_runs, is_labeled, + number_of_replicates = 1L, + noise_standard_deviation = 0.1, + censored_fraction = 0, seed = 1) { + set.seed(seed) + input <- data.table::CJ( + FEATURE = paste0("F", seq_len(number_of_features)), + RUN = paste0("R", seq_len(number_of_runs)), + LABEL = if (is_labeled) c("H", "L") else "L", + REPLICATE = seq_len(number_of_replicates) ) - dt_l <- data.table::data.table( - FEATURE = factor("F1"), - RUN = factor("R1"), - LABEL = "L", - newABUNDANCE = 14, - cen = 1L - ) - dt <- data.table::rbindlist(list(dt_h, dt_l)) - ref_vals <- ifelse(dt$LABEL == "L", as.character(dt$RUN), "0") - dt[["ref_covariate"]] <- factor(ref_vals, levels = c("0", runs)) - dt + input[, FEATURE := factor(FEATURE)] + input[, RUN := factor(RUN)] + input[, newABUNDANCE := 10 + as.integer(FEATURE) + as.integer(RUN) * 0.5 + + ifelse(LABEL == "L", 4, 0) + (REPLICATE - 1) * 0.1 + + rnorm(.N, sd = noise_standard_deviation)] + input[, REPLICATE := NULL] + input[, cen := 1L] + if (censored_fraction > 0) { + censoring_threshold <- stats::quantile(input$newABUNDANCE, + censored_fraction) + input[newABUNDANCE < censoring_threshold, + `:=`(cen = 0L, newABUNDANCE = censoring_threshold)] + } + if (is_labeled) add_reference_covariate(input) + input } -make_surv_unlabeled_single <- function() { - data.table::data.table( - FEATURE = factor(rep("F1", 15)), - RUN = factor(rep(paste0("R", 1:5), each = 3)), - LABEL = "L", - newABUNDANCE = seq(10, by = 0.5, length.out = 15), +make_labeled_input_with_too_few_observations <- function() { + input <- data.table::data.table( + FEATURE = factor(c(paste0("F", 1:8), "F1")), + RUN = factor(c(rep_len(paste0("R", 1:3), 8), "R1")), + LABEL = c(rep("H", 8), "L"), + newABUNDANCE = c(seq(10, by = 0.5, length.out = 8), 14), cen = 1L ) + add_reference_covariate(input) } -make_surv_unlabeled_multi_welldetermined <- function() { - dt <- data.table::CJ( - FEATURE = paste0("F", 1:3), - RUN = paste0("R", 1:5) +coefficient_names <- function(fit) names(coef(fit)) + +predictor_cases <- list( + "labeled single-feature" = list( + input = make_survival_input(number_of_features = 1, number_of_runs = 3, + is_labeled = TRUE, number_of_replicates = 3), + has_reference_covariate = TRUE, has_feature = FALSE), + "labeled multi-feature, enough observations" = list( + input = make_survival_input(number_of_features = 3, number_of_runs = 4, + is_labeled = TRUE), + has_reference_covariate = TRUE, has_feature = TRUE), + "labeled multi-feature, too few observations" = list( + input = make_labeled_input_with_too_few_observations(), + has_reference_covariate = TRUE, has_feature = FALSE), + "unlabeled single-feature" = list( + input = make_survival_input(number_of_features = 1, number_of_runs = 5, + is_labeled = FALSE, number_of_replicates = 3), + has_reference_covariate = FALSE, has_feature = FALSE), + "unlabeled multi-feature, enough observations" = list( + input = make_survival_input(number_of_features = 3, number_of_runs = 5, + is_labeled = FALSE), + has_reference_covariate = FALSE, has_feature = TRUE) +) + +for (case_name in names(predictor_cases)) { + case <- predictor_cases[[case_name]] + cholesky_coefficient_names <- coefficient_names( + MSstats:::.fitSurvival(case$input, number_of_iterations)) + expect_equal( + any(grepl("ref_covariate", cholesky_coefficient_names)), + case$has_reference_covariate, + info = paste(".fitSurvival", case_name, + ": reference covariate in coefficients") + ) + expect_equal( + any(grepl("^FEATURE", cholesky_coefficient_names)), case$has_feature, + info = paste(".fitSurvival", case_name, ": FEATURE in coefficients") + ) + expect_equal( + sort(coefficient_names( + MSstats:::.fitSurvivalCG(case$input, number_of_iterations))), + sort(cholesky_coefficient_names), + info = paste(case_name, ": the conjugate gradient solver must select", + "the same predictors as the Cholesky solver") ) - dt[, FEATURE := factor(FEATURE)] - dt[, RUN := factor(RUN)] - dt[, LABEL := "L"] - dt[, newABUNDANCE := seq(10, by = 0.5, length.out = .N)] - dt[, cen := 1L] - dt } -coef_names <- function(fit) names(coef(fit)) -surv_labeled_single <- MSstats:::.fitSurvival(make_surv_labeled_single(), 90) +check_solvers_agree <- function(input, label, use_jacobi_preconditioner) { + fit_cholesky <- MSstats:::.fitSurvival(input, number_of_iterations) + fit_conjugate_gradient <- MSstats:::.fitSurvivalCG( + input, number_of_iterations, + use_jacobi_preconditioner = use_jacobi_preconditioner) + summarize_fit <- function(fit) { + predictions <- predict(fit, newdata = input, se.fit = TRUE) + list( + coefficients = fit$coefficients[names(fit_cholesky$coefficients)], + scale = fit$scale, + `predicted values` = predictions$fit, + `prediction standard errors` = predictions$se.fit + ) + } + expected <- summarize_fit(fit_cholesky) + actual <- summarize_fit(fit_conjugate_gradient) + for (quantity in names(expected)) { + expect_equal( + actual[[quantity]], expected[[quantity]], + tolerance = 1e-4, check.attributes = FALSE, + info = paste(label, ":", quantity, + "should match the Cholesky solver") + ) + } +} -expect_true( - any(grepl("ref_covariate", coef_names(surv_labeled_single))), - info = ".fitSurvival labeled single-feature: ref_covariate must appear in coefficients" -) -expect_false( - any(grepl("^FEATURE", coef_names(surv_labeled_single))), - info = ".fitSurvival labeled single-feature: FEATURE must not appear (only one feature)" -) +make_noisy_censored_input <- function(is_labeled, seed) { + make_survival_input(number_of_features = 3, number_of_runs = 4, + is_labeled = is_labeled, + noise_standard_deviation = 0.7, + censored_fraction = 0.2, seed = seed) +} -surv_labeled_multi_wd <- MSstats:::.fitSurvival(make_surv_labeled_multi_welldetermined(), 90) -expect_true( - any(grepl("ref_covariate", coef_names(surv_labeled_multi_wd))), - info = ".fitSurvival labeled multi well-determined: ref_covariate must appear in coefficients" -) -expect_true( - any(grepl("^FEATURE", coef_names(surv_labeled_multi_wd))), - info = ".fitSurvival labeled multi well-determined: FEATURE must appear in coefficients" +noisy_inputs <- list( + labeled = make_noisy_censored_input(is_labeled = TRUE, seed = 1), + unlabeled = make_noisy_censored_input(is_labeled = FALSE, seed = 2) ) +for (input_name in names(noisy_inputs)) { + for (use_jacobi_preconditioner in c(FALSE, TRUE)) { + check_solvers_agree( + noisy_inputs[[input_name]], + label = paste0(input_name, ", noisy, censored", + if (use_jacobi_preconditioner) + ", with Jacobi preconditioner"), + use_jacobi_preconditioner = use_jacobi_preconditioner + ) + } +} -surv_labeled_under <- MSstats:::.fitSurvival(make_surv_labeled_underdetermined(), 90) +noisy_input <- make_noisy_censored_input(is_labeled = FALSE, seed = 3) +fit_cholesky <- MSstats:::.fitAFTModel(noisy_input, number_of_iterations, + "cholesky") +expect_inherits( + fit_cholesky, "survreg", + info = ".fitAFTModel(aft_solver = 'cholesky') should return a survreg fit" +) expect_true( - any(grepl("ref_covariate", coef_names(surv_labeled_under))), - info = ".fitSurvival labeled underdetermined: ref_covariate must appear in fallback coefficients" + is.null(fit_cholesky$cg_diagnostics), + info = "the Cholesky solver should not attach cg_diagnostics" ) -expect_false( - any(grepl("^FEATURE", coef_names(surv_labeled_under))), - info = ".fitSurvival labeled underdetermined: FEATURE must not appear in fallback formula" +for (aft_solver in c("cg", "pcg")) { + fit <- MSstats:::.fitAFTModel(noisy_input, number_of_iterations, aft_solver) + expect_false( + is.null(fit$cg_diagnostics), + info = paste0(".fitAFTModel(aft_solver = '", aft_solver, + "') should attach cg_diagnostics") + ) +} +expect_error( + MSstats:::.fitAFTModel(noisy_input, number_of_iterations, "cgp"), + pattern = "aft_solver", + info = ".fitAFTModel should reject an unsupported aft_solver instead of falling back to Cholesky" ) -surv_unlabeled_single <- MSstats:::.fitSurvival(make_surv_unlabeled_single(), 90) -expect_false( - any(grepl("ref_covariate", coef_names(surv_unlabeled_single))), - info = ".fitSurvival unlabeled single-feature: ref_covariate must not appear" +expect_silent( + MSstats:::.fitSurvivalCG(noisy_input, number_of_iterations, verbose = FALSE) +) +expect_message( + MSstats:::.fitSurvivalCG(noisy_input, number_of_iterations, verbose = TRUE), + pattern = "\\[AFT-CG\\] starting fit", + info = "verbose = TRUE should report the problem size at the start of the fit" ) -expect_false( - any(grepl("^FEATURE", coef_names(surv_unlabeled_single))), - info = ".fitSurvival unlabeled single-feature: FEATURE must not appear (only one feature)" +expect_message( + MSstats:::.fitSurvivalCG(noisy_input, number_of_iterations, verbose = TRUE), + pattern = "\\[AFT-CG\\] finished", + info = "verbose = TRUE should report a summary once fitting finishes" ) -surv_unlabeled_multi_wd <- MSstats:::.fitSurvival(make_surv_unlabeled_multi_welldetermined(), 90) -expect_false( - any(grepl("ref_covariate", coef_names(surv_unlabeled_multi_wd))), - info = ".fitSurvival unlabeled multi well-determined: ref_covariate must not appear" + +fit_with_diagnostics <- MSstats:::.fitSurvivalCG(noisy_input, + number_of_iterations) +expect_equal( + nrow(fit_with_diagnostics$cg_diagnostics), fit_with_diagnostics$iter, + info = "cg_diagnostics should have one row per Newton-Raphson iteration taken" ) expect_true( - any(grepl("^FEATURE", coef_names(surv_unlabeled_multi_wd))), - info = ".fitSurvival unlabeled multi well-determined: FEATURE must appear in coefficients" + all(fit_with_diagnostics$cg_diagnostics$cg_iterations >= 0), + info = paste("the conjugate gradient iteration count should be", + "non-negative for every Newton iteration") ) diff --git a/man/MSstatsSummarizeSingleLinear.Rd b/man/MSstatsSummarizeSingleLinear.Rd index 8595bc16..ed60f8c4 100644 --- a/man/MSstatsSummarizeSingleLinear.Rd +++ b/man/MSstatsSummarizeSingleLinear.Rd @@ -10,7 +10,9 @@ MSstatsSummarizeSingleLinear( censored_symbol, remove50missing, aft_iterations = 90, - equal_variances = TRUE + equal_variances = TRUE, + aft_solver = "cholesky", + aft_verbose = FALSE ) } \arguments{ @@ -25,6 +27,15 @@ MSstatsSummarizeSingleLinear( \item{aft_iterations}{number of iterations for AFT model fitting} \item{equal_variances}{if TRUE, observation are assumed to be homoskedastic} + +\item{aft_solver}{Which linear solve to use for the AFT imputation +model's Newton-Raphson step: "cholesky" (default, via +\code{survival::survreg}), "cg" (conjugate gradient), or "pcg" +(conjugate gradient with a Jacobi/inverse-diagonal preconditioner).} + +\item{aft_verbose}{If \code{TRUE}, log AFT fitting diagnostics for +every protein fit. See \code{.fitSurvival}'s and +\code{.fitSurvivalCG}'s \code{verbose}.} } \value{ list with protein-level data diff --git a/man/MSstatsSummarizeSingleTMP.Rd b/man/MSstatsSummarizeSingleTMP.Rd index cd115723..82046bbd 100644 --- a/man/MSstatsSummarizeSingleTMP.Rd +++ b/man/MSstatsSummarizeSingleTMP.Rd @@ -9,7 +9,9 @@ MSstatsSummarizeSingleTMP( impute, censored_symbol, remove50missing, - aft_iterations = 90 + aft_iterations = 90, + aft_solver = "cholesky", + aft_verbose = FALSE ) } \arguments{ @@ -22,6 +24,16 @@ MSstatsSummarizeSingleTMP( \item{remove50missing}{only for method = "TMP"; drops proteins missing >=50\% per peptide in every run} \item{aft_iterations}{number of iterations for AFT model fitting} + +\item{aft_solver}{only used when impute = TRUE; linear solve +used in the AFT imputation model's Newton-Raphson step: "cholesky" +(default, via \code{survival::survreg}), "cg" (conjugate gradient), or +"pcg" (conjugate gradient with a Jacobi/inverse-diagonal preconditioner). +"cg"/"pcg" are experimental.} + +\item{aft_verbose}{if \code{TRUE}, \code{message()} AFT fitting diagnostics +(problem size, elapsed time, and for "cg"/"pcg" per-Newton-iteration +conjugate-gradient counts) for every protein fit. Default \code{FALSE}.} } \value{ list of two data.tables: one with fitted survival model, diff --git a/man/MSstatsSummarizeWithMultipleCores.Rd b/man/MSstatsSummarizeWithMultipleCores.Rd index b95b1e81..cc33e43d 100644 --- a/man/MSstatsSummarizeWithMultipleCores.Rd +++ b/man/MSstatsSummarizeWithMultipleCores.Rd @@ -13,6 +13,8 @@ MSstatsSummarizeWithMultipleCores( equal_variance, numberOfCores = 1L, aft_iterations = 90L, + aft_solver = "cholesky", + aft_verbose = FALSE, verbose = FALSE, BPPARAM = NULL, track_memory = FALSE, @@ -36,6 +38,16 @@ MSstatsSummarizeWithMultipleCores( \item{aft_iterations}{number of AFT model iterations} +\item{aft_solver}{only used when impute = TRUE; linear solve +used in the AFT imputation model's Newton-Raphson step: "cholesky" +(default, via \code{survival::survreg}), "cg" (conjugate gradient), or +"pcg" (conjugate gradient with a Jacobi/inverse-diagonal preconditioner). +"cg"/"pcg" are experimental.} + +\item{aft_verbose}{if \code{TRUE}, \code{message()} AFT fitting diagnostics +(problem size, elapsed time, and for "cg"/"pcg" per-Newton-iteration +conjugate-gradient counts) for every protein fit. Default \code{FALSE}.} + \item{verbose}{whether to print verbose output} \item{BPPARAM}{optional \code{BiocParallelParam} instance} diff --git a/man/MSstatsSummarizeWithSingleCore.Rd b/man/MSstatsSummarizeWithSingleCore.Rd index fc711d42..c004bded 100644 --- a/man/MSstatsSummarizeWithSingleCore.Rd +++ b/man/MSstatsSummarizeWithSingleCore.Rd @@ -11,7 +11,9 @@ MSstatsSummarizeWithSingleCore( censored_symbol, remove50missing, equal_variance, - aft_iterations = 90 + aft_iterations = 90, + aft_solver = "cholesky", + aft_verbose = FALSE ) } \arguments{ @@ -28,6 +30,16 @@ MSstatsSummarizeWithSingleCore( \item{equal_variance}{only for method = "linear"; assume equal variance among feature intensities} \item{aft_iterations}{Number of iterations for AFT model fitting. Default is 90.} + +\item{aft_solver}{only used when impute = TRUE; linear solve +used in the AFT imputation model's Newton-Raphson step: "cholesky" +(default, via \code{survival::survreg}), "cg" (conjugate gradient), or +"pcg" (conjugate gradient with a Jacobi/inverse-diagonal preconditioner). +"cg"/"pcg" are experimental.} + +\item{aft_verbose}{if \code{TRUE}, \code{message()} AFT fitting diagnostics +(problem size, elapsed time, and for "cg"/"pcg" per-Newton-iteration +conjugate-gradient counts) for every protein fit. Default \code{FALSE}.} } \value{ list of length one with run-level data. diff --git a/man/dataProcess.Rd b/man/dataProcess.Rd index 1ed4fe03..99c6497d 100644 --- a/man/dataProcess.Rd +++ b/man/dataProcess.Rd @@ -25,7 +25,9 @@ dataProcess( verbose = TRUE, log_file_path = NULL, numberOfCores = 1, - aft_iterations = 90 + aft_iterations = 90, + aft_solver = "cholesky", + aft_verbose = FALSE ) } \arguments{ @@ -121,6 +123,22 @@ a logfile named `MSstats_dataProcess_log_progress.log` is created to track progress. Only works for Linux & Mac OS. Default is 1.} \item{aft_iterations}{Number of iterations for AFT model fitting. Default is 90.} + +\item{aft_solver}{Which linear solve to use for the AFT imputation +model's Newton-Raphson step: "cholesky" (default) delegates to +\code{survival::survreg}, which solves it via Cholesky factorization. +"cg" solves the same Newton step with a vendored conjugate-gradient +routine instead; "pcg" is the same conjugate-gradient routine with a +Jacobi (inverse-diagonal) preconditioner, which can reduce the number +of conjugate-gradient iterations needed. "cg"/"pcg" are experimental +alternatives, currently opt-in only.} + +\item{aft_verbose}{If \code{TRUE}, \code{message()} diagnostics for +every protein fit: problem size and elapsed fitting time for all +solvers, plus per-Newton-iteration conjugate-gradient iteration counts +and timing when \code{aft_solver} is "cg" or "pcg" - useful for +evaluating solver time complexity, but produces one block of output +per protein, so leave at the default \code{FALSE} for routine runs.} } \value{ A list containing: