This vignette compares fastglm with reference implementations across sample sizes and model classes:

library(fastglm)
suppressPackageStartupMessages({
    library(MASS)
    library(pscl)
    library(logistf)
    library(brglm2)
    library(speedglm)
    library(Matrix)
    library(bigmemory)
    library(microbenchmark)
})

bench_unit <- "milliseconds"
.fmt <- function(mb) {
    s <- summary(mb, unit = bench_unit)
    s <- s[, c("expr", "min", "median", "mean", "max")]
    names(s) <- c("expr", "min (ms)", "median (ms)", "mean (ms)", "max (ms)")
    s
}
# Summary of one comparison, as a table.
.fmt_table <- function(mb, caption = NULL) {
    knitr::kable(.fmt(mb), digits = 2, row.names = FALSE, caption = caption)
}
# Median times with one row per method and one column per sample size.
.timing_table <- function(tim) {
    methods <- unique(tim$expr)
    sizes <- sort(unique(tim$n))
    medians <- tapply(tim$median, list(tim$expr, tim$n), identity)[methods, as.character(sizes), drop = FALSE]
    labels <- paste0("n = ", format(sizes, big.mark = ",", scientific = FALSE, trim = TRUE))
    knitr::kable(data.frame(Method = methods, medians, check.names = FALSE, row.names = NULL),
                 col.names = c("Method", labels), digits = 2,
                 caption = "Median time in milliseconds")
}
# Paired medians for fastglm and brglm2, by family and sample size.
.pair_table <- function(tim) {
    keys <- unique(tim[, c("family", "n")])
    median_of <- function(method, i) tim$median[tim$expr == method & tim$family == keys$family[i] & tim$n == keys$n[i]]
    fast <- vapply(seq_len(nrow(keys)), function(i) median_of("fastglm", i), 0)
    reference <- vapply(seq_len(nrow(keys)), function(i) median_of("brglm2", i), 0)
    knitr::kable(data.frame(Family = keys$family, n = format(keys$n, big.mark = ",", scientific = FALSE, trim = TRUE),
                            fastglm = fast, brglm2 = reference, Ratio = reference / fast),
                 col.names = c("Family", "n", "fastglm (ms)", "brglm2 (ms)", "Time ratio"),
                 digits = c(0, 0, 2, 2, 1), row.names = FALSE)
}
# Short scientific notation for numbers placed in the text.
.sci <- function(x) formatC(x, format = "e", digits = 1)
.collect <- function(mb, n) {
    s <- summary(mb, unit = bench_unit)
    data.frame(n = n, expr = as.character(s$expr), median = s$median,
               stringsAsFactors = FALSE)
}
# One style for every figure: fastglm is dominant, its solver variants share a
# color and differ by linetype, and comparison methods stay visible but thinner.
library(ggplot2)
.method_colors <- c(fastglm = "#0B5CAD", fastglm_qr = "#0B5CAD", fastglm_chol = "#0B5CAD",
    fastglm_nb = "#0B5CAD", fastglm_hurdle = "#0B5CAD", fastglm_zi = "#0B5CAD",
    glm.fit = "#D1661F", glm.nb = "#D1661F", brglm2 = "#D1661F", pscl_hurdle = "#D1661F",
    pscl_zi = "#D1661F", speedglm = "#1F9E8F", glm2_fit2 = "#9A62C4")
