f_lmer: Fit a linear mixed model with 'lme4::lmer()' including...

In rfriend: Provides Batch Functions and Visualisation for Basic Statistical Procedures

Fit a linear mixed model with lme4::lmer() including assumption checks, diagnostics, R-squared and post hoc tests.

Description

Fits a linear mixed-effects model using lme4::lmer() (with p-values from lmerTest) and produces a fully-formatted report containing the fixed-effects table, random-effects variance components, model-fit indices (AIC, BIC, logLik, marginal & conditional R^2), residual and BLUP diagnostics, convergence / singular-fit warnings, and post hoc comparisons (emmeans) on factor fixed effects. Results can be returned to the console or written to 'pdf', 'Word' or 'Excel'.

Usage

f_lmer(
  formula,
  data = NULL,
  REML = TRUE,
  ddf = "Satterthwaite",
  alpha = 0.05,
  adjust = "sidak",
  diagnostic_plots = TRUE,
  effect_plot = TRUE,
  contrast_plots = FALSE,
  post_hoc = TRUE,
  intro_text = TRUE,
  output_type = "default",
  save_as = NULL,
  save_in_wdir = FALSE,
  close_generated_files = FALSE,
  open_generated_files = interactive(),
  norm_plots = lifecycle::deprecated(),
  ...
)

Arguments

