Numerical dependencies

The default correlation solver handles regular-grid AR(1) and exchangeable working correlation directly. SuperGauss is suggested rather than required and is used only when corr.solver = "supergauss" is requested explicitly:

install.packages("SuperGauss")

Irregularly sampled continuous-time AR(1) correlations use an exact tridiagonal Markov precision implemented inside fastFGEE from the formulas of Allevius (2018). Symmetric positive-definite coefficient-space solves use registered Rcpp/LAPACK routines. The archived irregulAR1 and sanic packages are no longer dependencies.

What "auto" does

For Gaussian identity-link models, "auto" uses "fastk_staged". The exact held-out Gaussian loss can be written in terms of per-fold coefficient-space sufficient statistics, so once those are built the objective no longer depends on the number of functional observations. Tuning takes about the same time whether the data have \(10^5\) or \(2 \times 10^7\) scalar observations.

For supported non-Gaussian families – binomial with a logit link, Poisson, Gamma and negative binomial with a log link, and beta regression with a logit link – "auto" uses "fastk_grad_fast". There is no finite sufficient-statistic reduction for the exact non-Gaussian held-out loss, so every objective evaluation touches the data and the cost of tuning is dominated by how many evaluations the search needs.

Negative binomial and beta need their nuisance parameter

Gamma is the easy case: its dispersion enters the log-likelihood multiplicatively, plus terms free of the linear predictor, so it cancels out of the criterion and Gamma tuning never needed it. Negative binomial and beta are different – theta appears as \(\theta(\mu - y)/(\mu + \theta)\) and the precision inside \(\psi(\mu\phi)\) – so the criterion cannot be evaluated without a value.

The value used is the one that formed the working variance, not a fresh estimate. That is a correctness requirement rather than a convenience: the held-out loss scores per-fold coefficients that were themselves computed from \(\bar W\) and \(\bar d\) under that variance, so evaluating them under a different nuisance would no longer be a held-out predictive loss for the model being fitted.

Practically, that means supplying it through the family so it survives into the initial fit:

fit_nb <- fgee(
  formula = Y ~ X1 + X2, data = dat, cluster = "ID",
  family = mgcv::nb(),                 # or mgcv::nb(theta = 3) to fix it
  time = "time", corr_long = "exchangeable", corr_fn = "ar1"
)

fit_beta <- fgee(
  formula = Y ~ X1 + X2, data = dat, cluster = "ID",
  family = mgcv::betar(),
  time = "time", corr_long = "exchangeable", corr_fn = "ar1"
)

MASS::negative.binomial() will not work: mgcv rejects a plain family object that lacks variance derivatives, inside refund::pffr(), before any fastFGEE code runs. mgcv::nb(theta = ) is the way to fix theta.

One caveat to carry. Whatever mgcv::nb() or mgcv::betar() reports is estimated by REML under an independence likelihood, which is misspecified for correlated functional data, and a flexible smoother can absorb variability that would otherwise appear as overdispersion. The package therefore records which nuisance value was used and updates it only at the designated final working-state calculation; it does not re-estimate the nuisance parameter for every candidate smoothing vector.

Dedicated selector simulations support retaining fastK as the non-Gaussian default. Relative to fastk_grad_fast, fastK was never materially worse in the examined tiers, whereas pure sandwich_qreml had coefficient-RMSE ratios as large as 1.5653 for a concurrent negative-binomial functional-covariate design and 1.4785 for the corresponding beta design under functional-independence working-correlation misspecification. qreml_fastk reached the same fastK answers but was 1.16–1.31 times slower, and fastk_staged showed no advantage off Gaussian identity-link models. These results motivate the current automatic dispatch; they do not imply that fastK is uniformly optimal outside the validated simulation regimes.

Why the non-Gaussian default changed in 0.3.0

Before 0.3.0 the non-Gaussian default was "fastk_grad", which began by scoring a dense 55-point common-scale ray, a qREML candidate and two axis probes per smoothing parameter, and then ran four independent L-BFGS-B optimisations. With three smoothing parameters that is roughly 62 exact objective evaluations spent on choosing a starting point, against 4 to 46 spent in the optimiser. The start search, not the optimisation, was the expense.

"fastk_grad_fast" optimises the same criterion from a cheaper start:

1. a coarse common-scale ray, \(\lambda = 10^{c}\lambda_{\text{fREML}}\) for \(c \in \{-3, \dots, 3\}\), extended outward while its best point sits on an endpoint and then refined by one-dimensional search;
2. a sparse anisotropic stage around that point, letting the smoothing parameters separate – a \(3^q\) grid for small \(q\) and a coordinate sweep otherwise, so the cost grows linearly rather than exponentially in \(q\);
3. a single continuous polish, with a second start only when the optimiser fails to converge, stops on a boundary, leaves a large projected gradient, or a runner-up candidate is close in criterion but far away in \(\log\lambda\).