.method_lines <- function(methods) ifelse(methods == "fastglm_qr", "dashed", "solid")
.method_widths <- function(methods) ifelse(grepl("^fastglm", methods), 1.2, 0.6)
.benchmark_theme <- function() {
    theme_bw(base_size = 12) +
        theme(panel.grid.minor = element_blank(), legend.position = "bottom",
              strip.background = element_blank(), strip.text = element_text(hjust = 0))
}
.size_axis <- function() {
    scale_x_log10(labels = function(x) format(x, big.mark = ",", scientific = FALSE, trim = TRUE))
}
.plot_scaling <- function(df, title, facet = NULL) {
    methods <- unique(as.character(df$expr))
    plot <- ggplot(df, aes(n, median, colour = expr, linetype = expr, linewidth = expr, group = expr)) +
        geom_line() + geom_point(size = 2.2) +
        .size_axis() + scale_y_log10() +
        scale_colour_manual(values = .method_colors[methods], breaks = methods) +
        scale_linetype_manual(values = setNames(.method_lines(methods), methods), breaks = methods) +
        scale_linewidth_manual(values = setNames(.method_widths(methods), methods), breaks = methods) +
        labs(x = "Number of observations, log scale",
             y = paste0("Median time (", bench_unit, "), log scale"),
             title = title, colour = NULL, linetype = NULL, linewidth = NULL) +
        .benchmark_theme()
    if (!is.null(facet)) plot <- plot + facet_wrap(facet, scales = "free_y")
    plot
}
sessionInfo()
#> R version 4.5.1 (2025-06-13)
#> Platform: aarch64-apple-darwin20
#> Running under: macOS Sequoia 15.7.3
#> 
#> Matrix products: default
#> BLAS:   /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib 
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
#> 
#> locale:
#> [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
#> 
#> time zone: America/Chicago
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#>  [1] ggplot2_4.0.1        microbenchmark_1.5.0 bigmemory_4.6.4     
#>  [4] speedglm_0.3-5       biglm_0.9-3          DBI_1.2.3           
#>  [7] Matrix_1.7-5         brglm2_1.1.1         logistf_1.26.1      
#> [10] pscl_1.5.9           MASS_7.3-66          fastglm_0.1.3       
#> 
#> loaded via a namespace (and not attached):
#>  [1] gtable_0.3.6           shape_1.4.6.1          xfun_0.57             
#>  [4] formula.tools_1.7.1    lattice_0.22-7         numDeriv_2016.8-1.1   
#>  [7] vctrs_0.7.2            tools_4.5.1            Rdpack_2.6.4          
#> [10] generics_0.1.4         enrichwith_0.5.0       sandwich_3.1-1        
#> [13] tibble_3.3.0           pan_1.9                pkgconfig_2.0.3       
#> [16] jomo_2.7-6             RColorBrewer_1.1-3     S7_0.2.1              
#> [19] uuid_1.2-1             lifecycle_1.0.5        farver_2.1.2          
#> [22] compiler_4.5.1         statmod_1.5.1          codetools_0.2-20      
#> [25] glmnet_4.1-10          mice_3.18.0            pillar_1.11.1         
#> [28] nloptr_2.2.1           tidyr_1.3.1            reformulas_0.4.1      
#> [31] iterators_1.0.14       rpart_4.1.24           boot_1.3-31           
#> [34] multcomp_1.4-28        foreach_1.5.2          mitml_0.4-5           
#> [37] nlme_3.1-168           tidyselect_1.2.1       mvtnorm_1.3-7         
#> [40] dplyr_1.1.4            purrr_1.2.1            splines_4.5.1         
#> [43] operator.tools_1.6.3.1 grid_4.5.1             cli_3.6.6             
#> [46] magrittr_2.0.4         survival_3.8-3         broom_1.0.12          
#> [49] TH.data_1.1-3          bigmemory.sri_0.1.8    withr_3.0.2           
#> [52] scales_1.4.0           backports_1.5.0        estimability_1.5.1    
#> [55] emmeans_1.11.2         nnet_7.3-20            lme4_1.1-37           
#> [58] zoo_1.8-14             evaluate_1.0.5         knitr_1.51            
#> [61] rbibutils_2.3          mgcv_1.9-3             nleqslv_3.3.5         
#> [64] rlang_1.3.0            Rcpp_1.1.1-1.1         xtable_1.8-4          
#> [67] glue_1.8.0             minqa_1.2.8            R6_2.6.1

Standard GLMs

Logistic regression with an increasing number of observations, holding p = 30 fixed. Comparing the default fastglm path (method = 0, pivoted QR), the LLT-Cholesky path (method = 2), stats::glm.fit(), glm2::glm.fit2(), and speedglm::speedglm.wfit().

set.seed(1)
p <- 30

run_one_logit <- function(n) {
    X <- matrix(rnorm(n * p), n, p)
    beta <- rnorm(p) * 0.05
    y <- rbinom(n, 1, plogis(X %*% beta))
    microbenchmark(
        fastglm_qr   = fastglm(X, y, family = binomial(),            method = 0),
        fastglm_chol = fastglm(X, y, family = binomial(),            method = 2),
        glm.fit      = glm.fit(X, y, family = binomial()),
        glm2_fit2    = glm2::glm.fit2(X, y, family = binomial()),
        speedglm     = speedglm::speedglm.wfit(y = y, X = X,
                                               family = binomial(),
                                               intercept = FALSE),
        times = 5L
    )
}

logit_tim <- list()
for (n in c(1e3, 1e4, 1e5)) {
    mb <- run_one_logit(n)
    logit_tim[[length(logit_tim) + 1]] <- .collect(mb, n)
}
logit_tim <- do.call(rbind, logit_tim)
.timing_table(logit_tim)

Median time in milliseconds
Method	n = 1,000	n = 10,000	n = 100,000
fastglm_qr	1.26	11.30	120.40
fastglm_chol	0.67	4.66	47.89
glm.fit	2.75	29.02	343.30
glm2_fit2	2.89	29.29	402.48
speedglm	2.03	20.63	255.29

.plot_scaling(logit_tim, "Logistic regression: median time vs n (p = 30)")

$plot of chunk unnamed-chunk-2$

plot of chunk unnamed-chunk-2

The same comparison for Poisson regression with a log link:

run_one_poisson <- function(n) {
    X <- matrix(rnorm(n * p), n, p)
    beta <- rnorm(p) * 0.05
    y <- rpois(n, exp(X %*% beta + 1))
    microbenchmark(
        fastglm_qr   = fastglm(X, y, family = poisson(),            method = 0),
        fastglm_chol = fastglm(X, y, family = poisson(),            method = 2),
        glm.fit      = glm.fit(X, y, family = poisson()),
        glm2_fit2    = glm2::glm.fit2(X, y, family = poisson()),
        speedglm     = speedglm::speedglm.wfit(y = y, X = X,
                                               family = poisson(),
                                               intercept = FALSE),
        times = 5L
    )
}

pois_tim <- list()
for (n in c(1e3, 1e4, 1e5)) {
    mb <- run_one_poisson(n)
    pois_tim[[length(pois_tim) + 1]] <- .collect(mb, n)
}
pois_tim <- do.call(rbind, pois_tim)
.timing_table(pois_tim)

