5  Generalized Additive Model (GAM) based on Quantiles

ImportantDisclaimer

These packages (Note 1) are a one-person project undergoing rapid evolution. Backward compatibility (per Hadley Wickham) is provided as a courtesy rather than a guarantee.

Until further notice, these packages should

  • not be used as a basis for research grant applications,
  • not be cited as an actively maintained tool in a peer-reviewed manuscript,
  • not be used to support or fulfill requirements for pursuing an academic degree.

In addition, work primarily based on these packages (Note 1) should not be presented at academic conferences or similar scholarly venues.

Furthermore, a person’s ability to use these packages (Note 1) does not necessarily imply an understanding of their underlying mechanisms. Accordingly, demonstration of their use alone should not be considered sufficient evidence of expertise, nor should it be credited as a basis for academic promotion or advancement.

These statements do not apply to the contributors (Tip 1) to these packages (Note 1) with respect to their specific contributions.

These statements do not apply when the maintainer of these packages (Note 1), Tingting Zhan, is credited as the first author, the lead author, and/or the corresponding author in a peer-reviewed manuscript, or as the Principal Investigator or Co-Principal Investigator in a research grant application and/or a final research progress report.

These statements are advisory in nature and do not modify or restrict the rights granted under the GNU General Public License https://www.r-project.org/Licenses/.

Note

The examples in Chapter 5 require

library(hyper.gam)
library(survival)

5.1 Compute Aggregated Quantiles

We use the data example Ki67 from package groupedHyperframe (v0.4.0, GPL-2) in this non-spatial application.

Listing 5.1: Data: a hyper data frame Ki67q 
Ki67q = groupedHyperframe::Ki67 |>
  within.data.frame(expr = {
    x = y = NULL # remove x- and y-coords for non-spatial application
  }) |>
  aggregate2hyper(by = ~ patientID/tissueID) |>
  within(expr = {
    logKi67.q = quantile(logKi67, probs = seq.int(from = .01, to = .99, by = .01))
  }) |>
  aggregate() |>
  within(expr = {
    logKi67.q = logKi67.q |> do.call(what = pmean)
  })
# Hypercolumn(s) logKi67 created!
# Variable(s) tissueID removed from aggregation

The returned object Ki67q is a hyper data frame with

  • a numeric-hypercolumn of the aggregated sample quantiles logKi67.q per patientID. These quantiles are calculated at a pre-specified grid of probabilities \(\{p_k, k=1,\cdots,K \} \in [0,1]\). Note that the aggregation must be performed at the level of biologically independent clusters, e.g., ~patientID, to produce independent quantile predictors.
  • Metadata, including the outcome of interest, e.g., progression free survival PFS, Her2, HR, etc.
A hyper data frame Ki67q: aggregated quantiles 
Ki67q |>
  head()
# Hyperframe:
#   Tstage  PFS adj_rad adj_chemo histology  Her2   HR  node  race age patientID   logKi67 logKi67.q
# 1      2 100+   FALSE     FALSE         3  TRUE TRUE  TRUE White  66   PT00037 (anylist) (numeric)
# 2      1   22   FALSE     FALSE         3 FALSE TRUE FALSE Black  42   PT00039 (anylist) (numeric)
# 3      1  99+   FALSE        NA         3 FALSE TRUE FALSE White  60   PT00040 (anylist) (numeric)
# 4      1  99+   FALSE      TRUE         3 FALSE TRUE  TRUE White  53   PT00042 (anylist) (numeric)
# 5      1  112    TRUE      TRUE         3 FALSE TRUE  TRUE White  52   PT00054 (anylist) (numeric)
# 6      4   12    TRUE     FALSE         2  TRUE TRUE  TRUE Black  51   PT00059 (anylist) (numeric)

Readers are encouraged to learn more about the hyper data frame (hyperframe) from package spatstat.geom (Baddeley et al. 2015; Baddeley and Turner 2005) and each function in this pipeline from,

Sections in Earlier Chapters
Function Purpose Where
aggregate2hyper() to aggregate a data frame into a (grouped) hyper data frame Chapter 2
within(expr = quantile) to calculate the quantiles of each numeric-vector in the numeric-hypercolumn Section 3.3.3
aggregate() |> within(expr = pmean) to aggregate the quantiles over multiple tissues per patient by point-wise means Section 3.4

