---
title: "fastFGEE: Functional Generalized Estimating Equations"
author: "Gabriel Loewinger"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{fastFGEE: Functional Generalized Estimating Equations}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  warning = FALSE,
  message = FALSE,
  fig.align = "center",
  fig.width = 8.5,
  fig.height = 5.2,
  out.width = "100%"
)
library(fastFGEE)
```

# Installation

Install the CRAN release with:

```{r, eval = FALSE}
install.packages("fastFGEE")
```

The development version can be installed from GitHub:

```{r, eval = FALSE}
remotes::install_github("gloewing/fastFGEE", build_vignettes = TRUE)
```

# Quick start

`fastFGEE` fits marginal regression models for longitudinal functional outcomes.
Each row of the input data is one repeated observation, while the response is a
curve stored in a matrix-valued column.

## Loading example data

The package includes a small simulated binary dataset, `d`.

```{r data-preview}
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)
)

head(dat[c("ID", "X1", "X2", "time")])
head(Y[, 1:4])
```

`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.

```{r basic-fit, results = "hide"}
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:

```{r basic-fit-output}
fit_1step$sp.method
fit_1step$lambda
```

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.

```{r basic-plot}
fgee.plot(
  fit_1step,
  xlab = "Functional domain",
  title_names = c("Intercept", "X1", "X2")
)
```

The plotting data can also be returned for custom figures:

```{r, eval = FALSE}
plot_data <- fgee.plot(fit_1step, return = TRUE)
head(plot_data[[1]])
```

# Choosing a working correlation

There are two directions of dependence in a longitudinal functional dataset:

- `corr_long` describes association across repeated observations within a cluster;
- `corr_fn` describes association across points of the functional outcome.

The available choices are:

| Argument | Choices |
|---|---|
| `corr_long` | `"independent"`, `"exchangeable"`, `"ar1"` |
| `corr_fn` | `"independent"`, `"exchangeable"`, `"ar1"`, `"fpca"` |

Other structured working correlations, including stationary AR($p$) and Matérn
models on regular grids, are natural extensions. They are discussed under
[Numerical dependencies](#numerical-dependencies) below.

A working correlation is a device for improving efficiency. It does not need to
be the true covariance model for the marginal mean model to be useful.

For longitudinal dependence, `"exchangeable"` is useful when correlations are
roughly similar across visits. `"ar1"` is useful when observations closer in
time are expected to be more strongly related. For the functional direction,
`"ar1"` gives a simple local correlation model, while `"fpca"` provides a more
flexible low-rank covariance estimate.

## 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.

| `corr_long` | `corr_fn` | Interpretation |
|---|---|---|
| `"ar1"` or `"exchangeable"` | `"independent"` | model repeated-observation correlation only |
| `"independent"` | `"ar1"` or `"exchangeable"` | model functional-domain correlation only |
| `"independent"` | `"fpca"` | FPCA working covariance in the functional direction |
| non-independent | non-independent | separable / Kronecker working covariance |

For example:

```{r, eval = FALSE}
# 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.

# Variance estimators and confidence intervals

The covariance estimator and interval calibration are controlled separately.

`var.type` chooses the covariance estimate for the spline coefficients:

| `var.type` | Use |
|---|---|
| `"sandwich"` | default cluster-robust sandwich covariance |
| `"fastboot"` | fast cluster-bootstrap covariance; preferred over `"boot"` |
| `"boot"` | cluster bootstrap retained for compatibility |

`joint.CI = "wild"` uses the studentized wild cluster bootstrap to construct
pointwise intervals and an optional simultaneous band over the whole coefficient
function. This is the recommended setting for routine use.

For example, the first few interval limits for the `X1` coefficient in the fit
above are stored directly in the fitted object:

```{r ci-output}
head(fit_1step$crit$ci$ci_pointwise[[2]], 3)
head(fit_1step$crit$ci$ci_joint[[2]], 3)
```

```{r, eval = FALSE}
fit_sw <- fgee(
  Y ~ X1 + X2, data = dat, cluster = "ID",
  family = binomial(), time = "time",
  corr_long = "ar1", corr_fn = "ar1",
  var.type = "sandwich",
  joint.CI = "wild"
)

fit_fb <- fgee(
  Y ~ X1 + X2, data = dat, cluster = "ID",
  family = binomial(), time = "time",
  corr_long = "ar1", corr_fn = "ar1",
  var.type = "fastboot",
  boot.samps = 2000,
  joint.CI = "wild"
)
```

A pointwise interval answers a different question from a simultaneous band. The
simultaneous band is calibrated to cover the whole coefficient function in the
resampling procedure. As with other cluster-robust methods, finite-sample
performance can be poor when there are very few independent clusters or when a
small number of clusters have high leverage.

Set `joint.CI = FALSE` when confidence intervals are not needed.

# Starting from your own `refund::pffr()` fit

`fgee()` normally fits the initial `refund::pffr()` model for you. Supplying
`pffr.mod` is useful when you want more control over that initial fit, when the
functional grid is irregular, or when you want to compare the initial fit with
the one-step fGEE update.