Median time in milliseconds
Method	n = 1,000	n = 10,000	n = 100,000
fastglm_qr	1.21	10.93	114.10
fastglm_chol	0.63	4.27	43.23
glm.fit	2.79	27.77	353.56
glm2_fit2	2.89	28.33	358.98
speedglm	2.13	19.80	205.14

.plot_scaling(pois_tim, "Poisson regression: median time vs n (p = 30)")

plot of chunk unnamed-chunk-3

These timings compare two ways to solve the weighted regression in each fitting iteration: QR uses the predictor matrix directly, whereas Cholesky works with its smaller crossproduct. Their costs depend on the numbers of rows and columns, the response family, and the numerical libraries used. The observed ratios apply to the displayed cases; attributing a difference specifically to parallel computation would require a separate comparison with threading controlled.

Sparse and big.matrix paths

A sparse design matrix from one-hot encoding a 100-level factor against a continuous covariate. The column count is p = 200; fastglm on the dense matrix stores all 200 columns, while the sparse path operates only on the nonzero entries:

set.seed(2)
n   <- 5e4
g   <- factor(sample.int(100, n, replace = TRUE))
xn  <- rnorm(n)
Xd  <- model.matrix(~ g * xn)
Xs  <- as(Xd, "CsparseMatrix")
betas <- rnorm(ncol(Xd)) / 4
y   <- rbinom(n, 1, plogis(Xd %*% betas))

mb_sparse <- microbenchmark(
    fastglm_dense  = fastglm(Xd, y, family = binomial(), method = 2),
    fastglm_sparse = fastglm(Xs, y, family = binomial(), method = 2),
    times = 5L
)
.fmt_table(mb_sparse)

expr	min (ms)	median (ms)	mean (ms)	max (ms)
fastglm_dense	283.81	285.47	285.93	288.70
fastglm_sparse	22.06	23.94	24.13	27.88

The design has 200 columns, and 2% of its entries are nonzero.

A big.matrix comparison on a large dense design, with n = 5e5 and p = 20. The file-backed big.matrix reads blocks of rows and gives the same estimates:

set.seed(3)
n  <- 5e5
p  <- 20
X  <- matrix(rnorm(n * p), n, p)
Xb <- as.big.matrix(X)
y  <- rbinom(n, 1, plogis(X %*% rnorm(p) * 0.05))

mb_big <- microbenchmark(
    dense       = fastglm(X,  y, family = binomial(), method = 2),
    big.matrix  = fastglm(Xb, y, family = binomial(), method = 2),
    times = 3L
)
.fmt_table(mb_big)

expr	min (ms)	median (ms)	mean (ms)	max (ms)
dense	171.44	173.77	199.10	252.08
big.matrix	214.23	218.48	217.31	219.22


f1 <- fastglm(X,  y, family = binomial(), method = 2)
f2 <- fastglm(Xb, y, family = binomial(), method = 2)

The largest coefficient difference between the two fits is 1.2e-15. Reading a file-backed matrix in blocks adds overhead but can reduce the need to copy a large predictor matrix. Another option is to supply blocks through fastglm_source() and fit them with fastglm_streaming(), as described in the large-data vignette. These choices have different memory requirements and should be compared on the intended data size.

Negative-binomial regression

fastglm_nb() against MASS::glm.nb() across a range of sample sizes. The data-generating process is NB(mu, theta = 2) with three covariates and an intercept:

set.seed(4)
run_one_nb <- function(n) {
    X  <- cbind(1, matrix(rnorm(n * 3), n, 3))
    mu <- exp(X %*% c(0.5, 0.4, -0.2, 0.3))
    y  <- MASS::rnegbin(n, mu = mu, theta = 2)
    df <- data.frame(y = y, x1 = X[, 2], x2 = X[, 3], x3 = X[, 4])
    microbenchmark(
        fastglm_nb = fastglm_nb(X, y),
        glm.nb     = MASS::glm.nb(y ~ x1 + x2 + x3, data = df),
        times = 5L
    )
}

nb_tim <- list()
for (n in c(1e3, 1e4, 1e5, 5e5)) {
    mb <- run_one_nb(n)
    nb_tim[[length(nb_tim) + 1]] <- .collect(mb, n)
}
nb_tim <- do.call(rbind, nb_tim)
.timing_table(nb_tim)

Median time in milliseconds
Method	n = 1,000	n = 10,000	n = 100,000	n = 500,000
fastglm_nb	1.62	13.21	111.02	577.99
glm.nb	7.71	63.77	540.85	2215.27

.plot_scaling(nb_tim, "Negative-binomial joint MLE: median time vs n")

plot of chunk unnamed-chunk-6

For the largest case, both functions target the same maximum likelihood fit, and numerical agreement also depends on convergence and stopping tolerances. Dense fastglm fits perform both coefficient and theta updates in C++, which reduces repeated calls between R and compiled code.

set.seed(99)
n  <- 5e5
X  <- cbind(1, matrix(rnorm(n * 3), n, 3))
mu <- exp(X %*% c(0.5, 0.4, -0.2, 0.3))
y  <- MASS::rnegbin(n, mu = mu, theta = 2)
df <- data.frame(y = y, x1 = X[, 2], x2 = X[, 3], x3 = X[, 4])
f_nb <- fastglm_nb(X, y)
m_nb <- MASS::glm.nb(y ~ x1 + x2 + x3, data = df)

