Speed matters most when a workflow repeats the same fit many times. Parametric bootstrap is the canonical case: a 1,000-replicate bootstrap means roughly 1,000 refits, plus the original fit and summary work. The same logic applies to cluster bootstrap, simulation studies, multi-seed reanalysis, and anything whose inner loop is “fit the same model again with new data”. A small per-refit speedup compounds into the difference between an overnight run and a coffee break.
mixeff is roughly two to five times faster than
lme4 on the parity benchmark — enough to make bootstrap a
routine option rather than a special-occasion calculation. This vignette
shows the figures from a single run of that benchmark and then documents
the harness used to make them. The numbers themselves are machine- and
model-dependent, so the scripts also write CSV files that you can
re-plot and rerun on the designs you actually care about.
What does the benchmark show?
The benchmark sweeps five mixed-model designs — varying rows, varying
subject levels, adding random slopes, crossing subject and item, and
combining crossed grouping with a random slope — across a log-spaced
grid of scale values. At each cell, mixeff::lmm() and
lme4::lmer() fit the same model three times. The summary
table records median elapsed seconds and speedup over lme4
per cell.
benchmark_path <- system.file(
"extdata", "lme4-scaling-summary.csv",
package = "mixeff"
)
benchmark <- read.csv(benchmark_path, stringsAsFactors = FALSE)
benchmark$scenario_label <- factor(
benchmark$scenario,
levels = c("rows", "groups", "slopes", "crossed", "crossed_slope"),
labels = c(
"vary rows (1 | subject)",
"vary subject levels (1 | subject)",
"add random slope (1 + x | subject)",
"crossed (1 | subject) + (1 | item)",
"crossed + slope (1 + x | subject) + (1 | item)"
)
)
head(benchmark[, c("scenario", "scale_value", "engine",
"median_sec", "speedup_vs_lme4")])
#> scenario scale_value engine median_sec speedup_vs_lme4
#> 1 crossed 10 lme4 0.008 1.000000
#> 2 crossed 10 mixeff 0.002 4.000000
#> 3 crossed 15 lme4 0.010 1.000000
#> 4 crossed 15 mixeff 0.003 3.333333
#> 5 crossed 20 lme4 0.011 1.000000
#> 6 crossed 20 mixeff 0.004 2.750000The first figure plots the median fit time against design scale, on
log-log axes, with one panel per scenario. mixeff runs
below lme4 everywhere on the grid; the gap widens as design
complexity grows.
library(ggplot2)
ggplot(benchmark,
aes(scale_value, median_sec * 1000, colour = engine,
group = engine)) +
geom_line() +
geom_point() +
facet_wrap(~ scenario_label, scales = "free_x", ncol = 2) +
scale_x_log10() +
scale_y_log10() +
scale_colour_manual(values = c(lme4 = "#888888", mixeff = "#3366cc")) +
labs(
x = "design scale (rows, subject levels, etc.)",
y = "median fit time (ms, log scale)",
colour = NULL
) +
theme_minimal(base_size = 11) +
theme(legend.position = "top",
panel.grid.minor = element_blank())
The second figure folds the same data into one panel: the speedup of
mixeff over lme4 at each cell. The dashed
reference line at 1× would mean the two engines are equally fast.
mixeff_rows <- subset(benchmark, engine == "mixeff")
ggplot(mixeff_rows,
aes(scale_value, speedup_vs_lme4,
colour = scenario_label, group = scenario_label)) +
geom_hline(yintercept = 1, linetype = 2, colour = "grey50") +
geom_line() +
geom_point() +
scale_x_log10() +
scale_y_continuous(breaks = seq(0, 8, by = 1)) +
labs(
x = "design scale",
y = "speedup over lme4 (×)",
colour = NULL
) +
theme_minimal(base_size = 11) +
theme(legend.position = "bottom",
legend.direction = "vertical",
panel.grid.minor = element_blank()) +
guides(colour = guide_legend(ncol = 1))
Read the speedup figure with the caveats it deserves. Three
replications per cell give a stable median but a wide envelope; the
absolute timings are sub-15 ms throughout, so an extra millisecond of
system noise moves the ratio visibly. The qualitative pattern is robust:
mixeff is faster everywhere, and the speedup grows with the
design’s complexity — exactly the regime where bootstrap workflows spend
most of their time.
What should you measure?
There are two separate speed questions.
First, measure model fitting. This asks how quickly each engine solves the same mixed model as the number of rows, groups, slopes, or crossed random effects changes.
Second, measure bootstrap routes. This asks how long the user-facing inference verbs take when they repeatedly refit models.
data.frame(
bootstrap_replicates = c(50L, 200L, 1000L),
approximate_refits = c(50L, 200L, 1000L)
)
#> bootstrap_replicates approximate_refits
#> 1 50 50
#> 2 200 200
#> 3 1000 1000The arithmetic is simple, but it is the core product point: small
per-refit speed differences compound quickly when nsim is
large.
How do you benchmark fitting?
Use the scaling benchmark when you want to compare
mixeff::lmm() with lme4::lmer() over several
synthetic LMM designs.
system2(
"Rscript",
c(
"inst/benchmarks/lme4-scaling.R",
"--reps=5",
"--scenarios=rows,groups,slopes,crossed",
"--no-plots"
)
)That script writes one row per engine, design, scale value, and
timing repetition, plus a summary table with median elapsed seconds and
speed relative to lme4.
data.frame(
file = c(
"benchmarks/lme4-scaling/lme4-scaling-raw.csv",
"benchmarks/lme4-scaling/lme4-scaling-summary.csv"
),
contents = c(
"one timing row per engine/design/repetition",
"median seconds, fits per second, and relative speed"
)
)
#> file
#> 1 benchmarks/lme4-scaling/lme4-scaling-raw.csv
#> 2 benchmarks/lme4-scaling/lme4-scaling-summary.csv
#> contents
#> 1 one timing row per engine/design/repetition
#> 2 median seconds, fits per second, and relative speedHow do you benchmark bootstrap?
Use the bootstrap benchmark when you want to time the inference routes that make the audit-then-bootstrap story practical.
The benchmark includes
test_effect(method = "bootstrap"),
confint(method = "bootstrap"),
compare(method = "bootstrap"), and
lme4::bootMer() as an R-engine baseline. If
pbkrtest is installed, it also times
pbkrtest::PBmodcomp().
data.frame(
route = c(
"mixeff_test_effect_bootstrap",
"mixeff_confint_bootstrap",
"mixeff_compare_bootstrap_lrt",
"lme4_bootMer_fixef_distribution",
"pbkrtest_PBmodcomp"
),
target = c(
"term-level fixed-effect p-value",
"fixed-effect confidence interval",
"model comparison p-value",
"fixed-effect bootstrap distribution",
"parametric bootstrap model comparison"
)
)
#> route target
#> 1 mixeff_test_effect_bootstrap term-level fixed-effect p-value
#> 2 mixeff_confint_bootstrap fixed-effect confidence interval
#> 3 mixeff_compare_bootstrap_lrt model comparison p-value
#> 4 lme4_bootMer_fixef_distribution fixed-effect bootstrap distribution
#> 5 pbkrtest_PBmodcomp parametric bootstrap model comparisonWhat does inference method choice cost?
Speed is only one part of the inference decision. The route table also needs calibration evidence: how often a method rejects under a null, how often it detects a small alternative, and whether interval routes cover that small alternative.
The lightweight simulation scaffold records those three quantities over four fixtures: an interior fit, a boundary-like fit, a reduced-rank fit, and a small-group fit. The checked-in summary is produced by the default fast mode, which excludes bootstrap routes so the script runs in seconds. Slow mode adds bootstrap and bootstrap-LRT routes for local evidence gathering.
simulation_path <- system.file(
"extdata", "inference-method-simulation-summary.csv",
package = "mixeff"
)
if (!nzchar(simulation_path)) {
simulation_file <- file.path("inst", "extdata",
"inference-method-simulation-summary.csv")
simulation_path <- if (file.exists(simulation_file)) {
simulation_file
} else {
file.path("..", simulation_file)
}
}
simulation <- read.csv(simulation_path, stringsAsFactors = FALSE)
head(simulation)
#> method fixture mode n_reps type_I_error power_at_alt
#> 1 asymptotic_wald_z interior fast 2 0 1
#> 2 satterthwaite interior fast 2 0 1
#> 3 kenward_roger interior fast 2 NA NA
#> 4 bootstrap interior fast 0 NA NA
#> 5 profile interior fast 2 NA NA
#> 6 bootstrap_lrt interior fast 0 NA NA
#> coverage_at_alt
#> 1 1
#> 2 1
#> 3 NA
#> 4 NA
#> 5 NA
#> 6 NA
system2(
"Rscript",
c(
"inst/benchmarks/inference-method-simulation.R",
"--mode=slow",
"--reps=100",
"--nsim=199"
)
)The plot below shows whichever metrics are available for each route.
P-value routes contribute Type I error and power; interval routes
contribute coverage. Routes not wired in the current package version
remain in the CSV with n_reps = 0, so downstream reports
can distinguish “not run in fast mode” from “not part of the
contract”.
library(ggplot2)
simulation_long <- stats::reshape(
simulation,
varying = c("type_I_error", "power_at_alt", "coverage_at_alt"),
v.names = "value",
timevar = "metric",
times = c("Type I error", "Power at small alternative", "Coverage at small alternative"),
direction = "long"
)
simulation_long <- subset(simulation_long, n_reps > 0 & is.finite(value))
simulation_long$metric <- factor(
simulation_long$metric,
levels = c("Type I error", "Power at small alternative",
"Coverage at small alternative")
)
ggplot(simulation_long, aes(method, value, fill = method)) +
geom_col(width = 0.7, show.legend = FALSE) +
facet_grid(metric ~ fixture) +
coord_cartesian(ylim = c(0, 1)) +
labs(x = NULL, y = "empirical proportion") +
theme_minimal(base_size = 10) +
theme(
axis.text.x = element_text(angle = 45, hjust = 1),
panel.grid.minor = element_blank()
)
How do you read the result?
For a local speed claim, read the summary CSV and compute ratios from
elapsed seconds. Values above 1 mean the baseline route took longer than
the mixeff route.
summary <- read.csv("benchmarks/bootstrap-inference/bootstrap-inference-summary.csv")
mixeff_ci <- subset(summary, route == "mixeff_confint_bootstrap")
lme4_boot <- subset(summary, route == "lme4_bootMer_fixef_distribution")
data.frame(
comparison = "lme4 bootMer / mixeff bootstrap CI",
speedup = lme4_boot$median_sec / mixeff_ci$median_sec
)Treat that ratio as local evidence, not a universal constant. It
depends on the formula, optimizer, random-effects structure,
nsim, and machine. The point of the benchmark scripts is to
make the claim reproducible for the model class you are actually
fitting.
What should go in a report?
When you report benchmark results, include the model formula, number
of rows, number of grouping levels, random-effects structure,
nsim, repetitions, and the package versions. The raw CSV
has the timing rows; the summary CSV is for tables and figures.
data.frame(
field = c("formula", "nobs", "groups", "random_terms", "nsim", "reps", "versions"),
reason = c(
"defines the fixed and random structure being refit",
"controls the amount of data processed per fit",
"controls the random-effect dimension",
"controls covariance structure and optimizer work",
"controls bootstrap cost",
"shows timing stability",
"makes the run reproducible"
)
)
#> field reason
#> 1 formula defines the fixed and random structure being refit
#> 2 nobs controls the amount of data processed per fit
#> 3 groups controls the random-effect dimension
#> 4 random_terms controls covariance structure and optimizer work
#> 5 nsim controls bootstrap cost
#> 6 reps shows timing stability
#> 7 versions makes the run reproducible