---
title: 'gamFactory: GAM & Additive Modeling'
url: https://www.emergentmind.com/topics/gamfactory
type: topic
---

# gamFactory: GAM & Additive Modeling

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) [2005.10092] [2511.19234]. 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 [2005.10092]. 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) [2511.19234]. A separate paper also uses “gamFactory” to denote a factor-graph-based accelerator concept for autonomous machine computing, presented through LiDAR-Inertial Odometry (LIO) [2209.02207]. This usage is distinct from the R package.

## 1. Definition and scope

In "Additive stacking for disaggregate electricity demand forecasting" [2005.10092], 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 [2005.10092].

In "Integrating Complex Covariate Transformations in Generalized Additive Models" [2511.19234], 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 [2511.19234].

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 [2005.10092] [2511.19234]. 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 [2005.10092] 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 [2005.10092]. 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 [2005.10092].

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 [2005.10092]. 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 [2005.10092].

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
$$
\eta_{ki} = \sum_{j\in I_k} f_{kj}(\bm x_i), \qquad k=2,\dots,K,
$$
where the $f_{kj}$ terms can be parametric linear effects, smooth effects based on spline bases, or random effects [2005.10092]. 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 $\gamma_{ki}^{cu}$ [2005.10092]. The latter are defined as normalized geometric means of predictive densities over lag windows, thereby encoding how well each expert has performed recently [2005.10092].

## 3. Estimation strategy and computational formulation

The fitting strategy in [2005.10092] is based on a penalized posterior:
$$
\log p(\bm \beta \mid \bm y, \bm \lambda) = \sum_{i=1}^N \log \sum_{k=1}^K \alpha_{ki}(\bm \beta)\, p_k(y_i\mid \bm x_i) - \frac{1}{2} \sum_{g=1}^G \lambda_g \bm\beta^\top {\bf S}_g \bm\beta,
$$
up to a constant. Coefficients are estimated by MAP, while the smoothing parameters $\bm\lambda$ are selected by maximizing the Laplace approximate marginal likelihood
$$
\tilde{\mathcal V}(\bm\lambda) = \mathcal L(\hat{\bm\beta}) +\frac{1}{2}\log|{\bf S}^{\bm\lambda}|_{+} -\frac{1}{2}\log|\bm{\mathcal H}| +\frac{M_p}{2}\log(2\pi),
$$
where $\bm{\mathcal H}$ is the negative Hessian at the posterior mode $\hat{\bm\beta}$ [2005.10092].

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 [2005.10092].

The workflow across households is deliberately heterogeneous. Experts $M_1$, $M_2$, and $M_3$ are fit separately for each household, whereas $M_4$ 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 [2005.10092]. To avoid overfitting, the stacking model is trained on out-of-sample expert forecasts generated by a rolling forecasting scheme [2005.10092].

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 [2005.10092].

## 4. Integrated covariate transformations within GAMs

The later gamFactory framework in [2511.19234] shifts the emphasis from combining experts to transforming predictors inside the model. In a standard GAM or GAMLSS-type formulation, each distributional parameter $\theta_{ji}$ is linked to a linear predictor through
$$
g_j(\theta_{ji})=\eta_{ji}.
$$
The extension allows predictors to contain both ordinary smooth terms and nested smooth effects:
$$
\eta_{ji} = \mathbf Z^0_{ji}\gamma_{j0} + \sum_{k=1}^{K_j} f_{jk}\!\left(\mathbf x_i^{S_{jk}}\right) + \sum_{u=1}^{U_j} s_{ju}\!\left(\tilde s_{ju}\!\left(\mathbf x_i^{\tilde S_{ju}}\right)\right), \qquad j=1,\dots,m.
$$
Here, $f_{jk}(\cdot)$ are standard GAM smooths, $\tilde s_{ju}(\cdot)$ are covariate transformations, and $s_{ju}(\cdot)$ are outer smooth effects applied to the transformed covariate [2511.19234].

The crucial definitional point is that a transformation is a scalar-valued function
$$
\tilde s:\mathbb R^p \to \mathbb R
$$
parameterized by $\mathbf a$, with the requirement that it be sufficiently differentiable; the paper uses fourth-order differentiability to support the fitting machinery [2511.19234]. The input may nevertheless be high-dimensional and complex, including time series, spatial fields, lagged histories, or other high-dimensional objects [2511.19234].