The largest coefficient difference is 2.2e-14, the difference in theta is 9.7e-12, and the difference in log likelihood is 2.3e-09.

Firth bias-reduced GLMs

Firth bias reduction (Firth, 1993; Kosmidis and Firth, 2009) is supported for all standard GLM families via firth = TRUE. We benchmark against brglm2::brglmFit(type = "AS_mean"), which implements the same AS_mean adjustment, and against logistf::logistf() for the binomial-logit special case.

First, the canonical logistic benchmark on the sex2 dataset and a larger simulated example:

data(sex2, package = "logistf")
X_sex2 <- model.matrix(~ age + oc + vic + vicl + vis + dia, data = sex2)
y_sex2 <- sex2$case

mb_firth_logistf <- microbenchmark(
    fastglm = fastglm(X_sex2, y_sex2, family = binomial(), firth = TRUE),
    logistf = logistf(case ~ age + oc + vic + vicl + vis + dia, data = sex2),
    brglm2  = glm(case ~ age + oc + vic + vicl + vis + dia, data = sex2,
                  family = binomial(), method = "brglmFit", type = "AS_mean"),
    times = 10L
)
.fmt_table(mb_firth_logistf, caption = paste0("Heinze-Schemper sex2 data, n = ", nrow(sex2)))

Heinze-Schemper sex2 data, n = 239
expr	min (ms)	median (ms)	mean (ms)	max (ms)
fastglm	0.63	0.66	0.68	0.82
logistf	7.18	7.30	7.38	8.26
brglm2	2.57	2.65	2.94	5.52

Next, a comprehensive comparison of fastglm(firth = TRUE) against brglm2 across all supported families and links. We fix n = 10000 and p = 5 to see the per-family speedup:

set.seed(123)
n <- 10000;  p_firth <- 5
X_firth <- cbind(1, matrix(rnorm(n * p_firth), n, p_firth))
eta_firth <- X_firth %*% c(0.5, 0.3, -0.2, 0.1, 0.05, -0.1)

firth_families <- list(
    list(fam = binomial("logit"),       label = "binomial logit",
         y = rbinom(n, 1, plogis(eta_firth))),
    list(fam = binomial("probit"),      label = "binomial probit",
         y = rbinom(n, 1, plogis(eta_firth))),
    list(fam = binomial("cloglog"),     label = "binomial cloglog",
         y = rbinom(n, 1, plogis(eta_firth))),
    list(fam = poisson("log"),          label = "poisson log",
         y = rpois(n, exp(eta_firth))),
    list(fam = poisson("sqrt"),         label = "poisson sqrt",
         y = rpois(n, exp(eta_firth))),
    list(fam = Gamma("log"),            label = "Gamma log",
         y = rgamma(n, shape = 5, rate = 5 / exp(eta_firth))),
    list(fam = gaussian("identity"),    label = "gaussian identity",
         y = as.numeric(eta_firth) + rnorm(n, 0, 0.5)),
    list(fam = gaussian("log"),         label = "gaussian log",
         y = pmax(as.numeric(eta_firth) + rnorm(n, 0, 0.5), 0.01))
)

xnames <- paste0("x", seq_len(p_firth))
df_firth <- data.frame(X_firth[, -1])
names(df_firth) <- xnames
fml_firth <- as.formula(paste("y ~", paste(xnames, collapse = " + ")))

firth_all_tim <- list()
for (fi in firth_families) {
    df_firth$y <- fi$y
    mb <- microbenchmark(
        fastglm = fastglm(X_firth, fi$y, family = fi$fam, firth = TRUE),
        brglm2  = glm(fml_firth, data = df_firth, family = fi$fam,
                      method = "brglmFit", type = "AS_mean"),
        times = 10L
    )
    firth_all_tim[[length(firth_all_tim) + 1]] <-
        transform(.collect(mb, n), family = fi$label)
}
firth_all_tim <- do.call(rbind, firth_all_tim)
.pair_table(firth_all_tim)

Family	n	fastglm (ms)	brglm2 (ms)	Time ratio
binomial logit	10,000	4.76	14.17	3.0
binomial probit	10,000	7.43	15.68	2.1
binomial cloglog	10,000	8.22	17.39	2.1
poisson log	10,000	5.88	123.80	21.0
poisson sqrt	10,000	3.04	123.02	40.5
Gamma log	10,000	3.37	105.35	31.2
gaussian identity	10,000	1.46	21.27	14.6
gaussian log	10,000	7.40	29.01	3.9

The median time for each family, with each pair of methods joined by a line:

firth_all_tim$family <- reorder(firth_all_tim$family, firth_all_tim$median, max)
ggplot(firth_all_tim, aes(median, family, colour = expr, group = family)) +
    geom_line(colour = "grey70") + geom_point(size = 2.4) +
    scale_x_log10() +
    scale_colour_manual(values = .method_colors[c("fastglm", "brglm2")], breaks = c("fastglm", "brglm2")) +
    labs(x = paste0("Median time (", bench_unit, "), log scale"), y = NULL, colour = NULL,
         title = paste0("Firth bias-reduced GLMs: fastglm and brglm2 (n = ", n, ")")) +
    .benchmark_theme()