formula	A two-sided formula passed to lme4::lmer(), e.g. y ~ treatment + time + (1 | subject) or y ~ treatment * time + (1 + time | subject). The right-hand side must contain at least one random-effects term in the (varying | grouping) syntax. See Details for a guide to reading the random-effects syntax in study-design terms. More than one response variable can be supplied on the left-hand side using + (e.g. y1 + y2 ~ treatment + (1 | subject)). A separate model is then fit for each response, sharing the same right-hand side, and a multiple-testing warning is added to the report. See the Multiple Testing Across Response Variables section.
data	A data frame containing the variables in the model.
REML	Logical. If TRUE (default), the model is fit with restricted maximum likelihood, which gives less biased variance component estimates and is the appropriate choice for inference on fixed effects with Kenward-Roger or Satterthwaite degrees of freedom. Set to FALSE only when comparing nested models that differ in their fixed-effects structure.
ddf	Character. Method for computing denominator degrees of freedom for fixed-effects p-values. One of: "Satterthwaite" (default) Fast and accurate for most designs. Provided by lmerTest. "Kenward-Roger" Considered the gold standard, especially for small samples and unbalanced designs. Slower, and requires the pbkrtest package. "lme4" No p-values; only t-statistics are reported. Equivalent to a plain lme4::lmer() summary.
alpha	Numeric. Significance level for the fixed-effects table and the post hoc tests. Default is 0.05.
adjust	Character. Method used to adjust p-values for multiple pairwise comparisons in the post hoc step (passed to emmeans::emmeans()). One of "sidak" (default), "tukey", "bonferroni", "fdr", "none".
diagnostic_plots	Logical. If TRUE (default), diagnostic plots (residuals vs fitted, Q-Q of level-1 residuals, Q-Q of random-effect BLUPs, scale-location) are included in the output.
effect_plot	Logical. If TRUE (default), an estimated marginal means plot (estimate \pm 95% CI, with jittered raw data and compact-letter-display labels) is added after each categorical fixed-effect term's post hoc table. For a significant interaction between categorical factors, interaction plots are drawn instead, matching f_aov: a two-way interaction uses an x-axis factor and a colour trace and is shown in both orientations; three- and four-way interactions add facet panels for the remaining factor(s), with one plot per choice of x-axis factor. Interactions of order greater than four (which would need an illegible nested facet grid) are skipped with a note, but the interaction cell-means post hoc table still reports every combination. For a significant numeric x categorical interaction (the covariate slope differs across factor levels), f_lmer additionally draws a slope plot - a scatter of the raw data coloured by the factor with one model-fitted regression line per level and a confidence band. This goes beyond f_aov, which holds covariates at their mean and does not plot slopes. Numeric x numeric interactions and higher-order numeric/categorical mixes have no standard 2-D plot and are skipped with a note pointing to the coefficient-table slopes and emmeans::emtrends(). All effect and interaction plots are ggplot2 objects and are stored in the returned object (e.g. out$y1$effect_plot_treatment, out$y1$interaction_plot_a_b_1, out$y1$interaction_plot_a_b_c_1, out$y1$interaction_plot_dose_treatment_1 for a slope plot) so they can be retrieved and customised afterwards (themes, colours, axis labels, etc.).
contrast_plots	Logical. If TRUE, a contrast forest plot is added for each categorical post hoc term: one row per pairwise comparison, showing the estimated difference between two levels with its confidence interval and a reference line at zero. A CI that excludes zero indicates a significant difference; because the interval is on the difference itself, this "excludes zero" reading is exact (it is the same information the compact-letter display encodes, but it also shows the direction and magnitude of each difference). Default FALSE because the number of pairwise contrasts grows quickly with the number of factor levels (k levels give k(k-1)/2 contrasts); turn it on when you want the detailed pairwise picture. No cap is applied - if you enable it for a many-level factor you will get a tall figure, which is your choice. Main-effect and interaction contrast plots are kept separate: main-effect plots are stored as out$y1$contrast_plot_<term> (e.g. contrast_plot_treatment) and rendered with the main-effect post hoc tables, while interaction cell-contrast plots are stored as out$y1$interaction_contrast_plot_<term> (e.g. interaction_contrast_plot_a_b) and rendered with the interaction cell-means tables. Contrast CIs use the same adjust method as the post hoc p-values, so figure and table agree.
post_hoc	Logical. If TRUE (default), runs emmeans::emmeans() pairwise comparisons only when the linear mixed model finds a significant fixed-effect term (ANOVA p-value below alpha). The post hoc is performed on each significant factor fixed-effect term separately. Numeric covariates are skipped because their slope is already reported in the fixed-effects coefficient table; pairwise contrasts are not meaningful for a continuous predictor. If no fixed-effect term is significant, no post hoc is run. When a significant interaction between categorical factors is present, an additional cell-means post hoc table is produced for that interaction (estimated mean for every factor-level combination, compared simultaneously, with a compact letter display and pairwise contrasts), matching f_aov. It is stored alongside the main-effect results under the interaction term name, e.g. out$y1$post_hoc[["a:b"]]. In addition, when a main-effect term takes part in a significant interaction, a caution note is printed directly above that term's marginal-means table (and the term's heading is annotated), warning that the marginal means average over the interacting factor and can hide or reverse the real pattern; the interaction cell-means table and plot(s) should be read instead.
intro_text	Logical. If TRUE (default), prepends an explanation of LMM assumptions and the random-effects syntax linked to study design.
output_type	Character. Output format. One of: "default": returns the f_lmer object; auto-prints if unassigned. "console": forces immediate printing. "pdf", "word", "excel": writes a file. "rmd": stores the raw markdown string in the returned object for embedding in an R Markdown chunk with {r, echo=FALSE, results='asis'}.
save_as	Character. Output file path. See f_aov for the resolution rules; default is file.path(tempdir(), "<dataname>_lmer_output.<ext>").
save_in_wdir	Logical. If TRUE, save in the working directory instead of tempdir(). Default FALSE.
close_generated_files	Logical. Closes any open Word or Excel files before writing. Cross-platform (Windows taskkill, macOS / Linux pkill). Default FALSE. WARNING: save your work first.
open_generated_files	Logical. Whether to open the generated output files after creation. Defaults to TRUE in an interactive R session and FALSE otherwise (e.g. in scripts or automated pipelines). Set to TRUE or FALSE to override this behaviour explicitly.
norm_plots	'r lifecycle::badge("deprecated")' Deprecated in version 4.0.0. Use diagnostic_plots instead. If supplied, its value is passed through to diagnostic_plots with a warning.
...	Additional arguments forwarded to lmer. The arguments subset, na.action, and weights are handled specially: when supplied, they are applied via model.frame once before the per-response loop, so every response in a multi-response call is fitted on the identical row set. Other lmer() arguments (e.g. control = lmerControl(...), contrasts, offset) are forwarded unchanged on every fit.

