All modelling methods at a glance
4 November 2026
Package
prolfqua 0.99.8
Contents
1 Purpose
build_contrast_analysis() selects the model with its method argument. This vignette runs every method listed by list_facades() on one small simulated dataset and shows one line of output per method. Plots and a discussion of the methods are in the article Comparing modelling methods on one dataset.
Figure 1: build_contrast_analysis() is a factory: the method key selects the model, and every model rolls out with the same interface
2 Data
We simulate peptide intensities for 20 proteins in three groups and compare group A with the control.
library(prolfqua)
sim <- sim_lfq_data_peptide_config(Nprot = 20, seed = 42)
lfq_peptide <- LFQData$new(sim$data, sim$config)
lfq_peptide <- lfq_peptide$get_Transformer()$log2()$lfq
lfq_protein <- lfq_peptide$get_Aggregator()$aggregate()
lfq_limpa <- AggregateLimpa$new(lfq_peptide, "protein")$aggregate()
contrasts <- c(A_vs_Ctrl = "group_A - group_Ctrl")Each method declares the input it needs:
same: protein-level input, herelfq_proteinaggregated with median polish.limpatakeslfq_limpa, aggregated byAggregateLimpa, which adds the per-protein standard errors it needs.nested: peptide-level input,lfq_peptide. The method models the peptides and reports protein-level contrasts.
3 Fit all methods
facades <- list_facades()
input_for <- function(method) {
if (method == "limpa") {
return(lfq_limpa)
}
if (facades[[method]]$needs == "nested") lfq_peptide else lfq_protein
}
model_for <- function(method) {
if (method == "lmer_nested") "~ group_ + (1 | peptide_Id) + (1 | sampleName)" else "~ group_"
}
results <- lapply(names(facades), function(method) {
fa <- build_contrast_analysis(input_for(method), model_for(method), contrasts, method = method)
fa$get_contrasts()
})
names(results) <- names(facades)4 One line per method
about <- c(
lm = "linear model per protein",
lm_impute = "lm; groups without observations imputed at the limit of detection (LOD)",
rlm = "robust linear model (MASS::rlm)",
rfit = "rank-based regression (Rfit)",
rfit_impute = "rfit with LOD imputation",
limma = "limma, empirical-Bayes moderated variances",
limma_impute = "limma with LOD imputation",
limma_voom = "limma with vooma precision weights",
limma_voom_impute = "limma_voom with LOD imputation",
deqms = "limma with DEqMS peptide-count variance moderation",
deqms_voom = "deqms with vooma precision weights",
limpa = "limma on limpa DPC-quantified proteins",
limpa_nested = "limpa, DPC quantification from peptides",
firth = "Firth logistic regression, protein detection",
firth_nested = "Firth logistic regression, peptide detection",
binomial_nested = "quasibinomial model, peptide detection counts",
binomial_firth_nested = "binomial_nested with Firth bias reduction",
hurdle_nested = "lm_impute intensity + binomial_firth_nested detection",
lmer_nested = "mixed model, peptide and sample random effects",
ropeca_nested = "peptide fold changes summarised per protein (ROPECA)"
)
overview <- data.frame(
method = names(results),
input = vapply(facades, function(x) x$needs, character(1)),
proteins = vapply(results, nrow, integer(1)),
FDR_lt_0.1 = vapply(results, function(x) sum(x$FDR < 0.1, na.rm = TRUE), integer(1)),
model = about[names(results)]
)
knitr::kable(overview, row.names = FALSE)| method | input | proteins | FDR_lt_0.1 | model |
|---|---|---|---|---|
| binomial_firth_nested | nested | 20 | 0 | binomial_nested with Firth bias reduction |
| binomial_nested | nested | 20 | 0 | quasibinomial model, peptide detection counts |
| deqms | same | 20 | 14 | limma with DEqMS peptide-count variance moderation |
| deqms_voom | same | 20 | 13 | deqms with vooma precision weights |
| firth | same | 20 | 0 | Firth logistic regression, protein detection |
| firth_nested | nested | 20 | 0 | Firth logistic regression, peptide detection |
| hurdle_nested | nested | 20 | 14 | lm_impute intensity + binomial_firth_nested detection |
| limma | same | 20 | 14 | limma, empirical-Bayes moderated variances |
| limma_impute | same | 20 | 14 | limma with LOD imputation |
| limma_voom | same | 20 | 14 | limma with vooma precision weights |
| limma_voom_impute | same | 20 | 14 | limma_voom with LOD imputation |
| limpa | same | 20 | 16 | limma on limpa DPC-quantified proteins |
| limpa_nested | nested | 20 | 16 | limpa, DPC quantification from peptides |
| lm | same | 20 | 14 | linear model per protein |
| lm_impute | same | 20 | 14 | lm; groups without observations imputed at the limit of detection (LOD) |
| lmer_nested | nested | 12 | 5 | mixed model, peptide and sample random effects |
| rfit | same | 20 | 14 | rank-based regression (Rfit) |
| rfit_impute | same | 20 | 14 | rfit with LOD imputation |
| rlm | same | 20 | 14 | robust linear model (MASS::rlm) |
| ropeca_nested | nested | 20 | 14 | peptide fold changes summarised per protein (ROPECA) |
The detection models (firth, firth_nested, binomial_nested, binomial_firth_nested) test differences in how often a protein or its peptides are observed. The simulated data has few missing values, so they find little.
5 One protein, all methods
All methods return the same columns. Here is the protein with the smallest lm p-value. For the detection models diff is a log odds ratio, not a log2 fold change. This protein is detected in every sample, so the detection models report diff 0 and a p-value of 1 or NaN.
top <- results$lm$protein_Id[which.min(results$lm$p.value)]
one <- lapply(results, function(x) x[x$protein_Id == top, c("diff", "p.value", "FDR")])
one <- dplyr::bind_rows(one, .id = "method")
knitr::kable(one, digits = 3)| method | diff | p.value | FDR |
|---|---|---|---|
| binomial_firth_nested | 0.000 | 1.000 | 1.000 |
| binomial_nested | 0.000 | 1.000 | 1.000 |
| deqms | -0.469 | 0.000 | 0.000 |
| deqms_voom | -0.469 | 0.000 | 0.000 |
| firth | 0.000 | 1.000 | 1.000 |
| firth_nested | 0.000 | 1.000 | 1.000 |
| hurdle_nested | -0.469 | 0.000 | 0.000 |
| limma | -0.469 | 0.000 | 0.000 |
| limma_impute | -0.469 | 0.000 | 0.000 |
| limma_voom | -0.469 | 0.000 | 0.000 |
| limma_voom_impute | -0.469 | 0.000 | 0.000 |
| limpa | -0.272 | 0.000 | 0.000 |
| limpa_nested | -0.272 | 0.000 | 0.000 |
| lm | -0.469 | 0.000 | 0.000 |
| lm_impute | -0.469 | 0.000 | 0.000 |
| lmer_nested | -0.276 | 0.001 | 0.003 |
| rfit | -0.463 | 0.000 | 0.000 |
| rfit_impute | -0.463 | 0.000 | 0.000 |
| rlm | -0.469 | 0.000 | 0.000 |
| ropeca_nested | -0.495 | 0.000 | 0.000 |
6 Session info
sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.5 LTS
#>
#> Matrix products: default
#> BLAS: /home/biocbuild/bbs-3.24-bioc/R/lib/libRblas.so
#> LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0 LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
#> [3] LC_TIME=en_GB LC_COLLATE=C
#> [5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
#> [7] LC_PAPER=en_US.UTF-8 LC_NAME=C
#> [9] LC_ADDRESS=C LC_TELEPHONE=C
#> [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: America/New_York
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] ggplot2_4.0.3 prolfqua_0.99.8 dplyr_1.2.1 BiocStyle_2.41.0
#>
#> loaded via a namespace (and not attached):
#> [1] RColorBrewer_1.1-3 jsonlite_2.0.0 shape_1.4.6.1
#> [4] magrittr_2.0.5 TH.data_1.1-5 estimability_2.0.0
#> [7] jomo_2.7-6 magick_2.9.1 farver_2.1.2
#> [10] logistf_1.26.1 nloptr_2.2.1 rmarkdown_2.32
#> [13] GlobalOptions_0.1.4 vctrs_0.7.3 Cairo_1.7-0
#> [16] minqa_1.2.8 tinytex_0.61 htmltools_0.5.9
#> [19] forcats_1.0.1 progress_1.2.3 broom_1.0.13
#> [22] mitml_0.4-5 sass_0.4.10 KernSmooth_2.23-27
#> [25] bslib_0.12.0 htmlwidgets_1.6.4 plyr_1.8.9
#> [28] sandwich_3.1-3 emmeans_2.0.4 plotly_4.12.1
#> [31] zoo_1.9-1 cachem_1.1.0 mime_0.13
#> [34] lifecycle_1.0.5 iterators_1.0.14 pkgconfig_2.0.3
#> [37] Matrix_1.7-6 R6_2.6.1 fastmap_1.2.0
#> [40] rbibutils_2.4.1 shiny_1.14.0 clue_0.3-68
#> [43] digest_0.6.39 numDeriv_2016.8-1.1 colorspace_2.1-4
#> [46] S4Vectors_0.51.10 crosstalk_1.2.2 labeling_0.4.3
#> [49] httr_1.4.9 mgcv_1.9-5 compiler_4.6.1
#> [52] withr_3.0.3 doParallel_1.0.17 S7_0.2.2
#> [55] backports_1.5.1 UpSetR_1.4.1 pan_2.0
#> [58] MASS_7.3-66 rjson_0.2.23 tools_4.6.1
#> [61] otel_0.2.0 httpuv_1.6.17 nnet_7.3-21
#> [64] glue_1.8.1 nlme_3.1-172 promises_1.5.0
#> [67] grid_4.6.1 cluster_2.1.8.3 generics_0.1.4
#> [70] operator.tools_1.6.3.1 gtable_0.3.6 formula.tools_1.7.1
#> [73] tidyr_1.3.2 data.table_1.18.6.1 hms_1.1.4
#> [76] utf8_1.2.6 BiocGenerics_0.59.12 ggrepel_0.9.8
#> [79] foreach_1.5.2 pillar_1.11.1 limma_3.99.4
#> [82] later_1.4.8 circlize_0.4.18 splines_4.6.1
#> [85] lattice_0.23-1 survival_3.8-12 tidyselect_1.2.1
#> [88] ComplexHeatmap_2.29.0 knitr_1.52 reformulas_0.4.4
#> [91] gridExtra_2.3.1 bookdown_0.48 IRanges_2.47.5
#> [94] stats4_4.6.1 xfun_0.62 statmod_1.5.2
#> [97] matrixStats_1.5.0 stringi_1.8.9 yaml_2.3.12
#> [100] boot_1.3-32 evaluate_1.0.5 codetools_0.2-20
#> [103] Rfit_0.27.0 tibble_3.3.1 BiocManager_1.30.27
#> [106] cli_3.6.6 rpart_4.1.27 xtable_1.8-8
#> [109] Rdpack_2.6.6 jquerylib_0.1.4 dichromat_2.0-1
#> [112] Rcpp_1.1.2 coda_0.19-4.1 png_0.1-9
#> [115] parallel_4.6.1 prettyunits_1.2.0 lme4_2.1-0
#> [118] glmnet_5.2 viridisLite_0.4.3 limpa_1.5.2
#> [121] mvtnorm_1.4-2 lmerTest_3.2-1 scales_1.4.0
#> [124] purrr_1.2.2 crayon_1.5.3 GetoptLong_1.1.1
#> [127] rlang_1.3.0 multcomp_1.4-32 mice_3.19.0