plot of chunk unnamed-chunk-10

Scaling with sample size for a representative subset of families:

set.seed(123)
firth_scaling_fams <- list(
    list(fam = binomial("logit"),    label = "binomial logit"),
    list(fam = poisson("log"),       label = "poisson log"),
    list(fam = Gamma("log"),         label = "Gamma log"),
    list(fam = gaussian("identity"), label = "gaussian identity")
)

firth_scale_tim <- list()
for (n_val in c(1e3, 1e4, 5e4)) {
    X_s <- cbind(1, matrix(rnorm(n_val * p_firth), n_val, p_firth))
    eta_s <- X_s %*% c(0.5, 0.3, -0.2, 0.1, 0.05, -0.1)
    df_s <- data.frame(X_s[, -1])
    names(df_s) <- xnames

    for (fi in firth_scaling_fams) {
        y_s <- switch(fi$fam$family,
            "binomial" = rbinom(n_val, 1, plogis(eta_s)),
            "poisson"  = rpois(n_val, exp(eta_s)),
            "Gamma"    = rgamma(n_val, shape = 5, rate = 5 / exp(eta_s)),
            "gaussian" = as.numeric(eta_s) + rnorm(n_val, 0, 0.5)
        )
        df_s$y <- y_s
        mb <- microbenchmark(
            fastglm = fastglm(X_s, y_s, family = fi$fam, firth = TRUE),
            brglm2  = glm(fml_firth, data = df_s, family = fi$fam,
                          method = "brglmFit", type = "AS_mean"),
            times = 5L
        )
        firth_scale_tim[[length(firth_scale_tim) + 1]] <-
            transform(.collect(mb, n_val), family = fi$label)
    }
}
firth_scale_tim <- do.call(rbind, firth_scale_tim)
.pair_table(firth_scale_tim)

Family	n	fastglm (ms)	brglm2 (ms)	Time ratio
binomial logit	1,000	1.00	2.45	2.5
poisson log	1,000	0.93	13.25	14.2
Gamma log	1,000	0.71	12.19	17.2
gaussian identity	1,000	1.15	3.70	3.2
binomial logit	10,000	5.30	14.39	2.7
poisson log	10,000	5.76	123.61	21.4
Gamma log	10,000	3.52	99.90	28.4
gaussian identity	10,000	1.40	21.82	15.6
binomial logit	50,000	21.92	69.31	3.2
poisson log	50,000	22.24	605.09	27.2
Gamma log	50,000	16.00	517.47	32.3
gaussian identity	50,000	5.40	109.83	20.4


.plot_scaling(firth_scale_tim, "Firth bias reduction: median time by family", facet = ~family)

plot of chunk unnamed-chunk-11

Firth across all decomposition methods

Firth bias reduction works with all six decomposition methods. This example checks that they give the same coefficients and compares their timings:

set.seed(123)
n_meth <- 10000;  p_meth <- 10
X_meth <- cbind(1, matrix(rnorm(n_meth * (p_meth - 1)), n_meth, p_meth - 1))
y_meth <- rbinom(n_meth, 1, plogis(X_meth %*% rnorm(p_meth) * 0.1))

ref_coef <- coef(fastglm(X_meth, y_meth, family = binomial(), method = 2, firth = TRUE))
method_names <- c("ColPivQR", "QR", "LLT", "LDLT", "FullPivQR", "SVD")
compared_methods <- c(0L, 1L, 3L, 4L, 5L)
method_differences <- vapply(compared_methods, function(m) {
    max(abs(coef(fastglm(X_meth, y_meth, family = binomial(), method = m, firth = TRUE)) - ref_coef))
}, 0)
knitr::kable(data.frame(Method = method_names[compared_methods + 1], Difference = .sci(method_differences)),
             col.names = c("Method", "Largest coefficient difference from LLT"), row.names = FALSE)

Method	Largest coefficient difference from LLT
ColPivQR	5.8e-16
QR	6.9e-16
LDLT	8.3e-17
FullPivQR	6.7e-16
SVD	3.5e-16


mb_methods <- microbenchmark(
    ColPivQR = fastglm(X_meth, y_meth, family = binomial(), method = 0, firth = TRUE),
    QR       = fastglm(X_meth, y_meth, family = binomial(), method = 1, firth = TRUE),
    LLT      = fastglm(X_meth, y_meth, family = binomial(), method = 2, firth = TRUE),
    LDLT     = fastglm(X_meth, y_meth, family = binomial(), method = 3, firth = TRUE),
    FullPivQR= fastglm(X_meth, y_meth, family = binomial(), method = 4, firth = TRUE),
    SVD      = fastglm(X_meth, y_meth, family = binomial(), method = 5, firth = TRUE),
    times = 5L
)
.fmt_table(mb_methods, caption = paste0("Firth logistic regression, n = ", n_meth, ", p = ", p_meth))

Firth logistic regression, n = 10000, p = 10
expr	min (ms)	median (ms)	mean (ms)	max (ms)
ColPivQR	6.89	8.73	8.33	9.49
QR	6.82	7.35	7.65	8.71
LLT	4.88	5.28	5.83	7.23
LDLT	5.74	7.30	7.00	8.10
FullPivQR	3020.98	3033.14	3052.33	3107.26
SVD	9.71	11.24	11.05	12.71

