# BigGAM


GAM for large datasets using discretized fitting.


Usage

``` python
BigGAM(
    formula,
    family=None,
    *,
    n_discrete=200,
)
```


[BigGAM](BigGAM.md#whittaker.BigGAM) is a drop-in subclass of `~whittaker.gam.GAM` for datasets too large to fit comfortably with the standard dense design matrix (roughly `n > 1_000_000`). It uses the `bam` approach of Wood, Li, & Shaddick (2017): each covariate is rounded onto a grid of at most [n_discrete](BigGAM.md#whittaker.BigGAM.n_discrete) representative values, and the smooth basis is evaluated only once per unique (combination of) discretized value(s) rather than once per observation. An index array records, for every observation, which unique bin it fell into (see `DiscretizedBlock.indices` in `build_discretized_model_matrix`).

Fitting still runs penalized iteratively reweighted least squares (P-IRLS), exactly as in `~whittaker.gam.GAM.fit`, alternating an inner coefficient update with an outer smoothing parameter selection. The difference is purely computational: instead of forming the full `n x p` design matrix `X` and computing `X'WX` and `X'Wz` directly, `bam_fit` (in `whittaker.fitting.bam`) accumulates these quantities per bin -- for each smooth's discretized block, observation weights are aggregated into per-bin totals via `numpy.bincount` over the bin indices, and the resulting `d x d` cross-products of unique basis rows ([d](TPRS.md#whittaker.TPRS.d) = number of unique bins) are scattered into the correct `p x p` block of `X'WX` (see `_compute_XtWX`). Because [d](TPRS.md#whittaker.TPRS.d) can be orders of magnitude smaller than `n`, this reduces the memory needed for the cross-product step from `O(n p)` to `O(d p)`, and the resulting fit closely approximates the exact (non-discretized) GAM fit on the same data.

[BigGAM](BigGAM.md#whittaker.BigGAM) is a drop-in subclass of [GAM](GAM.md#whittaker.GAM): `predict()`, `summary()`, [plot()](GAM.md#whittaker.GAM.plot), and [check()](check.md#whittaker.check) all work the same way as for [GAM](GAM.md#whittaker.GAM). The one internal difference is that `self._model_matrix.X` is an empty array (the dense design matrix is never materialized) so operations that would otherwise reconstruct per-term columns (e.g. `smooth_tests()`) instead re-expand each smooth's columns on demand from its `DiscretizedBlock` via `_expand_block_columns`.


## Parameters


`formula: str | Formula`  
Model formula as a string (e.g. `"y ~ s(x1) + s(x2) + x3"`), or an already-parsed [Formula](Formula.md#whittaker.Formula) object. Same syntax as `~whittaker.gam.GAM`.

`family: Family | None = None`  
Response distribution family, e.g. [Gaussian()](Gaussian.md#whittaker.Gaussian), [Binomial()](Binomial.md#whittaker.Binomial), [Poisson()](Poisson.md#whittaker.Poisson), [Gamma()](Gamma.md#whittaker.Gamma), or [Tweedie()](Tweedie.md#whittaker.Tweedie). Defaults to [Gaussian()](Gaussian.md#whittaker.Gaussian).

`n_discrete: int = ``200`  
Maximum number of unique representative values per covariate (or per combination of covariates, for multi-dimensional smooths). Defaults to `200`.


## Notes

[n_discrete](BigGAM.md#whittaker.BigGAM.n_discrete) controls the accuracy/memory tradeoff directly. Larger values give a discretized grid that more finely resolves each covariate's range, so the fit approaches the exact (non-discretized) GAM fit at the cost of more unique bins [d](TPRS.md#whittaker.TPRS.d) and therefore more memory and computation in the `X'WX` accumulation step. Smaller values reduce memory and speed up fitting, but coarsen the covariate resolution: because all observations within a bin share the same basis row, this can slightly bias smooths with high curvature or steep local features, since fine-scale variation within a bin is averaged away. The default of `200` is usually more than enough resolution for typical smooth terms; it rarely needs to be increased unless a covariate has an unusually large number of important local features. The benefit of discretization (versus plain [GAM](GAM.md#whittaker.GAM)) is only realized once `n` is much larger than [n_discrete](BigGAM.md#whittaker.BigGAM.n_discrete), i.e. for large datasets -- for small or moderate `n`, [GAM](GAM.md#whittaker.GAM) is simpler and just as fast.


## Examples


``` python
import numpy as np
import whittaker as wt
from whittaker.bam import BigGAM

rng = np.random.default_rng(0)
n = 5_000
x1 = rng.uniform(0, 1, n)
x2 = rng.uniform(0, 1, n)
y = np.sin(2 * np.pi * x1) + x2**2 + rng.normal(scale=0.2, size=n)

model = BigGAM("y ~ s(x1) + s(x2)", n_discrete=100).fit({"x1": x1, "x2": x2, "y": y})
print(model.summary())
```


    GAM fit summary
    ============================================================
    Formula:    y ~ s(x1) + s(x2)
    Family:     Gaussian(link='identity')
    Inference:  REML
    Observations: 0
    Coefficients: 19

    Parametric coefficients:
      Term                       Estimate    Std.Err    t value    p-value
      ------------------------ ---------- ---------- ---------- ----------
      (Intercept)                  0.3376     0.0000      0.000          1

    Approximate significance of smooth terms:
      Term                        EDF Ref.df     Chi.sq    p-value
      ------------------------ ------ ------ ---------- ----------
      s(x1)                      8.91      9  62254.971    < 1e-16
      s(x2)                      6.26      7  11059.123    < 1e-16

    Total EDF:  16.17
    Scale est:  0.040068
    Deviance:   199.6941
    Null dev:   3148.4427
    Dev. expl:  93.7%
    GCV score:  0.040198
    AIC:        -1880.28
    BIC:        -1774.91


This example uses a modest `n` for speed, but [BigGAM](BigGAM.md#whittaker.BigGAM)'s memory and speed advantage over plain [GAM](GAM.md#whittaker.GAM) really shows up once `n` reaches into the millions, where materializing the full design matrix would be impractical.


## Attributes

| Name | Description |
|----|----|
| [n_discrete](#n_discrete) | Number of discretization grid points per covariate. |

------------------------------------------------------------------------


### n_discrete


Number of discretization grid points per covariate.


`n_discrete: int`


This is the [n_discrete](BigGAM.md#whittaker.BigGAM.n_discrete) value passed to `__init__`: the maximum number of unique representative values each covariate (or combination of covariates, for multi-dimensional smooths) is rounded onto before the smooth basis is evaluated. It bounds the number of unique rows [d](TPRS.md#whittaker.TPRS.d) used when accumulating `X'WX` during fitting (see the class docstring), and therefore controls the tradeoff between fit accuracy and memory/speed.


## Methods

| Name | Description |
|----|----|
| [fit()](#fit) | Fit the BigGAM using discretized P-IRLS. |
| [smooth_tests()](#smooth_tests) | Approximate significance tests for each smooth term. |

------------------------------------------------------------------------


### fit()


Fit the BigGAM using discretized P-IRLS.


Usage

``` python
fit(
    data,
    *,
    smoothing_params=None,
    method="fREML",
    weights=None,
    select=False,
    vi_options=None,
    mcmc_options=None,
)
```


Builds a `DiscretizedModelMatrix` via `build_discretized_model_matrix` and fits it with `bam_fit`, which accumulates `X'WX` and `X'Wz` from per-bin basis blocks instead of the full design matrix (see the class docstring for how the discretization and cross-product accumulation work).


#### Parameters


`data: dict[str, numpy.ndarray] or InputData`  
Column-oriented data as `{name: 1-D array}` (or any `InputData`-compatible object). All columns referenced by the formula must be present and of equal length.

`smoothing_params: list of float = None`  
Fixed smoothing parameters `lambda_j`, one per smooth term, in formula order. If `None` (the default), smoothing parameters are selected automatically according to [method](CausalGAM.md#whittaker.CausalGAM.method).

`method: str = ``"fREML"`  
Criterion used to select smoothing parameters when `smoothing_params` is `None`. One of `"fREML"` (default, fast discretized REML), `"REML"`, `"ML"`, or `"GCV"`. See `~whittaker.gam.GAM.fit` for a description of each criterion.

`weights: numpy.ndarray = None`  
Observation (prior) weights, shape `(n,)`. Must be strictly positive.

`select: bool = ``False`  
If `True`, add an extra penalty on each smooth's null space so the term can be shrunk to exactly zero (double-penalty selection). Defaults to `False`.


#### Returns


`BigGAM`  
Returns `self` for method chaining, e.g. `model = BigGAM("y ~ s(x)").fit(data)`.


------------------------------------------------------------------------


### smooth_tests()


Approximate significance tests for each smooth term.


Usage

``` python
smooth_tests()
```


This is the [BigGAM](BigGAM.md#whittaker.BigGAM) counterpart to `~whittaker.gam.GAM.smooth_tests`, adapted to the discretized fitting path: because `self._model_matrix.X` is never materialized, the per-term design columns needed for the test are reconstructed on demand from each smooth's `DiscretizedBlock` via `_expand_block_columns` (or read directly from `parametric_cols` for non-smooth columns) rather than sliced out of a dense `X`. For each smooth term, the coefficient sub-vector, its covariance sub-block, and its expanded design columns are passed to `~whittaker.fitting.inference._smooth_test` to obtain an approximate chi-squared test of whether the term is uniformly zero.


#### Returns


`list of SmoothTestResult`  
One result per smooth term (or per `by=` level), each with `term_label`, [stat](PPCResult.md#whittaker.PPCResult.stat), `edf`, `ref_df`, and [p_value](PPCResult.md#whittaker.PPCResult.p_value) attributes.