gamFactory: GAM & Additive Modeling

• gamFactory is a computational framework for GAMs that integrates probabilistic additive stacking and complex covariate transformations.
• It employs additive multinomial-logit structures to model ensemble weights that depend on covariates, using MAP estimation and LAML for smoothing.
• The framework jointly estimates transformation and regression parameters, enabling efficient uncertainty quantification for large-scale forecasting tasks.

Searching arXiv for recent and relevant papers on “gamFactory” to ground the article. gamFactory is an R package and modelling framework associated, in the provided literature, with two closely related lines of work in generalized additive modelling: probabilistic additive stacking for disaggregate electricity demand forecasting and the integration of complex covariate transformations directly into multi-parameter Generalized Additive Models (GAMs) (Capezza et al., 2020, Collarin et al., 24 Nov 2025). In the first setting, gamFactory implements additive stacking, a probabilistic mixture-of-experts framework in which ensemble weights vary with covariates through an additive multinomial-logit structure (Capezza et al., 2020). In the second, it implements methods for embedding scalar-valued covariate transformations of potentially high-dimensional and complex covariates within the model itself, with joint estimation of transformation parameters and regression coefficients via maximum a posteriori (MAP) methods and smoothing parameter selection by Laplace approximate marginal likelihood (LAML) (Collarin et al., 24 Nov 2025). A separate paper also uses “gamFactory” to denote a factor-graph-based accelerator concept for autonomous machine computing, presented through LiDAR-Inertial Odometry (LIO) (Hao et al., 2022). This usage is distinct from the R package.

1. Definition and scope

In "Additive stacking for disaggregate electricity demand forecasting" (Capezza et al., 2020), gamFactory is presented as the implementation of a new ensemble method for probabilistic forecasting in which expert forecasts are aggregated using covariate-dependent weights modelled by linear combinations of parametric, smooth, or random effects. The package is described as providing the machinery to fit probabilistic additive stacking models, estimate MAP coefficients, select smoothing parameters via LAML, handle nonlinear multinomial-logit weight models, and work with large-scale data efficiently (Capezza et al., 2020).

In "Integrating Complex Covariate Transformations in Generalized Additive Models" (Collarin et al., 24 Nov 2025), gamFactory is presented as the implementation of a framework for integrating interpretable covariate transformations within multi-parameter GAMs. In that formulation, transformations are treated as integral model components rather than as external preprocessing, and their parameters are estimated jointly with the rest of the model (Collarin et al., 24 Nov 2025).

These two uses are methodologically aligned. Both papers describe penalized-likelihood or Bayesian-prior views of estimation, MAP optimization for coefficient estimation, and empirical Bayes selection of smoothing parameters through a Laplace approximation to the marginal likelihood (Capezza et al., 2020, Collarin et al., 24 Nov 2025). This suggests a unifying characterization of gamFactory as a computational environment for GAM-style models in which either the mixture weights or the covariates themselves are modelled through flexible structured components.

2. Additive stacking for probabilistic forecasting

The forecasting formulation in (Capezza et al., 2020) targets household-level electricity demand, where the object of interest is the conditional predictive density

$p(y_i^c \mid \bm x_i^c).$

The paper defines four probabilistic experts $M_1,\dots,M_4$: LastMonth, GaulssInd, Dynamic, and GaulssCommon, each capturing different demand dynamics (Capezza et al., 2020). For household $c$, expert $k$ yields a predictive density $p_k^c(y_i^c \mid \bm x_i^c)$, and the stacked model combines them as

$$p^c(y_i^c \mid \bm x_i^c) = \sum_{k=1}^K \alpha_k(\bm x_i^c)\, p_k^c(y_i^c \mid \bm x_i^c),$$

with simplex constraints on the weights (Capezza et al., 2020).

The paper explicitly frames this as an extension of Breiman-style regression stacking. Two generalizations are central. First, the method performs probabilistic stacking by combining full densities rather than point predictions. Second, the weights $\alpha_k$ vary with covariates rather than remaining constant (Capezza et al., 2020). The probabilistic stacking objective is

$$\hat{\bm \alpha} = \arg\max_{\bm \alpha} \sum_{i=1}^N \log \sum_{k=1}^K \alpha_k \, p_k(y_i\mid \bm x_i),$$

again subject to positivity and simplex constraints (Capezza et al., 2020).

The additive structure enters through a multinomial-logit parameterization:

$$\alpha_{ki} = \frac{\exp(\eta_{ki})}{\sum_{a=1}^{K}\exp(\eta_{ai})}, \qquad k=2,\dots,K,$$

with $\eta_{1i}=0$ for identifiability, and

$M_1,\dots,M_4$0

where the $M_1,\dots,M_4$1 terms can be parametric linear effects, smooth effects based on spline bases, or random effects (Capezza et al., 2020). In this way, the “additive” in additive stacking refers specifically to the additive predictors governing the expert weights.

The paper emphasizes that the weight model can depend on contextual variables that the experts themselves may not use. In the Irish household study, the linear predictors incorporate day of week, historical average demand, historical standard deviation, days out of home, time of day, day of year, and recent relative expert performance summaries $M_1,\dots,M_4$2 (Capezza et al., 2020). The latter are defined as normalized geometric means of predictive densities over lag windows, thereby encoding how well each expert has performed recently (Capezza et al., 2020).

3. Estimation strategy and computational formulation

The fitting strategy in (Capezza et al., 2020) is based on a penalized posterior:

$M_1,\dots,M_4$3

up to a constant. Coefficients are estimated by MAP, while the smoothing parameters $M_1,\dots,M_4$4 are selected by maximizing the Laplace approximate marginal likelihood

$M_1,\dots,M_4$5

where $M_1,\dots,M_4$6 is the negative Hessian at the posterior mode $M_1,\dots,M_4$7 (Capezza et al., 2020).

This MAP-plus-LAML scheme is presented as a major reason the framework is practical for large-scale forecasting. The motivating application involves a dataset with over 30 million observations, and the paper states that gamFactory was developed in part to remain computationally feasible at that scale (Capezza et al., 2020).

The workflow across households is deliberately heterogeneous. Experts $M_1,\dots,M_4$8, $M_1,\dots,M_4$9, and $c$0 are fit separately for each household, whereas $c$1 is fit jointly to all households. The stacking weights are then estimated using the full dataset across all households, allowing the model to borrow information globally while still producing household-specific weights through covariates (Capezza et al., 2020). To avoid overfitting, the stacking model is trained on out-of-sample expert forecasts generated by a rolling forecasting scheme (Capezza et al., 2020).

The paper also emphasizes interpretability. Although ALE plots are produced using mgcViz, they are described as part of the broader modelling workflow supported by the additive stacking approach implemented in gamFactory (Capezza et al., 2020).

4. Integrated covariate transformations within GAMs

The later gamFactory framework in (Collarin et al., 24 Nov 2025) shifts the emphasis from combining experts to transforming predictors inside the model. In a standard GAM or GAMLSS-type formulation, each distributional parameter $c$2 is linked to a linear predictor through

$c$3

The extension allows predictors to contain both ordinary smooth terms and nested smooth effects:

$c$4

Here, $c$5 are standard GAM smooths, $c$6 are covariate transformations, and $c$7 are outer smooth effects applied to the transformed covariate (Collarin et al., 24 Nov 2025).

The crucial definitional point is that a transformation is a scalar-valued function

$c$8

parameterized by $c$9, with the requirement that it be sufficiently differentiable; the paper uses fourth-order differentiability to support the fitting machinery (Collarin et al., 24 Nov 2025). The input may nevertheless be high-dimensional and complex, including time series, spatial fields, lagged histories, or other high-dimensional objects (Collarin et al., 24 Nov 2025).

The paper gives three main classes of transformations. The first is adaptive exponential smoothing:

$k$0

with $k$1 and $k$2, where $k$3 is the logistic function (Collarin et al., 24 Nov 2025). The second is multivariate kernel smoothing for spatially misaligned covariates:

$k$4

The third is linear combinations or single-index effects:

$k$5

The paper notes that $k$6 can itself be a model matrix, basis expansion, or lagged-history vector (Collarin et al., 24 Nov 2025).

This formulation replaces external feature engineering with embedded model components. The transformation is not precomputed and fixed; its parameters become part of the global coefficient vector and are estimated jointly with the rest of the model (Collarin et al., 24 Nov 2025).

5. Joint estimation, uncertainty, and implicit differentiation

For the integrated-transformation framework, the posterior log-density is written as

$k$7

up to a constant (Collarin et al., 24 Nov 2025). For fixed $k$8, MAP estimation is done with Newton’s method (Collarin et al., 24 Nov 2025).