All six methods agree to machine precision. The LLT and LDLT Cholesky methods are the fastest because they only decompose the p x p cross-product, while the QR and SVD methods work with the full n x p matrix.

Firth on sparse and streaming designs

This example compares several ways to supply the predictor matrix: ordinary storage, sparse storage of nonzero entries, a big.matrix, and a function that reads blocks of rows as needed. Bias reduction uses leverage values, which describe the relative influence of predictor rows in a weighted regression. Some block-based calculations reuse these values from the preceding coefficient update. Agreement with the direct calculation therefore needs to be checked after convergence.

set.seed(123)
n_sp <- 5000;  p_sp <- 6
X_sp <- cbind(1, matrix(rnorm(n_sp * (p_sp - 1)), n_sp, p_sp - 1))
Xs_sp <- as(X_sp, "CsparseMatrix")
Xb_sp <- as.big.matrix(X_sp)
y_sp <- rbinom(n_sp, 1, plogis(X_sp %*% c(-0.3, 0.4, -0.2, 0.1, 0.05, -0.15)))

chunk_size_sp <- 500
chunks_sp <- function(k) {
    idx <- ((k - 1) * chunk_size_sp + 1):(k * chunk_size_sp)
    list(X = X_sp[idx, , drop = FALSE], y = y_sp[idx])
}

mb_backends <- microbenchmark(
    dense      = fastglm(X_sp,  y_sp, family = binomial(), method = 2, firth = TRUE),
    sparse     = fastglm(Xs_sp, y_sp, family = binomial(), method = 2, firth = TRUE),
    big.matrix = fastglm(Xb_sp, y_sp, family = binomial(), method = 2, firth = TRUE),
    streaming  = fastglm_streaming(chunks_sp, n_chunks = n_sp / chunk_size_sp,
                                   family = binomial(), method = 2, firth = TRUE),
    times = 5L
)
.fmt_table(mb_backends, caption = paste0("Firth logistic regression, n = ", n_sp, ", p = ", p_sp))

Firth logistic regression, n = 5000, p = 6
expr	min (ms)	median (ms)	mean (ms)	max (ms)
dense	2.56	2.63	2.63	2.75
sparse	6.29	6.70	6.88	7.72
big.matrix	2.46	3.27	3.14	3.44
streaming	5.01	5.05	5.12	5.37


# Verify agreement
ref_sp <- coef(fastglm(X_sp, y_sp, family = binomial(), method = 2, firth = TRUE))
sp_c   <- coef(fastglm(Xs_sp, y_sp, family = binomial(), method = 2, firth = TRUE))
bm_c   <- coef(fastglm(Xb_sp, y_sp, family = binomial(), method = 2, firth = TRUE))
st_c   <- coef(fastglm_streaming(chunks_sp, n_chunks = n_sp / chunk_size_sp,
                                  family = binomial(), method = 2, firth = TRUE))
knitr::kable(data.frame(Representation = c("sparse", "big.matrix", "streaming"),
                        Difference = .sci(c(max(abs(sp_c - ref_sp)), max(abs(bm_c - ref_sp)), max(abs(st_c - ref_sp))))),
             col.names = c("Representation", "Largest coefficient difference from dense"), row.names = FALSE)

Representation	Largest coefficient difference from dense
sparse	2.1e-10
big.matrix	2.1e-10
streaming	1.1e-10

In this logistic example, the reported coefficient differences across data representations are below $10^{-7}$. Memory and runtime comparisons remain specific to the displayed design. For other families, the mean-bias adjustment is defined by adjusted estimating equations and need not correspond to a penalized likelihood.

Hurdle models

fastglm_hurdle against pscl::hurdle across a range of sample sizes, Poisson count component:

set.seed(6)
run_one_hurdle <- function(n) {
    x1 <- rnorm(n);  x2 <- rnorm(n)
    lam    <- exp(0.7 + 0.4 * x1 - 0.3 * x2)
    is_pos <- rbinom(n, 1, plogis(-0.4 + 0.5 * x1 + 0.2 * x2))
    yt     <- integer(n)
    for (i in seq_len(n)) {
        repeat { v <- rpois(1, lam[i]); if (v > 0) { yt[i] <- v; break } }
    }
    y  <- ifelse(is_pos == 1, yt, 0L)
    df <- data.frame(y = y, x1 = x1, x2 = x2)
    microbenchmark(
        fastglm_hurdle = fastglm_hurdle(y ~ x1 + x2, data = df, dist = "poisson"),
        pscl_hurdle    = pscl::hurdle (y ~ x1 + x2, data = df, dist = "poisson"),
        times = 5L
    )
}

hurdle_tim <- list()
for (n in c(1e3, 1e4, 5e4, 1e5)) {
    mb <- run_one_hurdle(n)
    hurdle_tim[[length(hurdle_tim) + 1]] <- .collect(mb, n)
}
hurdle_tim <- do.call(rbind, hurdle_tim)
.timing_table(hurdle_tim)

Median time in milliseconds
Method	n = 1,000	n = 10,000	n = 50,000	n = 100,000
fastglm_hurdle	0.61	3.55	18.22	37.49
pscl_hurdle	4.95	43.69	290.87	512.86

.plot_scaling(hurdle_tim, "Hurdle Poisson: median time vs n")

plot of chunk unnamed-chunk-14