Details

What is a linear mixed model?
A linear mixed model (LMM) extends ordinary regression / ANOVA by allowing two kinds of effects:

• Fixed effects - factors you actively manipulated or whose specific levels you care about (treatment, dose, time, genotype). Reported as estimates with confidence intervals.

• Random effects - grouping structure that creates non-independence in your data but whose levels are a random sample from a larger population (subjects measured repeatedly, plots within fields, observers, batches). Reported as variance components.

Use an LMM whenever observations share something that makes them more alike than two random observations from the dataset. Ignoring such grouping (running a plain aov or lm) is pseudoreplication, i.e. treating non-independent observations as if they were independent: standard errors shrink, p-values shrink, false positives explode.

Vocabulary.
Before going further, a few terms used throughout the report:

• Subject - the experimental unit that is measured repeatedly (a person, animal, pot and plot, cell line); in lme4 syntax it is the grouping factor on the right of the |, e.g. (1 | subject).

• Within-subject factor - a predictor whose levels vary within the same subject (time in a longitudinal study, treatment in a cross-over study).

• Between-subject factor - a predictor whose levels vary across subjects but are constant within a subject (sex, genotype, treatment arm in a parallel-groups trial). Both within- and between-subject factors are fixed effects.

• BLUP - Best Linear Unbiased Predictor. The model's estimate of the random-effect value for each subject (e.g. how much a particular subject deviates from the population intercept). BLUPs are checked for normality just like residuals.

• ICC - intraclass correlation coefficient. The share of total variance attributable to between-group differences. ICC = 0 means the grouping factor is irrelevant; ICC = 1 means observations within a group are identical.

• REML - restricted maximum likelihood. The default fitting method for variance components; gives less biased estimates than ordinary maximum likelihood.

• Satterthwaite / Kenward-Roger - methods to approximate the denominator degrees of freedom for fixed-effect p-values, since there is no exact df in an LMM.

Reading the (1 | group) syntax.
Every random-effects term has the form ( <varying> | <group> ). The bar reads as "varies by". The grouping factor on the right is what creates the non-independence. The left side is what is allowed to differ between groups. Common patterns:

• (1 | subject) - random intercept per subject (each subject has its own baseline). Repeated measures, longitudinal data.

• (1 | field) - randomised block design or multi-site trial; one intercept per block.

• (1 | field/plot) - plot nested in field; equivalent to (1|field) + (1|field:plot). Split-plot or hierarchical sampling.

• (1 + time | subject) - random intercept and random slope of time per subject. Subjects differ both in baseline and in how fast they change. Growth curves.

• (1 | subject) + (1 | observer) - crossed random effects: every observer can rate every subject. Inter-rater designs.

Rule of thumb: if you can answer "if I duplicated this experiment, would I draw new levels of this factor?" with yes, it belongs on the right of a |. If you would re-use the exact same levels (e.g. control vs treated) it is a fixed effect.

When to use a linear mixed model.
The most common reason is a repeated-measures design, in which the same experimental units are measured on more than one occasion or under more than one treatment. Compared with a between-groups design analysed by plain ANOVA this gives two real advantages: fewer experimental units are needed (each subject acts as its own control, removing between-subject variation from the comparison) and individual differences cannot bias the treatment groups (in a cross-over design every subject receives every treatment). Two canonical examples:

• Longitudinal study - same subjects measured at several time points: y ~ time + (1 | subject). If subjects also differ in how fast they change, add a random slope: y ~ time + (1 + time | subject).

• Cross-over design - every subject receives every treatment in sequence: y ~ treatment + (1 | subject). If carry-over between periods is a concern, add period as a fixed effect.

LMMs also apply to non-repeated structures that still create non-independence: randomised block designs, split-plot trials, multi-site studies, inter-rater designs.

