Draft FPC

[strengejacke]

library(glmmTMB)
library(Matrix)
library(lme4)

# 1. Fit your glmmTMB model
model <- lme4::lmer(Reaction ~ Days + (Days | Subject), data = sleepstudy)

vcovFPC <- function(object, popsize2 = NULL, popsize1 = NULL, KR = FALSE) {
  # Obtained finite-population-adjusted standard errors for fixed effect
  # estimates for a fitted multilevel model
  #
  # Args:
  # object: an R object of class merMod as resulting from lmer()
  # popsize2: population size at level-2; if NULL, an infinite level-2
  # population is assumed
  # popsize1: population size at level-1; if NULL, an infinite level-1
  # population is assumed
  # KR: Whether Kenward-Roger approximation of standard errors should be used,
  # which is recommended for small number of clusters and average cluster size.
  # Default to FALSE.
  #
  # Returns:
  # The variance-covariance matrix of the fixed effect estimates, as
  # returned by vcov()
  if (!inherits(object, "merMod")) {
    stop("Wrong input: Not a fitted model from lmer() with class merMod")
  }
  if (length(object@flist) != 1) {
    stop("Wrong input: Only models with two levels are supported")
  }
  if (is.null(popsize1) & is.null(popsize2)) {
    message("No FPC specified; return results from lme4::vcov.merMod()")
    return(vcov(object))
  }
  PR <- object@pp
  N <- unname(object@devcomp$dims["n"])
  nclus <- unname(ngrps(object))
  if (isTRUE(popsize2 > nclus)) {
    fpc2 <- 1 - nclus / popsize2
  } else {
    fpc2 <- 1
    message("No FPC needed at level-2")
  }
  if (isTRUE(popsize1 > N)) {
    fpc1 <- 1 - N / popsize1
  } else {
    fpc1 <- 1
    message("No FPC needed at level-1")
  }
  if (fpc1 == 1 & fpc2 == 1) {
    message("Return results from lme4::vcov.merMod()")
    return(vcov(object))
  }
  A <- PR$Lambdat %*% PR$Zt
  Astar <- A * sqrt(fpc2)
  X <- PR$X
  Astar_X <- Astar %*% X
  D <- Matrix::Diagonal(nrow(Astar), fpc1) + tcrossprod(Astar)
  Fisher_I <- (crossprod(X) - crossprod(solve(t(chol(D)), Astar_X))) / fpc1
  Phi <- solve(Fisher_I) * sigma(object)^2
  Phi <- as(Phi, "dpoMatrix")
  nmsX <- colnames(X)
  dimnames(Phi) <- list(nmsX, nmsX)
  if (!KR) {
    return(Phi)
  } else {
    if (!require("pbkrtest")) {
      stop(
        "Please install the `pbkrtest` package for the use of Kenward-Roger correction!"
      )
    } else {
      SigmaG <- pbkrtest::get_SigmaG(object, details = 0)
      vcov_kr <- pbkrtest:::vcovAdj16_internal(Phi, SigmaG, X, details = 0)
      vcov_kr <- as(Phi, "dpoMatrix")
      return(vcov_kr)
    }
  }
}

parameters::model_parameters(model, effects = "fixed")
#> # Fixed Effects
#> 
#> Parameter   | Coefficient |   SE |           95% CI | t(174) |      p
#> ---------------------------------------------------------------------
#> (Intercept) |      251.41 | 6.82 | [237.94, 264.87] |  36.84 | < .001
#> Days        |       10.47 | 1.55 | [  7.42,  13.52] |   6.77 | < .001
#> 
#> Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
#>   using a Wald t-distribution approximation.
parameters::model_parameters(
  model,
  vcov = vcovFPC(model, 20),
  effects = "fixed"
)
#> No FPC needed at level-1
#> # Fixed Effects
#> 
#> Parameter   | Coefficient |   SE |           95% CI | t(174) |      p
#> ---------------------------------------------------------------------
#> (Intercept) |      251.41 | 4.00 | [243.52, 259.29] |  62.91 | < .001
#> Days        |       10.47 | 0.80 | [  8.89,  12.04] |  13.13 | < .001
#> 
#> Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
#>   using a Wald t-distribution approximation.
parameters::model_parameters(
  model,
  vcov = "fpc",
  vcov_args = list(cluster_size = 20),
  effects = "fixed"
)
#> # Fixed Effects
#> 
#> Parameter   | Coefficient |   SE |           95% CI | t(174) |      p
#> ---------------------------------------------------------------------
#> (Intercept) |      251.41 | 4.00 | [243.52, 259.29] |  62.91 | < .001
#> Days        |       10.47 | 0.80 | [  8.89,  12.04] |  13.13 | < .001
#> 
#> Uncertainty intervals (equal-tailed) and p-values (two-tailed) computed
#>   using a Wald t-distribution approximation.

^{Created on 2026-05-19 with reprex v2.1.1}

[gemini-code-assist]

Code Review

This pull request implements Finite Population Correction (FPC) for mixed models and refactors the Kenward-Roger degrees of freedom calculation for better readability and performance. The code review identified a logic error in the FPC factor calculation where population sizes equal to sample sizes were not handled correctly. Additionally, feedback was provided to maintain sparse matrix structures to prevent excessive memory usage in large models and to ensure that all external dependencies, specifically the 'Matrix' package, are verified before use.

[gemini-code-assist]