NB hurdle (with simultaneous theta MLE) on a single moderately large example:

set.seed(7)
n <- 1e4
x1 <- rnorm(n);  x2 <- rnorm(n)
lam <- exp(0.8 + 0.4 * x1 - 0.3 * x2)
is_pos <- rbinom(n, 1, plogis(-0.4 + 0.5 * x1))
yt <- integer(n)
for (i in seq_len(n)) {
    repeat {
        v <- MASS::rnegbin(1, mu = lam[i], theta = 1.5)
        if (v > 0) { yt[i] <- v; break }
    }
}
y  <- ifelse(is_pos == 1, yt, 0L)
df <- data.frame(y = y, x1 = x1, x2 = x2)

mb_hurdle_nb <- microbenchmark(
    fastglm_hurdle = fastglm_hurdle(y ~ x1 + x2, data = df, dist = "negbin"),
    pscl_hurdle    = pscl::hurdle (y ~ x1 + x2, data = df, dist = "negbin"),
    times = 3L
)
.fmt_table(mb_hurdle_nb)

expr	min (ms)	median (ms)	mean (ms)	max (ms)
fastglm_hurdle	31.55	32.07	33.09	35.65
pscl_hurdle	50.60	50.97	55.02	63.49

Zero-inflated models

Zero-inflated Poisson against pscl::zeroinfl() across sample sizes:

set.seed(8)
run_one_zi <- function(n) {
    x1 <- rnorm(n);  x2 <- rnorm(n)
    eta_c <- 0.7 + 0.4 * x1 - 0.3 * x2
    eta_z <- -0.4 + 0.5 * x1 + 0.2 * x2
    z     <- rbinom(n, 1, plogis(eta_z))
    y     <- ifelse(z == 1, 0L, rpois(n, exp(eta_c)))
    df    <- data.frame(y = y, x1 = x1, x2 = x2)
    microbenchmark(
        fastglm_zi = fastglm_zi(y ~ x1 + x2, data = df, dist = "poisson"),
        pscl_zi    = pscl::zeroinfl(y ~ x1 + x2, data = df, dist = "poisson"),
        times = 5L
    )
}

zi_tim <- list()
for (n in c(1e3, 1e4, 5e4, 1e5)) {
    mb <- run_one_zi(n)
    zi_tim[[length(zi_tim) + 1]] <- .collect(mb, n)
}
zi_tim <- do.call(rbind, zi_tim)
.timing_table(zi_tim)

Median time in milliseconds
Method	n = 1,000	n = 10,000	n = 50,000	n = 100,000
fastglm_zi	2.47	25.25	125.92	238.41
pscl_zi	9.19	109.32	546.18	1087.10

.plot_scaling(zi_tim, "Zero-inflated Poisson: median time vs n")

plot of chunk unnamed-chunk-16

Zero-inflated NB on a single moderately large example:

set.seed(9)
n  <- 1e4
x1 <- rnorm(n);  x2 <- rnorm(n)
eta_c <- 0.8 + 0.4 * x1 - 0.3 * x2; lam <- exp(eta_c)
eta_z <- -0.4 + 0.5 * x1 + 0.2 * x2
z <- rbinom(n, 1, plogis(eta_z))
y <- ifelse(z == 1, 0L, MASS::rnegbin(n, mu = lam, theta = 2.0))
df <- data.frame(y = y, x1 = x1, x2 = x2)

mb_zi_nb <- microbenchmark(
    fastglm_zi = fastglm_zi(y ~ x1 + x2, data = df, dist = "negbin",
                            em.tol = 1e-8, em.maxit = 200L),
    pscl_zi    = pscl::zeroinfl(y ~ x1 + x2, data = df, dist = "negbin"),
    times = 3L
)
.fmt_table(mb_zi_nb)

expr	min (ms)	median (ms)	mean (ms)	max (ms)
fastglm_zi	73.99	78.61	80.56	89.08
pscl_zi	144.19	155.03	152.10	157.09

Summary

A compact summary of the speed advantage fastglm delivers, expressed as the ratio of reference-implementation median time to fastglm median time. Larger is better. The reference for each model class is the canonical R implementation; for the standard GLMs we report the ratio against the fastest among glm.fit, glm2, and speedglm so the comparison is conservative.

.speedup <- function(df, fast_label, ref_labels) {
    out <- lapply(split(df, df$n), function(sub) {
        fast <- sub$median[sub$expr == fast_label]
        if (length(fast) == 0) return(NULL)
        ref  <- min(sub$median[sub$expr %in% ref_labels])
        data.frame(n = sub$n[1], speedup = ref / fast)
    })
    do.call(rbind, out)
}

su_logit  <- .speedup(logit_tim,  "fastglm_chol",
                      c("glm.fit", "glm2_fit2", "speedglm"))
su_pois   <- .speedup(pois_tim,   "fastglm_chol",
                      c("glm.fit", "glm2_fit2", "speedglm"))
su_nb     <- .speedup(nb_tim,     "fastglm_nb",     "glm.nb")
firth_avg_tim <- aggregate(median ~ n + expr, data = firth_scale_tim, FUN = mean)
su_firth  <- .speedup(firth_avg_tim, "fastglm",  "brglm2")
su_hurdle <- .speedup(hurdle_tim, "fastglm_hurdle", "pscl_hurdle")
su_zi     <- .speedup(zi_tim,     "fastglm_zi",     "pscl_zi")