`pffr.mod` supplies the aligned functional response and design information from
`pffr()`. The `data` argument should still be the original wide data set used to
define clusters and repeated observations.

## Common-grid example

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

```{r, eval = FALSE}
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.

```{r, eval = FALSE}
# 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.

```{r, eval = FALSE}
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.

# What the one-step fit does

The public estimator follows the same basic sequence for every supported
family:

1. fit the initial function-on-scalar model with `refund::pffr()` (unless
   `pffr.mod` is supplied);
2. estimate the working correlation from the initial residuals;
3. choose the smoothing parameters (fast cluster cross-validation by default);
4. take one penalized GEE update; and
5. rebuild the working quantities at the updated mean for variance estimation
   and confidence intervals.

This is the one-step estimator described in the paper. It avoids repeatedly
refitting the full estimating equation while still using the working covariance
to improve efficiency.

# Advanced options

Most analyses can use the defaults above. The options in this section are
mainly useful for large datasets, reproducibility work, or method development.

## Reducing the fitted object size

Large fitted objects can be reduced after fitting with:

```{r, eval = FALSE}
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:

| `sp.method` | Description |
|---|---|
| `"auto"` | recommended family-specific default |
| `"fastk_staged"` | staged fast cluster CV; Gaussian identity default |
| `"fastk_grad_fast"` | exact fastK criterion with a shorter adaptive search; non-Gaussian default |
| `"fastk_grad"` | earlier analytic-gradient fastK search |
| `"sandwich_qreml"` | experimental coefficient-space working restricted quasi-likelihood |
| `"qreml_fastk"` | qREML start followed by exact fastK optimization |

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.

### Why the non-Gaussian default changed in 0.2.0

Earlier builds used `"fastk_grad"` for non-Gaussian models. Much of its run time
was spent evaluating many starting values before the continuous optimizer did
very much work. `"fastk_grad_fast"` minimizes the same fastK criterion but uses
a shorter adaptive start search and usually one continuous optimization. In the
selector simulations it gave essentially the same estimation and coverage as
`"fastk_grad"` with less tuning time. Pure qREML could be more efficient when
the working correlation was well specified, but it was noticeably worse under
some functional-correlation misspecification settings, so it remains opt-in.

## 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:

```{r, eval = FALSE}
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)`.

```{r, eval = FALSE}
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:

```{r, eval = FALSE}
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"
)
```

### Other structured working correlations

The current public interface keeps the correlation choices deliberately small,
but several other structures can be handled efficiently without forming and
inverting a dense covariance matrix. On a regularly spaced one-dimensional
grid, a stationary correlation that depends only on lag has a Toeplitz
correlation matrix. Examples include stationary AR($p$) models and Matérn
kernels. Matérn $\nu=1/2$ is the exponential kernel; Matérn $\nu=3/2$ and
$\nu=5/2$ are common smoother alternatives. On an equally spaced grid these
structures are Toeplitz and can, in principle, use fast Toeplitz solvers such as
those in `SuperGauss` rather than dense inversion. Irregular spacing generally
breaks the Toeplitz structure and requires a different operator.

FPCA-based working correlations are a different type of extension. They are not
generally Toeplitz, but their low-rank structure can still be exploited with
low-rank matrix identities. `corr_fn = "fpca"` is the current example; richer
longitudinal or functional FPCA working correlations could be added in the
future.

AR($p$), Matérn, and additional FPCA-based structures are not currently public
`corr_long` or `corr_fn` options. If one of these structures would be useful for
your application, please contact the package author. Use cases are helpful for
prioritizing which correlation models to add next. For background, see Ling and
Lysy (2022) for fast stationary Toeplitz calculations in `SuperGauss`, and
Rasmussen and Williams (2006) for Matérn covariance functions.

The archived `irregulAR1` and `sanic` packages are not required by
`fastFGEE` 0.2.0.

# What changed in fastFGEE 0.2.0

The main user-facing changes are straightforward:

- public fitting uses the one-step estimator;
- `corr_long` and `corr_fn` make the two working-correlation directions explicit;
- `var.type` controls the variance estimator;
- `sp.method = "auto"` uses staged fastK for Gaussian identity models and
  `fastk_grad_fast` for supported non-Gaussian models;
- negative binomial and beta regression are supported through `mgcv::nb()` and
  `mgcv::betar()`;
- `fastk.memory = "lowmem"` is deprecated; and
- internal Rcpp/LAPACK and irregular-AR(1) routines remove the old `sanic` and
  `irregulAR1` package dependencies.

# Notes for users of older versions

The working-correlation directions are specified separately: `corr_long` acts
across repeated observations within a cluster and `corr_fn` acts along the
functional domain. The older `cov.type` interface is replaced by these two
arguments, and the older `sandwich` argument is replaced by `var.type`.

The public `fgee()` interface in 0.2.0 exposes only the validated one-step
estimator. Arguments such as `exact`, `gee.fit`, `max.iter`, `tune.method`, and
`working.engine` are not accepted by `fgee()`.

Historical exact-GLS, pffr-only, legacy-engine, and fully iterated
implementations remain unexported for regression testing and future method
development.