Probe evaluations drop from about 62 to about 34, isolated tuning is roughly three times faster, and complete fits are 1.2 to 1.5 times faster, with the larger gains at larger problems.

"fastk_grad" is still available and unchanged:

fit_prev <- fgee(
  formula = Y ~ X1 + X2, data = dat, cluster = "ID",
  family = binomial(), time = "time",
  corr_long = "exchangeable", corr_fn = "ar1",
  sp.method = "fastk_grad"
)

Evidence for the change

Changing a default selector changes which \(\lambda\) is returned, so the two selectors were compared in simulation rather than on timing alone.

Data were generated through a Gaussian copula: correlated latent normals are mapped to uniforms and then through the inverse marginal distribution function, so the response is exactly Bernoulli, Poisson or Gamma at the intended mean while remaining correlated within a cluster. The marginal coefficient functions that a GEE targets are therefore known in closed form. This matters: adding a correlated error to the linear predictor instead would make the marginal mean an attenuated version of the coefficient functions under a non-linear link, and comparing estimates against them would show bias that is an artefact of the simulation rather than a property of the estimator.

With 300 replicates per configuration, four cluster and visit configurations, and both selectors seeing identical data at each replicate:

• no comparison of slope root-mean-square error was significant after adjustment for multiplicity;
• pointwise and simultaneous coverage of the studentised wild cluster bootstrap intervals differed by at most 0.008, against a Monte Carlo standard error of 0.013;
• the same held under three working-correlation misspecification scenarios: independent longitudinally, independent functionally, and both.

A caution about reading \(\lambda\)

The fastK criterion is very flat in places. Two selectors can return smoothing parameters differing by orders of magnitude while agreeing on the fitted coefficient functions to five significant figures. When comparing selectors, compare the criterion value, the fitted coefficients, the effective degrees of freedom, or coefficient error – not \(\lambda\).

The optional compiled kernel

Each non-Gaussian objective evaluation computes, for every fold, a linear predictor, a loss and its derivative, and a crossproduct for the score contribution. In R that allocates several vectors as long as the fold for every fold of every evaluation. A compiled kernel does the same arithmetic in one pass with no intermediate allocation, making objective evaluations two to three times faster.

The kernel is an accelerator, not a change of method: it reproduces the R objective to about \(10^{-16}\) in both the score and the gradient, so turning it off cannot move a fitted model. It is used automatically when available, and every entry point falls back to the R implementation otherwise. To disable it:

fit_no_kernel <- fgee(
  formula = Y ~ X1 + X2, data = dat, cluster = "ID",
  family = binomial(), time = "time",
  corr_long = "exchangeable", corr_fn = "ar1",
  fastk.kernel = FALSE
)

# or for the whole session
options(fastFGEE.kernel = FALSE)

The compiled working-correlation inverse

A second compiled kernel sits under getD()/getW(). When both correlations are independence, exchangeable, or regular-grid AR(1), the closed-form precision operators are applied in one allocation and two strided passes, rather than by reshaping each cluster block, permuting it with aperm() and permuting it back. On the largest benchmark cell that took the working-statistics build from 2.27 s to 1.07 s and its share of a whole fit from 24% to 13%.

The same caveat as the fastK kernel applies with one addition. The kernel reproduces the R operator to about \(10^{-16}\), so it cannot meaningfully change a fitted model — but the fastK criterion is flat in some directions, and a perturbation that small can move a selected lambda on a plateau by orders of magnitude. Compare fit$beta rather than fit$lambda when checking a run with the kernel on against one with it off:

options(fastFGEE.corr.kernel = FALSE)

FPCA functional correlation and an explicitly requested corr.solver = "supergauss" remain separate paths. Irregular-grid AR(1) now uses the package’s compiled exact tridiagonal Markov precision, avoiding both a dense inverse and the archived irregulAR1 dependency.

Workspace memory for non-Gaussian tuning

fastk.memory = "balanced" is the recommended default. It is the only current layout that builds prep$X_eval, which is required by the compiled fastK loss/gradient kernel. The "speed" and "lowmem" layouts therefore use the R evaluation path.

Fresh-process measurements found no regime in which "lowmem" reduced peak resident memory; it retained the long working table and was substantially slower. It is deprecated and will be removed in a future release. "speed" remains available for cases where a modest reduction in retained workspace is worth slower R evaluation.

Peak memory must be assessed in a fresh R process because the long data, model matrix, package dependencies, and transient working-statistics construction usually dominate the retained fastK workspace.

Block diagonal versus Kronecker product working covariance

When exactly one direction is non-independent, the working covariance is block diagonal. This means the package models dependence in only one direction:

• corr_long != "independence" and corr_fn = "independence" gives a block diagonal covariance across longitudinal observations,
• corr_long = "independence" and corr_fn != "independence" gives a block diagonal covariance across functional observations.