su_all <- rbind(
    transform(su_logit,  model = "logit"),
    transform(su_pois,   model = "poisson"),
    transform(su_nb,     model = "neg-binomial"),
    transform(su_firth,  model = "Firth (avg)"),
    transform(su_hurdle, model = "hurdle"),
    transform(su_zi,     model = "zero-inflated")
)

su_all$model <- factor(su_all$model, levels = unique(su_all$model))
ggplot(su_all, aes(n, speedup, colour = model, shape = model, group = model)) +
    geom_hline(yintercept = 1, linetype = 2, colour = "grey60") +
    geom_line(linewidth = 0.8) + geom_point(size = 2.4) +
    .size_axis() + scale_y_log10() +
    scale_colour_manual(values = c("#0B5CAD", "#D1661F", "#1F9E8F", "#9A62C4", "#9C8200", "#D0628E")) +
    labs(x = "Number of observations, log scale",
         y = "Reference time / fastglm time, log scale",
         title = "fastglm time ratio against each reference, by model class",
         colour = NULL, shape = NULL) +
    .benchmark_theme()

plot of chunk unnamed-chunk-18

Speed advantages depend on the design, so these comparisons do not establish a uniform speed or convergence advantage. Any timing comparison should first check coefficients, fitted means, uncertainty calculations, and convergence.

Ordinal models

These comparisons use packages that fit the same cumulative-link models. Each design has 10 predictors, and each time includes coefficient covariance. A time is compared only when that fit converged and agreed with fastglm on the log likelihood and common slopes.

plot of chunk ordinal-times

Link	Method	Seconds	Time ratio
logit	`fastglm`	0.067	1.0
logit	`rms::orm`	0.197	2.9
logit	`ordinal::clm`	0.690	10.3
logit	`MASS::polr`	1.405	20.9
logit	`VGAM::vglm`	1.762	26.2
logit	`WeightIt::ordinal_weightit`	3.049	45.4
probit	`fastglm`	0.101	1.0
probit	`rms::orm`	0.219	2.2
probit	`ordinal::clm`	0.688	6.8
probit	`MASS::polr`	1.800	17.9
probit	`VGAM::vglm`	2.007	20.0
probit	`WeightIt::ordinal_weightit`	6.429	63.9
cloglog	`fastglm`	0.100	1.0
cloglog	`ordinal::clm`	0.763	7.7
cloglog	`MASS::polr`	1.477	14.8
cloglog	`VGAM::vglm`	6.579	NA
cloglog	`WeightIt::ordinal_weightit`	5.922	59.5

The table shows 100,000 observations and 5 categories. The ratio is the comparison time divided by the fastglm time. rms::orm was timed for logit and probit only. For cloglog, VGAM::vglm warned about intersecting predictors and returned a different fit, so it has no ratio.

The threshold-specific panel has no comparison at 1,000 observations, because the fastglm fit reported nonconvergence for that design: the best fit with ordered cumulative probabilities at every row lies on the ordering boundary. ordinal::clm reported convergence with a higher log likelihood, but its cumulative predictors cross at 3 of the 1,000 rows. At 100,000 observations, ordinal::clm reported that it had not converged, so it is omitted there as well.

With 20,000 observations, rms::orm overtakes fastglm between 30 and 300 categories:

plot of chunk ordinal-categories

Categories	Method	Seconds	Time ratio
30	`fastglm`	0.017	1.00
30	`rms::orm`	0.044	2.56
30	`ordinal::clm`	0.649	38.06
30	`MASS::polr`	0.891	52.25
30	`VGAM::vglm`	5.443	319.08
30	`WeightIt::ordinal_weightit`	2.578	151.13
300	`fastglm`	0.098	1.00
300	`rms::orm`	0.049	0.50
300	`ordinal::clm`	41.008	418.83
300	`MASS::polr`	29.603	302.35
300	`WeightIt::ordinal_weightit`	47.976	490.00
1000	`fastglm`	0.316	1.00
1000	`rms::orm`	0.085	0.27
1000	`ordinal::clm`	483.794	1530.82
1000	`MASS::polr`	277.030	876.58
1000	`WeightIt::ordinal_weightit`	330.964	1047.24

rms::orm is designed for outcomes with very many distinct values and is the faster choice in that range. VGAM::vglm was not run above 30 categories because it would need about 160 GB of memory at 1,000.

None of these packages fits mean bias-reduced cumulative-link models or reads ordinal data in blocks:

plot of chunk ordinal-streaming

Observations	Fit	Seconds	Largest coefficient difference from dense fit
100,000	Dense ML	0.068	0.0e+00
100,000	Dense mean bias reduction	0.154	0.0e+00
100,000	Streamed ML	0.157	2.2e-16
100,000	Streamed mean bias reduction	0.362	1.1e-16
1,000,000	Dense ML	0.673	0.0e+00
1,000,000	Dense mean bias reduction	1.536	0.0e+00
1,000,000	Streamed ML	1.570	1.1e-16
1,000,000	Streamed mean bias reduction	3.540	3.0e-17

Streamed fits use blocks of 10,000 rows. The last column compares each streamed fit with the dense fit of the same estimator.

Benchmark Study for ‘fastglm’

Jared Huling

2026-11-06