5.2 Estimate Integrand Surface

Linear quantile index (QI) (Equation 5.1) is a predictor in a functional generalized linear model (James 2002) for outcomes from the exponential family of distributions, or a linear functional Cox model (Gellar et al. 2015) for survival outcomes,

\[ \text{QI}_{i}=\int_{0}^{1} \beta(p)Q_i(p)dp \tag{5.1}\]

where \(Q_i(p)\) is the (aggregated) sample quantiles logKi67.q for the \(i\)-th subject, and \(\beta(p)\) is the unknown coefficient function to be estimated. Listing 5.2 fits a generalized additive model (gam) with integrated linear spline-based smoothness estimation using the function s() (Wood 2003, v1.9.4, GPL (>= 2)). This is a scalar-on-function model (Reiss et al. 2017) that predicts a scalar outcome (e.g., progression free survival time PFS[,1L]) using the aggregated quantiles function as a functional predictor.

Listing 5.2: Linear functional Cox model for survival outcome (Gellar et al. 2015) (Listing 5.1) 
m0 = hyper_gam(PFS ~ logKi67.q, data = Ki67q)

Nonlinear quantile index (nlQI) (Equation 5.2) is a predictor in the functional generalized additive model (McLean et al. 2014) for outcomes from the exponential family of distributions, or an additive functional Cox model (Cui et al. 2021) for survival outcomes.

\[ \text{nlQI}_{i}= \int_{0}^{1} F\big(p, Q_i(p)\big)dp \tag{5.2}\]

where \(F(\cdot,\cdot)\) is an unknown bivariate twice differentiable function. Listing 5.3 fits a generalized additive model (gam) with tensor product interaction estimation using the function ti() (Wood 2006, v1.9.4, GPL (>= 2)).

Listing 5.3: Additive functional Cox model for survival outcome (Cui et al. 2021) (Listing 5.1) 
m1 = hyper_gam(PFS ~ logKi67.q, data = Ki67q, nonlinear = TRUE)

The fitted functional model m0 and m1 have the S3 class 'hyper_gam' (Chapter 27).

5.2.1 Visualization

The function integrandSurface() creates an interactive htmlwidget (Vaidyanathan et al. 2023, v1.6.4, MIT + file LICENSE) visualization of the estimated integrand surfaces for the linear (Equation 5.1) or nonlinear quantile index (Equation 5.2) using package plotly (Sievert 2020, v4.12.0, MIT + file LICENSE). The integrand surfaces, defined on \(p\in[0,1]\) and \(q\in\text{range}\big\{Q_i(p), i=1,\cdots,n\big\}\), are

\[ \hat{S}(p,q) = \begin{cases} \hat{\beta}(p)\cdot q\\ \hat{F}(p,q) \end{cases} \tag{5.3}\]

Also in this interactive visualization are

  • the contour lines on the integrand surfaces (Equation 5.3), as well as their projections along the \(s\)-axis, i.e., onto the \((p,q)\)-plane (a.k.a., the “floor”);
  • the estimated linear integrand paths \(\hat{\beta}(p)Q_i(p)\) or the nonlinear integrand paths \(\hat{F}(p, Q_i(p))\) on the integrand surfaces (Equation 5.3);
  • the sample quantiles \(Q_i(p)\), i.e., the projections of the estimated linear or nonlinear integrand path along the \(s\)-axis, i.e., onto the \((p,q)\)-plane (a.k.a., the “floor”);
  • the projections of the estimated linear or nonlinear integrand path along the \(q\)-axis, i.e., onto the \((p,s)\)-plane (a.k.a., the “backwall”), so that the area under each projected path is equal to the estimated linear (Equation 5.1) or nonlinear quantile index (Equation 5.2).

Figure 5.1 is a collage (via package htmltools, Cheng et al. 2025, v0.5.9, GPL (>= 2)) of the interactive htmlwidget visualizations of the linear and nonlinear integrand surfaces, contours, integrand paths and their projections to the “floor” and “backwall”. Listing 5.4 uses n=21L to reduce the htmlwidget objects size, in order to comply with Quarto Pub file size limit.