For the parametric block-diagonal cases ("ar1" or "exchangeable"), the package estimates a different correlation parameter for each block across the other dimension:

• longitudinal block-diagonal fits estimate rho_long(s) across the functional domain,
• functional block-diagonal fits estimate rho_fn(t) across longitudinal time.

If rho.smooth = TRUE, those block-specific correlation curves are smoothed across the corresponding index.

When both directions are non-independent, the package uses a separable / Kronecker product working covariance. In this case the correlation parameters are pooled and fixed within direction:

• one pooled longitudinal correlation parameter (or longitudinal factor), and
• one pooled functional correlation parameter (or functional factor).

So the main difference is:

• block diagonal: correlation varies by block,
• Kronecker / separable: correlation is fixed within direction and shared across blocks.

For Kronecker fits, each cluster must be observed on a complete longitudinal-by-functional grid.

Longitudinal block-diagonal example

fit_long_block <- fgee(
  formula = Y ~ X1 + X2,
  data = dat,
  cluster = "ID",
  family = "binomial",
  time = "time",
  corr_long = "ar1",
  corr_fn = "independence",
  rho.smooth = TRUE,
  joint.CI = "wild",
  var.type = "sandwich"
)

This fit models within-cluster dependence across repeated observations, and allows the longitudinal correlation estimate to vary over the functional domain.

Functional block-diagonal example

fit_fn_block <- fgee(
  formula = Y ~ X1 + X2,
  data = dat,
  cluster = "ID",
  family = "binomial",
  time = "time",
  corr_long = "independence",
  corr_fn = "exchangeable",
  rho.smooth = TRUE,
  joint.CI = "wild",
  var.type = "sandwich"
)

This fit models dependence along the functional domain while treating repeated observations as working independent.

Separable / Kronecker example

fit_sep <- fgee(
  formula = Y ~ X1 + X2,
  data = dat,
  cluster = "ID",
  family = "binomial",
  time = "time",
  corr_long = "ar1",
  corr_fn = "ar1",
  joint.CI = "wild",
  var.type = "sandwich"
)

This fit uses a separable working covariance with one longitudinal AR1 factor and one functional AR1 factor.

FPCA-based functional covariance

fastFGEE also allows an FPCA-based covariance in the functional direction by setting corr_fn = "fpca".

fit_fpca_block <- fgee(
  formula = Y ~ X1 + X2,
  data = dat,
  cluster = "ID",
  family = "binomial",
  time = "time",
  corr_long = "independence",
  corr_fn = "fpca",
  joint.CI = "wild",
  var.type = "sandwich"
)

fit_fpca_sep <- fgee(
  formula = Y ~ X1 + X2,
  data = dat,
  cluster = "ID",
  family = "binomial",
  time = "time",
  corr_long = "ar1",
  corr_fn = "fpca",
  joint.CI = "wild",
  var.type = "sandwich"
)

Example: irregular functional grid

The following code shows the workflow using a user-supplied pffr() fit.

# Start from the wide data object used above
Y_wide <- as.matrix(dat$Y)
colnames(Y_wide) <- paste0("Y_", seq_len(ncol(Y_wide)))

dat_wide <- data.frame(
  ID = dat$ID,
  X1 = dat$X1,
  X2 = dat$X2,
  time = dat$time,
  Y_wide
)

# Convert the matrix outcome to long format
# (shown here with tidyr for readability)
dat_long <- tidyr::pivot_longer(
  dat_wide,
  cols = tidyselect::starts_with("Y_"),
  names_to = "yindex",
  names_prefix = "Y_",
  values_to = "Y",
  values_drop_na = FALSE
)

dat_long$yindex <- as.integer(dat_long$yindex)
dat_long$time <- as.numeric(dat_long$time)

# Construct the ydata object expected by refund::pffr()
Y.mat <- data.frame(
  .obs = seq_len(nrow(dat_long)),
  .index = dat_long$yindex,
  .value = dat_long$Y
)

fit_pffr <- refund::pffr(
  formula = Y ~ X1 + X2,
  algorithm = "bam",
  family = binomial(),
  discrete = TRUE,
  yind = Y.mat$.index,
  ydata = Y.mat,
  bs.yindex = list(bs = "bs", k = 11),
  data = dat_long
)

Now update that initial fit with fgee(). Notice that pffr.mod receives the fitted pffr() object, but data is still the original wide data frame.

fit_from_pffr <- fgee(
  formula = Y ~ X1 + X2,
  pffr.mod = fit_pffr,
  data = dat,
  cluster = "ID",
  family = "binomial",
  time = "time",
  corr_long = "exchangeable",
  corr_fn = "independence",
  joint.CI = "wild",
  var.type = "sandwich"
)

fgee.plot(fit_from_pffr)