Assumptions of a linear mixed model:

1. Linearity in the parameters of the fixed-effects part.

2. Independence of observations conditional on the random effects. If structure remains (e.g. temporal autocorrelation), more random effects or a correlation structure are needed.

3. Normality of level-1 residuals (Q-Q plot of residuals(m)).

4. Normality of the random-effect BLUPs (Q-Q plot of ranef(m)). This is the assumption most users forget.

5. Homoscedasticity: residual variance roughly constant across fitted values and across grouping levels.

6. At least ~5 levels of each grouping factor; with 3-4 levels it is usually better to treat the factor as fixed.

After the diagnostics the report adds a Recommendations for Heteroscedasticity and/or non-normal residuals section, but only when a flagged violation survives a follow-up check rather than on a raw significant p-value alone. A significant Levene's test triggers it only when the Scale-Location panel shows a corroborating variance trend across the fitted scale (a by-group Levene test on its own over-fires); a significant Shapiro-Wilk test on the level-1 residuals triggers it only when the rejection reflects genuine skew rather than a few outliers or heavy but symmetric tails (judged by the tail-trimmed Q-Q correlation and the residual skewness); a significant Shapiro-Wilk test on the random-effect BLUPs triggers it directly. The section gives concrete next steps (generalised mixed model, transformation).

Convergence and singular fits.
f_lmer surfaces lme4 convergence warnings and the "boundary (singular) fit" message prominently in the output. A singular fit usually means the random-effects structure is too complex for the data (often a random slope with too few levels) - simplify the model before interpreting results.

This function requires Pandoc (>= 1.12.3) for pdf, word and rmd output. See f_aov for installation notes.

Value

An object of class f_lmer: a named list containing the fitted lmerModLmerTest model, the ANOVA-style fixed-effects table (fixed_effects; the displayed table reports NumDF, DenDF, F and p only, while the full lmerTest table including the non-additive Sum Sq / Mean Sq columns is kept in fixed_effects_full), the variance components and ICC (plus a per-grouping-factor ICC breakdown in icc_by_group when the model has two or more intercept grouping factors), the R^2 values, the observed descriptives table (raw-data n, mean, sd, se, min, Q1, median, Q3, max grouped by the categorical fixed-effect predictors), post hoc results (if any; per main-effect terms plus, for a significant categorical interaction, a cell-means entry keyed by the interaction term name), diagnostic plots, and convergence diagnostics. When more than one response variable is supplied on the left-hand side, these elements are nested one level deep under each response name, e.g. out$y1$fixed_effects, out$y2$fixed_effects. When output_type = "rmd" the markdown string is stored in $rmd.

Multiple Testing Across Response Variables

When several response variables are analysed in a single call (e.g. y1 + y2 + y3 ~ treatment + (1 | subject)), each linear mixed model is an independent null-hypothesis test at level alpha. The post hoc adjustments (adjust = "sidak", "tukey", etc.) only control the family-wise error rate within one model (across pairwise contrasts for that response). They do not protect against the inflation of Type I error across the set of responses.

Practical implication: With k independent response variables all tested at \alpha = 0.05, the probability of obtaining at least one false positive is 1 - (1 - 0.05)^k, which reaches ~40% for k = 10.

When this matters: The risk is highest in exploratory studies where many responses are screened simultaneously without a clear a priori hypothesis for each one. It is less of a concern when each response is a pre-specified primary outcome with its own biological rationale.

Possible remedies:

• Bonferroni correction across responses: use alpha = 0.05 / k where k is the number of response variables. Conservative but simple.

• False Discovery Rate (FDR): apply p.adjust(p_values, method = "fdr") to the vector of per-response fixed-effect p-values after the fact.

• Multivariate model: if the responses are correlated and you want a single omnibus test, fit a joint multivariate mixed model (e.g. MCMCglmm, brms) before interpreting individual responses.

• Pre-registration: declare primary vs. exploratory responses before data collection to justify differential correction thresholds.

Author(s)