Listing 5.4: Figure: visualize linear and nonlinear quantile index (Listing 5.2, Listing 5.3)
Code 
scene = list(
  xaxis = list(title = 'Probability (p)', tickformat = '.0%', color = 'dodgerblue'), 
  yaxis = list(title = 'Quantile (q)', color = 'deeppink'),
  zaxis = list(title = 'Integrand (s)', color = 'darkolivegreen')
)
htmltools::tagList(
  m0 |> integrandSurface(n = 21L) |> plotly::layout(scene = scene), 
  m1 |> integrandSurface(n = 21L) |> plotly::layout(scene = scene)
) |> 
  htmltools::browsable()
Figure 5.1: Linear (top) and nonlinear (bottom) integrand surfaces, contours, integrand paths and their projections to the “floor” and “backwall”

Static illustrations of the estimated integrand surfaces, e.g., the perspective and contour plots (Section 27.2), are produced by calling the S3 generic functions graphics::persp() and graphics::contour() in package graphics (R version 4.5.3 (2026-03-11)). These figures are suppressed to reduce the file size of this vignette.

Listing 5.5: Static figures using package graphics
Code 
m0 |> persp()
m0 |> contour()
m1 |> persp()
m1 |> contour()

5.3 Compute Quantile Index Predictor

Linear and nonlinear quantile indices are the predictors in the functional models (Equation 5.1) and (Equation 5.2), respectively. Let’s consider a conventional scenario that we first fit a hyper_gam model to the training data set, then compute the quantile index predictors in the training and/or test data set using the training model.

First, we partition the 622 patients in hyper data frame Ki67q into a training data set with 498 patients and a test data set with 124 patients, i.e., a 80% versus 20% partition.

set.seed(16); id = Ki67q |> nrow() |> seq_len() |> caret::createDataPartition(p = .8)
Ki67q_0 = Ki67q[id[[1L]],] # training set
Ki67q_1 = Ki67q[-id[[1L]],] # test set

Next, we fit a functional generalized additive model to the the training data set Ki67q_0,

Listing 5.6: Functional generalized additive model 
m1a = hyper_gam(PFS ~ logKi67.q, nonlinear = TRUE, data = Ki67q_0)

We can, but we should not, use the quantile index predictors (Section 27.3) of the training data set for downstream analysis, because these quantile index predictors are optimized on the training data set and the results would be optimistically biased.

Optimistically biased!! 
Ki67q_0[,c('PFS', 'age', 'race')] |> 
  as.data.frame() |> # invokes spatstat.geom::as.data.frame.hyperframe()
  data.frame(nlQI = predict(m1a, newdata = Ki67q_0)) |>
  coxph(formula = PFS ~ age + nlQI, data = _)
# Call:
# coxph(formula = PFS ~ age + nlQI, data = data.frame(as.data.frame(Ki67q_0[, 
#     c("PFS", "age", "race")]), nlQI = predict(m1a, newdata = Ki67q_0)))
# 
#           coef exp(coef)  se(coef)      z       p
# age  -0.023320  0.976950  0.008321 -2.803 0.00507
# nlQI  1.118713  3.060911  0.343152  3.260 0.00111
# 
# Likelihood ratio test=23.24  on 2 df, p=8.993e-06
# n= 498, number of events= 99

Instead, we should use the quantile index predictors computed in the test data set for downstream analysis,

Ki67q_1[,c('PFS', 'age', 'race')] |> 
  as.data.frame() |> # invokes spatstat.geom::as.data.frame.hyperframe()
  data.frame(nlQI = predict(m1a, newdata = Ki67q_1)) |>
  coxph(formula = PFS ~ age + nlQI, data = _)
# Call:
# coxph(formula = PFS ~ age + nlQI, data = data.frame(as.data.frame(Ki67q_1[, 
#     c("PFS", "age", "race")]), nlQI = predict(m1a, newdata = Ki67q_1)))
# 
#          coef exp(coef) se(coef)      z     p
# age  -0.01636   0.98377  0.01862 -0.879 0.380
# nlQI  1.21050   3.35516  0.75858  1.596 0.111
# 
# Likelihood ratio test=3.42  on 2 df, p=0.1805
# n= 124, number of events= 19