Joint uncertainty is quantified by approximating the posterior of $k$9 with a Gaussian centered at the MAP estimate:

$p_k^c(y_i^c \mid \bm x_i^c)$0

where

$p_k^c(y_i^c \mid \bm x_i^c)$1

This yields joint uncertainty for standard GAM coefficients, smoothing coefficients, transformation parameters, and nested-effect parameters (Collarin et al., 24 Nov 2025). For a nested effect $p_k^c(y_i^c \mid \bm x_i^c)$2, the paper uses the Delta method to propagate uncertainty (Collarin et al., 24 Nov 2025).

Smoothing parameter selection again proceeds by empirical Bayes through a Laplace approximation to the marginal likelihood:

$p_k^c(y_i^c \mid \bm x_i^c)$3

with approximation

$p_k^c(y_i^c \mid \bm x_i^c)$4

The optimization is carried out over $p_k^c(y_i^c \mid \bm x_i^c)$5 to enforce positivity (Collarin et al., 24 Nov 2025).

A major technical feature of the paper is efficient derivative computation for models in which transformation parameters affect smooth basis evaluations. The difficulty is that, for nested effects, changes in inner parameters alter the transformed covariate and hence the basis used by the outer smooth. The paper extends implicit differentiation ideas from Wood et al. to derive efficient derivatives for such nested effects (Collarin et al., 24 Nov 2025). It presents a modular derivative structure for gradients and Hessian blocks, which is what makes MAP and LAML optimization computationally feasible for these richer models (Collarin et al., 24 Nov 2025).

The paper also identifies a scaling and identifiability issue: as the transformation parameters change, the scale and range of the transformed covariate may change, leading to poor knot placement, basis support mismatch, undersmoothing, and inefficiency. Its solution is to standardize the transformation and add a penalty that controls the sample variance of the standardized transformed covariate within estimation (Collarin et al., 24 Nov 2025).

6. Applications, related usages, and nomenclature

The additive-stacking paper develops gamFactory in the context of household electricity demand forecasting from smart-meter data, where accurate disaggregate forecasts are described as crucial for downstream tasks like battery scheduling (Capezza et al., 2020). The integrated-transformation paper demonstrates the later framework through applications to forecasting electricity net-demand in Great Britain and modelling house prices in London (Collarin et al., 24 Nov 2025). In the electricity application, the model jointly estimates linear effects, ordinary smooth effects, lagged demand single-index effects, spatial or region-weighted weather effects, and exponentially smoothed temperature effects with learned smoothing rates (Collarin et al., 24 Nov 2025). In the London house-price application, it combines standard covariate effects, a smooth spatial surface, and a nested kernel-smoothed neighborhood price index whose decay parameter is estimated jointly with the response model (Collarin et al., 24 Nov 2025).

A separate usage of the name appears in "Factor Graph Accelerator for LiDAR-Inertial Odometry" (Hao et al., 2022), where gamFactory denotes a factor-graph-based accelerator concept for autonomous machine computing rather than an R package. That work treats factor graph inference as a first-class computing primitive and presents a hardware architecture with factor, partial-QR decomposition, and back-substitution blocks for LiDAR, IMU, and GPS fusion (Hao et al., 2022). The nomenclature is therefore overloaded across distinct research contexts.

The term should also be distinguished from several similarly named but unrelated systems. "GFA: Exploratory Analysis of Multiple Data Sources with Group Factor Analysis" describes the GFA R package for Bayesian group factor analysis of multiple co-occurring data sources (Leppäaho et al., 2016). "GameFactory: Creating New Games with Generative Interactive Videos" concerns action-controlled scene-generalizable game video generation using pre-trained video diffusion models (Yu et al., 14 Jan 2025). "GameGPT: Multi-agent Collaborative Framework for Game Development" concerns LLM-based automation of game development workflows (Chen et al., 2023). These systems share lexical similarity with gamFactory but address different methodological domains.

Taken together, the GAM-oriented literature presents gamFactory as a computational framework for extending generalized additive modelling beyond standard fixed-covariate formulations. One strand models the expert-combination weights themselves as additive structures in a probabilistic stacking ensemble (Capezza et al., 2020). The other embeds transformations of complex covariates directly inside the model, with joint estimation, approximate Bayesian uncertainty quantification, and empirical Bayes smoothing selection (Collarin et al., 24 Nov 2025). The common methodological theme is the elevation of flexible structured components—whether weights or transformations—from external heuristics to estimable parts of the model.