The paper gives three main classes of transformations. The first is adaptive exponential smoothing:
$$
\tilde s_i = \omega_i \tilde s_{i-1} + (1-\omega_i)x_i,\qquad i\ge 1,
$$
with $\tilde s_0=x_0$ and $\omega_i = \phi(\tilde{\mathbf x}_i^\top \mathbf a)$, where $\phi$ is the logistic function [2511.19234]. The second is multivariate kernel smoothing for spatially misaligned covariates:
$$
\tilde s(\mathbf x_i) = \frac{\sum_{j\in\mathcal N_i} K_{\mathbf a}(\mathbf x_i,\mathbf x_j) z_j} {\sum_{q\in\mathcal N_i} K_{\mathbf a}(\mathbf x_i,\mathbf x_q)}.
$$
The third is linear combinations or single-index effects:
$$
\tilde s(\mathbf x_i)=\mathbf x_i^\top \mathbf a.
$$
The paper notes that $\mathbf x_i$ can itself be a model matrix, basis expansion, or lagged-history vector [2511.19234].

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 [2511.19234].

## 5. Joint estimation, uncertainty, and implicit differentiation

For the integrated-transformation framework, the posterior log-density is written as
$$
\mathcal L(\boldsymbol\zeta) = \log p(\boldsymbol\zeta\mid \mathbf y,\boldsymbol\lambda) = \sum_{i=1}^n \ell_i(\boldsymbol\zeta\mid y_i) -\frac12\sum_{g=1}^G \lambda_g \boldsymbol\zeta^\top \mathbf S_g \boldsymbol\zeta,
$$
up to a constant [2511.19234]. For fixed $\boldsymbol\lambda$, MAP estimation is done with Newton’s method [2511.19234].

Joint uncertainty is quantified by approximating the posterior of $\boldsymbol\zeta$ with a Gaussian centered at the MAP estimate:
$$
\boldsymbol\zeta \mid \mathbf y,\boldsymbol\lambda \approx \mathcal N(\hat{\boldsymbol\zeta},\mathbf V_{\boldsymbol\zeta}),
$$
where
$$
\mathbf V_{\boldsymbol\zeta} = (\hat{\boldsymbol{\mathcal I}}+\mathbf S_{\boldsymbol\lambda})^{-1}.
$$
This yields joint uncertainty for standard GAM coefficients, smoothing coefficients, transformation parameters, and nested-effect parameters [2511.19234]. For a nested effect $s(\tilde s(\mathbf x))$, the paper uses the Delta method to propagate uncertainty [2511.19234].

Smoothing parameter selection again proceeds by empirical Bayes through a Laplace approximation to the marginal likelihood:
$$
\mathcal V(\boldsymbol\lambda) = \log p(\mathbf y\mid \boldsymbol\lambda) = \log \int p(\mathbf y\mid \boldsymbol\zeta)p(\boldsymbol\zeta\mid \boldsymbol\lambda)\,d\boldsymbol\zeta,
$$
with approximation
$$
\tilde{\mathcal V}(\boldsymbol\lambda) = \mathcal L(\hat{\boldsymbol\zeta}) +\frac12 \log |\mathbf S^{\boldsymbol\lambda}|_+ -\frac12 \log |\mathbf H| +\frac{M_p}{2}\log(2\pi).
$$
The optimization is carried out over $\rho_g = \log \lambda_g$ to enforce positivity [2511.19234].

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 [2511.19234]. 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 [2511.19234].

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 [2511.19234].

## 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 [2005.10092]. The integrated-transformation paper demonstrates the later framework through applications to forecasting electricity net-demand in Great Britain and modelling house prices in London [2511.19234]. 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 [2511.19234]. 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 [2511.19234].

A separate usage of the name appears in "Factor Graph Accelerator for LiDAR-Inertial Odometry" [2209.02207], 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 [2209.02207]. 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 [1611.01534]. "GameFactory: Creating New Games with Generative Interactive Videos" concerns action-controlled scene-generalizable game video generation using pre-trained video diffusion models [2501.08325]. "GameGPT: Multi-agent Collaborative Framework for Game Development" concerns LLM-based automation of game development workflows [2310.08067]. 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 [2005.10092]. The other embeds transformations of complex covariates directly inside the model, with joint estimation, approximate Bayesian uncertainty quantification, and empirical Bayes smoothing selection [2511.19234]. The common methodological theme is the elevation of flexible structured components—whether weights or transformations—from external heuristics to estimable parts of the model.

Source: https://www.emergentmind.com/topics/gamfactory