Loading example data

The package includes a small simulated binary dataset, d.

data("d", package = "fastFGEE")
dat <- d
Y <- as.matrix(dat$Y)

c(
  rows = nrow(dat),
  clusters = length(unique(dat$ID)),
  functional_grid_points = ncol(Y)
)
#>                   rows               clusters functional_grid_points 
#>                    300                     30                    100

head(dat[c("ID", "X1", "X2", "time")])
#>   ID        X1         X2 time
#> 1  1 0.2199248 -1.2810522    1
#> 2  1 0.2199248 -3.6745461    2
#> 3  1 0.2199248 -3.0682195    3
#> 4  1 0.2199248 -1.3994682    4
#> 5  1 0.2199248 -0.8282399    5
#> 6  1 0.2199248 -0.4335023    6
head(Y[, 1:4])
#>      Y_1 Y_2 Y_3 Y_4
#> [1,]   1   0   1   1
#> [2,]   1   1   1   1
#> [3,]   0   1   0   1
#> [4,]   1   1   1   1
#> [5,]   0   1   1   1
#> [6,]   1   1   1   1

ID identifies independent clusters, time orders repeated observations, and Y contains the functional response on a common grid. Scalar predictors can be cluster-level or observation-level variables.

Fit a one-step fGEE

A typical fit specifies the mean model, cluster variable, response family, and working correlation in the longitudinal and functional directions.

fit_1step <- fgee(
  formula = Y ~ X1 + X2,
  data = dat,
  cluster = "ID",
  family = binomial(link = "logit"),
  time = "time",
  corr_long = "ar1",
  corr_fn = "independent",
  rho.smooth = TRUE,
  var.type = "sandwich",
  joint.CI = "wild",
  verbose.tuning = FALSE
)

The selected tuning method and smoothing parameters are stored in the fitted object:

fit_1step$sp.method
#> [1] "fastk_grad_fast"
fit_1step$lambda
#> [1]  3.466562  2.923348 15.437922

All public fgee() fits use the one-step estimator. There is no iteration setting to choose for the final coefficient fit.

Plot the coefficient functions

fgee.plot() plots the estimated functional coefficients. When interval estimates are available, the plot also shows the pointwise intervals and the simultaneous whole-curve band.

fgee.plot(
  fit_1step,
  xlab = "Functional domain",
  title_names = c("Intercept", "X1", "X2")
)

The plotting data can also be returned for custom figures:

plot_data <- fgee.plot(fit_1step, return = TRUE)
head(plot_data[[1]])

Common correlation patterns

When only one direction is modeled as correlated, the working covariance is block diagonal. When both directions are correlated, fastFGEE uses a separable (Kronecker product) working covariance.

For example:

# Longitudinal correlation only
fit_long <- fgee(
  Y ~ X1 + X2, data = dat, cluster = "ID",
  family = binomial(), time = "time",
  corr_long = "ar1", corr_fn = "independent"
)

# Correlation in both directions
fit_sep <- fgee(
  Y ~ X1 + X2, data = dat, cluster = "ID",
  family = binomial(), time = "time",
  corr_long = "ar1", corr_fn = "ar1"
)

# Flexible functional covariance
fit_fpca <- fgee(
  Y ~ X1 + X2, data = dat, cluster = "ID",
  family = binomial(), time = "time",
  corr_long = "independent", corr_fn = "fpca"
)

Kronecker fits require a complete longitudinal-by-functional grid within each cluster.

Smoothing or pooling correlation estimates

With a block-diagonal working covariance, the correlation parameter can vary over the other index. rho.smooth = TRUE smooths those pointwise estimates. rho.pool can instead use one pooled parameter when that is more appropriate. The default, rho.pool = "fn", pools a functional-only correlation estimate; "long" pools a longitudinal-only estimate, "both" allows either, and "none" keeps the pointwise behavior used in older versions.

Common-grid example

For data observed on a common grid, you can fit pffr() directly and pass the result to fgee():

fit_pffr <- refund::pffr(
  Y ~ X1 + X2,
  data = dat,
  family = binomial(),
  algorithm = "bam",
  method = "fREML",
  discrete = TRUE,
  bs.yindex = list(bs = "bs", k = 11, m = c(2, 1))
)

fit_from_pffr <- fgee(
  Y ~ X1 + X2,
  pffr.mod = fit_pffr,
  data = dat,
  cluster = "ID",
  family = binomial(),
  time = "time",
  corr_long = "exchangeable",
  corr_fn = "independent"
)

fgee.plot(fit_from_pffr)

Irregular functional grid

If the functional response is observed on an irregular grid, it is usually easier to build the pffr() input explicitly and then pass the fitted object to fgee(). The important pieces are the response value and its functional-domain location. The grid values do not have to be equally spaced.

The example below starts from the same wide data layout used above and shows the full conversion. In your own data, replace s_grid with the actual functional domain locations.

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