Sander H. van Delden [e-mail]

Examples

# sleepstudy: reaction time vs days of sleep deprivation,
# repeated measures within Subject (ships with lme4).
data(sleepstudy, package = "lme4")

# 1) Random intercept per subject - the simplest mixed model.
#    Each subject has its own baseline reaction time; the fixed
#    effect of Days is the average slope across subjects.
#    With output_type = "default" (the default), the result auto-
#    prints if not assigned, so no print() call is needed.
f_lmer_out <- f_lmer(Reaction ~ Days + (1 | Subject),
                     data = sleepstudy)


# Re-print the stored result and show the diagnostic plots.
print(f_lmer_out)
plot(f_lmer_out)

# 2) Random intercept AND random slope of Days per subject,
#    fitted with Kenward-Roger denominator df, saved to MS Word.
f_lmer(Reaction ~ Days + (1 + Days | Subject),
       data = sleepstudy,
       ddf  = "Kenward-Roger",
       output_type = "word"
       )

# 3) A factor fixed effect triggers a post hoc test.
#    Bin Days into three sleep-deprivation phases so that the
#    fixed effect is categorical and emmeans pairwise comparisons
#    with a compact letter display are produced automatically.
sleepstudy$Phase <- cut(sleepstudy$Days,
                        breaks = c(-Inf, 2, 6, Inf),
                        labels = c("early", "mid", "late"))
f_lmer(Reaction ~ Phase + (1 | Subject),
       data = sleepstudy,
       adjust = "tukey")

# 4) A minimal report: suppress the intro text and the diagnostic
#    plots, and save it directly to MS Word. Useful when embedding
#    many models in one document or when you only need the tables.
f_lmer(Reaction ~ Days + (1 | Subject),
       data = sleepstudy,
       intro_text = FALSE,
       diagnostic_plots = FALSE,
       output_type = "word"
       )

# 5) Get the raw markdown back for embedding in an R Markdown
#    document. Use it inside a chunk with results = 'asis'.
f_lmer_rmd_out <- f_lmer(Reaction ~ Days + (1 | Subject),
                         data = sleepstudy,
                         output_type = "rmd")
cat(f_lmer_rmd_out$rmd)

# 6) Two response variables analysed in one call. A separate model
#    is fit for each, sharing the same right-hand side. The results
#    are nested under each response name.
sleepstudy$Reaction2 <- sleepstudy$Reaction + rnorm(nrow(sleepstudy), 0, 5)
multi_out <- f_lmer(Reaction + Reaction2 ~ Days + (1 | Subject),
                    data = sleepstudy,
                    intro_text = FALSE,
                    diagnostic_plots = FALSE)
multi_out$Reaction$fixed_effects
multi_out$Reaction2$fixed_effects

# 7) Blocks, treatment and time together: a randomized complete
#    block design with repeated measures. The bundled plant_trial
#    dataset has five field blocks (block), one plant per treatment
#    per block measured at three time points (plant_id), four
#    treatments and a time factor.
data(plant_trial)

# 'treatment' and 'time_weeks' are crossed fixed effects (we care about their
# main effects and their interaction). 'block' and 'plant' are sources
# of unwanted variation we want to account for, not estimate, so
# they are random. 'plant_id' is nested in block - written compactly
# as (1 | block/plant_id), which expands to
# (1 | block) + (1 | block:plant_id). The random plant intercept is
# what makes this a repeated-measures model: the three time points
# on one plant share that plant's level.
f_lmer(height_cm ~ treatment * time_weeks + (1 | block/plant_id),
       data = plant_trial)

# The same design fitted with Kenward-Roger denominator df (the
# gold standard for small, balanced designs like this one), with
# the contrast forest plots turned on and the report saved to Word.
f_lmer(height_cm ~ treatment * time_weeks + (1 | block/plant_id),
       data = plant_trial,
       ddf = "Kenward-Roger",
       contrast_plots = TRUE,
       output_type = "word"
       )

rfriend documentation built on July 7, 2026, 1:06 a.m.