The logic for converting population sizes to finite population correction (FPC) factors should use >= instead of >. If the population size provided equals the sample size ($N=n$), the correction factor should be 0 ($1 - n/n$). Currently, if $N=n$, the condition fpc1 > n fails and the raw population size is used in the variance calculation, leading to incorrect results. Additionally, the default value of 1 (set on lines 211 and 214) should be excluded from this conversion to avoid treating a sample size of 1 as a population size of 1.

  # adjust population size factor when it's larger than or equal to actual sample or group size
  if (!is.null(vcov_args$popsize_level1) && fpc1 >= n) {
    fpc1 <- 1 - n / fpc1
  }
  if (!is.null(vcov_args$popsize_level2) && fpc2 >= n_grplevel) {
    fpc2 <- 1 - n_grplevel / fpc2
  }

[gemini-code-assist]

Similar to the previous point, converting the sparse block-diagonal Cholesky factor to a dense matrix using as.matrix() is unnecessary and potentially memory-intensive. Keeping it sparse allows the subsequent multiplication on line 234 to remain efficient and reduces the overall memory footprint of the FPC calculation.

    t(L_matrix / stats::sigma(model))

[mattansb]

Should also add a method for lm like models:

vcovFPC.lm <- function(object, popsize = NULL, ..., verbose = TRUE) {
  N <- insight::n_obs(object)
  V <- stats::vcov(object)
  if (popsize <= N) {
    if (verbose) {
      insight::format_alert("No FPC needed.")
    }
    return(V)
  }

  fpc <- (popsize - N) / (popsize - 1)
  return(V / fpc)
}

[strengejacke]

@gemini-code-assist review

[Copilot]

Pull request overview

This PR adds finite population correction (FPC) support to insight::get_varcov() (and related docs/tests), including optional Kenward–Roger adjustment, with the intent of enabling FPC-adjusted fixed-effect uncertainty estimates for (at least) two-level mixed models.

Changes:

• Add new "fpc" option to get_varcov() and route it through .get_varcov_sandwich() to a new vcovFPC() implementation.
• Introduce vcovFPC() S3 generic + methods (mixed models, lm/glm, glmmTMB alias) and add tests validating results against a reference implementation.
• Update documentation (Rd + roxygen + NEWS) for "fpc" and minor KR wording fixes.

Reviewed changes

Copilot reviewed 9 out of 13 changed files in this pull request and generated 7 comments.

Show a summary per file

File	Description
R/get_varcov.R	Documents new vcov = "fpc" option and adds reference.
R/get_varcov_sandwich.R	Adds "fpc" shortcut mapping to vcovFPC() and tweaks error formatting.
R/get_varcov_fpc.R	New implementation of vcovFPC() (FPC + optional KR) and helper for transposed Cholesky extraction.
tests/testthat/test-vcov_fpc.R	New tests for FPC behavior and error cases.
R/get_df.R	Minor doc wording change (“Kenward-Roger”).
R/get_df_kenward-rogers.R	Refactor/commenting improvements in KR ddf internals.
NEWS.md	Announces new "fpc" support in get_varcov().
NAMESPACE	Exports vcovFPC() and registers S3 methods.
man/get_varcov.Rd	Mirrors roxygen doc update for "fpc" + reference.
man/get_predicted.Rd	Mirrors vcov option docs update (includes "fpc").
man/get_predicted_ci.Rd	Mirrors vcov option docs update (includes "fpc").
man/get_df.Rd	Mirrors doc wording change (“Kenward-Roger”).
DESCRIPTION	Bumps package version.

Files not reviewed (4)

• man/get_df.Rd: Language not supported
• man/get_predicted.Rd: Language not supported
• man/get_predicted_ci.Rd: Language not supported
• man/get_varcov.Rd: Language not supported

---

💡 Add Copilot custom instructions for smarter, more guided reviews. Learn how to get started.

[gemini-code-assist]

Code Review

This pull request implements Finite Population Correction (FPC) for variance-covariance matrix estimation, adding support for lme4, glmmTMB, and basic linear models. The code review identified several critical issues, including a mathematical error where the FPC factor was used to divide rather than multiply the variance, and logic errors in population size boundary checks. High-priority feedback was also provided regarding performance, specifically advising against converting sparse matrices to dense formats to prevent memory exhaustion. Additionally, the reviewer noted that the glmmTMB implementation incorrectly relies on lme4-specific functions and suggested using internal package helpers for consistent matrix processing.

[mattansb]

@gemini-code-assist Does the glmmTMB implementation support fixed effects only models? Does it check for 2 levels?

[gemini-code-assist]

Hi @mattansb, thanks for the questions.

Regarding the glmmTMB implementation in vcovFPC:

1. Fixed effects only models: No, the current implementation does not support fixed-effects-only models. The function explicitly checks is_mixed_model(model), which will return FALSE for models without random effects, resulting in an error.

2. Level check: Yes, it does check for the number of levels. Specifically, it uses length(find_random(model)$random) to ensure that the model has exactly one grouping factor, which corresponds to a two-level model. If the model has more or fewer than one grouping factor, it will throw an error stating that only two-level models are currently supported.

[mattansb]

@strengejacke Then we should allow glmmTMB to also use the "regular" FPC when it only has fixed effects.

[strengejacke]

Ok, and what about the comments on your lm-method?

[mattansb]

DANIEL YOU'RE TOO FAST!

[strengejacke]

Maybe you're just too slow... 😎