# Functional-domain locations. These may be irregularly spaced.
s_grid <- attr(dat$Y, "yindex")
if (is.null(s_grid)) {
  s_grid <- seq_len(ncol(Y_wide))
}

stopifnot(length(s_grid) == ncol(Y_wide))

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

# Convert the matrix response to long form.
# tidyr is used here only to make the reshaping easy to read.
dat_long <- tidyr::pivot_longer(
  dat_wide,
  cols = tidyselect::starts_with("Y_"),
  names_to = "yindex_col",
  names_prefix = "Y_",
  values_to = "Y",
  values_drop_na = FALSE
)

# Map each response column back to its functional-domain location.
dat_long$yindex_col <- as.integer(dat_long$yindex_col)
dat_long$yindex <- s_grid[dat_long$yindex_col]
dat_long$time <- as.numeric(dat_long$time)

head(dat_long[c("ID", "time", "yindex", "Y")])

# 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_irregular <- 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 pass that initial fit to fgee(). Notice that pffr.mod receives the fitted pffr() object, while data is still the original wide data frame containing the cluster and longitudinal variables.

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

fgee.plot(fit_from_irregular_pffr)

This route is also useful when the initial pffr() fit needs a custom basis, custom smoothing settings, or other preprocessing that you want to control directly.

Reducing the fitted object size

Large fitted objects can be reduced after fitting with:

fit_small <- fgee(
  Y ~ X1 + X2,
  data = dat,
  cluster = "ID",
  family = binomial(),
  time = "time",
  corr_long = "exchangeable",
  corr_fn = "independent",
  joint.CI = FALSE,
  keep.data = FALSE,
  keep.initial.fit = FALSE,
  keep.working.stats = FALSE
)

working.retain = "auto" keeps the smallest set of cluster working statistics needed by the requested variance and CI procedure. The other choices are "aggregate", "scores", and "full".

Smoothing-parameter tuning

sp.method = "auto" is recommended. It uses "fastk_staged" for Gaussian identity-link models and "fastk_grad_fast" for the supported non-Gaussian families.

Available smoothing selectors include:

The last two are useful for research and diagnostics. sandwich_qreml is not the automatic default because it can be sensitive to working-correlation misspecification. qreml_fastk is retained for reproducibility; in current simulations it reached the same fastK solutions by a slower route. In the latest selector simulations, its run time was about 1.16–1.31 times that of fastk_grad_fast. Pure sandwich_qreml had worst-case coefficient-RMSE ratios of 1.565 for negative binomial and 1.479 for beta under functional-independence working-correlation misspecification. fastk_staged remains useful for Gaussian identity models but has not shown an advantage for the non-Gaussian families.

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.

Optional compiled fastK kernel

For supported non-Gaussian families, the balanced workspace can evaluate the fastK loss and gradient with the package’s compiled kernel. This changes only the computation, not the tuning criterion. It is used automatically when available and can be disabled for reference comparisons:

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

Negative binomial and beta families

Negative binomial and beta models each have one extra parameter that matters to the working variance and to the non-Gaussian fastK loss: the NB2 size parameter and the beta precision parameter. Use the mgcv family functions so those values are carried through the initial pffr() fit. In the examples below, nb_dat denotes count-response data and beta_dat denotes proportion-response data with values strictly inside (0, 1).

fit_nb <- fgee(
  Y ~ X1 + X2, data = nb_dat, cluster = "ID",
  family = mgcv::nb(),
  time = "time", corr_long = "exchangeable", corr_fn = "ar1"
)

# Fix theta if you want to supply it rather than estimate it initially
fit_nb_fixed <- fgee(
  Y ~ X1 + X2, data = nb_dat, cluster = "ID",
  family = mgcv::nb(theta = 3),
  time = "time", corr_long = "exchangeable", corr_fn = "ar1"
)

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

fastFGEE keeps the initial nuisance value fixed while smoothing parameters are tuned, then updates it at the final working-state calculation. It does not re-estimate the nuisance parameter for every candidate smoothing vector. MASS::negative.binomial() is not supported here; use mgcv::nb() instead. For beta regression, responses should lie strictly between 0 and 1. Boundary adjustments are reported rather than silently ignored.

Numerical dependencies

You normally do not need to select a numerical backend. Regular-grid AR(1) and exchangeable working correlations use direct package operators. Irregularly sampled continuous-time AR(1) uses an internal tridiagonal precision implementation based on Allévius (2018), and symmetric positive-definite solves use registered Rcpp/LAPACK routines.

SuperGauss is optional. Install it only if you want to request the Toeplitz backend explicitly:

install.packages("SuperGauss")

fit_sg <- fgee(
  Y ~ X1 + X2, data = dat, cluster = "ID",
  family = binomial(), time = "time",
  corr_long = "ar1", corr_fn = "independent",
  corr.solver = "supergauss"
)